2013-team4/code/matrix

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

{{{
矩阵快速幂
}}}
{{{
typedef int DATA;
const int MAXN = 202;

struct mat{
    int n;
    DATA v[MAXN][MAXN];
    mat(int n, bool one = false){
        this->n = n;
        memset(v, 0, sizeof(v));
        if(one) for(int i = 0; i < n; i++) v[i][i] = 1;
    }
    mat operator*(const mat &a) const{
        mat ret(n);
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                for(int k = 0; k < n; k++)
                    ret.v[i][j] += v[i][k] * a.v[k][j];
        return ret;
    }
    mat operator^(int p) const{
        mat ret(n, true), a = (*this);
        while(p > 0){
            if(p & 1)
                ret = ret * a;
            p >>= 1;
            a = a * a;
        }
        return ret;
    }
};
}}}
{{{
列主元法求逆矩阵
}}}
{{{
#include <vector>
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cassert>
using namespace std;
typedef long double Double;
typedef vector<Double> VD;
typedef vector<VD> VVD;
#define ONE 1
#define ZERO 0
#define REP(i,n) for(int i=0; i<int(n); i++)
#define eps 1e-10
inline int sgn(double x){
    if(fabs(x)<eps)return 0;
    else return x>0?1:-1;
}

class Matrix{
public:
    VVD data;
    VD &operator[](size_t i){return data[i];}
    Matrix(size_t n, size_t m, int type=0){
        data=VVD(n,VD(m,0));
        switch(type){
            case ONE: {
                int nn=min(n,m);
                REP(i,nn)data[i][i]=1;
                break;
            }
            case ZERO: break;
            default: assert(type==ONE||type==ZERO);
        }
    }
    Matrix(const VVD &d):data(d){}
    Matrix(const Matrix &A):data(A.data){}

    size_t size(){return data.size();}

    void print(void){
        for(size_t i=0; i<data.size(); i++){
            for(size_t j=0; j<data[i].size(); j++)
                printf("%.2lf ",data[i][j]);
            printf("\n");
        }
    }

    bool isSquare(){
        if(data.size()==0)return true;
        else if(data.size()==data[0].size())return true;
        else return false;
    }

    //return value means whether this method is working properly
    bool transpose(void){
        if(!isSquare())return false;
        for(size_t i=0; i<data.size(); i++){
            for(size_t j=i+1; j<data.size(); j++)
                swap(data[i][j],data[j][i]);
        }
        return true;
    }

    void swapRow(size_t i, size_t j){swap(data[i],data[j]);}
    void mulRow(size_t i, Double times){REP(j,data[i].size())data[i][j]*=times;}
    //add times*Row i to row j
    void addMulRow(Double times, size_t i, size_t j)
    {REP(k,data[i].size())data[j][k]+=data[i][k]*times;}

    bool invert(Matrix &M){
        if(!isSquare())return false;
        int n=data.size();
        M=Matrix(n, n, ONE);
        REP(k, n){//kth column
            int row=k; Double maxv=fabs(data[k][k]);
            for(size_t i=k+1; i<n; i++)//Find the maximum element in the kth column
                if(fabs(data[i][k])>maxv)maxv=data[row=i][k];
            if(sgn(maxv)==0)return false;
            swapRow(k, row);
            M.swapRow(k, row);
            maxv=data[k][k];
            for(size_t i=0; i<n; i++){
                if(i==k)continue;
                Double t=-data[i][k]/maxv;
                addMulRow(t, k, i);
                M.addMulRow(t, k, i);
                assert(sgn(data[i][k])==0);
            }
            mulRow(k, 1/maxv);
            M.mulRow(k, 1/maxv);
        }
        return true;
    }
};

int main(){
    int n;
    scanf("%d",&n);
    Matrix A(n,n);
    REP(i,n)REP(j,n)
        cin>>A[i][j];
    A.transpose();
    Matrix M(0,0);
    if(A.invert(M)){
        /*A.print();
        printf("-------\n");
        M.print();*/
        bool flag=false;
        REP(i,M.size()){
            REP(j,M[i].size()){
                if(sgn(round(M[i][j])-M[i][j])!=0)flag=true;
            }
        }
        if(flag)printf("Power of magic saves lives\n");
        else printf("Death\n");
    } else {
        printf("Power of magic saves lives\n");
    }
}
}}}

矩阵快速幂

