#define mat_type double

int stepen(int a, int b) {
int r = 1;
for (int l = 1; l <= b; l++)
r *= a;

return r;
}

class Vektor {
private:
int n; //razmer vektora
mat_type *data;

public:
Vektor(int in);
Vektor(int in, float a[]);
Vektor(Vektor *a);
~Vektor();

int getSize();
mat_type get(int x);

void put(int x, mat_type d);
void print();
};

Vektor::Vektor(int in) {
n = in;

if ((data = (mat_type *) malloc(n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

mat_type *temp = data;

for(int i = 0; i < n; i++) {
*temp = 0;
temp++;
}
}

Vektor::Vektor(int in, float a[]) {
n = in;

if ((data = (mat_type *) malloc(n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

mat_type *temp = data;

for(int i = 0; i < n; i++) {
*temp = a[i];
temp++;
}
}

Vektor::Vektor(Vektor *a) {
n = a->getSize();

if ((data = (mat_type *) malloc(n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

for(int i = 1; i <= n; i++) {
*data = a->get(i);
data++;
}

data -= n;
}

Vektor::~Vektor() {
free(data);
}

int Vektor::getSize() {
return n;
}

mat_type Vektor::get(int x) {
return *(data + x - 1);
}

void Vektor::put(int x, mat_type d) {
if ((x < 1) || (x > n)) exit(1);

*(data + x - 1) = d;
}

void Vektor::print() {
cout << "(";
for (int i = 1; i <= n; i++) {
cout << *data;
if (i != n) cout << ", ";
data++;
}
cout << ")";
data -= n;
}

class Matrix {
private:
int n; //razmer matrici
mat_type *data;

public:
Matrix(int in); //input matrix size for get mem
Matrix(Matrix *a);
~Matrix();

void put(int x, int y, mat_type el);
int getSize();
void trans();

mat_type get(int x, int y);
mat_type minor(int x, int y); //absolyutni sdvig
mat_type norma();

void swapCol(int y, int y1);
void swapRow(int x, int x1);

float det();
int getMax(int y); //Index max elem v stroke y
Matrix *copy();
Matrix *obrat();

void print();
};

Matrix::Matrix(int in) {
n = in;

if ((data = (mat_type *) malloc(n * n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

mat_type *temp = data;

for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++) {
*temp = 0;
temp++;
}

}

Matrix::Matrix(Matrix *a) {
n = a->getSize();

if ((data = (mat_type *) malloc(n * n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

for (int j = 1; j<=n; j++)
for (int i = 1; i<=n; i++) {
*data = a->get(i, j);
data++;
}

data -= n * n;
}

Matrix::~Matrix() {
free(data);
}

void Matrix::put(int x, int y, mat_type el) {
*(data + (y - 1) * n + (x - 1)) = el;
}

mat_type Matrix::get(int x, int y) {
return *(data + (y - 1) * n + (x - 1));
}

float Matrix::det() {
if (n == 1) return get(1, 1);

mat_type d = 0;
int y;

Matrix new_mat(n - 1);

for(int l = 1; l <= n; l++) {

for(int j = 1; j <= n-1; j++)
for(int i = 1; i <= n-1; i++)
if (i < l) new_mat.put(i, j, get(i, j+1));
else new_mat.put(i, j, get(i+1, j+1));

d += stepen(-1, l + 1) * get(l, 1) * new_mat.det();
}

return d;
}

int Matrix::getMax(int y) {
int max = 1;
mat_type *temp = data + (y - 1) * n;

for(int j = 1; j <= n; j++) {
if (get(max, y) < *temp) max = j;
temp++;
}

return max;
}

mat_type Matrix::minor(int x, int y) {

int sx, sy;
Matrix temp(n - 1);

sx = 0;
sy = 0;
for(int j = 1; j <= n-1; j++) {
if (j >= y) sy = 1;
sx = 0;
for(int i = 1; i <= n-1; i++) {
if (i >= x) sx = 1;
temp.put(i, j, get(i + sx, j + sy));
}
}

return temp.det();
}

void Matrix::swapCol(int x, int x1) {
mat_type buf;
mat_type *temp = data + (x - 1);
mat_type *temp1 = data + (x1 - 1);

for(int j = 0; j < n; j++) {
buf = *temp;
*temp = *temp1;
*temp1 = buf;
temp += n;
temp1 += n;
}
}

void Matrix::swapRow(int y, int y1) {
mat_type buf;
mat_type *temp = data + (y - 1) * n;
mat_type *temp1 = data + (y1 - 1) * n;

for(int i = 0; i < n; i++) {
buf = *temp;
*temp = *temp1;
*temp1 = buf;
temp++;
temp1++;
}
}

int Matrix::getSize() {
return n;
}

void Matrix::trans() {
mat_type d;

for (int i = 2; i <= n; i++)
for (int j = i-1; j >= 1; j--)
if (i != j) {
d = get(i, j);
put(i, j, get(j, i));
put(j, i, d);
}

}

Matrix *Matrix::obrat() {
Matrix *S = new Matrix(n);
float d = det();

for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
S->put(i, j, minor(i, j) / d * stepen(-1, i+j));

S->trans();
return S;
}

void Matrix::print() {

mat_type *temp = data;

for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++) {
cout << *temp << "\t";
temp++;
}
cout << "\n";
}
}

mat_type Matrix::norma() {
mat_type r = 0;

for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
r += get(i, j) * get(i, j);

return sqrt(r);
}

Vektor *minVV(Vektor *a, Vektor *x) {
if (a->getSize() != x->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int n = a->getSize();

Vektor *b = new Vektor(n);

for (int j = 1; j <= n; j++) {
b->put(j, a->get(j) - x->get(j));
}

return b;
}

Vektor *mulMV(Matrix *a, Vektor *x) {
if (a->getSize() != x->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int n = a->getSize();
mat_type summ;

Vektor *b = new Vektor(n);

for (int j = 1; j <= n; j++) {
summ = 0;
for (int i = 1; i <= n; i++) {
summ += x->get(i) * a->get(i, j);
}
b->put(j, summ);
}

return b;
}

Matrix *mulMM(Matrix *a, Matrix *b) {
if (a->getSize() != b->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int n = a->getSize();

Matrix *c = new Matrix(n);
mat_type summ;

for (int j = 1; j <= n; j++)
for (int i = 1; i <= n; i++) {
summ = 0;

for (int k = 1; k <= n; k++)
summ += a->get(k, j) * b->get(i, k);

c->put(i, j, summ);
}

return c;
}

Vektor *metodKK(Matrix *A, Vektor *c) {

int n = A->getSize();
Vektor *b = new Vektor(c);
Matrix *S = new Matrix(n);

mat_type summ = 0;

for (int j = 1; j <= n; j++)
for (int i = 1; i <= n; i++) {
if (i == j) {
summ = 0;

for (int p = 1; p <= i-1; p++)
summ += S->get(p, i) * S->get(p, i);

S->put(i, i, sqrt(abs(A->get(i, i) - summ)));
}

if (j > i) {
summ = 0;

for (int p = 1; p <= i-1; p++)
summ += S->get(p, i) * S->get(p, j);

S->put(i, j, (A->get(i, j) - summ) / S->get(i, i));

}
if (i > j) S->put(i, j, 0);
}

Matrix *St = new Matrix(S);
St->trans();

for (int i = 1; i <=n; i++) {
b->put(i, b->get(i) / S->get(i, i));
S->put(i, i, 1);
for (int j = i+1; j <=n; j++) {
b->put(j, b->get(j) - b->get(i) * S->get(i, j));
S->put(i, j, 0);
}
}

for (i = n; i >= 1; i--) {
b->put(i, b->get(i) / St->get(i, i));
St->put(i, i, 1);

for (j = i-1; j >= 1; j--) {
b->put(j, b->get(j) - b->get(i) * St->get(i, j));
St->put(i, j, 0);
}
}

delete S;

return b;
}

Vektor *metodG(Matrix *r, Vektor *c) {
int n = r->getSize();

Vektor *b = new Vektor(c);
Matrix *a = new Matrix(r);

for(int j = 1; j <= n; j++) {

a->swapCol(j, a->getMax(j));

b->put(j, b->get(j) / a->get(j, j));

for(int i = n; i >= 1; i--)
a->put(i, j, a->get(i, j) / a->get(j, j));

for(int k = j + 1; k <= n; k++) {
float ak = a->get(j, k);

b->put(k, b->get(k) - ak * b->get(j));

for(int i = 1; i <= n; i++)
if (i > j)
a->put(i, k, a->get(i, k) - ak * a->get(i, j));
else
a->put(i, k, 0);
}
}

// obratnii hod

for (int i = n; i >=2; i--)
for (int j = i-1; j >= 1; j--) {
b->put(j, b->get(j) - b->get(i) * a->get(i, j));
a->put(i, j, 0);
}

delete a;

return b;
}

void debug(Matrix *a, Vektor *x) {

if (a->getSize() != x->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int y;

Vektor *b = mulMV(a, x);

Vektor *xg = metodG(a, b);
Vektor *ag = mulMV(a, xg);
Vektor *rg = minVV(b, ag);

Vektor *xk = metodKK(a, b);
Vektor *ak = mulMV(a, xk);
Vektor *rk = minVV(b, ak);

Matrix *a1 = a->obrat();
Matrix *c = mulMM(a, a1);
float cond = a->norma() + a1->norma();

cout << "Matrix A:\n"; a->print(); cout << "\n";
cout << "Matrix A^(-1): \n"; a1->print(); cout << "\n";
cout << "Matrix C = A * A^(-1): \n"; c ->print(); cout << "\n";

cout << "Obuslovlennost: " << cond << ", ";
if (cond > 10) cout << "ploho";
else cout << "horosho";
cout << " obuslovlena.\n";

cin >> y;

clrscr();
cout << "Matrix A:\n"; a->print(); cout << "\n";
cout << "Vektor x-tochnii: "; x->print(); cout << "\t";

cout << "b = A * x = "; b->print(); cout << "\n";

cout << "Metod Gaussa: "; xg->print(); cout << "\n";
cout << "\tVektor nevyazki: "; rg->print(); cout << "\n";

cout << "Metod kv.kornya: "; xk->print(); cout << "\n";
cout << "\tVektor nevyazki: "; rk->print(); cout << "\n";

delete b;
delete xg;
delete ak;
delete ag;
delete rk;
delete rg;
delete xk;
delete a1;
delete c;
}

int main() {
clrscr();

unsigned long StartMem = coreleft();

int n = 4;

// float v[3] = {531, -460, 193};
float v[4] = {1, 1, 1, 1};

Vektor *x = new Vektor(4, v);

Matrix *a1 = new Matrix(n);
Matrix *a2 = new Matrix(4);
Matrix *a3 = new Matrix(4);

for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++) {
a1->put(i, j, 1/(i+j-1));
}

/* a1->put(1, 1, 81); a1->put(2, 1, -45); a1->put(3, 1, 45);
a1->put(1, 2, -45); a1->put(2, 2, 50); a1->put(3, 2, -15);
a1->put(1, 3, 45); a1->put(2, 3, -15); a1->put(3, 3, 38);*/

a2->put(1, 1, 1.06); a2->put(2, 1, 2.14); a2->put(3, 1, 0.08); a2->put(4, 1, 1.17);
a2->put(1, 2, 3.13); a2->put(2, 2, 4.17); a2->put(3, 2, 0.10); a2->put(4, 2, 9.12);
a2->put(1, 3, 5.04); a2->put(2, 3, 2.20); a2->put(3, 3, 6.11); a2->put(4, 3, 7.19);
a2->put(1, 4, 1.14); a2->put(2, 4, 3.19); a2->put(3, 4, 2.22); a2->put(4, 4, 0.07);

a3->put(1, 1, 5.88); a3->put(2, 1, 0.04); a3->put(3, 1, 1.28); a3->put(4, 1, 0.82);
a3->put(1, 2, 0.04); a3->put(2, 2, 7.87); a3->put(3, 2, 0.12); a3->put(4, 2, 1.32);
a3->put(1, 3, 1.28); a3->put(2, 3, 0.12); a3->put(3, 3, 3.19); a3->put(4, 3, 0.82);
a3->put(1, 4, 0.82); a3->put(2, 4, 1.32); a3->put(3, 4, 0.82); a3->put(4, 4, 6.22);

debug(a3, x);

delete a1;
delete a2;
delete a3;
delete x;

cout << "\nLost memory: " << StartMem - coreleft();

return 0;
};

�����