Yet Another Number Sequence

题意：

$F_0 = 0, F_1 = 1, F_n = F_{n-1} + F_{n-2}$

求解$sum_{i=1}^n{ F_i i^K }  mod  10^9+7$。

解法：

记$S(n,m) = sum_{i=1}^n { F_i i^m}$

这样有：

$$S(2n,m)  = sum_{i=1}^n{F_i i^m} + sum_{i=1}^n{F_{i+n} (i+n)^m}$$

记$G = lgroup egin{matrix} 1 & 1 \ 1 & 0 end{matrix} group$

有$F_n = G_{1,1}$

将$S(n,m)$的含义换为对应求和的矩阵。

1.考虑从$S(n,m)$推到$S(2n,m)$

$$S(2n,m) = sum_{i=1}^n{G^n i^m} + G^n sum_{i=1}^n{G^{i} (i+n)^m}$$

$$S(2n,m) = sum_{i=1}^n{G^n i^m} + G^n sum_{i=1}^n{G^{i} sum_{r=0}^m{i^rm^{m-r}C_m^r}}$$

$$S(2n,m) = sum_{i=1}^n{G^n i^m} + G^n sum_{r=0}^m{m^{m-r}C_m^r S(n,r)}$$

这样$O(K)$完成单次转移。 

2.考虑从$S(n,m)$推到$S(n+1,m)$

$$S(n+1,m) = S(n,m) + (n+1)^m G^{n+1}$$

这样分治下去，总效率$O(K^2 logn)$

（注意本题中n过大，要先对于n取模再进行乘法）

#include <iostream>
#include <cstdio>
#include <cstring>

#define LL long long
#define P 1000000007LL

using namespace std;
    
struct MA
{
    LL a[2][2];
    
    void init()
    {
        memset(a,0,sizeof(a));
    }
    
    MA operator*(const MA &x)const
    {
        MA c;
        for(int i=0,j,k;i<2;i++)
            for(j=0;j<2;j++)
            {
                c.a[i][j]=0;
                for(k=0;k<2;k++)
                {
                    c.a[i][j]+=a[i][k]*x.a[k][j]%P;
                    if(c.a[i][j]>=P) c.a[i][j]-=P;
                }
            }
        return c;
    }
    
    MA operator*(const LL &x)const
    {
        MA c;
        for(int i=0,j;i<2;i++)
            for(j=0;j<2;j++)
                c.a[i][j] = a[i][j]*x%P;
        return c;
    }
    
    MA operator+(const MA &x)const
    {
        MA c;
        for(int i=0,j;i<2;i++)
            for(j=0;j<2;j++)
                c.a[i][j] = (a[i][j]+x.a[i][j])%P;
        return c;
    }
    
    void print()
    {
        puts("Matrix");
        for(int i=0,j;i<2;i++)
        {
            for(j=0;j<2;j++) cout<<a[i][j]<<' ';
            cout<<endl;
        }
    }
};

MA G,Gn;
MA S[2][50];
LL n,power[50],C[50][50];
int K,now;

void solve(LL n)
{
    if(n==1)
    {
        now=0;
        Gn = G;
        for(int i=0;i<=K;i++) S[now][i] = G;
        return;
    }
    solve(n>>1);
    now^=1;
    MA tmp;
    power[0]=1;
    for(int k=1;k<=K;k++) power[k] = power[k-1]*((n>>1LL)%P)%P;
    for(int k=0;k<=K;k++)
    {
        tmp.init();
        for(int r=0;r<=k;r++)
            tmp = tmp + ( S[now^1][r] * (C[k][r]*power[k-r]%P) );
        S[now][k] = S[now^1][k] + (Gn * tmp);
    }
    Gn = Gn*Gn;
    if(n&1)
    {
        Gn = Gn*G;
        power[0]=1LL;
        for(int k=0;k<=K;k++)
        {
            S[now][k] = S[now][k] + (Gn * power[k]);
            power[k+1] = power[k]*(n%P)%P;
        }
    }
}

int main()
{
    G.a[0][0]=1; G.a[0][1]=1;
    G.a[1][0]=1; G.a[1][1]=0;
    while(~scanf("%I64d%d",&n,&K))
    {
        C[0][0]=1;
        for(int i=1;i<=K;i++)
        {
            C[i][0]=1;
            for(int j=1;j<=i;j++)
            {
                C[i][j] = C[i-1][j-1]+C[i-1][j];
                if(C[i][j]>=P) C[i][j] -= P;
            }
        }
        solve(n);
        cout << S[now][K].a[0][0] << endl;
    }
    return 0;
}