bzoj 3657 斐波那契数列(fib.cpp/pas/c/in/out)
空间 512M 时限2s
【题目描述】
有n个大于1的正整数a1,a2,…,an，我们知道斐波那契数列的递推式是f(i)=f(i-1)+f(i-2)，现在我们修改这个递推式变为f(i)=f(i-1)+f(i-2)+r(i-1)，其中r(x)为a1,a2,…,an中为x的约数的个数。现在要求f(m) mod 19940417的值。注：初值f(1)=1,f(2)=1
输入格式：
第一行两个数n,m。
接下来一行n个正整数a1,a2,…,an。
输出格式：
输出一行仅一个数，f(m) mod 19940417的值。
样例输入：
3 7
2 2 3
样例输出：
33
数据范围：
30%的数据n<=1000,m<=1000
另外20%的数据 n=0,m<=109
100%的数据n<=100000,m<=109，2<=ai<=109
题解：
对于100%的数据，我们可以先考虑把fib[i]=fib[i-1]+fib[i-2] 的答案先用矩阵快速幂跑出来。然后依次输入ai，来看每个ai对fib[m]的影响，因为fib(i)=fib(i-1)+fib(i-2)+r(i-1)，所以每一个ai，在它k倍(k*ai<=M)的地方的斐波那契值都会产生+1的影响。我们考虑如果对于斐波那契数列的第i项我们对它加一个1并且继续进行后面的递推的话，那么第j项(j>i)的值就是fib[j]+fib[j-i+1]。所以实际上我们可以对于每个ai分别处理，对于ai，它会给最后的答案贡献fib[m mod ai]+fib[ai+(m mod ai)]+…保证[]内的值小于等于m。
但如果只是一个一个让答案加上fib[k*ai+(m%ai)]，还是会超时，肯定要用到矩阵快速幂来优化，假设我们让B为表示fib[m%ai]的矩阵，那么f[k*ai+(m%ai)]可以表示为B*A^k*ai，然后解决的就是SUM = (A + A^2 + A^3 + ... + A^B)%C的问题(讲解)。
1 #include <cstdio> 2 #include <algorithm> 3 using namespace std; 4 typedef long long LL; 5 const int mod=19940417; 6 struct mat { 7 int a,b,c,d; 8 }ZR,E,F,Ans; 9 int n,m; 10 mat pre[33],pw[33]; 11 mat operator*(mat X,mat Y) { 12 mat Z; 13 Z.a=((LL)X.a*Y.a+(LL)X.b*Y.c)%mod; 14 Z.b=((LL)X.a*Y.b+(LL)X.b*Y.d)%mod; 15 Z.c=((LL)X.c*Y.a+(LL)X.d*Y.c)%mod; 16 Z.d=((LL)X.c*Y.b+(LL)X.d*Y.d)%mod; 17 return Z; 18 } 19 mat operator+(mat X,mat Y) { 20 mat Z; 21 Z.a=(X.a+Y.a)%mod; 22 Z.b=(X.b+Y.b)%mod; 23 Z.c=(X.c+Y.c)%mod; 24 Z.d=(X.d+Y.d)%mod; 25 return Z; 26 } 27 mat fpm(mat a,int b) { 28 mat w=E; 29 while(b){ 30 if(b&1) w=w*a; 31 a=a*a; 32 b>>=1; 33 } 34 return w; 35 } 36 mat vsum(int n){ 37 if(n==0) return ZR; 38 if(n==1) return E; 39 int m=1,t=0; 40 while(m<=n) m<<=1,++t; 41 m>>=1,--t; 42 return pre[t]+pw[t]*vsum(n-m); 43 } 44 void prepare(mat A){//A矩阵是系数矩阵的ai次方 45 for(int i=0;i<=30;++i){ 46 if(i==0) pw[i]=A; 47 else pw[i]=pw[i-1]*pw[i-1]; 48 if(i==0) pre[i]=E;//单位矩阵 49 else pre[i]=pre[i-1]*(E+pw[i-1]); 50 } 51 } 52 mat solve(int d) { 53 if(d>=m) return ZR; 54 int k=(m-1)/d; 55 prepare(fpm(F,d)); 56 return fpm(F,m-1-k*d)*vsum(k); 57 } 58 int main() { 59 // freopen("fib.in" , "r", stdin); 60 // freopen("fib.out", "w", stdout); 61 scanf("%d%d",&n,&m); 62 if(m<=2){ 63 printf("1\n"); 64 return 0; 65 } 66 E.a=E.d=1; 67 F.a=F.b=F.c=1; 68 Ans=Ans+fpm(F,m-1);//先算出纯 fib序列 69 70 for(int i=1;i<=n;++i){// 71 int x; 72 scanf("%d",&x); 73 Ans=Ans+solve(x); 74 } 75 printf("%d\n",Ans.a); 76 }
