大步小步算法(BSGS)及其拓展
大步小步算法,即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; }