BZOJ 2480 && 3239 && 2995 高次不定方程（高次同余方程）

链接

BZOJ 2480

虽然是个三倍经验题（2333），但是只有上面这道（BZOJ2480）有 p = 1 的加强数据，推荐大家做这道。

题解

这是一道BSGS（Baby Step Giant Step）的棵题，要考虑的细节还是很多的……

先讲一下题解吧。

由于每个版本的题设的字母都不尽相同，所以……在下文中我使用的方程是

\[A^x\equiv B\pmod C \]

考虑暴力枚举，可以证明，\(A^x \bmod C\)是周期性的，且周期长度不超过\(C\)。那么暴力枚举\([0, C - 1]\)的整数就可以得到答案。

考虑优化——BSGS算法。

将所有要枚举的数分成\(n = \lfloor \sqrt C \rfloor\)块，每块有\(m\)个数。

如果我们按块枚举，在第\(i\)块中的倒数第\(y\)个位置找到了答案，则有：

\[A^{i * m - y} \equiv B \pmod C \]

可以把\(A^{-y}\)移到等式右边，得：

\[A^{i * m} \equiv A^y*B \pmod C \]

那么可以\(O(m)\)枚举\(y \in [1, m]\)，将得到的所有\(A^y*B \bmod C\)加入哈希表，然后\(O(n)\)枚举\(i\)得到\(A^{i * m} \bmod C\)，如果在哈希表中能找到它对应的\(y\)，则\(i * m - y\)就是一个答案。

听起来没有问题？
很好写？

其实还有一个问题……
上面的叙述其实是默认\(A, C\)互质的……
而当\(A, C\)不互质，\(A^{-1} \pmod C\)即A关于C的逆元是不存在的，所以不能吧\(A^{-y}\)看作一个数移到等式右边！

那怎么办呢……？
必须强行让\(A, C\)互质了！

下文为了方便，搞了个\(A^x\equiv B_1\pmod {C_1}\)的等价式子$$A^x + C_1 * y = B_1$$

把等式两边同时除以\(g_1 = gcd(A, C)\)：

\[\frac{A}{g_1} A^{x - 1} + C_2 * y = B_2 \]

其中\(C_2 = \frac{C_1}{g_1}, B_2 = \frac{B_1}{g_1}\)。

显然，如果\(B_1 \bmod g_1 \not= 0\)，则无解。

这样是不是就完成了呢？

不是！

此时\(A\)和\(C_2\)不一定互质，例如\(A = 12, C_1 = 9, C_2 = 3\)，\(12\)和\(3\)就不互质。

那么我们要不断地往下除\(gcd(A, C)\)，直到二者互质，这个操作的次数是\(O(\log n)\)的。

最后，我们得到的方程像这样：

\[\frac{A^{cnt}}{d_1d_2d_3...d_n}A^{x - nct} \equiv B_n \pmod {C_n} \]

设\(D = \frac{A^{cnt}}{d_1d_2d_3...d_n}\)，则把\(D\)移到等式右边：

\[A^{x - cnt} \equiv D^{-1}B_n \pmod {C_n} \]

此时\(A\)和\(C_n\)是互质的，解这个方程就好了。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <queue>
using namespace std;
typedef long long ll;
#define enter putchar('\n')
#define space putchar(' ')
template <class T>
void read(T &x){
    char c;
    bool op = 0;
    while(c = getchar(), c > '9' || c < '0')
	if(c == '-') op = 1;
    x = c - '0';
    while(c = getchar(), c >= '0' && c <= '9')
	x = x * 10 + c - '0';
    if(op) x = -x;
}
template <class T>
void write(T x){
    if(x < 0) putchar('-'), x = -x;
    if(x >= 10) write(x / 10);
    putchar('0' + x % 10);
}
const int P = 999979, N = 100005; //P是哈希表用的质数
ll A, B, C, cnt, D;
ll adj[P], nxt[N], num[N], val[N], stk[N], top;
ll gcd(ll a, ll b){
    return b ? gcd(b, a % b) : a;
}
void exgcd(ll a, ll b, ll &x, ll &y){
    if(!b) return (void)(x = 1, y = 0);
    exgcd(b, a % b, y, x);
    y -= a / b * x;
}
ll inv(ll a, ll p){ //求逆元
    ll x, y;
    exgcd(a, p, x, y);
    return (x % p + p) % p;
}
void clear(){ //高效率清空哈希表
    while(top) adj[stk[top--]] = 0;
}
ll find(ll x){ //找到哈希表中val最大（即最后添加）、num == x的元素的val
    for(ll i = adj[x % P]; i; i = nxt[i])
        if(num[i] == x) return val[i];
    return -1;
}
void insert(ll x, ll y){ //向哈希表插入一个元素，num = x，val = y
    stk[++top] = x % P;
    nxt[top] = adj[stk[top]];
    adj[stk[top]] = top;
    num[top] = x;
    val[top] = y;
}
bool check(){ //使AC互质，同时排除部分无解情况（B % gcd(A, C) != 0）
    if(C == 0 || (A == 0 && B != 0)) return 0;
    cnt = 0, D = 1;
    for(ll g = gcd(A, C); g != 1; g = gcd(A, C))
        if(B % g) return 0;
        else cnt++, B /= g, C /= g, D = D * A / g % C;
    B = B * inv(D, C) % C;
    return 1;
}
bool force(){ //为了排除解小于cnt的情况，先进行小范围暴力
    ll sum = 1 % C; //一个神犇给BZOJ2480加了一组 C == 1 的数据！
    for(int i = 0; i <= 30; i++){
        if(sum == B) return write(i), enter, 1;
        sum = sum * A % C;
    }
    return 0;
}
int main(){
    while(read(A), read(C), read(B), A + B + C != 0){
        A %= C, B %= C;
        if(force()) continue;
        else if(!check()) puts("No Solution");
        else{
            clear();
            bool solved = 0;
            ll sum = 1, n = sqrt(C), m = ceil((double)C / n); //分成n块，每块m个
            for(ll i = 1; i <= m; i++){ //将 pow(A, i) * B 插入哈希表
                sum = sum * A % C;
                insert(sum * B % C, i);
            }
            for(ll i = 1, tot = 1, y; i <= n && !solved; i++){ //检查每个 pow(A, i * m) 是否在哈希表中
                tot = tot * sum % C;
                if((y = find(tot)) != -1){
                    write(i * m - y + cnt), enter;
                    solved = 1;
                }
            }
            if(!solved) puts("No Solution");
        }
    }
    return 0;
}