HIT2255 Not Fibonacci 变形斐波那契前n项和

题目链接：http://acm.hit.edu.cn/hoj/problem/view?id=2255

题目大意：f[0]=a,f[1]=b,f[n]=p×f[n-1]+q×f[n-2]，计算∑f[i]mod 10^7（i from s to e）.

分析：对于斐波那契数列的前n项和，我们可以参考HIT2060.对于这题，我们知道：

|0 1|  |f(0)|
|q p|  |f(1)|
这个可以利用矩阵幂次求得任何一个F(n),可惜这里是求和,
那么令A=前面那个p,q的矩阵
则S(n) = f(0) + f(1) + ... + f(n) = f(0) + (A + A^2 + ... + A^n)*F所得矩阵的右上角的值
那么如何求得A + A^2 + ... + A^n呢
可以继续构造如下的分块矩阵，其中I是单位矩阵
|A I|
|0 I|
令R等于上面的矩阵，则
R^2 = |A^2   A*I + I|
      |0     I      |
R^3 = |A^3   A^2 * I + A * I + I|
      | 0            I          |
可以发现右上角即为I + A + A^2 + ... + A^n，多个I后面给减掉就可以了
这样我们同样可以再次利用矩阵幂次求得R^n
则S(n) = ((R^(n+1).m12 - I) * F).a12 + f(0)
可以很完美的解决这个问题.

实现代码如下：

#include <cstdio>
#include <iostream>
using namespace std;
const int mod=10000000;
typedef struct
{
    long long m[3][3];
}matrix;
matrix I={1,0,0,
          0,1,0,
          0,0,1};
matrix mul(matrix a,matrix b)
{
    int i,j,k;
    matrix c;
    for(i=0;i<3;i++)
      for(j=0;j<3;j++)
      {
          c.m[i][j]=0;
          for(k=0;k<3;k++)
            c.m[i][j]+=((a.m[i][k]*b.m[k][j])%mod+mod)%mod;
          c.m[i][j]=(c.m[i][j]%mod+mod)%mod;
      }
    return c;
}
matrix quick_mod(long long n,int p,int q)
{
    matrix P={1,p,q,
          0,p,q,
          0,1,0};
    matrix a=P,b=I;
    while(n>0)
    {
        if(n&1)
          b=mul(b,a);
        n=n>>1;
        a=mul(a,a);
    }
    return b;
}
int main()
{
    int a,b,p,q,a1,b1;
    int t;
    long long ans;
    while(scanf("%d",&t)!=-1)
    while(t--)
    {
        scanf("%d%d%d%d%d%d",&a1,&b1,&p,&q,&a,&b);
        if(a==0&&b==1) ans=a1+b1;
        else if(a==0&&b==0)  ans=a1;
        else if(a==0&&b>=2)
        {
            matrix temp2=quick_mod(b-1,p,q);
            ans=((temp2.m[0][0]*(a1+b1))%mod+(temp2.m[0][1]*b1)%mod+(temp2.m[0][2]*a1)%mod)%mod;
        }
        else if(a==1&&b==1)  ans=b1;
        else
        {
            long long head,tail;
            matrix temp1,temp2;
            temp1=quick_mod(a-1,p,q);
            temp2=quick_mod(b-1,p,q);
            head=((temp1.m[0][0]*(a1+b1))%mod+(temp1.m[0][1]*b1)%mod+(temp1.m[0][2]*a1)%mod)%mod;
            tail=((temp2.m[0][0]*(a1+b1))%mod+(temp2.m[0][1]*b1)%mod+(temp2.m[0][2]*a1)%mod)%mod;
            ans=((tail-head)%mod+(temp1.m[1][1]*b1)%mod+(temp1.m[1][2]*a1)%mod)%mod;
        }
        if(ans<0)  ans=(ans+mod)%mod;
        printf("%lld\n",ans);
    }
    return 0;
}