HDU 3306 Another kind of Fibonacci

题意：A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2) (N >= 2)；给定三个值N，X，Y求S(N):S(N) = A(0)2 +A(1)2+……+A(n)2。

思路：原来我们讲的斐波那契数列是： F(0) = 1, F(1) = 1, F(N) = F(N - 1) + F(N - 2)

这道题规定了另一种斐波那契数列形式：A(0) = 1 , A(1) = 1 , A(N) = X * A(N - 1) + Y * A(N - 2)

其实原理是一样的！
我们以前快速求Fibonacci数列f[n]=f[n-1]+f[n-2],f[1]=f[2]=1第n项的方法是！！？？！！
构造常系数矩阵！
（一）Fibonacci数列f[n]=f[n-1]+f[n-2],f[1]=f[2]=1的第n项的快速求法（不考虑高精度）.
解法：
考虑1×2的矩阵【f[n-2],f[n-1]】。根据fibonacci数列的递推关系，我们希望通过乘以一个2×2的矩阵，得到矩阵【f[n-1],f[n]】=【f[n-1],f[n-1]+f[n-2]】
很容易构造出这个2×2矩阵A，即：
０ １
１ １
所以，有【f[1],f[2]】×A＝【f[2],f[3]】
又因为矩阵乘法满足结合律，故有：
【f[1],f[2]】×A n-1=【f[n],f[n+1]】
这个矩阵的第一个元素即为所求。
至于如何快速求出A n-1，相信大家都会，即递归地：n为偶数时，An=(A n/2)2；n为奇数时，An=(A n/2)2*A。
问题（一）解决。
（二）数列f[n]=f[n-1]+f[n-2]+1,f[1]=f[2]=1的第n项的快速求法（不考虑高精度）.
解法：
仿照前例，考虑1×3的矩阵【f[n-2],f[n-1],1】，希望求得某3×3的矩阵A，使得此1×3的矩阵乘以A得到矩阵：【f[n-1],f[n],1】＝【f[n-1],f[n-1]+f[n-2]+1,1】
容易构造出这个3×3的矩阵A，即：
０ １ ０
１ １ ０
０ １ １
问题（二）解决。
（三）数列f[n]=f[n-1]+f[n-2]+n+1,f[1]=f[2]=1的第n项的快速求法（不考虑高精度）.
解法：
仿照前例，考虑1×4的矩阵【f[n-2],f[n-1],n,1】，希望求得某4×4的矩阵A，使得此1×4的矩阵乘以A得到矩阵：
【f[n-1],f[n],n+1,1】＝【f[n-1],f[n-1]+f[n-2]+n+1,n+1,1】
容易构造出这个4×4的矩阵A，即：
０ １ ０ ０
１ １ ０ ０
０ １ １ ０
０ １ １ １
问题（三）解决……
（四）数列f[n]=f[n-1]+f[n-2],f[1]=f[2]=1的前n项和s[n]的快速求法（不考虑高精度）.
解法：
虽然我们有S[n]=F[n+2]-1，但本文不考虑此方法，我们想要得到更一般的方法。
考虑（一）的矩阵A，容易发现我们要求【f[1],f[2]】×（A+A2+A3+…+AN-1）。很多人使用一种很数学的方法构造一个2r*2r（r是A的阶数，这里为2）的矩阵来计算，这种方法比较麻烦且很慢，这里不再介绍。下面考虑一种新方法。
仿照之前的思路，考虑1×3的矩阵【f[n-2],f[n-1],s[n-2]】，我们希望通过乘以一个3×3的矩阵A，得到1×3的矩阵：
【f[n-1],f[n],s[n-1]】=【f[n-1],f[n-1]+f[n-2],s[n-2]+f[n-1]】
容易得到这个3×3的矩阵是：
０ １ ０
１ １ １
０ ０ １
然后…………容易发现，这种方法的矩阵规模是(r+1)*(r+1)，比之前流行的方法好得多。
（五）数列f[n]=f[n-1]+f[n-2]+n+1,f[1]=f[2]=1的前n项和s[n]的快速求法（不考虑高精度）.
解法：
结合（三）（四），容易想到……
考虑1×5的矩阵【f[n-2],f[n-1],s[n-2],n,1】,
我们需要找到一个5×5的矩阵A，使得它乘以A得到如下1×5的矩阵：
【f[n-1],f[n],s[n-1],n+1,1】
=【f[n-1], f[n-1]+f[n-2]+n+1,s[n-2]+f[n-1],n+1,1】
容易构造出A为：
０ １ ０ ０ ０
１ １ １ ０ ０
０ ０ １ ０ ０
０ １ ０ １ ０
０ １ ０ １ １
然后……问题解决。
一般地，如果有f[n]=p*f[n-1]+q*f[n-2]+r*n+s
可以构造矩阵A为：
０ q 0 0 0
1 p 1 0 0
0 0 1 0 0
0 r 0 1 0
0 s 0 1 1
更一般的，对于f[n]=Sigma(a[n-i]*f[n-i])+Poly(n)，其中0<i<=某常数c, Poly (n)表示n的多项式，我们依然可以构造类似的矩阵A来解决问题。
设Degree(Poly(n))=d, 并规定Poly(n)=0时，d=-1，此时对应于常系数线性齐次递推关系。则本方法求前n项和的复杂度为：
((c+1)+(d+1))3*logns

