BSGS（baby step gaint step）+快速幂+exgcd逆元

一个神奇的算法，快速求出 A^x=B(mod C)，C是素数，求最小的非负x值若x有解，那么 0<=x < C，为什么？

由费马小定理可以得a^(c-1) = 1(mod c)
那么0~c-1必定是一个循环节（不一定是最小的）。既然是%c，那么B一定是0到c-1之间的一个数。最坏的条件下，a的c-1以内次方%c的余数各不相同，那么在0~C-1时一定存在一个x满足条件。
那么不是最坏情况时，就能找到一个更小的在0~c-1的循环节。

BSGS的算法是这样的：

首先取m=sqrt(c)向上取整。（为什么取sqrt(c)？我也不是很懂）

然后先预处理a的0到m次方。

a^x=b ( %c ) 设x=i*m+j; 即: i为x/m，j为x%m。 a^(i*m+j)=b; b * (a^(-m))^i =
a^j ( %c )

先枚举j，把右边存起来（Hash 下一步用二分查找，或哈希表）枚举i，如果左边的数值曾经存储过（b * (a^(-m))^i =
a^j）,则 x=i*m+j，得解。
快速幂，其实很好理解，对于求x^y mod p
我们可以把每次的乘积存下来，最后一起累乘，即
2^10
= 2^5 * 2^5
那么每次把y>>=1，如果y是奇数那么就把这次分出的答案乘进答案
xy mod p = z mod p
求最小的x，其实很简单，可以将方程变成ax+by=c
同时除以gcd（a， b），如果c除不尽那么无解
否则变成a0x+b0y=c0;->(a0x+b0y=1;x，y乘以c0)

下面给一个例题；
给出下列三种操作
1、给定 y、z、p，计算 y^z mod p 的值；
2、给定 y、z、p，计算满足 xy ≡z(mod p)的最小非负整数 x；
3、给定 y、z、p，计算满足 y^x ≡z(mod p)的最小非负整数 x.
三个操作都在上面讲过了下面给出代码

#include <bits/stdc++.h>
using namespace std;

#define hash Hash
#define rep(i, s, t) for(int i = s; i <= t; ++i)

typedef long long ll;

ll p;
ll calc1(ll x, ll y) {
    ll res = 1LL;
    for(; y; y>>=1) {if(y & 1) res = res*x % p; x = x*x % p;}
    return res;
}

ll gcd(ll a, ll b) {
    while(b) {
        ll t = a;
        a = b;
        b = t%b;
    }
    return a;
}

void exgcd(ll a, ll b, ll &x, ll &y) {
    if(!b) {
        x = 1; y = 0;
    }
    else {
        exgcd(b, a%b, y, x);
        y -= a/b*x;
    }
}

struct HASH_TABLE {
    const ll mod = 100007;
    ll h[100007+10][2];

    ll find(ll x) {
        ll t = x % mod;
        while(h[t][0] != x && h[t][0] != -1) t = (t+1) % mod;
        return t;
    }

    void insert(ll x, ll id) {
        ll t = find(x);
        if(h[t][0] == -1)
            h[t][0] = x, h[t][1] = id;
    }

    ll get(ll x) {
        ll t = find(x);
        return h[t][1];
    }

    void BSGS(ll a, ll b, ll p) {
        memset(h, -1, sizeof(h));
        ll m = ceil(sqrt(p)), t = 1;
        rep(i, 0, m-1) {
            insert(t, i);
            t = t*a%p;
        }
        ll d=1, x, y;
        ll ans = -1;
        rep(i, 0, m-1) {
            exgcd(d, p, x, y);
            x = ((x*b) % p + p)%p;
            y = get(x);
            //printf("%lld,%lld
",x,y);
            if(y ^ -1) {
                ans = i*m+y;
                break;
            }
            d = d*t % p;
        }
        if(ans == -1) puts("Orz, I cannot find x!");
        else printf("%lld
", ans);
    }
}T;

int main() {
#ifndef ONLINE_JUDGE
    freopen("calc.in", "r", stdin);
    freopen("calc.out", "w", stdout);
#endif
    int _, cmd;
    scanf("%d%d", &_, &cmd);
    while(_--) {
        ll y, z, a, b;
        scanf("%lld%lld%lld", &y, &z, &p);
        if(cmd == 1) cout << calc1(y, z) << endl;
        else if(cmd == 2) {
            ll Gcd = gcd(y, p);
            if(z % Gcd)
                puts("Orz, I cannot find x!");
            else {
                z %= p;
                y /= Gcd; p /= Gcd; z /= Gcd;
                exgcd(y, p, a, b);
                a = (a * z % p + p) % p;
                printf("%lld
", a);
            }
        }
        else {
            if(y%p==0) {
                if(z==1) printf("0
"); else puts("Orz, I cannot find x!");
                continue;
            }
            T.BSGS(y, z, p);
        }
    }
    return 0;
}