#include <stdio.h>
#include <stdlib.h>
//#include <i86.h>
#include <math.h>
//#include <conio.h>
class matrix
{
private:
    int row, col;
    int sparse;
    double * *y;
public:
    matrix(int n = 1, int m = 1, int s = 0);
    matrix(const matrix &);
    ~matrix();
    matrix operator + (const matrix & a) const;
    matrix operator - (const matrix & a) const;
    void print();
    double & operator()(int n, int m) const;
    double & operator()(int n) const;
    matrix operator = (const matrix & a);
    matrix operator * (const matrix & a) const;
    matrix operator | (const matrix & b) const;
    matrix operator + (); //unary operator for upper triangulation
    double det();
    double norm() const;
    matrix column(int);
    friend matrix operator * (double, matrix &);
    friend matrix operator * (matrix &, double);

    friend int eigJ(matrix &a, double eps, double *, matrix &);
};

//              end class matrix

matrix::matrix(int n, int m, int s) {
    row = n;
    col = m;
    int j, i;
    sparse = s;
    y = new double*[row];
    for (i = 0; i < n; i++) y[i] = new double[col];
    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            y[i][j] = 0.0;
}

matrix::matrix(const matrix &a) {
    row = a.row;
    col = a.col;
    int j, i;
    sparse = a.sparse;
    y = new double*[row];
    for (i = 0; i < row; i++) y[i] = new double[col];
    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            y[i][j] = a.y[i][j];
}

matrix::~matrix() {
    for (int i = 0; i < row; i++)delete y[i];
    delete y;
}

matrix matrix::operator = (const matrix &a){
    int j, i;
    for (i = 0; i < row; i++) delete y[i];
    delete y;
    row = a.row;
    col = a.col;
    y = new double*[row];
    for (i = 0; i < row; i++) y[i] = new double[col];
    for (i = 1; i <= row; i++)
        for (j = 1; j <= col; j++) {
            (*this)(i, j) = a(i, j);
        }
    return *this;
}

