矩阵快速运算_总结

矩阵乘法快速幂

struct matrix{

int mat[size][size];

matrix(){

int i,j;

for(i=0;i<d;i++)

for(j=0;j<d;j++)

mat[i][j]=0;

}

}; matrix mul(matrix a,matrix b){

matrix ret;

int i,j,k;

for(i=0;i<d;i++)

for(j=0;j<d;j++)

for(k=0;k<d;k++){

ret.mat[i][j]=(ret.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;

}

return ret;

}

matrix pow_mul(matrix a,int x){

matrix pw;

for(int i=0;i<d;i++)

pw.mat[i][i]=1;//注意初始化单位矩阵

while(x>0){

if(x&1)

pw=mul(pw,a);

a=mul(a,a);

x>>=1;

}

return pw;

}

#include<stdio.h>

#define N 2

#define mod 7

typedef struct Matrix{

int M[N][N];

}Matrix;

Matrix P;

Matrix I={1,0,0,1};

Matrix MatrixmulMatrix a,Matrix b){

Matrix ret;

int i,j,k;

for(i=0;i<N;i++)

for(j=0;j<N;j++){

ret.M[i][j]=0;

for(k=0;k<N;k++)

ret.M[i][j]+=(a.M[i][k]*b.M[k][j])%mod;

ret.M[i][j]%=mod;

}

return ret;

}

Matrix quiclpow(__int64 n)

{

Matrix m=p,b=I;

while(n>0){

if(n&1)

b=Matrixmul(b,m);

n>>=1;

m=Matrixmul(m,m);

}

return b;

}

strassen矩阵乘法算法

数学思想，分治法(二分查找，快速排序，大整数乘法strassen矩阵乘法算法这四种算法都运用了分治法的思想)和递归。

{ T(2)=b

  T(n)=7T(n/2)+an^2   n>2

#include<iostream>

using namespace std;

const int N=6;//define the size of the matrix

template<typename T>

void Strassen(int n, T a[][N],T b[][N],T c[][N]);

template<typename T>   

void Matrix_Multiply(T A[][N], T B[][N], T C[][N]) {  //Calculating A*B->C  

     for(int i=0; i<2; i++) {   

        for(int j=0; j<2; j++) {   

            C[i][j] = 0;         

           for(int t=0; t<2; t++) {   

               C[i][j] = C[i][j] + A[i][t]*B[t][j];           

            }     

         }           

      }   

}

template <typename T>   

void Matrix_Add(int n, T X[][N], T Y[][N], T Z[][N]) {   

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

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

            Z[i][j] = X[i][j] + Y[i][j];           

         }           

      }        

}  

template <typename T>   

void Matrix_Sub(int n, T X[][N], T Y[][N], T Z[][N]) {   

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

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

            Z[i][j] = X[i][j] - Y[i][j];           

         }           

      }        

}   

int main(){

int a[N][N],b[N][N],c[N][N];

for(int i=0; i<N; i++) {   

       for(int j=0; j<N; j++) {   

           a[i][j] = i * j;   

           b[i][j] = i * j;      

        }           

     }

Strassen(N,a,b,c);

return 0;

}

template<typename T>

void Strassen(int n,T a[][N],T b[][N],T c[][N]){

T a11[N][N],a12[N][N],a21[N][N],a22[N][N];

T b11[N][N],b12[N][N],b21[N][N],b22[N][N];//分块矩阵

T c11[N][N],c12[N][N],c21[N][N],c22[N][N];//

T m1[N][N],m2[N][N],m3[N][N],m4[N][N],m5[N][N],m6[N][N],m7[N][N];//关键矩阵

T aa[N][N],bb[N][N];

if(n==2)

Matrix_Multiply(a,b,c);

else{

for(int i=0;i<n/2;i++){

for(int j=0;j<n/2;j++){

  a11[i][j]=a[i][j]; 

               a12[i][j] = a[i][j+n/2];   

               a21[i][j] = a[i+n/2][j];   

               a22[i][j] = a[i+n/2][j+n/2];                    

               b11[i][j] = b[i][j];   

               b12[i][j] = b[i][j+n/2];   

               b21[i][j] = b[i+n/2][j];   

               b22[i][j] = b[i+n/2][j+n/2];       

            }           

         }              

        //Calculate M1 = (A0 + A3) × (B0 + B3)  

         Matrix_Add(n/2, a11, a22, aa);   

         Matrix_Add(n/2, b11, b22, bb);   

         Strassen(n/2, aa, bb, m1);           

        //Calculate M2 = (A2 + A3) × B0  

         Matrix_Add(n/2, a21, a22, aa);   

         Strassen(n/2, a, b11, m2);   

        //Calculate M3 = A0 × (B1 - B3)  

         Matrix_Sub(n/2, b12, b22, bb);   

         Strassen(n/2, a11, bb, m3);   

        //Calculate M4 = A3 × (B2 - B0)  

         Matrix_Sub(n/2, b21, b11, bb);   

         Strassen(n/2, a22, bb, m4);   

        //Calculate M5 = (A0 + A1) × B3  

         Matrix_Add(n/2, a11, a12, aa);   

         Strassen(n/2, aa, b22, m5);   

        //Calculate M6 = (A2 - A0) × (B0 + B1)  

         Matrix_Sub(n/2, a21, a11, a);   

         Matrix_Add(n/2, b11, b12, bb);   

         Strassen(n/2, aa, bb, m6);            

        //Calculate M7 = (A1 - A3) × (B2 + B3)  

         Matrix_Sub(n/2, a12, a22, aa);   

         Matrix_Add(n/2, b21, b22, bb);   

         Strassen(n/2, aa, bb, m7);   

        //Calculate C0 = M1 + M4 - M5 + M7  

         Matrix_Add(n/2, m1, m4, aa);   

         Matrix_Sub(n/2, m7, m5, bb);   

         Matrix_Add(n/2, aa, bb, c11);   

        //Calculate C1 = M3 + M5  

         Matrix_Add(n/2, m3, m5, c12);   

        //Calculate C2 = M2 + M4  

         Matrix_Add(n/2, m2, m4, c21);   

        //Calculate C3 = M1 - M2 + M3 + M6  

         Matrix_Sub(n/2, m1, m2, aa);   

         Matrix_Add(n/2, m3, m6, bb);   

         Matrix_Add(n/2, aa, bb, c22);   

        //Set the result to C[][N] 

        for(i=0; i<n/2; i++) {   

           for(int j=0; j<n/2; j++) {   

               c[i][j] = c11[i][j];   

               c[i][j+n/2] = c12[i][j];   

               c[i+n/2][j] = c21[i][j];   

               c[i+n/2][j+n/2] = c22[i][j];           

            }           

         }   

      }   

}  

C++ 类 矩阵存储及乘法

 struct Matrix{

int r, c;

int a[10000];

Matrix(int r,int c):r(r),c(c){

memset(a,0,sizeof(a));

}

const int* operator[](int i)const{return a+i*c}

int* operator[](int i){return a+i*c};

friend Matrix operator*(const Matrix& m1,Matrix& m2){

Matrix m(m1.r,m2.c);

for(int i=0;i<m1.r;i++)

for(int j=0;j<m1.c;j++)

for(int k=0;k<m2.c;k++)

m[i][k]+=m1[i][j]*m2[j][k];

return m;

}

friend Matrix operator^(const Matrix& m1,const int& n){

if(n==1)

return m1;

Matrix m=m1^(n>>1);

if(n&1)

return m*m*m1;

else

return m*m;

}

};