zoukankan      html  css  js  c++  java
  • 快速幂与矩阵—>快速矩阵幂

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

    还是用了二分的思想:

    下面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
    ");
              else if (n <= 4)
                      printf("%d
    ",f[n-2]);
                   else {
                          re = quickpow(n - 4);
                          printf("%d
    ",(((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 + A2 + A3 + … + Ak.

    Input

    The input contains exactly one test case. The first line of input contains three positive integersn (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follown lines each containing n 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 4
    0 1
    1 1

    Sample Output

    1 2
    2 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);
    				else
    					return 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("
    ");
    		}
    		
    						  
    		//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 66
    
    struct{
    	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]);
    			else
    				printf("%d ",ans.mat[i][j+n]-1);
    			printf("
    ");
    	}
    }
    
    int main(){
    	while(~scanf("%d%d%d",&n,&k,&m)){
    		Input();
    		Binary();
    		Output();
    	}
    	return 0;
    }


     

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

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

  • 相关阅读:
    HTC G12 HBOOT 2.0.0002 官解、降级精简教程!
    android call require api level
    刷机流程
    Android通过tcpdump抓包
    [原创]c#的线性表 Virus
    [原创]c++线性表的插入 Virus
    [原创]关于编程论坛封我ID的意见 Virus
    [原创]文件管理 Virus
    [原创]c#高级编程学习笔记(连载_委托) Virus
    [原创]c#中的内存管理 Virus
  • 原文地址:https://www.cnblogs.com/zswbky/p/5432058.html
Copyright © 2011-2022 走看看