HDU 3059 Fibonacci數列與矩陣求和 矩陣大小不固定

http://acm.hdu.edu.cn/showproblem.php?pid=3509

snowingsea is having Buge’s discrete mathematics lesson, Buge is now talking about the Fibonacci Number. As a bright student, snowingsea, of course, takes it as a piece of cake. He feels boring and soon comes over drowsy.
Buge，feels unhappy about him, he knocked at snowingsea’s head, says:”Go to solve the problem on the blackboard!”, snowingsea suddenly wakes up, sees the blackboard written :



snowingsea thinks a moment，and writes down:



snowingsea has a glance at Buge，Buge smiles without talking, he just makes a little modification on the original problem, then it becomes :



The modified problem makes snowingsea nervous, and he doesn't know how to solve it. By the way，Buge is famous for failing students, if snowingsea cannot solve it properly， Buge is very likely to fail snowingsea. But snowingsea has many ACM friends. So，snowingsea is asking the brilliant ACMers for help. Can you help him?

#include <stdio.h>
#include <string.h>
#include <iostream>
using namespace std;
typedef long long LL;
const int MAX=54;
LL MOD,n_size;//n_size是變化矩陣的規模
struct Matrix
{
    long long m[MAX][MAX];
};
Matrix P;
Matrix I;
LL quick_mod(LL m,LL n,LL k)
{
    if(n==0)
        return 1;
    int b=1;
    while(n>0)
    {
        if(n&1)
            b=(b*m)%k;
        n=n>>1;
        m=(m*m)%k;
    }
    return b;
}
Matrix matrixmul(Matrix a,Matrix b)
{
    int i,j,k;
    Matrix c;
    for(i=0; i<n_size; i++)
        for(int j=0; j<n_size; j++)
        {
            c.m[i][j]=0;
            for(k=0; k<n_size; k++)
                c.m[i][j]+=((a.m[i][k]%MOD)*(b.m[k][j]%MOD))%MOD;
            c.m[i][j]%=MOD;
        }
    return c;
}
Matrix quickpow(Matrix m, LL n)
{
    Matrix b=I;
    while(n>=1)
    {
        if(n&1)
            b=matrixmul(b,m);
        n=n>>1;
        m=matrixmul(m,m);
    }
    return b;
}
LL c[50][50];
int main()
{
    Matrix tmp;
    LL sum=0,temp1,temp2;
    LL f1,f2,a,b,k,n,m;
    memset(c,0,sizeof(c));
    for(int i=0; i<=49; i++)
    {
        c[i][0]=1;
        c[i][i]=1;
    }
    for(int i=1; i<=49; i++)
        for(int j=1; j<i; j++)
            c[i][j]=c[i-1][j]+c[i-1][j-1];
    int T;
    scanf("%d",&T);
    while(T--)
    {
        sum=0;
        scanf("%I64d %I64d %I64d %I64d %I64d %I64d %I64d",&f1,&f2,&a,&b,&k,&n,&MOD);
        if(k==0)
            printf("%I64d\n",n%MOD);
        if(k>=1)
        {
            if(n==1)
            {
                printf("%I64d\n",quick_mod(f1,k,MOD));
                continue;
            }
            if(n==2)
            {
                printf("%I64d\n",(quick_mod(f1,k,MOD)+quick_mod(f2,k,MOD))%MOD);
                continue;
            }
            n_size=k+2;
            memset(P.m,0,sizeof(P.m));
            memset(I.m,0,sizeof(I.m));
            for(int i=0; i<n_size; i++)
                I.m[i][i]=1;
            P.m[0][0]=1;
            P.m[0][n_size-1]=1;
            for(int u=1; u<n_size-1; u++)
                P.m[0][u]=0;
            for(int j=1; j<n_size; j++)
                for(int k=n_size-j,w=0; k<n_size; k++,w++)
                    P.m[j][k]=(((c[j-1][w]%MOD)*quick_mod(a,w,MOD))%MOD*quick_mod(b,j-1-w,MOD))%MOD;
            tmp=quickpow(P,n-1);
            sum=(sum+(tmp.m[0][0]%MOD)*quick_mod(f1,k,MOD))%MOD;
            for(int i=1; i<n_size; i++)
            {
                temp1=(quick_mod(f1,n_size-1-i,MOD)*quick_mod(f2,i-1,MOD))%MOD;
                temp2=(temp1*tmp.m[0][i]%MOD)%MOD;
                sum=(sum+temp2)%MOD;
            }
            printf("%I64d\n",sum%MOD);
        }
    }
    return 0;
}