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 \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \rgroup$

有$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取模再进行乘法）

 

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cstring>
  4 
  5 #define LL long long
  6 #define P 1000000007LL
  7 
  8 using namespace std;
  9     
 10 struct MA
 11 {
 12     LL a[2][2];
 13     
 14     void init()
 15     {
 16         memset(a,0,sizeof(a));
 17     }
 18     
 19     MA operator*(const MA &x)const
 20     {
 21         MA c;
 22         for(int i=0,j,k;i<2;i++)
 23             for(j=0;j<2;j++)
 24             {
 25                 c.a[i][j]=0;
 26                 for(k=0;k<2;k++)
 27                 {
 28                     c.a[i][j]+=a[i][k]*x.a[k][j]%P;
 29                     if(c.a[i][j]>=P) c.a[i][j]-=P;
 30                 }
 31             }
 32         return c;
 33     }
 34     
 35     MA operator*(const LL &x)const
 36     {
 37         MA c;
 38         for(int i=0,j;i<2;i++)
 39             for(j=0;j<2;j++)
 40                 c.a[i][j] = a[i][j]*x%P;
 41         return c;
 42     }
 43     
 44     MA operator+(const MA &x)const
 45     {
 46         MA c;
 47         for(int i=0,j;i<2;i++)
 48             for(j=0;j<2;j++)
 49                 c.a[i][j] = (a[i][j]+x.a[i][j])%P;
 50         return c;
 51     }
 52     
 53     void print()
 54     {
 55         puts("Matrix");
 56         for(int i=0,j;i<2;i++)
 57         {
 58             for(j=0;j<2;j++) cout<<a[i][j]<<' ';
 59             cout<<endl;
 60         }
 61     }
 62 };
 63 
 64 MA G,Gn;
 65 MA S[2][50];
 66 LL n,power[50],C[50][50];
 67 int K,now;
 68 
 69 void solve(LL n)
 70 {
 71     if(n==1)
 72     {
 73         now=0;
 74         Gn = G;
 75         for(int i=0;i<=K;i++) S[now][i] = G;
 76         return;
 77     }
 78     solve(n>>1);
 79     now^=1;
 80     MA tmp;
 81     power[0]=1;
 82     for(int k=1;k<=K;k++) power[k] = power[k-1]*((n>>1LL)%P)%P;
 83     for(int k=0;k<=K;k++)
 84     {
 85         tmp.init();
 86         for(int r=0;r<=k;r++)
 87             tmp = tmp + ( S[now^1][r] * (C[k][r]*power[k-r]%P) );
 88         S[now][k] = S[now^1][k] + (Gn * tmp);
 89     }
 90     Gn = Gn*Gn;
 91     if(n&1)
 92     {
 93         Gn = Gn*G;
 94         power[0]=1LL;
 95         for(int k=0;k<=K;k++)
 96         {
 97             S[now][k] = S[now][k] + (Gn * power[k]);
 98             power[k+1] = power[k]*(n%P)%P;
 99         }
100     }
101 }
102 
103 int main()
104 {
105     G.a[0][0]=1; G.a[0][1]=1;
106     G.a[1][0]=1; G.a[1][1]=0;
107     while(~scanf("%I64d%d",&n,&K))
108     {
109         C[0][0]=1;
110         for(int i=1;i<=K;i++)
111         {
112             C[i][0]=1;
113             for(int j=1;j<=i;j++)
114             {
115                 C[i][j] = C[i-1][j-1]+C[i-1][j];
116                 if(C[i][j]>=P) C[i][j] -= P;
117             }
118         }
119         solve(n);
120         cout << S[now][K].a[0][0] << endl;
121     }
122     return 0;
123 }
View Code

 

posted @ 2017-03-23 16:48  lawyer'  阅读(146)  评论(0)  编辑  收藏  举报