[BSGS算法]纯水斐波那契数列

学弟在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+];
int mod;
struct mat
{
int z[][];
mat(){memset(z,,sizeof(z));}
mat operator*(mat b)
{
mat c;int i,j,k;
for(i=;i<;++i)for(j=;j<;++j)
for(k=;k<;++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=;i<;++i)for(int j=;j<;++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>>=,t=t*t)if(x&)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[][]=f.z[][]=f.z[][]=;
fread(s,,MN,stdin);
for(l=,r=MN;l<=r;)
if(s[mid=l+r>>])l=mid+;
else i=mid,r=mid-;
while(s[--i]<''||s[i]>'')s[i]=;
while(s[i]>=''&&s[i]<='')--i;
for(l=i;s[l]<''||s[l]>'';)s[l--]=;
while(s[++i]>=''&&s[i]<='')mod=(mod<<)+(mod<<)+s[i]-'',s[i]=;
for(p=k=pow(K),i=;i<=K;++i)p=p*k,mp[p]=i;
for(p=f,i=;i<=K;++i,p=p*f)if(x=mp[p]){rp=x*K-i;break;}
for(x1=x2=x3=x4=,i=;s[i];i+=)
x1=(x1*+s[i-]-'')%rp,
x2=(x2*+s[i-]-'')%rp,
x3=(x3*+s[i-]-'')%rp,
x4=(x4*+s[i ]-'')%rp;
x=(x1*+x2*+x3*+x4)%rp;
for(i-=;s[i];++i)x=((x<<)+(x<<)+s[i]-'')%rp;
p.z[][]=p.z[][]=;if(x<)x+=rp;
printf("%d",(p*pow(x-)).z[][]);
}

正常版

#include<cstdio>
#include<cstring>
#include<map>
using namespace std;
#define ll long long
#define MN 30000000
#define K 115000
char s[MN+];
int mod;
struct mat
{
int z[][];
mat(){memset(z,,sizeof(z));}
mat operator*(mat b)
{
mat c;int i,j,k;
for(i=;i<;++i)for(j=;j<;++j)
for(k=;k<;++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=;i<;++i)for(int j=;j<;++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>>=,t=t*t)if(x&)r=r*t;
return r;
}
map<mat,int> mp;
int main()
{
int i;ll x,rp;
f.z[][]=f.z[][]=f.z[][]=;
scanf("%s%d",s,&mod);
for(p=k=pow(K),i=;i<=K;++i)p=p*k,mp[p]=i;
for(p=f,i=;i<=K;++i,p=p*f)if(x=mp[p]){rp=x*K-i;break;}
for(i=x=;s[i];++i)x=((x<<)+(x<<)+s[i]-'')%rp;
p.z[][]=p.z[][]=;if(x<)x+=rp;
printf("%d",(p*pow(x-)).z[][]);
}
上一篇:BSGS算法学习笔记


下一篇:浅谈Manacher算法与扩展KMP之间的联系