matrix matrix::operator + (const matrix & a) const {
    if ((row != a.row) || (col != a.col)) {
        //sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator+");
        exit(1);
    }
    int j;
    matrix t(row, col);
    for (int i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            t.y[i][j] = y[i][j] + a.y[i][j];
    return t;
}

matrix matrix::operator - (const matrix & a) const {
    if ((row != a.row) || (col != a.col)) {
        //sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator+");
        exit(1);
    }
    matrix t(row, col);
    int j;
    for (int i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            t.y[i][j] = y[i][j] - a.y[i][j];
    return t;
}

void matrix::print() {
    int j;
    for (int i = 0; i < row; i++) {
        printf("\n");
        for (j = 0; j < col; j++) printf(" %5.2f", y[i][j]);
    }
    printf("\n");
}

matrix matrix::operator*(const matrix &a) const {
    if (col != a.row) {
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator*");
        exit(1);
    }
    matrix t(row, a.col);
    int j, k;
    for (int i = 1; i <= t.row; i++)
        for (j = 1; j <= t.col; j++)
            for (k = 1; k <= col; k++) t(i, j) = t(i, j)+(*this)(i, k) * a(k, j);
    return (t);
}

matrix matrix::operator | (const matrix &b) const {
    if (col != b.row) {
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator|\n");
        exit(1);
    }
    if (b.col > 1) {
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("number columns for second argument must be 1 for operator|\n");
        exit(1);
    }
    matrix a1(row, col), c(row, col), x(row, 1), b1(row, 1);
    int i, j, h, n1, k;
    double v, r;
    a1 = *this;
    b1 = b;
    c = +a1;
    n1 = row - 1;
    //recalculation b1
    for (i = 1; i <= n1; i++)
        for (h = i + 1; h <= row; h++)
            b1(h, 1) = b1(h, 1) - c(h, i) * b1(i, 1);
    //solve triangular system
    int n = row, i1;
    x(n, 1) = b1(n, 1) / a1(n, n);
    i = n1;
l110:
    v = 0.0;
    i1 = i + 1;
    for (h = i1; h <= n; h++)v += a1(i, h) * x(h, 1);
    x(i, 1) = (b1(i, 1) - v) / a1(i, i);
    i--;
    if (i >= 1)goto l110;
    return x;
}
//**********************member function for triangulation********************

matrix matrix::operator + () {
    matrix c(row, col);
    int i, j, h, n1, k;
    double v, r;
    n1 = row - 1;
    //triangulation
    for (i = 1; i <= n1; i++) {
        h = i + 1;
        if ((*this)(i, i) != 0)goto l50;
        //find nonzero element
        for (k = h; k <= row; k++)if ((*this)(k, i) != 0)goto l30;
        //nonzero element not found
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("nonzero element not found in Gauss elemination");
        exit(1);

        //replace rows
l30:
        for (j = i; j <= row; j++) {
            v = (*this)(i, j);
            (*this)(i, j) = (*this)(k, j);
            (*this)(k, j) = v;
        }
        //erase a(h,i)
l50:
        if ((*this)(h, i) == 0.0)goto l80;
        r = -(*this)(h, i) / (*this)(i, i);
        c(h, i) = -r;
        for (j = i; j <= row; j++)(*this)(h, j) = (*this)(h, j) + r * (*this)(i, j);
l80:
        h++;
        if (h <= row)goto l50;
    }
    //end triangulation
    return c;
}
//************************determinant**********************************

double matrix::det() {
    matrix a1(row, col);
    double d;
    int i, j, h, n1, k, num = 0;
    double v, r;
    a1 = *this;
    n1 = row - 1;
    //triangulation
    for (i = 1; i <= n1; i++) {
        h = i + 1;
        if (a1(i, i) != 0)goto l50;
        //find nonzero element
        for (k = h; k <= row; k++)if (a1(k, i) != 0)goto l30;
        //nonzero element not found
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("nonzero element not found in Gauss elemination");
        exit(1);

        //replace rows
l30:
        for (j = i; j <= row; j++) {
            v = a1(i, j);
            a1(i, j) = a1(k, j);
            a1(k, j) = v;
        }
        num++;
        //erase a(h,i)
l50:
        if (a1(h, i) == 0.0)goto l80;
        r = -a1(h, i) / a1(i, i);
        for (j = i; j <= row; j++)a1(h, j) = a1(h, j) + r * a1(i, j);
l80:
        h++;
        if (h <= row)goto l50;
    }
    //end triangulation
    d = 1;
    for (i = 1; i <= row; i++)d *= a1(i, i);
    d *= pow(-1, num);
    return d;
}

//*********************** norm ***************************************

double matrix::norm() const {
    double a = 0;
    for (int i = 1; i <= row; i++)a = a + (*this)(i, 1)*(*this)(i, 1);
    return sqrt(a);
}

matrix matrix::column(int n) {
    if (n > col) {
        printf("number column exceed matrix dimension\n");
        exit(1);
    }
    matrix t(row, 1);
    for (int i = 1; i <= row; i++)t(i, 1) = (*this)(i, n);
    return t;
}

matrix operator*(double c, matrix &a) {
    matrix t(a.row, a.col);
    int j;
    for (int i = 1; i <= a.row; i++)
        for (j = 1; j <= a.col; j++)
            t(i, j) = c * a(i, j);
    return t;
}

matrix operator*(matrix &a, double c) {
    matrix t(a.row, a.col);
    int j;
    for (int i = 1; i <= a.row; i++)
        for (j = 1; j <= a.col; j++)
            t(i, j) = c * a(i, j);
    return t;
}

double & matrix::operator()(int n, int m) const {
    if ((n > row) || (m > col) || (n < 1) || (m < 1)) {
        //        sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("index exceeds matrix dimension");
        exit(1);
    }
    return y[n - 1][m - 1];
}

double & matrix::operator()(int n) const {
    if ((n > row * col) || (n < 1)) {
        //        sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("index exceeds matrix length");
        exit(1);
    }
    int n1, m1;
    m1 = n / row + 1;
    n1 = n % row;
    return y[n1 - 1][m1 - 1];
}

int eigJ(matrix &a, double eps, double *d, matrix &v) {
    int n = a.row, i, j, k, p, q;
    matrix u(n, n), ut(n, n), a1(n, n);
    double max, s, fi, t;
    for (i = 1; i <= n; i++) v(i, i) = u(i, i) = ut(i, i) = 1.0;
    a1 = a;
    t = 1;
    int n1 = n * (n - 1) / 2;
    while (t > eps) {
        max = fabs(a1(1, 2));
        p = 1;
        q = 2;
        for (i = 1; i <= n; i++)
            for (j = 1; j <= n; j++)if (i != j)
                    if (fabs(a1(i, j)) > max) {
                        p = i;
                        q = j;
                        max = fabs(a1(i, j));
                    }
        fi = 0.5 * atan(2 * a1(p, q) / (a1(p, p) - a1(q, q)));
        u(p, p) = u(q, q) = ut(p, p) = ut(q, q) = cos(fi);
        u(p, q) = ut(q, p) = -sin(fi);
        u(q, p) = ut(p, q) = sin(fi);
        a1 = ut * a1*u;
        v = v*u;
        t = 0.0;
        for (i = 1; i <= n; i++)
            for (j = i + 1; j <= n; j++) t += fabs(a1(i, j));
        t /= n1;
        u(p, p) = u(q, q) = ut(p, p) = ut(q, q) = 1.0;
        u(p, q) = ut(q, p) = 0.0;
        u(q, p) = ut(p, q) = 0.0;
    }
    for (i = 0; i < n; i++)d[i] = a1(i + 1, i + 1);
    return 1;
}

void main() {
    matrix a(3, 3), c(3, 3), a4(3, 3), b(4, 1), err;
    matrix a1(2, 2), a2(2, 2), a3(2, 2);
    matrix a6(4, 4), v6(4, 4), eigv6(4, 1), x6(4, 1);
    a1(1, 1) = a1(1, 2) = 2;
    a1(2, 1) = a1(2, 2) = 1.0;
    a2 = a1 * 2;
    double d6[4], eps = 0.1e-11;
    int i, j;
    printf("input matrix a1 and a2\n");
    a1.print();
    a2.print();
    printf("summa a1+a2\n");
    a3 = a1 + a2;
    a3.print();
    printf("press any key for continue\n");
    //getch();
    /* a(1,1)=1; a(1,2)=2;a(1,3)=3;
    a(2,1)=2;a(2,2)=3;a(2,3)=5;
    a(3,1)=3;a(3,2)=5;a(3,3)=7;
    */
    a6(1, 1) = 1;
    a6(1, 2) = 2;
    a6(1, 3) = 3;
    a6(1, 4) = 2;
    a6(2, 1) = 2;
    a6(2, 2) = 3;
    a6(2, 3) = 5;
    a6(2, 4) = 3;
    a6(3, 1) = 3;
    a6(3, 2) = 5;
    a6(3, 3) = 7;
    a6(3, 4) = 5;
    a6(4, 1) = 2;
    a6(4, 2) = 3;
    a6(4, 3) = 5;
    a6(4, 4) = 0.762;
    i = eigJ(a6, eps, d6, v6);
    for (i = 0; i < 4; i++)printf("d= %e\n", d6[i]);
    printf("error in eig\n");
    for (i = 1; i <= 4; i++) {
        eigv6 = v6.column(i);
        err = a6 * eigv6 - d6[i - 1] * eigv6;
        err.print();
    }
    b(1, 1) = 3;
    b(2, 1) = 5;
    b(3, 1) = 11;
    b(4, 1) = 0.657;
    //Solve system of linear equations a*x=b
    x6 = a6 | b;
    x6.print();
    printf("error for Gauss\n");
    err = a6 * x6 - b;
    err.print();
    //getch();
    printf("triangulation");
    c = +a6;
    a6.print();
    int nn = 500;
    matrix tst(nn, nn);
    for (i = 1; i <= nn; i++)for (j = 1; j <= nn; j++)tst(i, j) = ((double) rand()) / RAND_MAX;
    tst.operator + ();
    double det = 1.0;
    for (i = 1; i <= nn; i++)det *= tst(i, i);
    printf("det=%e\n", det);
}