[caioj 1488及hdu 1588]Gauss Fibonacci

有一個g數組，f數組。g[i]=k*i+b，k與i都是給定的常數。f[0]=0,f[1]=1,f[i]=f[i-1]+f[i-2] (i>=2)。現在，要你求f[g[0]]+f[g[1]]+…+f[g[n-1]]的值，但可能最後的值會很大，是以結果需要mod M。

這題看到了它的資料範圍，就令人很容易想到用矩陣乘法加快速幂，但是這是一道有難度的題，耗了我不少時間。它的難點就是如何構造矩陣，對于這個問題，網上有兩種做法，而本蒟蒻給出一種做法，是巧妙地運用到了矩陣乘法的結合律的，如果用這種方法就像是pku 3233和pku 3070的結合，大家可以先做一下這兩題，我的blog也詳解了這兩題。首先我們求出構造出Fibonacci數列的矩陣a，為{0,1}{1,1}。為什麼，因為a * {f1}{f2}=(f2}{f3(f1+f2)}（兩個相鄰括号之間有空行），然後初始序列s為{0}{1}， 因為a^x * {0}{1}={f[x]}{f[x+1]}。

這樣我們就可以用兩次快速幂來分别算出A矩陣（表示a^b),B矩陣（表示a^k)，而f[g[0]]+f[g[1]]+…+f[g[n-1]]可以看為f[b]+f[k+b]+…+f[k * (n-1)+b] ,則上式可為[A+a^(k+b)+…+a^(k^(n-1)+b)] * s，提公因式，則變成[a^b * (1+a^k+…+(a^k)^(n-1))] * s，轉換後的公式為[A*(1+B+B^2+…+B^(n-1))] * s。對于B+B^2+…+B^(n-1)可以用兩次快速幂做，詳見http://blog.csdn.net/lixuanjing2016/article/details/77235674

搞定B+B^2+…+B^(n-1)後，就可以算出轉換後的公式，這題就迎刃而解了。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
using namespace std;
struct node
{
    long long a[][];
    node()
    {
        memset(a,,sizeof(a));
    }
};
int M;
node jiafa(node a,node b)
{
    node c;
    for(int i=;i<=;i++)
    {
        for(int j=;j<=;j++)
        {
            c.a[i][j]=(c.a[i][j]+a.a[i][j]+b.a[i][j])%M;
        }
    }
    return c;
}
node chengfa1(node a,node b)
{
    node c;
    for(int i=;i<=;i++)
    {
        for(int j=;j<=;j++)
        {
            for(int k=;k<=;k++)
            {
                c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j])%M;
            }
        }
    }
    return c;
}
node chengfa2(node a,node b)
{
    node c;
    for(int i=;i<=;i++)
    {
        for(int k=;k<=;k++)
        {
            c.a[i][]=(c.a[i][]+a.a[i][k]*b.a[k][])%M;
        }
    }
    return c;
}
node power(node a,int x)
{
    node ans;
    ans.a[][]=ans.a[][]=;
    while(x>)
    {
        if(x%==)ans=chengfa1(a,ans);
        a=chengfa1(a,a);
        x/=;
    }
    return ans;
}
node getsum(node a,int k)
{
    if(k==)return a;
    node t;
    if(k%==)
    {
        t=getsum(a,k/);
        t=jiafa(t,chengfa1(t,power(a,k/)));
    }
    else
    {
        t=getsum(a,k/);
        t=jiafa(t,chengfa1(t,power(a,k/)));
        t=jiafa(t,power(a,k));
    }
    return t;
}
int main()
{
    int k,b,n,x;
    while(scanf("%d%d%d%d",&k,&b,&n,&M)!=EOF)
    {
        n--;
        node f,ans,ss,per1,per2,s;
        f.a[][]=;f.a[][]=;
        ss.a[][]=ss.a[][]=;
        s.a[][]=s.a[][]=s.a[][]=;
        per1=power(s,k);
        per2=power(s,b);
        ans=chengfa1(per2,jiafa(getsum(per1,n),ss));
        //機關矩陣ss就相當于1的存在
        f=chengfa2(ans,f);
        printf("%lld\n",f.a[][]);
    }
    return ;
}