engage/math/Mmn_tpl.h

using namespace std;

#include <cstdio>

#include "Mmn.h"

/****************************************************************
 * die Klasse "Mmn"                                           *
 ****************************************************************/
//------ Kontruktoren -------------------------------------------
template <class K>
Mmn<K>::Mmn(K* a, unsigned m, unsigned n) : m(m), n(n), Rg(0), Det_F(false) {
    _a=new (K*)[m];

    for(unsigned i=0; i<m; i++) {
        _a[i]=new K[n];
        for(unsigned j=0; j<n; j++)
            _a[i][j]=*(a+n*i+j);
    }
}

template <class K>
Mmn<K>::Mmn(unsigned m) : m(m), n(m), Rg(0), Det_F(false) {
    _a=new (K*)[m];

    for(unsigned i=0; i<m; i++) {
        _a[i]=new K[n];
        for(unsigned j=0; j<n; j++)
            _a[i][j]=(j==i)?1:0;
    }
}

template <class K>
Mmn<K>::Mmn(unsigned m, unsigned n) : m(m), n(n), Rg(0), Det_F(false) {
    _a=new (K*)[m];
    for(unsigned i=0; i<m; i++) {
        _a[i]=new K[n];
        for(unsigned j=0; j<n; j++)
            _a[i][j]=0;
    }
}

template <class K>
Mmn<K>::Mmn(const Mmn<K>& A) 
   : m(A.m), n(A.n), Rg(A.Rg), Det(A.Det), Det_F(A.Det_F) {
    _a=new (K*)[m];

    for(unsigned i=0; i<m; i++) {
        _a[i]=new K[n];
        for(unsigned j=0; j<n; j++)
            _a[i][j]=A._a[i][j];
    }
}

template <class K>
Mmn<K>::~Mmn() {
    freeall();
}
//---------------------------------------------------------------

//------ Operatoren ---------------------------------------------
template <class K>
const Mmn<K>& Mmn<K>::operator=(const Mmn<K>& A) {
    if(this==&A)
        return *this;

    if(!_a || m != A.m || n != A.n) {
        freeall();

        _a=new (K*)[A.m];
        for(unsigned i=0; i<A.m; i++)
            _a[i]=new K[A.n];
    }

    m=A.m;
    n=A.n;
    Rg=A.Rg;
    Det=A.Det;
    Det_F=A.Det_F;

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            _a[i][j]=A._a[i][j];

    return *this;
}

template <class K>
Mmn<K> Mmn<K>::operator*(const K& a) const {
    Mmn<K> A(m,n);

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            A._a[i][j]=_a[i][j]*a;

    return A;
}

template <class K>
Mmn<K> Mmn<K>::operator%(const Mmn<K>& A) const {
    if(n != A.m)
        return null();

    Mmn<K> B(m,A.n);

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<A.n; j++) {
            B._a[i][j]=0;

            for(unsigned ji=0; ji<n; ji++)
                B._a[i][j]+=_a[i][ji]*A._a[ji][j];
        }

    return B;
}

template <class K>
Mmn<K> Mmn<K>::operator+(const Mmn<K>& A) const {
    if(n != A.n || m != A.m)
        return null();

    Mmn<K> B(m,n);

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            B._a[i][j]=_a[i][j]+A._a[i][j];

    return B;
}

template <class K>
Mmn<K> Mmn<K>::operator-(const Mmn<K>& A) const {
    Mmn<K> B(m,n);

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            B._a[i][j]=-A._a[i][j];

    return *this+B;
}
//---------------------------------------------------------------

//------ Hilfsfunktionen ----------------------------------------
template <class K>
void Mmn<K>::freeall(void) {
    if(_a) {
        for(unsigned i=0; i<m; i++) {
            delete [] _a[i];
        }
        delete [] _a;
        _a=NULL;
    }

    m=n=Rg=0;
    Det_F=false;
}
//---------------------------------------------------------------

//------ Determinantenberechnung --------------------------------
template <class K>
K Mmn<K>::sarrus(void) const {
    return _a[0][0]*_a[1][1]*_a[2][2]\
          +_a[0][1]*_a[1][2]*_a[2][0]\
          +_a[0][2]*_a[1][0]*_a[2][1]\
          -_a[0][1]*_a[1][0]*_a[2][2]\
          -_a[0][0]*_a[1][2]*_a[2][1]\
          -_a[0][2]*_a[1][1]*_a[2][0];
}

