[HDU 3306]矩阵乘法、神一般的构造**

【题目大意】

f(0)=1,f(1)=1

给出递推式f(n)=x*f(n-1)+y*f(n-2)

令s(n)=f(0)^2+f(1)^2+……+f(n)^2

【算法分析】矩阵乘法，但是感觉这简直是神来之笔。。。

构造矩阵[S(n)   A(n-1)^2     A(n)*A(n-1)        A(n)^2]

转移的话乘这个矩阵

|   1           0             0             0            |

|   Y^2        0             0            Y^2         |

|   2*X*Y     0             Y            2*X*Y     |

|    X^2       1              X            X^2         |

【其它】WA了N遍，发现init矩阵的时候也要取余，因为有平方，X,Y<=2^31的话，乘起来就到int64了，然后矩阵乘法那里也可能爆，所以必须先取余

2073791 2010-02-07 14:46:22 Accepted 3306 687MS 220K 1660 B G++ cwj

【CODE】

#include #include const __int64 mod=10007;
__int64 c[4][4],a[4][4],t[4][4],p[4][4],n,x,y;

void init(){
    a[0][0]=2;
    a[0][1]=1;
    a[0][2]=1;
    a[0][3]=1;
    c[0][0]=1;    c[0][1]=0;     c[0][2]=0;      c[0][3]=0;
    c[1][0]=y*y; c[1][1]=0;     c[1][2]=0;      c[1][3]=y*y;
    c[2][0]=2*x*y;c[2][1]=0;     c[2][2]=y;      c[2][3]=2*x*y;
    c[3][0]=x*x; c[3][1]=1;     c[3][2]=x;      c[3][3]=x*x;
    for (int i=0;i<4;i++)
      for (int j=0;j<4;j++)
        c[i][j]%=mod;
}

void make(__int64 k){
    if (k==1){
        memcpy(p,c,sizeof(p));
        return;
    }   
    make(k/2);
    memset(t,0,sizeof(t));
    for (__int64 i=0;i<4;i++)
      for (__int64 j=0;j<4;j++)
        for (__int64 k=0;k<4;k++)
          t[i][j]=(t[i][j]+p[i][k]*p[k][j])%mod;
    memcpy(p,t,sizeof(t));
    if (k&1){
        memset(t,0,sizeof(t));
        for (__int64 i=0;i<4;i++)
          for (__int64 j=0;j<4;j++)
            for (__int64 k=0;k<4;k++)
              t[i][j]=(t[i][j]+p[i][k]*c[k][j])%mod;
        memcpy(p,t,sizeof(t));
    }   
}   

int main(){
    while (scanf("%I64d%I64d%I64d",&n,&x,&y)!=EOF){
        if (n==1){
            printf("2n");
            continue;
        }   
        if (n==0){
            printf("1n");
            continue;
        }   
        n–;
        init();       
        make(n);
        memset(t,0,sizeof(t));
        for (__int64 i=0;i<1;i++)
          for (__int64 j=0;j<4;j++)
            for (__int64 k=0;k<4;k++)
              t[i][j]=(t[i][j]+a[i][k]*p[k][j])%mod;
        printf("%I64dn",t[0][0]);
    }   
}