快速幂与矩阵—>快速矩阵幂

来源：互联网发布：sas 优化板块编辑：程序博客网时间：2024/05/24 07:04

矩阵的幂其实就是多次相乘，但是用最普通的方法求往往是超时的这就用到了快速幂的思想

还是用了二分的思想：

下面Down了网友的模板（可以结合者后面的具体题目理解~）

//HOJ 3493/*===================================*/|| 快速幂（quickpow）模板 || P 为等比，I 为单位矩阵|| MAX 要初始化！！！！||/*===================================*//*****************************************************/#include <cstdio>const int MAX = 3;typedef  struct{        int  m[MAX][MAX];}  Matrix;Matrix P = {5,-7,4,            1,0,0,            0,1,0,           };Matrix I = {1,0,0,            0,1,0,            0,0,1,           };           Matrix matrixmul(Matrix a,Matrix b) //矩阵乘法{       int i,j,k;       Matrix c;       for (i = 0 ; i < MAX; i++)           for (j = 0; j < MAX;j++)             {                 c.m[i][j] = 0;                 for (k = 0; k < MAX; k++)                     c.m[i][j] += (a.m[i][k] * b.m[k][j])%9997;                 c.m[i][j] %= 9997;             }       return c;}          Matrix quickpow(long long n){       Matrix m = P, b = I;       while (n >= 1)       {             if (n & 1)                b = matrixmul(b,m);             n = n >> 1;             m = matrixmul(m,m);       }       return b;}               /*************************************/int main(){    Matrix re;    int f[3] = {2,6,19};    long long n;    while (scanf("%I64d",&n) && n != 0)    {          if (n == 1)             printf("1\n");          else if (n <= 4)                  printf("%d\n",f[n-2]);               else {                      re = quickpow(n - 4);                      printf("%d\n",(((re.m[0][0]*f[2])                              + (re.m[0][1]*f[1]) + (re.m[0][2]*f[0])) %9997 + 9997) % 9997);                      }    }    return 0;}

下面举一个用快速矩阵幂求解的题目：

Matrix Power Series

题目描述：

Description

Given a n × n matrix A and a positive integer k, find the sumS =A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integersn (n ≤ 30),k (k ≤ 10⁹) and m (m < 10⁴). Then follown lines each containingn nonnegative integers below 32,768, givingA’s elements in row-major order.

Output

Output the elements of S modulo m in the same way as A is given.

Sample Input

2 2 40 11 1

Sample Output

1 22 3

矩阵快速幂。首先我们知道 A^x 可以用矩阵快速幂求出来（具体可见poj 3070）。其次可以对k进行二分，每次将规模减半，分k为奇偶两种情况，如当k = 6和k = 7时有：

k = 6 有： S(6) = (1 + A^3) * (A + A^2 + A^3) = (1 + A^3) * S(3)。

k = 7 有： S(7) = A + (A + A^4) * (A + A^2 + A^3) = A + (A + A^4) * S(3)。

所以有：

奇数: F[n]=F[n-1]+A^n

偶数: F[n]=F[n/2]+F[n/2]*An/2

（是不是和求数的快速幂很像？？？）

ps：对矩阵定义成结构体Matrix，求S时用递归，程序会比较直观，好写一点。当然定义成数组，然后再进行一些预处理，效率会更高些。

#include<iostream>#include <cstdio>using namespace std;const int MAX = 32;//int count=0;struct Matrix   //用于相同的类型的数组的开辟，比用单独的名称好{    int v[MAX][MAX];};int n, k, M;Matrix mtAdd(Matrix A, Matrix B)        // 求矩阵 A + B{int i, j;Matrix C;for(i=0; i<n; i++)        for(j=0; j<n; j++)C.v[i][j]=(A.v[i][j] + B.v[i][j]) % M;//count++;return C;}Matrix mtMul(Matrix A, Matrix B)        // 求矩阵 A*B{int i, j, k;Matrix C;for(i = 0; i < n; i ++)for(j = 0; j < n; j ++){            C.v[i][j] = 0;            for(k = 0; k < n; k ++)C.v[i][j] = (A.v[i][k] * B.v[k][j] + C.v[i][j]) % M;//count++;}return C;}Matrix mtPow(Matrix A, int k)           // 求矩阵 A^k{if(k == 0) {memset(A.v, 0, sizeof(A.v));for(int i = 0; i < n; i ++)A.v[i][i]=1;//count++;return A;}//count++;if(k == 1) return A;Matrix C = mtPow(A, k / 2);if(k % 2 == 0)return mtMul(C, C);elsereturn mtMul(mtMul(C, C), A);}Matrix mtCal(Matrix A, int k)           // 求S (k) = A + A2 + A3 + … + Ak{   //count++;if(k == 1) return A;Matrix B = mtPow(A, (k+1) / 2);Matrix C = mtCal(A, k / 2);if(k % 2 == 0)return mtMul( mtAdd(mtPow(A, 0), B), C );   // 如S(6) = ( 1 （单位矩阵） + A^3 ) * S(3)。 else return mtAdd( A, mtMul(mtAdd(A, B), C) );   // 如S(7) = A + (A + A^4) * S(3)}int main(){int i, j;Matrix A;//输入的原矩阵cin>>n>>k>>M;for(i=0; i<n; i++)for(j=0;j<n;j++)scanf("%d",&A.v[i][j]);A=mtCal(A, k);for(i=0; i<n; i++){for(j=0; j<n; j++){printf("%d ",A.v[i][j]);}printf("\n");}  //cout<<count;return 0;}

