事先声明,本博客代码主要模仿accepoc,且仅针对一般如本博主一样的蒟蒻。
这道题不得不说数据良心,给了75分的水分,但剩下25分真心很难得到,因此我们就来讲一讲这剩下的25分。
首先,有数据可知他无心炸long long,因此高精度什么的倒是不用,但10^18的数据范围明显O(n)递推使不靠谱的,又因为本题是建立在斐波那契数列之上,考虑矩阵快速幂优化。
首先先科普一下,斐波那契数列在mod一个数后会形成一个大循环,最大不超过6K(我不会证,有神犇路过望不吝赐教)。因此我们需要一个vi[i]数组记录斐波那契数列在%k下第一次出现i对应的是第几个斐波那契数(一会有用)。至于什么时候跳出循环嘛,第i个数和i-1个数为1时就跳出。
1 scanf("%lld%lld%lld",&n,&k,&p); 2 fib[1]=fib[2]=1; 3 for(int i=3;;i++) 4 { 5 fib[i]=fib[i-1]+fib[i-2]; 6 fib[i]%=k; 7 if(!vi[fib[i]]) vi[fib[i]]=i; 8 if(fib[i]==1&&fib[i-1]==1)break; 9 }
然后说本题最主要的部分,也就是最后25分,让我们以k=7时举例,那么新数列就为
1,1,2,3,5,0,
5,5,3,0,
3,3,6,2,0,
2,2,4,6,3,2,5,0,5,5,3,0,
3,3,6,2,0,
⋯
为0的地方就是原本%k为1的地方,不难看出这里有一个新的循环节,大家可以再举举别的例子,可以发现大部分数都是这样对,大部分。比如8貌似就不行。
其次,我们还可以注意到每一行除以第一个数都是一个斐波那契数列(当然是没有模过之前),因此,设该行长度为len,行首元素为x,那么x*f[len]%k==1,大家想到了什么,逆元!
逆元是个好东西,它可以帮助我们判断是否有循环节,如果x无逆元,说明无循环节,矩阵乘搞到头,那么,如果x有逆元呢?还记得之前的vi数组吗,对,由vi即可判断它到底是该行第几个,如果vi不存在,还是按上面处理。如果存在,直接推出来len,然后进行矩阵快速幂,因为mod1要自减,因此矩阵有两个,为了方便,我们称第一个为nol,第二个为de
上面比较常见,下面一个就是在行末才用的到的矩阵,来消掉1。
然后就是fin[i]数组了,用来标记以i为开头的数列是否已经被结算过,和res[i]搭配着用,res[i]就是以i为开头的数列其中nol和de相乘所得的最终结果,这样在利用循环节的时候就可以很方便的使用了。我们利用行首元素进行循环因为每个行首元素对应且仅对应唯一的数列,因此先通过之前处理出循环节所对应的矩阵求出一个包含整个循环节的矩阵,直接快速幂乱搞,把n%循环节总行数,将flag打上标记,用之前处理出的单行再次进行矩阵乘,直到乘完为止。
1 for(long long t=1;n;) //枚举各行开头 2 { 3 if(!inv[t]) inv[t]=exgcd(t,k,1); 4 if(inv[t]==-1) 5 { 6 ans=ans*ksm(nol,n); 7 break; 8 } 9 else 10 { 11 if(!fin[t]||flag) 12 { 13 fin[t]=1; 14 if(!vi[inv[t]]) 15 { 16 ans=ans*ksm(nol,n); 17 break; 18 } 19 len[t]=vi[inv[t]]; 20 if(n>=len[t]) 21 { 22 n-=len[t]; 23 res[t]=ksm(nol,len[t])*de; 24 ans=ans*res[t]; 25 (t*=fib[len[t]-1])%=k; 26 } 27 else 28 { 29 ans=ans*ksm(nol,n);//对剩下残余的n直接处理 30 break; 31 } 32 } 33 else 34 { 35 long long js=0; 36 node ret(3,3); 37 ret.a[0][0]=ret.a[1][1]=ret.a[2][2]=1; 38 for(long long i=t*fib[len[t]-1]%k;i!=t;(i*=fib[len[i]-1])%=k) 39 { 40 js+=len[i]; 41 ret=ret*res[i]; 42 } 43 js+=len[t]; 44 ret=res[t]*ret; 45 ans=ans*ksm(ret,n/js); 46 n%=js; 47 flag=1; 48 } 49 } 50 }
最后在说几个小坑,第一,尽管k在int范围内,但并不意味着我们就可以放心大胆的使用Int了毕竟在中间计算时会溢出,第二,矩阵乘并不满足交换律,虽然不小心把某些矩阵乘位置搞反并不会全WA,但当怎么看都觉得代码对的时候检查一下这个也是很有必要的。
对于这道题毕竟是NOI难度(虽然在cogs上只有两星半),如果实在无法理解上面的讲解可以自己先抄一下代码,根据代码去理解。
1 #include<iostream> 2 #include<cstdlib> 3 #include<cstdio> 4 #include<cstring> 5 #include<queue> 6 #include<algorithm> 7 #include<cmath> 8 using namespace std; 9 long long k,p; 10 long long fib[6000005]; 11 int vi[1000005]; 12 long long exgcd(long long a,long long b,long long c){ //求逆元 13 if(a==0)return -1; 14 else if(c%a==0) return c/a; 15 long long t=exgcd(b%a,a,((-c%a)+a)%a); 16 if(t==-1)return -1; 17 return (t*b+c)/a; 18 } 19 struct node{ 20 int n,m; 21 long long a[3][3]; 22 node(){} 23 node(int x,int y){ 24 n=x,m=y; 25 memset(a,0,sizeof(a)); 26 } 27 node operator*(const node &b) 28 { 29 node ans(n,b.m); 30 for(int i=0;i<n;i++) 31 { 32 for(int j=0;j<b.m;j++) 33 { 34 for(int k=0;k<m;k++) 35 { 36 (ans.a[i][j]+=(a[i][k]*b.a[k][j])%p)%=p; 37 } 38 } 39 } 40 return ans; 41 } 42 }res[1000005]; 43 node ksm(node a,long long n){ 44 long long x=n; 45 node ans(a.n,a.m); 46 ans.a[1][1]=ans.a[0][0]=ans.a[2][2]=1; 47 while(x) 48 { 49 if((x&1)) ans=ans*a; 50 a=a*a; 51 x>>=1; 52 } 53 return ans; 54 } 55 long long n,inv[1000005],len[1000005]; 56 bool fin[1000005]; 57 node ans,nol,de; 58 void solve(){ 59 nol.n=nol.m=de.n=de.m=3; 60 bool flag=0; 61 ans.n=1,ans.m=3; 62 ans.a[0][0]=ans.a[0][2]=1; 63 nol.a[0][1]=nol.a[1][1]=nol.a[1][0]=nol.a[2][2]=1;//对nol和ni矩阵初始化。 64 de.a[0][0]=de.a[1][1]=de.a[2][2]=1; 65 de.a[2][1]=-1; 66 for(long long t=1;n;) //枚举各行开头 67 { 68 if(!inv[t]) inv[t]=exgcd(t,k,1); 69 if(inv[t]==-1) 70 { 71 ans=ans*ksm(nol,n); 72 break; 73 } 74 else 75 { 76 if(!fin[t]||flag) 77 { 78 fin[t]=1; 79 if(!vi[inv[t]]) 80 { 81 ans=ans*ksm(nol,n); 82 break; 83 } 84 len[t]=vi[inv[t]]; 85 if(n>=len[t]) 86 { 87 n-=len[t]; 88 res[t]=ksm(nol,len[t])*de; 89 ans=ans*res[t]; 90 (t*=fib[len[t]-1])%=k; 91 } 92 else 93 { 94 ans=ans*ksm(nol,n);//对剩下残余的n直接处理 95 break; 96 } 97 } 98 else 99 { 100 long long js=0; 101 node ret(3,3); 102 ret.a[0][0]=ret.a[1][1]=ret.a[2][2]=1; 103 for(long long i=t*fib[len[t]-1]%k;i!=t;(i*=fib[len[i]-1])%=k)//直接跳转下一行开头 104 { 105 js+=len[i]; 106 ret=ret*res[i]; 107 } 108 js+=len[t]; 109 ret=res[t]*ret; 110 ans=ans*ksm(ret,n/js); 111 n%=js; 112 flag=1; 113 } 114 } 115 } 116 } 117 int main(){ 118 scanf("%lld%lld%lld",&n,&k,&p); 119 fib[1]=fib[2]=1; 120 for(int i=3;;i++) 121 { 122 fib[i]=fib[i-1]+fib[i-2]; 123 fib[i]%=k; 124 if(!vi[fib[i]]) vi[fib[i]]=i; 125 if(fib[i]==1&&fib[i-1]==1)break; 126 } 127 solve(); 128 printf("%lld ",(ans.a[0][1]+p)%p); 129 //while(1); 130 return 0; 131 } 132
关于这道题打法有很多,个人代码是看着上方博客打出来的,希望读者不要拘泥于这一种打法。