扩展BSGS算法

先来康康这个问题。
给定三个正整数$A,B,P$，求关于$x$的方程的最小正整数解：
$$A^x\equiv B(modP)$$
我们令$x=a\lceil \sqrt{P}\rceil -b,a\in N,b\in N$
由欧拉定理：
$$A^{\varphi (P)}\equiv 1(modP)$$
对同余方程变形：
$$A^{\varphi (P)+x}\equiv A^x(modP)⟺A^x\equiv A^x(modP)$$
所以说，我们只需要考虑$[1,\varphi (P)]$之间的解即可。
$\because$当$i$为质数时，$\varphi (i)=i-1$
$\therefore i-1\le \varphi (i)$
这样，我们只需考虑$[1,P-1]$之间的解。
只需令$a\in [1,\lceil \sqrt{P}\rceil ],b\in [0,\lceil \sqrt{P}\rceil )$，即可使$a\lceil \sqrt{P}\rceil -b$能取遍$[1,P-1]$中的每一个整数。
然后，对原方程变形：
$$A^{a\lceil \sqrt{P}\rceil }\equiv B\cdot A^b(modP)$$
由于方程左右两端最多有$\lceil \sqrt{P}\rceil$个值，所以我们可以把一边的值都算出来，存在一个哈希表中，然后再拿另一边的值去查找是否有相等的值。若有，则返回；否则无解。
上面的方法就是今天想介绍的BSGS（拔山盖世）算法。

---

但是，我们注意到一个细节，就是在这个地方：
$$A^{a\lceil \sqrt{P}\rceil -b}\equiv B(modP)$$
我们发现，这个式子成立，需要保证$A^b$在模$P$意义下的逆元存在，但是，当$gcd(A,P)\ne 1$时，它没有逆元。
欧几里得的定理告诉我们：
$$A^x\equiv B(modP)⟺A\ast A^{x-1}+kP=B，k\in \mathrm{Z}$$
我们可以多次将系数$A$,$B$,$P$除以$gcd(A,P)$。显然，当$B\nmid gcd(A,P)$时，原方程无解。
假设上面的过程反复进行了$d$次，系数$A$,$B$,$P$分别变为$A^{\prime }$,$B^{\prime }$,$P^{\prime }$，
则有：
$$A^{\prime }\ast A^{x-d-1}+k{P}^{\prime }=B^{\prime }⟺A^{x-d-1}\equiv B^{\prime }\ast A^{\prime -1}(modp)$$
然后套用上面的方法即可。这就是扩展BSGS（更加的拔山盖世）（大雾）算法。

---

下面是一道例题：
【BZOJ 2480 & SPOJ 3105】Mod
已知数a,p,b，求满足$a^x\equiv b(modp)$的最小自然数x。
一道裸题：

#include<bits/stdc++.h>
#define ll long long
using namespace std;
map<int, int> h;
int gcd(int a, int b) {return b ? gcd(b, a % b) : a;}
inline ll ksm(ll a, ll b, ll mod)
{
    ll ret = 1, h = a;
    while(b)
    {
        if(b & 1)
            ret *= h, ret %= mod;
        h *= h, h %= mod, b >>= 1;
    }
    return ret;
}
inline int exbsgs(int a, int b, int p)
{
    a %= p, b %= p;
    if(b == 1) return 0; 										//特判2
    int g, d = 0;
    ll c = 1;
    while((g = gcd(a, p)) > 1)
    {
        if(b % g) return -1;
        b /= g, p /= g, c = c * (a / g) % p, ++d;
        if(c == b) return d; 						//特判3（其实等价于特判2）
    }
    h.clear();
    int t = int(sqrt(p * 1.0) + 1);
    ll base = b;									//计算右半部分的值
    for(int i = 0; i < t; i++)
        h[base] = i, base = base * a % p;
    base = ksm(a, t, p); 							//用左半边的值验证
    ll now = c;
    for(int i = 1; i <= t + 1; i++)
    {
        now = now * base % p;
        if(h.count(now))
            return i * t - h[now] + d;				//返回解
    }
    return -1;										//无解
}
int main()
{
    int a, b, p;
    while(~scanf("%d%d%d", &a, &p, &b) && a && b && p)
    {
        if(p == 1) {puts("0"); continue;} 			//特判1
        int ans = exbsgs(a, b, p);
        if(ans == -1) puts("No Solution");
        else printf("%d\n", ans);
    }
}