template <class K>
K Mmn<K>::laplace(void) const {
    K detA=0;
    Mmn<K> A;

    // analytisches optimieren des Laplace-Algorithmus, durch finden der
    // Zeile/Spalte mit den meisten Nullen.
    //-------------------------------------------------------------------
    unsigned nulls_in_line[m];
    for(unsigned i=0; i<m; nulls_in_line[i++]=0);
    unsigned nulls_in_col[n];
    for(unsigned j=0; j<m; nulls_in_col[j++]=0);
    unsigned max_nulls=0;
    int max_nulls_line=0;
    int max_nulls_col=0;

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            if(_a[i][j]==(K)0) {
                nulls_in_line[i]++;
                nulls_in_col[j]++;

                if(nulls_in_line[i]>max_nulls)
                    max_nulls=nulls_in_line[i];
                if(nulls_in_col[i]>max_nulls)
                    max_nulls=nulls_in_col[i];

                // ist eine komplette Zeile oder Spalte 0, so koennen wir
                // direkt abbrechen und 0 zurueckgeben.
                if(max_nulls==m)
                    return (K)0;

                max_nulls_line=(max_nulls==nulls_in_line[i])?i:max_nulls_line;
                max_nulls_col=(max_nulls==nulls_in_col[j])?j:max_nulls_col;
            }

    if(nulls_in_line[max_nulls_line] > nulls_in_col[max_nulls_col]) {
        A=T(); 
        max_nulls_col=max_nulls_line;
    }
    else
        A=*this;
    //-------------------------------------------------------------------

    for(unsigned i=0; i<A.m; i++)
        if(A._a[i][0]!=(K)0)
            detA+=((K)(int)pow((double)-1,(double)(i+max_nulls_col)))
                  *A._a[i][max_nulls_col]
                  *A.Aij(i,max_nulls_col).det();

    return detA;
}

template <class K>
K Mmn<K>::det(void) const {
    K Det=(K)0;

    if(m!=n) {
        return Det;
    }
    else {
        switch(m) {
            case 1: Det=_a[0][0];
                    break;
            case 2: Det=_a[0][0]*_a[1][1]-_a[0][1]*_a[1][0];
                    break;
            case 3: Det=sarrus();
                    break;
            default: Det=laplace();
        }

        return Det;
    }
}

template <class K>
K Mmn<K>::det(void) {
    if(m!=n) {
        Det=(K)0;
        Det_F=false;

        return Det;
    }
    else {
        switch(m) {
            case 1: Det=_a[0][0];
                    break;
            case 2: Det=_a[0][0]*_a[1][1]-_a[0][1]*_a[1][0];
                    break;
            case 3: Det=sarrus();
                    break;
            default: Det=laplace();
        }

        Det_F=true;
        return Det;
    }
}
//---------------------------------------------------------------

//------ Matrizenmanipulationen ---------------------------------
template <class K>
Mmn<K> Mmn<K>::Aij(unsigned i, unsigned j) const {
    Mmn<K> A(m-1,n-1);
    unsigned _i, _l, _j, _k;

    for(_i=0, _l=0; _i<m; _i++, _l++) {
        if(_i==i) _i++;
        if(_i>=m) break;

        for(_j=0, _k=0; _j<n; _j++, _k++) {
            if(_j==j) _j++;
            if(_j>=n) break;

            A._a[_l][_k]=_a[_i][_j];
        }
    }

    return A;
}

template <class K>
Mmn<K> Mmn<K>::Ai(unsigned i) const {
    Mmn<K> A(m-1,n);
    unsigned _i, l;

    for(_i=0, l=0; i<m; _i++, l++) {
        if(_i==i) _i++;
        if(_i>=m) break;

        for(unsigned j=0; j<n; j++)
            A._a[l][j]=_a[_i][j];
    }

    return A;
}

template <class K>
Mmn<K> Mmn<K>::Aj(unsigned j) const {
    Mmn<K> A(m,n-1);
    unsigned _j, k;

    for(unsigned i=0; i<m; i++)
        for(_j=0, k=0; _j<n; _j++, k++) {
            if(_j==j) _j++;
            if(_j>=n) break;

            A._a[i][k]=_a[i][_j];
        }

    return A;
}

template <class K>
Mmn<K> Mmn<K>::T(void) const {
    Mmn<K> A(n,m);

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            A._a[j][i]=_a[i][j];

    return A;
}

template <class K>
Mmn<K> Mmn<K>::Ad_det(void) const {
    Mmn<K> A(n,m);

    for(unsigned i=0; i<m; i++)
        for(unsigned j=0; j<n; j++)
            A._a[j][i]=((K)(int)pow((double)-1,(double)(i+j)))*Aij(i,j).det();

    return A;
}

