BSGS(baby-step giant-step)及扩展

orz hhw posted @ 2015年7月14日 19:11 in 算法学习 with tags 模板 数论 算法学习 BSGS , 894 阅读

网上一直搜BSGS搜不出来,这几天BZOJ清理水题,发现SDOI2011计算器很多人A,我却不会第三个操作,于是发现了这个算法的全称baby-step giant-step,总算找到了一些资料,那些定理估计我也搞不懂,把模板放到这里

求最小的x使得a^x=b(mod p)的。如果p是质数,或者说p与a互质,那么:任意取一自然数m(m<=p)。可令x=k*m-r,其中(m>r>=0)。代入方程:a^(k*m)=b*(a^r)(mod p)。将r从0循环到m-1,计算出所有b*(a^r)(mod p),然后把每一个取值所对应的r存在一个哈希表里。如果p比较大就开map,如果p比较小就直接用数组把。并且这个取值和r不是一一对应的,但是这并无大碍。然后计算a^m,并且将k从1循环到(p/m),对于每一个k计算出a^(k*m),然后在哈希表中寻找有没有对应的r存在。如果存在那么就找到了一个解x=k*m-r。如果用map做哈希表:效率O((m+n/m)logm)。在m=[sqrt(n)]时效率最高。

简单来说,步骤为:①令m=ceil(√p)

②枚举A^0,A^1,...,A^(m-1),将这些数插入哈希表
③令D=A^m 枚举D^0,D^1,...,D^m
④对于枚举到的每一个D^i,利用EXGCD求出D^i*x≡B (mod C)的唯一非负整数解,去Hash表中寻找x是否出现过。
⑤如果找到x=A^j,那么答案为i*m+j
⑥如果枚举结束仍未找到x,则无解。
 
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
typedef long long LL;
long long v,n,m,ans,i,j,e,a,b;
map<LL,LL>hash;
LL pow_mod(LL a,LL b,LL n)  {
    LL s=1;  
    for(;b;b>>=1){
        if(b&1)s=(s*a)%n;
        a=(a*a)%n;
    }  
    return s;  
}
LL bsgs(LL a,LL b,LL n){
	m=ceil(sqrt(n+0.5));//x=i*m+j
	v=pow_mod(a,n-m-1,n);//a^m*v=1(mod n)  
	hash[1]=m;e=1;
	for(i=1;i<m;i++){//建哈希表,保存x^0,x^1,...,x^m-1
		e=(e*a)%n;
		if(!hash[e])hash[e]=i;
	}
	for(i=0;i<m;i++){//每次增加m次方,遍历所有1<=x<=n
		if(hash[b]){
			LL num=hash[b];
			if(num==m)num=0;
			hash.clear();//需要清空,不然会爆内存
			return i*m+num;
		}
		b=(b*v)%n;
	}
	return -1;
}
int main(){
	while(scanf("%lld%lld%lld",&n,&a,&b)!=EOF){
		LL ans=bsgs(a,b,n);
		if(ans==-1)puts("no solution");else printf("%lld\n",ans);
	}
}

SDOI2011计算器

#include<cstdio>
#include<map>
#include<cmath>
using namespace std;
typedef long long LL;
map<LL,LL>hash;
LL n,k,a,b,p,x,y,ans,t;
LL pow(LL a,LL b,LL p){//快速幂 
	LL sum=1;
	for(;b;b>>=1){
		if(b&1)sum=(sum*a)%p;
		a=(a*a)%p;
	}
	return sum;
}
void exgcd(LL a,LL b,LL &x,LL &y){//扩展欧几里得 
	if(!b){x=1;y=0;ans=a;return;}
	exgcd(b,a%b,x,y);
	LL t=x;x=y;y=t-a/b*y;
}
LL bsgs(LL a,LL b,LL p){//BSGS
	a%=p;b%=p;
	if(!a&&!b){return 1;}//特判 
	if(!a)return -1;//x=i*m+j
	LL m=ceil(sqrt(p)),v=pow(a,p-m-1,p),num,i,e=1;//a^m*v=1(mod n)
	hash.clear();//需要清空,不然会爆内存
	hash[1]=m;
	for(i=1;i<m;i++){//建哈希表,保存x^0,x^1,...,x^m-1
		e=(e*a)%p;
		if(!hash[e])hash[e]=i;
	}
	for(i=0;i<m;i++){
		if(hash[b]){
			num=hash[b];
			if(num==m)num=0;
			return i*m+num;
		}
		b=(b*v)%p;
	}
	return -1;
}
int main(){
	scanf("%lld%lld",&n,&k);
	while(n--){
		scanf("%lld%lld%lld",&a,&b,&p);
		if(k==1){
			printf("%lld\n",pow(a,b,p));
		}else if(k==2){
			exgcd(a,p,x,y);
			if(b%ans){puts("Orz, I cannot find x!");continue;}
			x=x*b/ans;t=p/ans;x=(x%t+t)%t;
			printf("%lld\n",x);
		}else if(k==3){
			ans=bsgs(a,b,p);
			if(ans==-1)puts("Orz, I cannot find x!");else printf("%lld\n",ans);
		}
	}
}

