【POJ】3233 Matrix Power Series

【算法】二分+矩阵快速幂

【题意】给定矩阵A和整数k,MOD，求A^0+A^1+A^2+...+A^k。

【题解】

定义题目要求的答案为f(n)，即：

$$f_n=sum_{i=0}^{n}A^i$$

当n为偶数时，可以拆成两半，后一半由前一半集体乘A(n/2)得到，即：

$$f_n=f_{frac{n}{2}}(A^{frac{n}{2}}+1)$$

当n为奇数时，直接递推：

$$f_n=f_{n-1}*A^n$$

复杂度O(n^3 log k)。

快速幂的单位矩阵是主对角线(左上到右下）全为1，其余全为0，不用memset就超时了，多用stdio.h。

#include<cstdio>
#include<algorithm>
#include<cstring>
#define ll long long
using namespace std;
const int maxn=40;
int n,MOD,kind;
struct Mat{ll a[maxn][maxn];}A;

Mat ch(Mat a,Mat b){
    Mat tmp;
    memset(tmp.a,0,sizeof(tmp.a));
    for(int k=1;k<=n;k++)
        for(int i=1;i<=n;i++)
            for(int j=1;j<=n;j++)
                tmp.a[i][j]=(tmp.a[i][j]+a.a[i][k]*b.a[k][j])%MOD;
    return tmp;
}
Mat pow(int k){
    Mat tmp=A,ans;
    memset(ans.a,0,sizeof(ans.a));
    for(int i=1;i<=n;i++)ans.a[i][i]=1;
    //快速幂初值为1（单位矩阵）！！！ 
    while(k>0){
        if(k&1)ans=ch(ans,tmp);
        tmp=ch(tmp,tmp);
        k>>=1;
    }
    return ans;
}
Mat plus(Mat a,Mat b){
    Mat tmp;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)tmp.a[i][j]=(a.a[i][j]+b.a[i][j])%MOD;
    return tmp;
}
Mat calc(int k){
    Mat tmp;
    if(k<=1)return A;
    if(k&1){
        tmp=plus(calc(k-1),pow(k));
    }
    else{
        Mat tmps=calc(k/2);
        tmp=plus(tmps,ch(tmps,pow(k/2)));
    }
    return tmp;
}
int main(){
    scanf("%d%d%d",&n,&kind,&MOD);
    for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)scanf("%lld",&A.a[i][j]);
    Mat ans=calc(kind);
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++)printf("%lld ",ans.a[i][j]%MOD);
        printf("
");
    }
    return 0;
}

还有一道题：HDU1588 Gauss Fibonacci

给定k,b,n,m，求：

$$ans=sum_{i=0}^{n-1}Fib(k*i+b) mod m$$

定义A^i表示Fib(i)的斐波那契矩阵（见Fibonacci，左下角项），那么：

$$sum=A^b imes sum_{i=0}^{n-1}(A^k)^i  mod m$$

后面将$A^k$视为整体后，就是本题的套路了。