template <class K>
Mmn<K> Mmn<K>::Ad_inv(void) const {
    // Das geht nur wenn die Determinante nicht 0 ist, da sonst die ganze
    // Ergebnismatrix automatisch 0 wird, was aber nicht unbedingt sein muß.
    // Dies gilt insbesondere Wenn die Matrize eine Null-Zeile oder Spalte
    // enthaelt.
    K d=Det_F?Det:det();

    if(d != (K)0)
        return d*gauss(INVERSE);
    else
        return Ad_det();
}

template <class K>
Mmn<K> Mmn<K>::inv_Ad(void) const {
    K Adet; 

    if((Adet=Det_F?Det:det())==(K)0) {
        return null();
    }
    else {
        return Ad_det()*((K)1/Adet);
    }
}

template <class K>
Mmn<K> Mmn<K>::gauss(gauss_typ T) const {
    Mmn<K> A=*this;
    Mmn<K> B(m);

    for(unsigned i=0; i<A.m; i++) {
        unsigned j;

        for(j=i; j<A.n; j++) {
            unsigned k;

            for(k=i; k<A.m && A._a[k][j]==(K)0; k++);
            if(k==A.m)
                continue;
            if(k!=i) {
                K* tmp=A._a[i];
                A._a[i]=A._a[k];
                A._a[k]=tmp;

                if(T==INVERSE) {
                    tmp=B._a[i];
                    B._a[i]=B._a[k];
                    B._a[k]=tmp;
                }
            }
            break;
        }

        if(j==A.m) break;

        if(A._a[i][j]!=(K)0) {
            K aii_inv=(K)1/A._a[i][j];
            A._a[i][j]=1;
            A.Rg++;
            B.Rg++;

            for(unsigned k=j+1; k<A.n; k++)
                A._a[i][k]*=aii_inv;
            if(T==INVERSE) {
                for(unsigned k=0; k<B.n; k++)
                    B._a[i][k]*=aii_inv;
            }

            for(unsigned k=0; k<A.m; k++) {
                if(k != i) {
                    K akj_neg=-A._a[k][j];

                    for(unsigned l=j; l<A.n; l++)
                        A._a[k][l]+=A._a[i][l]*akj_neg;
                    if(T==INVERSE) {
                        for(unsigned l=0; l<B.n; l++)
                            B._a[k][l]+=B._a[i][l]*akj_neg;
                    }
                }
            }
        }
    }

    if(T==TREPPENNORMALE)
        return A;
    else {
        // Es gibt nur eine Inverse wenn Rg(A)=m ist
        if(B.Rg==m)
            return B;
        else
            return null();
    }
}
//---------------------------------------------------------------

//------ lineare Gleichungssysteme ------------------------------
template <class K>
Mmn<K> Mmn<K>::solve(void) const {
    Mmn<K> A=gauss(TREPPENNORMALE);
    Mmn<K> B=Aj(n-1).gauss(TREPPENNORMALE);

    // Gibt es eine Loesung???
    if(B.Rg != A.Rg)
        return null();
    else {
        K l[B.n][B.n-B.Rg+1];
        K a[B.n][A.n];

        for(unsigned i=0, j=0; j<B.n && i<B.n; j++) {
            if(i>=A.m || A._a[i][j] != (K)1) {
                for(unsigned k=0; k<A.n; k++) {
                    if(k==j)
                        a[j][k]=(K)-1;
                    else
                        a[j][k]=(K)0;
                }
            }
            else {
                for(unsigned k=0; k<A.n; k++)
                    a[j][k]=A._a[i][k];
                i++;
            }
        }

        unsigned j=0;
        for (unsigned i=0; i<B.n; i++)
            l[i][j]=a[i][A.n-1];
        j++;

        for(unsigned k=0; k<B.n; k++) {
            if(a[k][k]==(K)-1) {
                for (unsigned i=0; i<B.n; i++)
                    l[i][j]=a[i][k];
                j++;
            }
        }

        return Mmn<K>((K*)l, B.n, B.n-B.Rg+1);
    }
}

template <class K>
Mmn<K> Mmn<K>::special_solve(void) const {
    Mmn<K> S=solve();
    Mmn<K> A(S.m,1);

    for(unsigned i=0; i<S.m; i++) {
        A._a[i][0]=(K)0;

        for(unsigned j=0; j<S.n; j++)
            A._a[i][0]+=S._a[i][j];
    }

    return A;
}

//---------------------------------------------------------------
/****************************************************************/