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

链接

BZOJ 2480

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

题解

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

先讲一下题解吧。

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

[A^xequiv Bpmod C ]

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

考虑优化——BSGS算法。

将所有要枚举的数分成(n = lfloor sqrt C floor)块，每块有(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 mod C)加入哈希表，然后(O(n))枚举(i)得到(A^{i * m} mod C)，如果在哈希表中能找到它对应的(y)，则(i * m - y)就是一个答案。

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

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

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

下文为了方便，搞了个(A^xequiv B_1pmod {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 mod g_1 ot= 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('
')
#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;
}