POJ 3233 Matrix Power Series（矩阵快速幂）

题面

Matrix Power Series
Time Limit: 3000MS Memory Limit: 131072K
Total Submissions: 19338 Accepted: 8161
Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’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

思路

这题暴力做肯定没救，我们需要优化。第一种做法是直接去递推矩阵关系，然后构造矩阵进行一个简化的快速幂，这样我们只需要操作一遍快速幂就好了。第二种做法的话，我们可以选择一个二分加速，求法类似于等比数列的二分求法，从数列中项开始提取公因式，不断递归就可以快速求得数列的和了，当然这里我们要把数列当成数字来处理，道理也是一样。所以，矩阵快速幂，包括一些数学题，当它的暴力做法很明显，而数据量又很大的时候，我们要尝试从数据量上去考虑如何加速。这题的话很容易想到二分，另外我们也需要对式子进行数学分析，看看是不是存在一些性质，可以让我们快速求解。

代码

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=65;
int n,mod,k;
struct Matrix {
    int maze[maxn][maxn];
    Matrix () {
        memset (maze,0,sizeof (maze));
    }
};
Matrix mul (Matrix a,Matrix b) {
    Matrix ans;
     for (int i=1;i<=n;i++) 
      for (int j=1;j<=n;j++) {
          for (int k=1;k<=n;k++) 
              (ans.maze[i][j]+=(a.maze[i][k]*b.maze[k][j])%mod)%=mod;
      }
      return ans;
}
Matrix fast_pow (Matrix a,int k) {
    Matrix ans;
    memset (ans.maze,0,sizeof (ans.maze));
    for (int i=1;i<=n;i++) ans.maze[i][i]=1;
    while (k) {
        if (k&1) ans=mul (ans,a);
        a=mul (a,a);
        k>>=1;
    }
    return ans;
}
int main () {
    while (cin>>n>>k>>mod) {
        Matrix a,b,c;
        for (int i=1;i<=n;i++)
         for (int j=1;j<=n;j++) {
             cin>>a.maze[i][j+n];
             b.maze[i+n][j+n]=a.maze[i][j+n];
         }  
        for (int i=1;i<=n;i++) {
           b.maze[i][i]=1;
           b.maze[i+n][i]=1;
        }
        n=n*2;
        c=mul (a,fast_pow (b,k));
        for (int i=1;i<=n/2;i++)
         for (int j=1;j<=n/2;j++) {
             if (j==n/2) cout<<c.maze[i][j]<<endl;
             else cout<<c.maze[i][j]<<" ";
         }
    }
    return 0;
}

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
const int MAXN = 35;
struct mat
{
    int a[MAXN][MAXN];
    mat()
    {
        memset(a,0,sizeof(a));
    }
};
mat m;
mat I;
int n,k,mod;
mat add(mat m1,mat m2)
{
    mat ans;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            ans.a[i][j]=(m1.a[i][j]+m2.a[i][j])%mod;
    return ans;
}
mat mul(mat m1,mat m2)
{
    mat ans;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
                for(int k=1;k<=n;k++)
                    ans.a[i][j]=(ans.a[i][j]+m1.a[i][k]*m2.a[k][j])%mod;
    return ans;
}
mat quickmul(mat m,int k)
{
    mat ans;
    for(int i=1;i<=n;i++) ans.a[i][i] = 1;
    while(k)
    {
        if(k&1) ans=mul(ans,m);
        m=mul(m,m);
        k>>=1;
    }
    return ans;
}
mat sum(int k)//等比求和
{
    if(k==1) return m;
    mat t=sum(k/2);
    mat ans;
    if(k&1)
    {
        mat cur=quickmul(m,k/2+1);
        ans=add(mul(add(I,cur),t),cur);
    }
    else
    {
        mat cur=quickmul(m,k/2);
        ans=mul(add(I,cur),t);
    }
    return ans;
}
void print(mat m)
{
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            printf("%d%c",m.a[i][j],j==n? '\n':' ');
}
int main()
{
    while(~scanf("%d%d%d",&n,&k,&mod))
    {
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
            {
                scanf("%d",&m.a[i][j]);
                I.a[i][j]=(i==j);
            }
        mat ans=sum(k);
        print(ans);
    }
    return 0;
}

原文地址：https://www.cnblogs.com/hhlya/p/13287240.html