扩展BSGS:可以处理模非质数的情况,要先gcd搞一搞,原理嘛并不知道

#include<cstdio>
#include<cmath>
#include<map>
#define LL __int64
LL gcd(LL a,LL b){return b==0?a:gcd(b,a%b);}
void gcd_mod(LL a,LL b,LL &d,LL &x,LL &y){
    if(!b){d=a;x=1;y=0;}
    else{gcd_mod(b,a%b,d,y,x);y-=x*(a/b);}
}
LL bsgs (LL a,LL b,LL n,LL c,LL d){
    if(b>=n)return -1;
    LL m,v,e=1,i,x,y,dd;
    m=ceil(sqrt(n+0.5));//x=i*m+j
    std::map<LL,LL>f;f[1]=m;
    for(i=1;i<m;i++){//建哈希表,保存a^0,a^1,...,a^m-1
        e=(e*a)%n;
        if(!f[e])f[e]=i;
    }
    e=(e*a)%n;//e=a^m
    for(i=0;i<m;i++){//每次增加m次方,遍历所有1<=f<=n
        gcd_mod(d,n,dd,x,y);//d*x+n*y=1-->(d*x)%n=1-->d*(x*b)%n==b
        x=(x*b%n+n)%n;
        if(f[x]){
            LL num=f[x];
            f.clear();//需要清空,不然会爆内存
            return c+i*m+(num==m?0:num);
        }
        d=(d*e)%n;
    }
    return -1;
}
int main(){
    LL a,b,n;//求a^x mod n=b 中x的最小值
    while(scanf("%I64d%I64d%I64d",&a,&n,&b)!=EOF){
        LL ans=0,c=0,d=1,t;
        while((t=gcd(a,n))!=1){
            if(b%t){ans=-1;break;}
            c++;n=n/t;b=b/t;d=d*a/t%n;
            if(d==b){ans=c;break;}//特判下是否成立。
        }
        if(ans!=0){if(ans==-1)printf("Orz,I can’t find D!\n");else printf("%I64d\n",ans);
        }else{
            ans=bsgs(a,b,n,c,d);
            if(ans==-1)printf("Orz,I can’t find D!\n");else printf("%I64d\n",ans);
        }
    }
}

BZOJ3122 先求出数列通项公式,Xi=aXi-1 +b,Xi+w=a(Xi-1+w),设Yi=Xi+w,则w=b/(a-1),求最小的i使得Y1*a^(i-1)=Yn,即a^(i-1)=Yn/Y1,套用BSGS解决。注意a=0和1时候的特判。a=1时要求逆元判断答案

#include<bits/stdc++.h>
using namespace std;typedef long long LL;LL v,m,i,e,w,T,p,a,b,x1,t;
LL pow(LL a,LL b,LL p){for(v=1,a%=p;b;a=a*a%p,b>>=1)if(b&1)v=v*a%p;return v;}
LL bsgs(LL a,LL b,LL p){
    LL m=ceil(sqrt(p)),v=pow(pow(a,m,p),p-2,p),i,e=1;map<LL,LL>mp;mp[1]=0;
    for(i=1;i<m;i++)if(!mp.count(e=e*a%p))mp[e]=i;
    for(i=0;i<m;b=b*v%p,i++)if(mp.count(b))return i*m+mp[b];
    return -2;
}
LL G1(){return b?((t-x1+p)%p*pow(b,p-2,p)%p+1):-1;}
LL G2(){return bsgs(a,(t*(1-a)-b)%p*pow((x1*(1-a)-b)%p,p-2,p)%p,p)+1;}
LL G(){
	if(x1==t)return 1;
	if(!a)return b==t?2:-1;
	return a==1?G1():G2();
}
int main(){
	for(scanf("%lld",&T);T--;printf("%lld\n",G()))scanf("%lld%lld%lld%lld%lld",&p,&a,&b,&x1,&t);		
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter