大步小步算法(BSGS)及其拓展

tonyfang posted @ 2015年8月11日 19:50 in Algorithm with tags c++ Algorithm , 1786 阅读

大步小步算法,即B(baby)S(step)G(giant)S(step)算法。

是指类似于$ A^{x} \equiv B (mod C) $的方程,已知A、B、C,求x的算法

一种朴素的方法,我们可以枚举$x$,从$0$到$C-2$(若$x=C-1$,那么$A^{C-1} \equiv 1 (mod C)$【费马小定理】)。

这样的复杂度为$O(C)$,对于C比较小的时候,仍是可以承受的,若在C比较大的情况下,我们就要使用BSGS算法了.

原始的BSGS只能解决$C$为质数的情况。或$gcd(A,C)=1$的情况。

$m=\lceil \sqrt{C} \rceil$,$x=i \times m+j$,然后$A^{x} = (A^{m}) ^{i} \times A^{j} $,其中$0 \leq i < m,0 \leq j < m$

那么我们枚举$i$即可,将时间复杂度降到$O(\sqrt{C})$

那么原式为$A^{(i+1)m} \times B^{-1} \equiv A^{m-j} (mod C)$,其中$B^{-1}$表示$B$的逆元

那么由于$m-j \leq m$,我们开一个数组哈希或map存即可,复杂度$O(1)$ 到 $O(logC)$

那么即可求出$A^{j}$,只需要哈希查表即可得到$j$。

那么由$x=i \times m+j$即可得到$x$。

拓展BSGS:

一开始的方程为 $A^x \times a + C \times b = B (a,b为整数) $

那么设$t=gcd(A,C)$,显然,$B$如果不能被$t$整除,就无解了。

下面来看个方程:

$ (\frac{A}{t})^{x'} \times a' + (\frac{C}{t}) \times b' = \frac {B}{t} $

两端乘$t$,得到$A^{x'+1} \times a' + C \times b' = B$

那么我们就有办法了,算$x'$,然后算$x=x'+1$即可。

继续迭代直至$gcd(A,C)=1$,假设做了$cnt$次,那么$x=x'+cnt$。

然后我们用普通BSGS求出$x'$即可。

HDU 2815 拓展的的BSGS,BSGS也可以快速求log: HERE

前面提交一直WA,原因是因为有个地方没取模……要用到逆元得写ex-gcd。

 

# include <stdio.h>
# include <string.h>
# include <map>
# include <math.h>
using namespace std;
typedef long long ull;
ull gcd(ull a,ull b) {
	return b?gcd(b,a%b):a;
}
inline void exgcd(ull a,ull b,ull &x,ull &y) {
	if(!b) x=1,y=0;
	else {
		exgcd(b,a%b,y,x);
		y-=a/b*x;
	}
}
inline ull inverse(ull a,ull b) {
	ull x,y;
	exgcd(a,b,x,y);
	if(x<0) x+=b;
	return x;
}
ull bsgs(ull A,ull B,ull C) {
    ull m,v,e=1,i;
    m=ceil(sqrt(C));
    map<ull,ull> hash;
    hash[1]=m;
    for (i=1;i<m;++i) {
        e=(e*A)%C;
        if(!hash[e]) hash[e]=i;
    } 
    e=e*A%C;
    e=inverse(e,C);
    for (i=0;i<m;++i) {
        if(hash[B]) {
            ull ret=hash[B];
            hash.clear();
            return i*m+(ret==m?0:ret);
        }
        B=(B*e)%C;
    }
    return -1;
}
ull extbsgs(ull a,ull b,ull c) {
	ull t,d=1,cnt=0; 
	while((t=gcd(a,c))!=1) {
		if(b%t) return -1;
		b/=t, c/=t;
		d=d*a/t%c;
		cnt++;
		if(d==b) return cnt;
	}
	b=b*inverse(d,c)%c;
	ull ret=bsgs(a,b,c);
	if(ret==-1) return -1;
	else return cnt+ret;
}
int main() {
	ull k,p,n,ans;
	while(~scanf("%I64d%I64d%I64d",&k,&p,&n)) {
		if (n>=p) {
			printf("Orz,I can’t find D!\n");
			continue;
		}
		else if (n==0) {
			printf("0\n"); 
			continue;
		}
		else {
			ans=extbsgs(k,n,p);
			if (ans==-1) printf("Orz,I can’t find D!\n");
			else printf("%I64d\n",ans);
		}
	}
	return 0;
}

POJ 2417 不拓展的BSGS(限定了$C$为质数):HERE   

# include <stdio.h>
# include <math.h>
# include <map>
using namespace std;
long long ksm(long long A,long long B,long long P) {
	if (B==0) return 1;
	if (B==1) return A%P;
	long long k=ksm(A,B>>1,P);
	k=(k*k)%P;
	if(B&1) return k*A%P;
	else return k; 
}
long long bsgs(long long A,long long B,long long C) {
	long long m,v,e=1,i;
	m=ceil(sqrt(C));
	v=ksm(A,C-1-m,C);
	map<long long,long long> hash;
	hash[1]=m;
	for (i=1;i<m;++i) {
		e=(e*A)%C;
		if(!hash[e]) hash[e]=i;
	} 
	for (i=0;i<m;++i) {
		if(hash[B]) {
			long long ret=hash[B];
			hash.clear();
			return i*m+(ret==m?0:ret);
		}
		B=(B*v)%C;
	}
	return -1;
}
int main() {
	long long A,B,C;
	while(scanf("%lld%lld%lld",&C,&A,&B)!=EOF) {
		long long ans=bsgs(A,B,C);
		if (ans==-1) printf("no solution\n");
		else printf("%lld\n",ans);
	}
	return 0;
}

登录 *


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