BSGS 算法，即 Baby Step,Giant Step 算法、拔山盖世算法。
计算 (a^x equiv b pmod p)。
(p)为质数时
特判掉 (a,p) 不互质的情况。
由于费马小定理 (x^{p-1} equiv 1 pmod p) 当 (p) 为质数，则要是暴力的话只需要枚举到 (p-1) 即可。
假设 (x=it-j)，其中 (t= lceil sqrt p
ceil,j in [0,t])，方程变为 (a^{it-j} equiv b pmod p)，即 (a^{it} equiv ba^j pmod p)。我们惊喜地发现，左右最多也就 (t) 个左右种可能的取值（这就是 (t) 为什么取那个值的原因），那我们枚举 (j)，把 (ba^j) 所对应的 (j) 都存起来，然后枚举 (i) 找有无对应即可。
保存的话用手写 hash 比 map 快很多的。
若 (ba^j) 冲突怎么办？答：存 (j) 大的。因为要想 (x) 尽量小，就要让 (j) 尽量大。
第一个找到的 (i) 和它对应的 (j) 就是答案。为什么？答：因为 (i) 变化 (1)，(x) 变化的幅度是 (t)。
时间复杂度 (mathrm{O}(sqrt{p}))。
计算器：
#include <iostream>
#include <cstdio>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
int T, k;
ll y, z, p;
map<int,int> d;
ll work1(ll a, ll b, ll p){
ll re=1;
while(b){
if(b&1) re = (re * a) % p;
a = (a * a) % p;
b >>= 1;
}
return re;
}
ll exgcd(ll a, ll b, ll &x, ll &y){
if(!b){
x = 1;
y = 0;
return a;
}
ll re=exgcd(b, a%b, x, y);
ll qwq=x;
x = y;
y = qwq - a / b * y;
return re;
}
void work2(){
ll u, v;
ll gcd=exgcd(y, p, u, v);
if(z%gcd) printf("Orz, I cannot find x!
");
else
printf("%lld
", ((u*z/gcd)%(p/gcd)+(p/gcd))%(p/gcd));
}
void work3(){
d.clear();
int m=sqrt(p);
if(m*m!=p) m++;
for(int i=0; i<=m; i++)
d[int((z*work1(y, i, p))%p)] = i;//当然这里也可以换掉快速幂，因为递推顺便就能求幂了
if(y%p==0){
if(z%p==1) printf("0
");
else if(z%p==0) printf("1
");
else printf("Orz, I cannot find x!
");
return ;
}
for(int i=1; i<=m; i++){
int tmp=work1(y, i*m, p);
if(d.find(tmp)!=d.end()){
printf("%lld
", (ll)i*m-d[tmp]);
return ;
}
}
printf("Orz, I cannot find x!
");
}
int main(){
cin>>T>>k;
while(T--){
scanf("%lld %lld %lld", &y, &z, &p);
if(k==1) printf("%lld
", work1(y, z, p));
if(k==2) work2();
if(k==3) work3();
}
return 0;
}
测试数据：
2 3
39 26 13
6 11 5
1
0
(p)无限制时
计算 (a^x equiv b pmod p)，最难过的事情就是 ((a,p)
ot = 1)。那我们就想办法消掉一些 (p) 的因子让 (a,p) 互素。记 (d=(a,p))。
倘若 (d
mid b)（不整除的语法是
mid），则 (b=1) 时方程有解，解为 (x=0)。
倘若 (d mid b) 且 (d=1)，则就是普通的 bsgs。
否则，
当然，还有更难过的是 ((a,p/d)) 可能还不为 (1)……那就继续搞下去，直到
且 ((a,frac{p}{prod_{i=1}^kd_i})=1)。
当然还有更更难过的是还没到互素的时候就有 (frac{p}{prod_{i=1}^kappa d_i}
mid b)……那唯一的可能就是 (x=0) 了。
然后发现如果真正的解 (x<k) 的话比较难过，那我们在累积公约数的时候可以顺便判断 (x in [0,k)) 时是否是解。
对于 (x geq k) 的情况，我们为了赏心悦目，把方程换元为
其中 ((a',p')=1)。
这样就是一个普通的 bsgs。最终的解是 (x=x'+k)。
Clever Y:
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
using namespace std;
typedef long long ll;
int a, b, p;
const int mod=32767;//这不是质数，不好
int gcd(int a, int b){
return !b?a:gcd(b,a%b);
}
int ksm(int a, int b, int p){
int re=1;
while(b){
if(b&1) re = ((ll)re * a) % p;
a = ((ll)a * a) % p;
b >>= 1;
}
return re;
}
struct HASHTABLE{
int nxt[200005], dai[200005], ret[200005], hea[200005], hmn;
void clear(){
memset(nxt, 0, sizeof(nxt));
memset(hea, 0, sizeof(hea));
hmn = 0;
}
void insert(int x, int y){
int tmp=x%mod;
for(int i=hea[tmp]; i; i=nxt[i])
if(dai[i]==x){
ret[i] = y;
return ;
}
nxt[++hmn] = hea[tmp];
dai[hmn] = x; ret[hmn] = y;
hea[tmp] = hmn;
}
int query(int x){
int tmp=x%mod;
for(int i=hea[tmp]; i; i=nxt[i])
if(dai[i]==x)
return ret[i];
return -1;
}
}hashTable;
int exbsgs(int a, int b, int p){
hashTable.clear();
a %= p; b %= p;
if(b==1) return 0;
int cnt=0, d=1, g;
do{
g = gcd(a, p);
if(b%g) return -1;
p /= g; b /= g;
d = ((ll)d * a / g) % p;
cnt++;
if(b==d) return cnt;
}while(g!=1);
int m=ceil(sqrt(p)), t=b;
for(int i=0; i<=m; i++){
hashTable.insert(t, i);
t = ((ll)t * a) % p;
}
g = ksm(a, m, p);
t = ((ll)d * g) % p;
for(int i=1; i<=m; i++){
int re=hashTable.query(t);
if(re>=0) return i*m-re+cnt;
t = ((ll)t * g) % p;
}
return -1;
}
int main(){
while(scanf("%d %d %d", &a, &p, &b)!=EOF){
if(a+b+p==0) break;
int re=exbsgs(a, b, p);
if(re==-1) printf("No Solution
");
else printf("%d
", re);
}
return 0;
}