广义斐波那契:,满足,求
acdreamers大佬的推导证明:广义Fibonacci数列找循环节
首先我们肯定知道,这是可以写成矩阵快速幂的形式,问题就在于,n很大的时候呢?
这就是有一个降幂公式,对于要模的数x,循环节就是,∏(pi-1)(pi+1) *x/∏pi
也就是把x的所有质因子pi,把(pi-1)*(pi+1)全乘起来,然后乘上x/pi
至于为什么,这个就不知道了。这个是由txdywjy灵机一闪,猜测出来的积性函数。
直接作为模板用上了。
所以牛客多校第5场的B可以直接降幂做,https://ac.nowcoder.com/acm/contest/885/B
1 #include<cstdio> 2 #include<cstring> 3 typedef long long ll; 4 const int N=1000000+10; 5 struct Node{ 6 int r,c; 7 ll jz[3][3]; 8 Node(){ 9 memset(jz,0,sizeof(jz)); 10 } 11 }; 12 ll x0,x1,a,b,n,md,mmd; 13 char s[N]; 14 int pn,pri[100000+10]; 15 void init(){ 16 ll lim; 17 for(int i=2;i<=100000;i++){ 18 if(!pri[i]) pri[pn++]=i; 19 for(int j=0;j<pn;j++){ 20 lim=1ll*i*pri[j]; 21 if(lim>100000) break; 22 pri[lim]=1; 23 if(i%pri[j]==0) break; 24 } 25 } 26 } 27 ll solve(ll x){ 28 ll ans1=1,ans2=1,xx=x; 29 for(int i=0;i<pn;i++){ 30 if(1ll*pri[i]*pri[i]>x) break; 31 if(x%pri[i]==0){ 32 ans1*=(pri[i]-1)*(pri[i]+1); 33 ans2*=pri[i]; 34 while(x%pri[i]==0) x/=pri[i]; 35 } 36 } 37 if(x>1){ 38 ans1*=(x-1)*(x+1); 39 ans2*=x; 40 } 41 return xx/ans2*ans1; 42 } 43 ll mul(ll x,ll y,ll z){ 44 x%=z; 45 y%=z; 46 ll ans=0; 47 while(y){ 48 if(y&1){ 49 ans+=x; 50 if(ans>=z) ans-=z; 51 } 52 x<<=1; 53 if(x>=z) x-=z; 54 y>>=1; 55 } 56 return ans; 57 } 58 Node Nmul(Node x,Node y,ll z){ 59 Node ans; 60 ans.r=x.r; 61 ans.c=y.c; 62 for(int i=0;i<x.r;i++) 63 for(int j=0;j<y.c;j++) 64 for(int k=0;k<x.c;k++){ 65 ans.jz[i][j]+=mul(x.jz[i][k],y.jz[k][j],z); 66 if(ans.jz[i][j]>=z) ans.jz[i][j]-=z; 67 } 68 return ans; 69 } 70 Node Npow(Node x,ll y,ll z){ 71 Node ans; 72 ans.r=x.r; 73 ans.c=x.c; 74 for(int i=0;i<ans.c;i++) ans.jz[i][i]=1; 75 while(y){ 76 if(y&1) ans=Nmul(ans,x,z); 77 x=Nmul(x,x,z); 78 y>>=1; 79 } 80 return ans; 81 } 82 int main(){ 83 init(); 84 while(~scanf("%lld%lld%lld%lld",&x0,&x1,&a,&b)){ 85 scanf("%s%lld",s,&md); 86 mmd=solve(md); 87 n=0; 88 int lens=strlen(s); 89 for(int i=0;i<lens;i++){ 90 n=mul(n,10,mmd)+s[i]-'0'; 91 if(n>=mmd) n-=mmd; 92 } 93 Node A,T; 94 A.r=2; 95 A.c=1; 96 A.jz[0][0]=x1; 97 A.jz[1][0]=x0; 98 T.r=2; 99 T.c=2; 100 T.jz[0][0]=a; 101 T.jz[0][1]=b; 102 T.jz[1][0]=1; 103 if(n>1){ 104 T=Npow(T,n-1,md); 105 A=Nmul(T,A,md); 106 } 107 printf("%lld ",A.jz[0][0]); 108 } 109 return 0; 110 }
但这题,n明确给出了,我们也可以用十进制的快速幂来处理指数来做。
1 #include<cstdio> 2 #include<cstring> 3 typedef long long ll; 4 const int N=1000000+10; 5 struct Node{ 6 int r,c; 7 ll jz[3][3]; 8 Node(){ 9 memset(jz,0,sizeof(jz)); 10 } 11 }; 12 ll x0,x1,a,b,n,md; 13 char s[N]; 14 Node Nmul(Node x,Node y,ll z){ 15 Node ans; 16 ans.r=x.r; 17 ans.c=y.c; 18 for(int i=0;i<x.r;i++) 19 for(int j=0;j<y.c;j++) 20 for(int k=0;k<x.c;k++){ 21 ans.jz[i][j]+=x.jz[i][k]*y.jz[k][j]%z; 22 if(ans.jz[i][j]>=z) ans.jz[i][j]-=z; 23 } 24 return ans; 25 } 26 Node Npow(Node x,ll y,ll z){ 27 Node ans; 28 ans.r=x.r; 29 ans.c=x.c; 30 for(int i=0;i<ans.c;i++) ans.jz[i][i]=1; 31 while(y){ 32 if(y&1) ans=Nmul(ans,x,z); 33 x=Nmul(x,x,z); 34 y>>=1; 35 } 36 return ans; 37 } 38 int main(){ 39 while(~scanf("%lld%lld%lld%lld",&x0,&x1,&a,&b)){ 40 scanf("%s%lld",s,&md); 41 Node A,T,B; 42 B.r=2;B.c=2; 43 for(int i=0;i<B.c;i++) B.jz[i][i]=1; 44 A.r=2;A.c=1; 45 A.jz[0][0]=x1; 46 A.jz[1][0]=x0; 47 T.r=2;T.c=2; 48 T.jz[0][0]=a; 49 T.jz[0][1]=b; 50 T.jz[1][0]=1; 51 n=0; 52 int lens=strlen(s); 53 for(int i=lens-1;i>=0;i--){ 54 if((s[i]-'0')!=0) B=Nmul(B,Npow(T,s[i]-'0',md),md); 55 T=Npow(T,10,md); 56 } 57 A=Nmul(B,A,md); 58 printf("%lld ",A.jz[1][0]); 59 } 60 return 0; 61 }
不过运行时间上,降幂是更快的,然后举一题需要降幂的例子。
Best Solver
题意:已知,给你整数x,和一个素数M,求[y]%M
对于这题,首先也是先写出矩阵快速幂的形式,这个很好写,就不多说了,
然后这里的指数是很长的,要求的,不能拆成十进制,所以需要降幂,因为m是质数,循环节就是(m-1)*(m+1)
最后还有个根号六要处理,引用别人的
首先(5+2sqrt(6))^n + (5-2sqrt(6))^n 是一个整数,且(5-2sqrt(6))^n < 1
所以 (5+2sqrt(6))^n + (5-2sqrt(6))^n = x + y sqrt(6) + x - ysqrt(6)
所以算出 (5+2sqrt(6))^n的整数部分x,答案就是2*x-1
1 #include<cstdio> 2 #include<cstring> 3 typedef long long ll; 4 struct Node{ 5 int r,c; 6 ll jz[3][3]; 7 Node(){ 8 memset(jz,0,sizeof(jz)); 9 } 10 }; 11 Node Nmul(Node x,Node y,ll z){ 12 Node ans; 13 ans.r=x.r; 14 ans.c=y.c; 15 for(int i=0;i<x.r;i++) 16 for(int j=0;j<y.c;j++) 17 for(int k=0;k<x.c;k++){ 18 ans.jz[i][j]+=x.jz[i][k]*y.jz[k][j]%z; 19 if(ans.jz[i][j]>=z) ans.jz[i][j]-=z; 20 } 21 return ans; 22 } 23 Node Npow(Node x,ll y,ll z){ 24 Node ans; 25 ans.r=x.r; 26 ans.c=x.c; 27 for(int i=0;i<ans.c;i++) ans.jz[i][i]=1; 28 while(y){ 29 if(y&1) ans=Nmul(ans,x,z); 30 x=Nmul(x,x,z); 31 y>>=1; 32 } 33 return ans; 34 } 35 ll poww(ll x,ll y,ll z){ 36 ll ans=1; 37 while(y){ 38 if(y&1) ans=ans*x%z; 39 x=x*x%z; 40 y>>=1; 41 } 42 return ans; 43 } 44 int main(){ 45 int t=1,C,md; 46 ll n,mmd; 47 scanf("%d",&C); 48 while(t<=C){ 49 scanf("%lld%d",&n,&md); 50 mmd=1ll*md*md-1; 51 Node A,T; 52 A.r=2;A.c=1; 53 A.jz[0][0]=5; 54 A.jz[1][0]=2; 55 T.r=2;T.c=2; 56 T.jz[0][0]=5; 57 T.jz[0][1]=12; 58 T.jz[1][0]=2; 59 T.jz[1][1]=5; 60 n=(poww(2,n,mmd)+1)%mmd; 61 if(n>1){ 62 T=Npow(T,n-1,md); 63 A=Nmul(T,A,md); 64 } 65 printf("Case #%d: %lld ",t++,(2ll*A.jz[0][0]-1+md)%md); 66 } 67 return 0; 68 }