hdu 1588（Fibonacci矩阵求和）

　　题目的大意就是求等差数列对应的Fibonacci数值的和，容易知道Fibonacci对应的矩阵为[1,1,1,0]，因为题目中f[0]=0，f[1]=1，所以推出最后结果f[n]=(A^n-1).a，所以 f(g(i))= f(k*i+b)= (A^(k*i+b-1)).a，i从 0取到 n-1，取出公因式 A^(b-1)（因为矩阵满足分配率），然后所求结果可化为 A^(b-1) * (A^0 + A^k + A^2k +....+ A^(n-1)k)，化到这里后难点就是求和了，一开始我尝试暴力求和（每个A^k可以用快速幂求出，logn级别），即O(n)的做法，结果TLE了，预料之中，这时我竟傻乎乎地套用等比数列求和公式，即(A^nk -A^0) /(A^k -E)，按比例放大再相减是没错，问题是不是简单的相除……总之思路应该是错的了，后来看别人的博客后才知道原来可以用二分来求和的，即 A^0 + A^k + A^2k +....+ A^(n-1)k = (A^0 + A^k + A^2k +...+ A^(n/2-1)k) *(E + A^(n/2)k)，递归来求和，处理好 n的奇偶性即可，但我还是调试了好久，就因为在一些细节问题上出错却检查不出来，下面附上代码：

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<cstdlib>
 4 using namespace std;
 5 typedef long long LL;
 6 LL mod= 1000000007;        //初不初始化都没问题，只是为了防止忘记读入时产生的异常退出
 7
 8 struct matrix{
 9     LL a,b,c,d;
10     matrix(LL a=0, LL b=0, LL c=0, LL d=0): a(a),b(b),c(c),d(d) {}
11     matrix operator +(const matrix &m2){
12         return matrix((a+m2.a)%mod,(b+m2.b)%mod,(c+m2.c)%mod,(d+m2.d)%mod);
13     }
14     matrix operator *(const matrix &m2){
15         return matrix((a*m2.a%mod+b*m2.c%mod)%mod, (a*m2.b%mod+b*m2.d%mod)%mod, (c*m2.a%mod+d*m2.c%mod)%mod, (c*m2.b%mod+d*m2.d%mod)%mod);
16     }
17     //一开始等比数列求和的思路要用到除法，后来才发现是错的，不过也不删了，就放在这吧
18     matrix operator /(const matrix &m2){
19         //二维求逆很好求的
20         matrix inv= matrix(m2.d,-m2.b,-m2.c,m2.a);
21         LL tmp= m2.a*m2.d-m2.b*m2.c;
22         return matrix((a*inv.a/tmp%mod+b*inv.c/tmp%mod)%mod, (a*inv.b/tmp%mod+b*inv.d/tmp%mod)%mod, (c*inv.a/tmp%mod+d*inv.c/tmp%mod)%mod, (c*inv.b/tmp%mod+d*inv.d/tmp%mod)%mod);
23     }
24 };
25 // A为 Fibonacci矩阵，E为单位矩阵，设为全局变量更方便一些
26 matrix A(1,1,1,0),E(1,0,0,1);
27
28 //简单的快速幂
29 matrix quick_mod(matrix m, LL b){
30     if(b==-1)    return matrix(0,1,1,-1);
31     //若指数为-1，返回矩阵 A^-1，相当于A^1的逆(算了好久T.T)
32     matrix res(E);        //res一开始为单位矩阵
33     while(b){
34         if(b&1)      res= res*m;
35         m= m*m;
36         b>>=1;
37     }
38     return res;
39 }
40
41 //二分法计算 A^0 + A^k + A^2k +....+ A^(n-1)k 的和
42 //即 sum =(A^0 + A^k + A^2k +...+ A^(n/2-1)k) *(E + A^(n/2)k)，分治的思想
43 matrix quick_sum(LL k, LL n){
44     if(n==1)    return E;    //若 n为1，即计算 A^0，此时返回的是 E！而不是 A！一开始没想到错在这里，卡了好久 T.T
45     if(n%2==0)    return quick_sum(k,n/2)*(E+quick_mod(A,(n/2)*k));
46     //else    return quick_sum(k,(n-1)/2)*(E+quick_mod(A,((n-1)/2)*k))*quick_mod(A,(n-1)*k);
47     else    return quick_sum(k,n-1) + quick_mod(A,(n-1)*k);
48     //这样的写法比起上一行的虽然会多一次调用函数的开销，但可读性增强，代码逻辑更清晰
49 }
50
51 int main()
52 {
53     //freopen("1588in.txt","r",stdin);
54     LL k,b,n;
55     while(~scanf("%I64d%I64d%I64d%I64d",&k,&b,&n,&mod)){
56         //最后的结果就是 A^b-1 *(A0 + A^k + A^2k +...+ A(n-1)k) 的 a
57         matrix ans= quick_mod(A,b-1) * quick_sum(k,n);
58         printf("%I64d\n",ans.a);
59     }
60     return 0;
61 }