typedef int DATA;
const int MAXN = 202;
struct mat{
    int n;
    DATA v[MAXN][MAXN];
    mat(int n, bool one = false){
        this->n = n;
        memset(v, 0, sizeof(v));
        if(one) for(int i = 0; i < n; i++) v[i][i] = 1;
    }
    mat operator*(const mat &a) const{
        mat ret(n);
        for(int i = 0; i < n; i++)
            for(int j = 0; j < n; j++)
                for(int k = 0; k < n; k++)
                    ret.v[i][j] += v[i][k] * a.v[k][j];
        return ret;
    }
    mat operator^(int p) const{
        mat ret(n, true), a = (*this);
        while(p > 0){
            if(p & 1)
                ret = ret * a;
            p >>= 1;
            a = a * a;
        }
        return ret;
    }
};

列主元法求逆矩阵

#include <vector>
#include <iostream>
#include <cstdio>
#include <cmath>
#include <algorithm>
#include <cassert>
using namespace std;
typedef long double Double;
typedef vector<Double> VD;
typedef vector<VD> VVD;
#define ONE 1
#define ZERO 0
#define REP(i,n) for(int i=0; i<int(n); i++)
#define eps 1e-10
inline int sgn(double x){
    if(fabs(x)<eps)return 0;
    else return x>0?1:-1;
}
class Matrix{
public:
    VVD data;
    VD &operator[](size_t i){return data[i];}
    Matrix(size_t n, size_t m, int type=0){
        data=VVD(n,VD(m,0));
        switch(type){
            case ONE: {
                int nn=min(n,m);
                REP(i,nn)data[i][i]=1;
                break;
            }
            case ZERO: break;
            default: assert(type==ONE||type==ZERO);
        }
    }
    Matrix(const VVD &d):data(d){}
    Matrix(const Matrix &A):data(A.data){}
    size_t size(){return data.size();}
    void print(void){
        for(size_t i=0; i<data.size(); i++){
            for(size_t j=0; j<data[i].size(); j++)
                printf("%.2lf ",data[i][j]);
            printf("\n");
        }
    }
    bool isSquare(){
        if(data.size()==0)return true;
        else if(data.size()==data[0].size())return true;
        else return false;
    }
    //return value means whether this method is working properly
    bool transpose(void){
        if(!isSquare())return false;
        for(size_t i=0; i<data.size(); i++){
            for(size_t j=i+1; j<data.size(); j++)
                swap(data[i][j],data[j][i]);
        }
        return true;
    }
    void swapRow(size_t i, size_t j){swap(data[i],data[j]);}
    void mulRow(size_t i, Double times){REP(j,data[i].size())data[i][j]*=times;}
    //add times*Row i to row j
    void addMulRow(Double times, size_t i, size_t j)
    {REP(k,data[i].size())data[j][k]+=data[i][k]*times;}
    bool invert(Matrix &M){
        if(!isSquare())return false;
        int n=data.size();
        M=Matrix(n, n, ONE);
        REP(k, n){//kth column
            int row=k; Double maxv=fabs(data[k][k]);
            for(size_t i=k+1; i<n; i++)//Find the maximum element in the kth column
                if(fabs(data[i][k])>maxv)maxv=data[row=i][k];
            if(sgn(maxv)==0)return false;
            swapRow(k, row);
            M.swapRow(k, row);
            maxv=data[k][k];
            for(size_t i=0; i<n; i++){
                if(i==k)continue;
                Double t=-data[i][k]/maxv;
                addMulRow(t, k, i);
                M.addMulRow(t, k, i);
                assert(sgn(data[i][k])==0);
            }
            mulRow(k, 1/maxv);
            M.mulRow(k, 1/maxv);
        }
        return true;
    }
};
int main(){
    int n;
    scanf("%d",&n);
    Matrix A(n,n);
    REP(i,n)REP(j,n)
        cin>>A[i][j];
    A.transpose();
    Matrix M(0,0);
    if(A.invert(M)){
        /*A.print();
        printf("-------\n");
        M.print();*/
        bool flag=false;
        REP(i,M.size()){
            REP(j,M[i].size()){
                if(sgn(round(M[i][j])-M[i][j])!=0)flag=true;
            }
        }
        if(flag)printf("Power of magic saves lives\n");
        else printf("Death\n");
    } else {
        printf("Power of magic saves lives\n");
    }
}