[POJ3233] Matrix Power Series && 线性递推关系矩阵加速解法

易知：

[A(n-1), A] * [ A, 0    = [A(n), A]

                      1, 1] 

可以用快速幂加速 

（如果不是练习矩阵加速有更好的做法 这里不给出）

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
using namespace std;
#define MAXN 30
#define MAXM 30
typedef int LL;
struct Maxtrix{  
    int n, m;  
    LL A[MAXN+10][MAXN+10];  
};
LL a[MAXN+10], f[MAXN+10];
int n, k;
LL p;
Maxtrix A, B, E;
Maxtrix ST[1][2], M[2][2], C[2][2];
void init(Maxtrix &X)
{
    X.n = X.m = 0;
    memset(X.A, 0, sizeof(X.A));
}
void Get_E()
{
    E.n = E.m = MAXN;
    for(int i = 0; i < E.n; i++) 
        E.A[i][i] = 1;
}
Maxtrix add(Maxtrix AA, Maxtrix BB)
{
    Maxtrix D; D = AA;
    for(int i = 0; i < AA.n; i++)
        for(int j = 0; j < AA.m; j++)
		{
			D.A[i][j] += BB.A[i][j];
			D.A[i][j] %= p;
		}
    return D;
}
Maxtrix mul(Maxtrix AA, Maxtrix BB)
{
    Maxtrix D;
    init(D);
    D.n = AA.n;
    D.m = BB.m;
    for(int i = 0; i < AA.n; i++)
        for(int j = 0; j < BB.n; j++)
            for(int k = 0; k < BB.m; k++)
            {
                D.A[i][k] += (AA.A[i][j] * BB.A[j][k]) % p;
                D.A[i][k] %= p;
            }
    return D;
}
Maxtrix pow_mod(Maxtrix &A, int k)
{
    Maxtrix ans = E, t = A;
    ans.n = A.n;
    ans.m = A.m;
    while(k)
    {
        if(k & 1)
            ans = mul(ans, t);
        t = mul(t, t);
        k >>= 1;
    }
    return ans;
}
void dou_mul(Maxtrix A[][2], Maxtrix BB[][2], Maxtrix C[][2])
{
    Maxtrix D[2][2];
    for(int i = 0; i < 2; i++) 
        for(int j = 0; j < 2; j++) 
            D[i][j] = B;
    for(int i = 0; i < 2; i++)
        for(int j = 0; j < 2; j++)
            for(int k = 0; k < 2; k++)
            {
                Maxtrix tmp;
                tmp = mul(A[i][j], BB[j][k]);
                D[i][k] = add(D[i][k], tmp);
            }
    memcpy(C, D, sizeof(D));
}
void dou_pow_mod(int k, Maxtrix Ans[][2])
{
    Maxtrix ans[2][2], t[2][2];
    for(int i = 0; i < 2; i++)  { ans[i][i] = E; ans[i][i ^ 1] = B; }
    memcpy(t, M, sizeof(M));
    while(k)
    {
        if(k & 1) dou_mul(ans, t, ans);
        dou_mul(t, t, t);
        k >>= 1;
    }
    memcpy(Ans, ans, sizeof(ans));
}
int main()
{
    Get_E();
    while(cin >> n >> k >> p && n+k+p)
    {
        init(A);
        init(B); B.n = B.m = n;
        E.n = E.m = A.n = A.m = n;
        for(int i = 0; i < n; i++) 
            for(int j = 0; j < n; j++)
                cin >> A.A[i][j];
        ST[0][0] = ST[0][1] = A;
        M[0][0] = A; M[0][1] = B; M[1][0] = E; M[1][1] = E;
        Maxtrix tmp[2][2];
        dou_pow_mod(k-1, tmp);
        dou_mul(ST, tmp, tmp);
        Maxtrix ANS = tmp[0][0];
        for(int i = 0; i < n; i++)
        {
            printf("%d", ANS.A[i][0]);
            for(int j = 1; j < n; j++)
                printf(" %d", ANS.A[i][j]);
            puts(" ");
        }
    }
}