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的奇偶性即可，但我还是调试了好久，就因为在一些细节问题上出错却检查不出来，下面附上代码：

#include<cstdio>
#include<cstring>
#include<cstdlib>
using namespace std;
typedef long long LL;
LL mod= 1000000007;        //初不初始化都没问题，只是为了防止忘记读入时产生的异常退出 

struct matrix{
    LL a,b,c,d;
    matrix(LL a=0, LL b=0, LL c=0, LL d=0): a(a),b(b),c(c),d(d) {}
    matrix operator +(const matrix &m2){
        return matrix((a+m2.a)%mod,(b+m2.b)%mod,(c+m2.c)%mod,(d+m2.d)%mod);
    }
    matrix operator *(const matrix &m2){
        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);
    }
    //一开始等比数列求和的思路要用到除法，后来才发现是错的，不过也不删了，就放在这吧 
    matrix operator /(const matrix &m2){
        //二维求逆很好求的 
        matrix inv= matrix(m2.d,-m2.b,-m2.c,m2.a);
        LL tmp= m2.a*m2.d-m2.b*m2.c;
        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);
    }
}; 
// A为 Fibonacci矩阵，E为单位矩阵，设为全局变量更方便一些 
matrix A(1,1,1,0),E(1,0,0,1);

//简单的快速幂 
matrix quick_mod(matrix m, LL b){
    if(b==-1)    return matrix(0,1,1,-1);
    //若指数为-1，返回矩阵 A^-1，相当于A^1的逆(算了好久T.T)
    matrix res(E);        //res一开始为单位矩阵 
    while(b){
        if(b&1)      res= res*m;
        m= m*m;
        b>>=1;
    }
    return res;
}

//二分法计算 A^0 + A^k + A^2k +....+ A^(n-1)k 的和
//即 sum =(A^0 + A^k + A^2k +...+ A^(n/2-1)k) *(E + A^(n/2)k)，分治的思想 
matrix quick_sum(LL k, LL n){
    if(n==1)    return E;    //若 n为1，即计算 A^0，此时返回的是 E！而不是 A！一开始没想到错在这里，卡了好久 T.T 
    if(n%2==0)    return quick_sum(k,n/2)*(E+quick_mod(A,(n/2)*k));
    //else    return quick_sum(k,(n-1)/2)*(E+quick_mod(A,((n-1)/2)*k))*quick_mod(A,(n-1)*k);
    else    return quick_sum(k,n-1) + quick_mod(A,(n-1)*k);
    //这样的写法比起上一行的虽然会多一次调用函数的开销，但可读性增强，代码逻辑更清晰 
}

int main()
{
    //freopen("1588in.txt","r",stdin);
    LL k,b,n;
    while(~scanf("%I64d%I64d%I64d%I64d",&k,&b,&n,&mod)){
        //最后的结果就是 A^b-1 *(A0 + A^k + A^2k +...+ A(n-1)k) 的 a
        matrix ans= quick_mod(A,b-1) * quick_sum(k,n);
        printf("%I64d
",ans.a);
    }
    return 0;
}