所以对于这道题，

我们可以根据(s[n-2], a[n-1]^2, a[n-1]*a[n-2], a[n-2]^2) * A = (s[n-1], a[n]^2, a[n]*a[n-1], a[n-1]^2)

能够求出关系矩阵

  　　  |  1  0  0  0    | 
A =    | 1 x^2 1  x   |
  　　  | 0  y^2  0 0  |
         | 0  2*x*y 0 y |

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=3306

 

#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <algorithm>
#include <iostream>
using namespace std;
const int N=5;
const int mod=10007;

typedef struct In{
    __int64 m[N][N];
}Matrix;
Matrix init,unit,S;
__int64 n,x,y;

void Init(){
    for(int i=1;i<N;i++){
        for(int j=1;j<N;j++){
            if(i==1&&j==1) init.m[i][j]=1;
            else if(i==2&&j==1) init.m[i][j]=1;
            else if(i==2&&j==2) init.m[i][j]=x*x%mod;
            else if(i==2&&j==3) init.m[i][j]=1;
            else if(i==2&&j==4) init.m[i][j]=x%mod;
            else if(i==3&&j==2)    init.m[i][j]=y*y%mod;
            else if(i==4&&j==2) init.m[i][j]=2*x*y%mod;
            else if(i==4&&j==4) init.m[i][j]=y%mod;
            else init.m[i][j]=0;
            unit.m[i][j]=(i==j);
            if(i==1) S.m[i][j]=1;
            else S.m[i][j]=0;
        }
    }
}

Matrix Mul(Matrix a,Matrix b){
    Matrix c;
    for(int i=1;i<N;i++)
        for(int j=1;j<N;j++){
            c.m[i][j]=0;
            for(int k=1;k<N;k++){
                c.m[i][j]+=a.m[i][k]*b.m[k][j];
                c.m[i][j]%=mod;
            }
        }
    return c;
}

Matrix Pow(Matrix a,Matrix b){
    while(n){
        if(n&1) b=Mul(a,b);
        a=Mul(a,a);
        n>>=1;
    }
    return b;
}

void Debug(Matrix a){
    for(int i=1;i<N;i++){
        for(int j=1;j<N;j++)
            printf("%d ",a.m[i][j]);
        puts("");
    }
}

int main(){
    
//    freopen("data.in","r",stdin);
//    freopen("data.out","w",stdout);
    
    while(scanf("%I64d%I64d%I64d",&n,&x,&y)!=EOF){
        Init();
        Matrix x=Pow(init,unit);
        x=Mul(S,x);
    //    Debug(x);
        printf("%I64d\n",x.m[1][1]%mod);
    }
    return 0;
}