[NOI2011]兔农 - jcvb
[NOI2011]兔农
题意:f[1]=f[2]=1;i>2时,f[i]=f[i-1]+f[i-2], if(f[i] % k == 1) f[n]--
求f[n]%p。
白送75分:同时记录f[i]%k和f[i]%p的值,暴力算即可。
100分算法:
以k=7为例,考虑f[i]%k组成的序列:
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,
…
把减1得0的位置标出,并以这些0为界分段,可以发现:
①每段开头必为相同两数,它恰是上一段的最末一位非0数;由于总共只有k-1种余数,所以不超过k段就会出现循环(如果有的话),比如上面k=7时的第3,4段就是循环节。
②记斐波那契数列为fib[i]。假如某段段首数字为x,那么这一段内第i个数即为x*fib[i] % k。若记这一段长度为len,则有x*fib[len]==1 (mod k)。
现在我们试图找到整个数列的循环结构:根据上式,①求x的逆元得到fib[len],②由fib[len]得知len,③用x*fib[len-1]%k算出下一段的段首,重复操作直到发现循环(或者发现这一段永远不终止)。
至于具体实现:①扩欧或者欧拉定理②预处理indfib[y]数组,表示斐波那契数列中模k余y的数第一次出现的下标③预处理fib[i]模k的值。有一个结论:斐波那契数列模k后一定是0,1,1开头的纯循环,而且这个循环节的长度<=6k(不知具体怎么证。。),所以只需暴力算fib数组并同时记录indfib[],发现循环即停止。
注意,假如第①步不存在逆元,或者第②步不存在符合的len,那么这一段将永远不会终止(比如k=8时就是这样),那么整个数列就不存在循环了(可以视作最后一段的长度为无穷大)。
接下来考虑如何用矩阵乘法计算f[n]%p。
[tex]\begin{pmatrix}
a_{i-1} & a_i & 1
\end{pmatrix}\begin{pmatrix}
0 &1 &0 \\
1&1 &0 \\
0&0 &1
\end{pmatrix}=\begin{pmatrix}
a_i &a_{i-1}+a_i &1
\end{pmatrix}[/tex]
[tex]\begin{pmatrix}
a_{i-1} & a_i & 1
\end{pmatrix}\begin{pmatrix}
1 &0 &0 \\
0&1 &0 \\
0&-1 &1
\end{pmatrix}=\begin{pmatrix}
a_{i-1} &a_i-1 &1
\end{pmatrix}[/tex]
分别记这两个3*3矩阵为A,B。令初始矩阵为[tex]\begin{pmatrix}
1 & 0& 1
\end{pmatrix}[/tex],通过对其不断右乘A和B便能实现累加、减1两种操作。对于分出的每一段算出一个矩阵A^len * B,表示这一段的“效果”。
接下来是喜闻乐见的分类讨论时间:假如整个数列是循环的,判断第n项是否在混循环的那部分里,若是则直接把前面几段乘起来,n所在这一段的零头直接用A的次幂算;若不是则先把混循环全部乘起来,然后把循环节全部乘起来,算出循环次数再快速幂,然后再像刚才一样算零头乘上去。若数列不循环倒方便些,也与上面类似,不多说了。
代码很丑陋
#include<iostream> #include<cstdio> #include<cstdlib> #include<algorithm> #include<cstring> using namespace std; typedef long long ll; ll n,k,p; ll fib[10000005],indfib[1000005]={0},inv[1000005]; ll exgcd(ll a,ll b,ll &lx,ll &ly){ ll x=0,y=1,q,t;lx=1,ly=0; while(b){ q=a/b; t=b;b=a%b;a=t; t=x;x=lx-q*x;lx=t; t=y;y=ly-q*y;ly=t; } return a; } void buildinv(){ ll tmp1,tmp2;inv[1]=1; for (ll i=2;i<k;i++) if(exgcd(i,k,tmp1,tmp2)!=1)inv[i]=0; else inv[i]=tmp1<0?tmp1+k:tmp1; } ll ord[1000005];ll tot=0; ll lensum[1000005]; ll vis[1000005]={0}; struct mat{ll a[3][3];void cl(){memset(a,0,sizeof(a));}void set1(){cl();a[0][0]=a[1][1]=a[2][2]=1;}}; mat mul(const mat &a,const mat &b){ mat c; for (int i=0;i<3;i++) for (int j=0;j<3;j++){ c.a[i][j]=0; for (int k=0;k<3;k++)c.a[i][j]+=a.a[i][k]*b.a[k][j]; c.a[i][j]%=p; } return c; } mat po(mat a,ll t){ mat c;c.set1(); while(t){ if(t&1)c=mul(c,a); a=mul(a,a); t>>=1; } return c; } mat s[1000005],A,B; int main() { cin>>n>>k>>p; fib[0]=0,fib[1]=fib[2]=1; for (ll i=3;;i++){ fib[i]=fib[i-1]+fib[i-2];if(fib[i]>=k)fib[i]-=k; if(!indfib[fib[i]])indfib[fib[i]]=i; if(fib[i]==1 && fib[i-1]==1)break; } buildinv(); ll cur=1;lensum[0]=0; A.cl();A.a[0][1]=A.a[1][0]=A.a[1][1]=A.a[2][2]=1; B.cl();B.a[0][0]=B.a[1][1]=B.a[2][2]=1;B.a[2][1]=-1; while(1){ if(vis[cur])break; if(!inv[cur] || !indfib[inv[cur]])break; tot++;vis[ord[tot]=cur]=tot; lensum[tot]=lensum[tot-1]+indfib[inv[cur]]; s[tot]=po(A,indfib[inv[cur]]);s[tot]=mul(s[tot],B); cur=1LL*cur*fib[indfib[inv[cur]]-1]%k; } mat N;N.set1(); if(vis[cur]){ int pre=vis[cur]-1,i;mat M; for (i=1;i<=pre && n>=lensum[i];i++)N=mul(N,s[i]); if(i<=pre){ M=po(A,n-lensum[i-1]); N=mul(N,M); }else{ ll cylen=lensum[tot]-lensum[pre]; n-=lensum[pre]; mat C;C.set1(); for (i=pre+1;i<=tot;i++)C=mul(C,s[i]); C=po(C,n/cylen);N=mul(N,C);n%=cylen; for (i=pre+1;n>=lensum[i]-lensum[pre];i++)N=mul(N,s[i]); M=po(A,n-(lensum[i-1]-lensum[pre])); N=mul(N,M); } }else{ mat M;int i; for (i=1;i<=tot && n>=lensum[i];i++)N=mul(N,s[i]); M=po(A,n-lensum[i-1]); N=mul(N,M); } ll ans=(N.a[0][1]+N.a[2][1])%p;if(ans<0)ans+=p; cout<<ans<<'\n'; return 0; }