UPDATE ON 2020.12.23 更新了原代码中的错误
BSGS算法
(BSGS(baby-step giant-step))，即大步小步算法。常用于求解离散对数问题。
形式化地说，该算法可以在 (O sqrt{p}) 的时间内求解
[a^x equiv b (mod p)
]
其中 (gcd(a,p)=1)
根据欧拉定理，当 (gcd(a,p)=1) 时，(a^{varphi(p)}mod p=1)
而 (a^0mod p=1)
所以从 (1 sim varphi(p)) 出现了一个循环节
因此，我们只需要求出 (a^1 sim a^{varphi(p)})
然后分别带入左右两边判断是否相等
一般的题目中 (p) 为质数，所以 (varphi(p)=p-1)
一般就用 (p) 替代
但是这样做时间复杂度太高，可以用分块的思想优化
令 (a^x=a^{km-t})，那么原式就变成了
[a^{km} equiv ba^t (mod p)
]
我们只需要分别求出 (t=1 sim m) 时 (ba^t) 的值，将其存到哈希表里
然后从 (1 sim p/m) 枚举 (k) ，判断哈希表中有没有和 (a^{km}) 相等的数
当 (m=sqrt{p}) 时最优
代码
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e5+5,mod=1e5+3;
int p,a,n,blo,h[maxn],tot=1;
struct asd{
int nxt,val,num;
}b[maxn];
void ad(int val,int id){
rg int now=val%mod;
for(rg int i=h[now];i!=-1;i=b[i].nxt){
rg int cs=b[i].val;
if(cs==val){
if(b[i].num<id) b[i].num=id;
return;
}
}
b[tot].val=val;
b[tot].nxt=h[now];
b[tot].num=id;
h[now]=tot++;
}
int cx(int val){
rg int now=val%mod;
for(rg int i=h[now];i!=-1;i=b[i].nxt){
rg int cs=b[i].val;
if(cs==val) return b[i].num;
}
return -1;
}
int main(){
memset(h,-1,sizeof(h));
p=read(),a=read(),n=read();
blo=sqrt(p)+1;//不要忘了加上一个1保证正确性
rg int now=n,cs=1,haha;
for(rg int i=1;i<=blo;i++){
now=1LL*now*a%p;
ad(now,i);
cs=1LL*cs*a%p;
}
now=1;
for(rg int i=1;i<=blo;i++){
now=1LL*now*cs%p;
haha=cx(now);
if(haha!=-1){
printf("%d
",i*blo-haha);
return 0;
}
}
printf("no solution
");
return 0;
}
扩展BSGS算法
求解
[a^x equiv b (mod p)
]
其中 (a,p) 不一定互质
这个时候就不能再找 (1 sim varphi(p)) 的循环节了，因为不满足欧拉定理的条件
我们把式子变一下，可以写成
[a imes a^{x-1}+kp=b
]
的形式
根据裴蜀定理，当 (gcd(a,p)|b) 时方程有解，否则无解
所以在有解的情况下，我们可以把方程左右两边同时除以 (gcd(a,p))，等号不变
此时原方程变为
[frac{a}{gcd(a,p)} imes a^{x-1}+k frac{p}{gcd(a,p)}=frac{b}{gcd(a,p)}
]
把 (frac{p}{gcd(a,p)}) 看成新的 (p)
把 (frac{b}{gcd(a,p)}) 看成新的 (b)
同时记录一下左边的系数之积 (frac{a}{gcd(a,p)})，并不断地提出 (a)
这样在最后一定会出现 (gcd(a,p)=1) 的情况，用普通的 (bsgs) 即可
不要忘了把提出去的 (a) 的系数加上
代码
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<cstring>
#define rg register
inline int read(){
rg int x=0,fh=1;
rg char ch=getchar();
while(ch<'0' || ch>'9'){
if(ch=='-') fh=-1;
ch=getchar();
}
while(ch>='0' && ch<='9'){
x=(x<<1)+(x<<3)+(ch^48);
ch=getchar();
}
return x*fh;
}
const int maxn=1e5+5,mod=1e5+3;
int h[maxn],tot=1;
struct asd{
int nxt,val,num;
}b[maxn];
void ad(int val,int id){
rg int now=val%mod;
for(rg int i=h[now];i!=-1;i=b[i].nxt){
rg int cs=b[i].val;
if(cs==val){
if(b[i].num<id) b[i].num=id;
return;
}
}
b[tot].val=val;
b[tot].nxt=h[now];
b[tot].num=id;
h[now]=tot++;
}
int cx(int val){
rg int now=val%mod;
for(rg int i=h[now];i!=-1;i=b[i].nxt){
rg int cs=b[i].val;
if(cs==val) return b[i].num;
}
return -1;
}
int bsgs(rg int a,rg int p,rg int n,rg int xs){
memset(h,-1,sizeof(h));
tot=1;
rg int blo=sqrt(p)+1;
rg int ac1=n,ac2=1,ac3=xs;
for(rg int i=1;i<=blo;i++){
ac1=1LL*ac1*a%p;
ad(ac1,i);
ac2=1LL*ac2*a%p;
}
for(rg int i=1;i<=blo;i++){
ac3=1LL*ac3*ac2%p;
ac1=cx(ac3);
if(ac1!=-1) return i*blo-ac1;
}
return -1;
}
int gcd(rg int aa,rg int bb){
return bb==0?aa:gcd(bb,aa%bb);
}
int main(){
rg int a,p,n,haha=1,now,ans,cnt;
rg bool jud=0;
while(1){
a=read(),p=read(),n=read();
if(a==0 && p==0 && n==0) break;
haha=1,jud=0,cnt=0;
a%=p,n%=p;
if(n==1){
printf("0
");
} else if(a==0 && n==0){
printf("0
");
} else if(a==0){
printf("No Solution
");
} else if(n==0){
while(1){
if(p==1){
printf("%d
",cnt);
break;
}
now=gcd(a,p);
if(now==1 && p!=-1){
printf("No Solution
");
break;
}
cnt++,p/=now;
}
} else {
while(1){
now=gcd(a,p);
if(n%now!=0){
jud=1;
break;
}
if(now==1){
ans=bsgs(a,p,n,haha);
if(ans==-1) jud=1;
ans+=cnt;
break;
}
cnt++,p/=now,n/=now;
haha=1LL*haha*(a/now)%p;
if(haha==n){
ans=cnt;
break;
}
}
if(jud){
printf("No Solution
");
} else {
printf("%d
",ans);
}
}
}
return 0;
}