BSGS学习记录/POJ 3696 The Luckiest number

BSGS

BSGS（baby-step giant-step），即大步小步算法。常用于求解离散对数问题。形式化地说，该算法可以在 (O(sqrt{p})) 的时间内求解

[a^x equiv b pmod p ]

其中 (aperp p)。方程的解 (x) 满足 (0 le x < p)。（在这里需要注意，只要 (aperp p) 就行了，不要求 (p) 是素数）

算法描述

令 (x = A left lceil sqrt p ight ceil - B)，其中 (0le A,B le left lceil sqrt p ight ceil)，则有 (a^{Aleft lceil sqrt p ight ceil -B} equiv b pmod p)，稍加变换，则有 (a^{Aleft lceil sqrt p ight ceil} equiv ba^B pmod p)。

我们已知的是 (a,b)，所以我们可以先算出等式右边的 (ba^B) 的所有取值，枚举 (B)，用 hash/map 存下来，然后逐一计算 (a^{Aleft lceil sqrt p ight ceil})，枚举 (A)，寻找是否有与之相等的 (ba^B)，从而我们可以得到所有的 (x)，(x=A left lceil sqrt p ight ceil - B)。

注意到 (A,B) 均小于 (left lceil sqrt p ight ceil)，所以时间复杂度为 (Thetaleft (sqrt p ight ))，用 map 则多一个 (log)。

上述介绍在CC BY-SA 4.0 和 SATA条款下，转自OIWIKI

P3846 [TJOI2007] 可爱的质数

模板

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
unordered_map<ll,ll> mp;

int main(){
	ll p,a,b;
	scanf("%lld%lld%lld",&p,&a,&b);
	ll sqp=ceil(sqrt(p));//square root of p
	ll rs=b,ap=1,ls=1;//rs/ls: right/left side
	//ap: sqp-th power of a
	for(int i=1;i<=sqp;++i){
		rs=(rs*a)%p;
		ap=(ap*a)%p;
		mp[rs]=i;
	}
	for(int i=1;i<=sqp;++i){
		ls=(ls*ap)%p;// ls=a^{i*sqrt{p}}
		if(mp[ls]){
			printf("%lld
",((i*sqp-mp[ls])%p+p)%p);
			return 0;
		}
	}
	puts("no solution");
	return 0;
}

POJ The Luckiest number

求 (underbrace{888 cdots88}_{ ext{x个}}equiv 0space (modspace L)最小的符合条件的正整数x).

即 ((10^x-1) cdotfrac{8}{9}equiv0space(modspace L))

(9n|8(10^x-1))

(ecause gcd(8,9)=1,gcd(frac{n}{gcd(L,8)},frac{L}{gcd(L,8)})=1)

( herefore 10^xequiv1(modspace frac{9n}{gcd(L,8)}))

注意到与此时x应严格大于0，在处理时，不能使 (x=A left lceil sqrt p ight ceil - B) （(a^{Aleft lceil sqrt p ight ceil} equiv ba^B pmod p)）的B不能为(left lceil sqrt p ight ceil)，第一次循环处理到(left lceil sqrt p ight ceil-1).

时间复杂度(O(sqrt L))

#include<cmath>
#include<cstdio>
#include<cstring>
#include<utility>
using namespace std;
typedef long long ll;
#define N 150000
ll hs[N],nxt[N],hd[N],val[N];
ll cnt;

void ins(ll x,ll v){
	ll p=x%N;
	hs[++cnt]=x;
	val[cnt]=v;
	nxt[cnt]=hd[p];
	hd[p]=cnt;
}

ll fd(ll x){
	ll p=x%N;
	for(ll i=hd[p];i!=-1;i=nxt[i]){
		if(hs[i]==x){
			return val[i];
		}
	}
	return -1;
}

ll __gcd(ll a,ll b){
	return b==0?a:__gcd(b,a%b);
}

ll mul(ll a,ll b,ll c){
	if(a<b){
		swap(a,b);
	}
	ll ret=0,tmp=a;
	while(b){
		if(b&1){
			ret+=tmp;
			if(ret>=c) ret-=c;	
		}
		tmp<<=1;
		if(tmp>=c) tmp-=c;
		b>>=1;
	}	
	return ret%c;
}

int main(){
	ll n;
	ll cs=0;
	while(scanf("%lld",&n),n){
		cnt=0;
		memset(hd,-1,sizeof hd);
		ll p=9*n/__gcd(n,8ll),a=10,b=1;
		if(__gcd(p,a)!=1){
			printf("Case %lld: 0
",++cs);
			continue;
		}
		ll sqp=ceil(sqrt((double)p));//square root of p
		ll rs=b,ap=1,ls=1;//rs/ls: right/left side
		//ap: sqp-th power of a
		for(ll i=1;i<sqp;++i){
			rs=(rs*a)%p;
			ap=(ap*a)%p;
			ins(rs,i);
		}
		if(sqp) ap=(ap*a)%p;
		bool f=0;
		for(ll i=1;i<=sqp;++i){
			ls=mul(ls,ap,p);// ls=a^{i*sqrt{p}}
			ll j;
			if((j=fd(ls))!=-1){
				f=1;
				printf("Case %lld: %lld
",++cs,((i*sqp-j)%p+p)%p);
				break;
			}
		}
		if(!f){
			printf("Case %lld: 0
",++cs);
		}
	}
	return 0;
}

另外本题还有另一种解法

易证:(aperp n, a^x equiv 1(modspace n))的最小整数解使(varphi(n))的约数，枚举约数即可。

#include<cmath>
#include<cstdio>
#include<cstring>
#include<utility>
#include<iostream>
using namespace std;
typedef long long ll;

ll __gcd(ll a,ll b){
	return b==0?a:__gcd(b,a%b);
}

ll mul(ll a,ll b,ll c){
	if(a<b){
		swap(a,b);
	}
	ll ret=0,tmp=a;
	while(b){
		if(b&1){
			ret+=tmp;
			if(ret>=c) ret-=c;	
		}
		tmp<<=1;
		if(tmp>=c) tmp-=c;
		b>>=1;
	}	
	return ret%c;
}

ll qpow(ll a,ll b,ll c){
	ll ret=1,tmp=a;
	while(b){
		if(b&1){
			ret=mul(ret,tmp,c);
		}
		tmp=mul(tmp,tmp,c);
		b>>=1;
	}
	return ret%c;
}

int main(){
	ll n;
	ll cs=0;
	while(scanf("%lld",&n),n){
		ll p=9*n/__gcd(n,8ll);
		if(__gcd(p,10)!=1){
			printf("Case %lld: 0
",++cs);
			continue;
		}
		ll fai=p,t=p,r=sqrt(p);
		for(ll i=2;i<=r;++i){
			if(t%i==0){
				fai=fai/i*(i-1);
				t/=i;
				while(t%i==0){
					t/=i;
				}	
			}
		}
		if(t>1){
			fai=fai/t*(t-1);
		}
		r=sqrt(fai);
		ll ans=1e18;
		bool f=0;
		for(ll i=1;i<=r;++i){
			if(fai%i==0){
				if(qpow(10,i,p)==1){
					ans=min(ans,i);
					f=1;
				}
				if(qpow(10,fai/i,p)==1){
					ans=min(ans,fai/i);
					f=1;
				}
			}
		}
		printf("Case %lld: %lld
",++cs,f?ans:0);
	}
	return 0;
}