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 }