poj-3233-Matrix Power Series

矩阵快速幂取模

点此查看原题链接

Description

Given a n × n matrix A and a positive integer k, find the sum $S = A + A^2 + A^3 + … + A^k.$

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

分析

一上来就使用单纯的快速幂，先求$A^k$，然后再 累加求和。但是k是$10^9$，严重超时。其实只要稍微研究一下还是可以发现规律的。
$s=A+A^2+…+A^k$
$s=A(1+A(1+A…（1+A))$
但是这样还是不够的，我们来依次列开看看：
$s1=A$
$s2=A+A^2$
$s3=A+A^2+A^3$
发现规律了吗？
$s2=A*(1+s1)$
$s3=A*(1+s2)$
是否发现了什么？这和快速幂求斐波那契数列是不是很相似？
于是我们可以构造矩阵

$$T= \begin{vmatrix} A & 1 \\ 0& 1\\ \end{vmatrix}$$
接下来看看。
$$T*T= \begin{vmatrix} A^2 & A +1\\ 0& 1\\ \end{vmatrix}$$
$$T*T*T= \begin{vmatrix} A ^3& A^2+A +1 \\ 0& 1\\ \end{vmatrix}$$
我们就可以求 $s$了，就是矩阵$T^{k+1}$的右上角，但是需要减去单位矩阵。
所以现在需要的就是构造矩阵 了。
构建数组cnt[2*n][2*n]，左上角用来存储矩阵A，右上角和右下角在对角线上预处理为零，就构造成功了。

最后就是贴上杂乱的代码了

#include <cstdio>
#include <iostream>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <string>
#include <queue>
#include <vector>
#include <stack>
#include <algorithm>
#include <set>
#include <map>
using namespace std;
const int N=61;
int ans[N][N];
int a[N][N];
int tmp[N][N];
int cnt[N][N];
int n,m,k;

void fast_mod(int x)
{
    memset(ans,0,sizeof(ans));
    for (int i=0;i<2*n;i++) //构建单位矩阵
        ans[i][i]=1;
    while(x)
    {
        if(x&1)
        {
            for (int i=0;i<2*n;i++)
                for (int j=0;j<2*n;j++)
                    tmp[i][j]=ans[i][j];
            memset(ans,0,sizeof(ans));
            for (int i=0;i<2*n;i++)
                for (int j=0;j<2*n;j++)
                    for (int kk=0;kk<2*n;kk++)
                        ans[i][j]=(ans[i][j]+(tmp[kk][j]*cnt[i][kk])%m)%m;
        }
        for (int i=0;i<2*n;i++)
            for (int j=0;j<2*n;j++)
                tmp[i][j]=cnt[i][j];
        memset(cnt,0,sizeof(cnt));
        for(int i=0;i<2*n;i++)
            for (int j=0;j<2*n;j++)
                for (int kk=0;kk<2*n;kk++)
                    cnt[i][j]=(cnt[i][j]+tmp[i][kk]*tmp[kk][j]%m)%m;
        x/=2;
    }
    return ;
}

int main ()
{
   while ( scanf ("%d%d%d",&n,&k,&m)!=EOF)
   {
        memset(a,0,sizeof(a));
        memset(cnt,0,sizeof(cnt));

        for (int i=0;i<n;i++)
        {
            for (int j=0;j<n;j++)
            {
                scanf ("%d",&a[i][j]);  //读入矩阵A
                cnt[i][j]=a[i][j];  //左上角填入矩阵A
            }
            cnt[i][i+n]=1;  //右上角对角线
            cnt[i+n][i+n]=1;    //右下角对角线
        }
        fast_mod(k+1);  //快速幂求解

        for (int i=0;i<n;i++)   //输出矩阵的右上角
        {
            for (int j=0;j<n-1;j++)
            {
                if(i!=j)
                    printf ("%d ",(ans[i][j+n]+m)%m);   //这里一定要先+m在%m，防止小于1
                else
                    printf ("%d ",(ans[i][j+n]-1+m)%m); //这里更加需要+m，减去单位矩阵
            }
            if(i!=n-1)
                printf ("%d\n",(ans[i][2*n-1]+m)%m);
            else
                printf ("%d\n",(ans[i][2*n-1]-1+m)%m);
        }
   }
   return 0;
}

再次提醒
最后输出时一定要+m再%m，我是绝不会告诉你，我就是这样WA了好多次。