还有一个，估计其实质是一样的但是现在还不是很懂，以后弄懂~~~

因为S可以看成S=A(I+A(I+A(I+...A(I+A)))) (I是单位矩阵)

拿k=3举例S=A(I+A(I+A))

那么我们想，可不可以构造一个矩阵T使得T*T（因为是k次幂）这样乘下去每次可以得到A*(A+I)

那么肯定T有个两个元素就是A与I

那么假设：T=A I
                     I I
那么T=T*T=A*A+I*I       A*I+I*I
                 A*I+I*I           I*I+I*I
这样存在一个I*（A+I）的式子，当T再乘以T的时候会出现A（A+I）

这个时候我们可以简化将T=A I

0 I

这样可以简化很多计算T*T=A*A A*I+I*I

0 I

那么容易得到T^(K+1)={A^(K+1) I+A+A^2+A^3+...+A^K}
0 I

这样我们只需要算T的k+1次幂就可以了

而k如此庞大所以需要二分来对T求k+1次幂

对T求快速幂：
首先我们知道A^19=A^16*A^2*A^1，因为19=B(10011)

那么就这样，拿A总是去和本身去乘，那么就可以取到 A A^2 A^4 A^8 A^16...

然后问A A^2 A^4 A^8 A^16 ...这么多项取与不取的抉择
其实就是一个求二进制的这样一个过程
19%2=1 那么A取进来
19/2=9
9%2=1 那么A^2取进来
9/2=4
4%2=0 那么A^4不用取进来
4/2=2
2%2=0 那么A^8不用去进来
2/2=1
1%2=1 那么A^16取进来
1/2=0 计算完毕

这个就是用模板做出来的~~~

抄了一点别人的代码~~~呵呵~~~

#include <iostream>#include <string.h>using namespace std;#define MAXV 66struct{int r,c;//c行r列int mat[MAXV][MAXV];///矩阵中再放入矩阵}Matrix;Matrix ans,cnt;int n,k,m;void Input(){int i,j;/*构造B矩阵*/memset(cnt.mat,0,sizeof(cnt.mat));memset(ans.mat,0,sizeof(ans.mat));for(i=0;i<n;i++){for(j=0;j<n;j++)scanf("%d",&cnt.mat[i][j]);}for(i=0;i<n;i++){cnt.mat[i+n][i+n]=cnt.mat[i][i+n]=1;ans.mat[i][i]=ans.mat[i+n][i+n]=1;}/*对矩阵B和B^（k+1）初始化*/cnt.c=cnt.r=2*n;ans.c=ans.r=2*n;}Matrix MatrixMul(Matrix x,Matrix y){//矩阵乘法Matrix t;int i,j,v;memset(t.mat,0,sizeof(t.mat));t.c=x.c;t.r=y.r;for(i=0;i<t.c;i++)for(j=0;j<t.r;j++){for(v=0;v<x.r;v++)t.mat[i][j]=t.mat[i][j]+((x.mat[i][v]*y.mat[v][j])%m);t.mat[i][j]=t.mat[i][j]%m;}return t;}void Binary(){//二分快速幂k++;while(k){if(k & 1)ans=MatrixMul(ans,cnt);cnt=MatrixMul(cnt,cnt);k=k>>1;}}void Output(){/*输出的时候要减去一个单位矩阵*/int i,j;for(i=0;i<n;i++){for(j=0;j<n;j++)if(i!=j)printf("%d ",ans.mat[i][j+n]);elseprintf("%d ",ans.mat[i][j+n]-1);printf("\n");}}int main(){while(~scanf("%d%d%d",&n,&k,&m)){Input();Binary();Output();}return 0;}

另外，本体还可以用类似于数列的前N项和来求解，这样的话用一次对幂求解一套公式就出来了，别说现在前N项和记不起来了~~~

A(1-A^q)/(1-q) (这里的1应该是单位向量了)

0 0