学弟在OJ上加了道“非水斐波那契数列”，求斐波那契第n项对1,000,000,007取模的值，n<=10^15，随便水过后我决定加一道升级版，说是升级版，其实也没什么变化，只不过改成n<=10^30000000，并对给定p取模，0<p<2^31。一样很水嘛大家说对不对。
下面来简单介绍一下BSGS算法，BSGS（Baby steps and giant steps），又称包身工树大步小步法，听上去非常高端，其实就是一个暴力搜索。比如我们有一个方程，a^x≡b (mod c)，x是未知数，怎么算？难道是什么神奇的数学公式？直觉和经验告诉我们这玩意儿并不好算，无奈的我们只好枚举x，欧拉dalao告诉我们a^φ(c)≡1 (mod c)，再怎么不济，我们O(c)也能算出答案……平时我们常见的模数一般都有十亿左右，明显要T啦，怎么办？
别急，我们来回顾一个经典问题，给你有n个数的集合，问他有多少个子集和为x。(n<=100,000，x<=10^18)
别想了，我也不会做，只好把数据改一改，n<=20？暴力就过啦。n<=50？此时有个常见的做法，我们枚举前n/2个数的和，开个map存下，再枚举后n/2个数，用x减掉枚举出的和再去map里找……复杂度O(2^(n/2))（可能还要多个log）。我们再考虑下子集积？发现并没有区别，枚举前n/2个数的积，存下，再枚举后n/2个数，用x乘上枚举出的积的逆元（假设取模）……是不是想到了什么？我们像快速幂那样把a^x表示成a^1,a^2,a^4,a^8……的积，解那个方程不就是做子集积吗？实际上不用这么麻烦，我们把x表示成只有两位的k进制数，我们就只要枚举两位了，令x=Ak+B，则a^(Ak+B)≡b (mod c)，移项得a^Ak≡b*a^-B (mod c)，我们先枚举A，算出等式左边，存入map里，再枚举B，算出等式右边，去map里找，就能算出答案啦。要保证这个k进制数只有两位，k自然大约是c^0.5。实际上我们可以稍微把k改大一点（或者A多枚举一点），把x改成x=Ak-B，这样有什么好处呢？你代回去再移一次项就会发现，我们不用算逆元啦！于是我们得到一个较为高效的求解指数方程的算法，复杂度大概是O(c^0.5)（视具体实现可能会多个log）。
好了，我们回到斐波那契数列，10^30000000明显就算我们用某klogklogn的算法也不是很好受（纯属口胡）。实际上很容易想到，斐波那契对一个数取模必然有循环节，因为每一项只和前两项有关，两项数的取值最多p^2种（p为模数），所以循环节至多p^2。事实上，某篇论文证明了，斐波那契数列对p取模的循环节不会超过6p（其实论文中还给出了对于不同的p，循环节的计算方法，不过复杂度应该和接下来要讲的BSGS做法相同(理论上论文做法更优，但我们肯定不想总依靠玄学的论文吧)，有兴趣的可以百度一下）。那么我们如何用BSGS算出循环节呢？假设斐波那契数列的递推矩阵为A，我们只要求出A^x≡A (mod p)的一个大于1的解就可以了，求解过程和整数方程并没有太大区别。求出循环节后把那巨长无比的数字串取模下，随便水过该题。
以下是加了一些奇怪的常数优化的代码……如果看的时候遇到一些不知道在干嘛的地方，跳过就好了。
#include<cstdio> #include<cstring> #include<map> using namespace std; #define ll long long #define MN 30000050 #define K 115000 char s[MN+5]; int mod; struct mat { int z[2][2]; mat(){memset(z,0,sizeof(z));} mat operator*(mat b) { mat c;int i,j,k; for(i=0;i<2;++i)for(j=0;j<2;++j) for(k=0;k<2;++k)c.z[i][j]=(c.z[i][j]+(ll)z[i][k]*b.z[k][j])%mod; return c; } friend bool operator<(mat a,mat b) { for(int i=0;i<2;++i)for(int j=0;j<2;++j) { if(a.z[i][j]<b.z[i][j])return true; if(a.z[i][j]>b.z[i][j])return false; } } }f,k,p; mat pow(ll x) { mat r=f,t=f; for(--x;x;x>>=1,t=t*t)if(x&1)r=r*t; return r; } map<mat,int> mp; int main() { int i,l,r,mid;ll x,x1,x2,x3,x4,rp; f.z[0][1]=f.z[1][0]=f.z[1][1]=1; fread(s,1,MN,stdin); for(l=1,r=MN;l<=r;) if(s[mid=l+r>>1])l=mid+1; else i=mid,r=mid-1; while(s[--i]<'0'||s[i]>'9')s[i]=0; while(s[i]>='0'&&s[i]<='9')--i; for(l=i;s[l]<'0'||s[l]>'9';)s[l--]=0; while(s[++i]>='0'&&s[i]<='9')mod=(mod<<3)+(mod<<1)+s[i]-'0',s[i]=0; for(p=k=pow(K),i=2;i<=K;++i)p=p*k,mp[p]=i; for(p=f,i=1;i<=K;++i,p=p*f)if(x=mp[p]){rp=x*K-i;break;} for(x1=x2=x3=x4=0,i=3;s[i];i+=4) x1=(x1*10000+s[i-3]-'0')%rp, x2=(x2*10000+s[i-2]-'0')%rp, x3=(x3*10000+s[i-1]-'0')%rp, x4=(x4*10000+s[i ]-'0')%rp; x=(x1*1000+x2*100+x3*10+x4)%rp; for(i-=3;s[i];++i)x=((x<<3)+(x<<1)+s[i]-'0')%rp; p.z[0][0]=p.z[0][1]=1;if(x<2)x+=rp; printf("%d",(p*pow(x-1)).z[0][0]); }
正常版
#include<cstdio> #include<cstring> #include<map> using namespace std; #define ll long long #define MN 30000000 #define K 115000 char s[MN+5]; int mod; struct mat { int z[2][2]; mat(){memset(z,0,sizeof(z));} mat operator*(mat b) { mat c;int i,j,k; for(i=0;i<2;++i)for(j=0;j<2;++j) for(k=0;k<2;++k)c.z[i][j]=(c.z[i][j]+(ll)z[i][k]*b.z[k][j])%mod; return c; } friend bool operator<(mat a,mat b) { for(int i=0;i<2;++i)for(int j=0;j<2;++j) { if(a.z[i][j]<b.z[i][j])return true; if(a.z[i][j]>b.z[i][j])return false; } } }f,k,p; mat pow(ll x) { mat r=f,t=f; for(--x;x;x>>=1,t=t*t)if(x&1)r=r*t; return r; } map<mat,int> mp; int main() { int i;ll x,rp; f.z[0][1]=f.z[1][0]=f.z[1][1]=1; scanf("%s%d",s,&mod); for(p=k=pow(K),i=2;i<=K;++i)p=p*k,mp[p]=i; for(p=f,i=1;i<=K;++i,p=p*f)if(x=mp[p]){rp=x*K-i;break;} for(i=x=0;s[i];++i)x=((x<<3)+(x<<1)+s[i]-'0')%rp; p.z[0][0]=p.z[0][1]=1;if(x<2)x+=rp; printf("%d",(p*pow(x-1)).z[0][0]); }