N次剩余

问题引入

求解$x^a\equiv b(modm)$。

离散对数

离散对数是基于同余运算和原根的对数运算，当模$m$有原根时，设$g$为其中一个原根，则当$x\equiv g^k(modm)$时，可以对两边求对数，这样问题变成了$lo{g}_{g}x\equiv k(mod\phi (m))$。

当$m$为质数时

设$m$为质数$p$，由于$p$是质数，那么$p$一定存在一个原根$g$，根据质数原根的性质，因此对于模$p$意义下的任何数$x(0\le x<p)$有且仅有一个数$i(0\le i<p-1)$满足$x\equiv g^i$。

方法一

令$x=g^c$（$g$是$p$的一个原根且$g,c$一定存在），问题转化为求解$(g^c{)}^a\equiv b(modp)$，然后可以得到：

$(g^a{)}^c\equiv b$

这样问题就变成了一个朴素$BSGS$问题，可以在$O(\sqrt{p})$求出$c$，这样得到了原方程的一个特解$x_0\equiv g^c(modp)$

方法二

令$x=g^c$，并且设$b=g^t$，那么得到$g^{ac}\equiv g^t(modp)$，两边同时取离散对数得到：

$ac\equiv t(mod\phi (p))$

这时可以通过$BSGS$求解$g^t\equiv b(modp)$得到$t$，那么就转化成了线性同余方程组的问题，这样最后得到的也是一个特解$x_0\equiv g^c(modp)$

求出所有解

已知$g^{\phi (p)}\equiv 1(modp)$，那么可以得到：

$x^a \equiv g^{ac+\varphi§t}\equiv b(mod~~p), \forall t \in \Z $

于是所有的解为：

$x\equiv g^{c+\frac{\phi (p)t}{a}}(modp),\forall t\in \mathbb{Z},a|\phi (p)t$

显然有$\frac{a}{gcd(a,\phi (p))}|t$，设$t=\frac{a}{gcd(a,\phi (p))}i$，那么得到原问题的所有解为：

$x\equiv g^{c+\frac{\phi (p)i}{gcd(a,\phi (p))}}(modp),\forall i\in \mathbb{Z}$

代码

以方法一为例，极限数据范围，求出的是$[0,p)$范围内的所有解，即一共$p$个上述$i$的取值范围是$[0,gcd(a,\phi (p))$。

bitset<maxn> vis;
unordered_map<ll, ll> mp;
ll prime[maxn], fac[maxn], ans[maxn];
int cnt, tot;

void euler() {
    vis.reset(), cnt = 0;
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) prime[++cnt] = i;
        for (int j = 1; j <= cnt && 1LL * i * prime[j] < maxn; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
        }
    }
}

void getFac(ll n) {
    tot = 0;
    for (int i = 1; prime[i] * prime[i] <= n; i++) {
        if (n % prime[i] == 0) {
            fac[++tot] = prime[i];
            while (n % prime[i] == 0) n /= prime[i];
        }
    }
    if (n > 1) fac[++tot] = n;
}

ll slow_mul(ll a, ll b, ll p) {
    ll ans = 0;
    a %= p;
    while (b) {
        if (b & 1) ans = (ans + a) % p;
        a = (a + a) % p;
        b >>= 1;
    }
    return ans;
}

ll qkp(ll x, ll n, ll p) {
    ll ans = 1;
    x %= p;
    while (n) {
        if (n & 1) ans = slow_mul(ans, x, p);
        x = slow_mul(x, x, p);
        n >>= 1;
    }
    return ans;
}


inline bool check(ll g, ll p) {
    for (int i = 1; i <= tot; i++) {
        ll t = (p - 1) / fac[i];
        if (qkp(g, t, p) == 1) return 0;
    }
    return 1;
}

ll getRoot(ll p) {
    getFac(p - 1);
    for (ll g = 2; g < p; g++) {
        if (check(g, p)) return g;
    }
}

inline ll gcd(ll a, ll b) {
    return b == 0 ? a : gcd(b, a % b);
}

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

ll getInv(ll a, ll p) {
    return qkp(a, p - 2, p);
}

ll bsgs(ll a, ll b, ll p) {
    mp.clear();
    a %= p, b %= p;
    if (b == 1 || p == 1) return 0;
    int m = sqrt(p) + 1;
    ll cur = b;
    for (int i = 0; i < m; i++) {
        if (!mp.count(cur)) mp[cur] = i;
        cur = slow_mul(cur, a, p);
    }
    ll x = qkp(a, m, p), ans = -1, s = 1;
    for (int i = 1; i <= m; i++) {
        s = slow_mul(s, x, p);
        if (mp.count(s)) {
            ans = 1LL * i * m - mp[s];
            break;
        }
    }
    return ans;
}

bool solve(ll a, ll b, ll p) {
    ll g = getRoot(p);
    ll t = bsgs(g, b, p), phi = p - 1;
    ll d = gcd(a, phi);
    if (t < 0 || t % d) return 0;
    ll x, y;
    exgcd(a, phi, x, y);
    t /= d, phi /= d;
    x = (x * t % phi + phi) % phi;
    ans[0] = qkp(g, x, p);
    for (int i = 1; i < d; i++) {
        x += phi;
        ans[i] = qkp(g, x, p);
    }
    sort(ans, ans + d);
    for (int i = 0; i < d; i++) printf("%lld\n", ans[i]);
    return 1;
}

当$m$为合数时

求出$[0,m)$范围内的解，模板中$mod$即为$m$

ll A, B, mod;

ll pow(ll x, ll y, ll mod = 0, ll ans = 1) {
    if (mod) {
        for (; y; y >>= 1, x = x * x % mod)
            if (y & 1) ans = ans * x % mod;
    } else {
        for (; y; y >>= 1, x = x * x)
            if (y & 1) ans = ans * x;
    }
    return ans;
}

struct factor {
    ll prime[20], expo[20], pk[20], tot;

    void factor_lleger(ll n) {
        tot = 0;
        for (ll i = 2; i * i <= n; ++i)
            if (n % i == 0) {
                prime[tot] = i, expo[tot] = 0, pk[tot] = 1;
                do ++expo[tot], pk[tot] *= i; while ((n /= i) % i == 0);
                ++tot;
            }
        if (n > 1) prime[tot] = n, expo[tot] = 1, pk[tot++] = n;
    }

    ll phi(ll id) const {
        return pk[id] / prime[id] * (prime[id] - 1);
    }
} mods, _p;

ll p_inverse(ll x, ll id) {
    assert(x % mods.prime[id] != 0);
    return pow(x, mods.phi(id) - 1, mods.pk[id]);
}

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;
}

ll inverse(ll x, ll mod) {
    assert(__gcd(x, mod) == 1);
    ll ret, tmp;
    exgcd(x, mod, ret, tmp), ret %= mod;
    return ret + (ret >> 31 & mod);
}

vector<ll> sol[20];

void solve_2(ll id, ll k) {
    ll mod = 1 << k;
    if (k == 0) {
        sol[id].emplace_back(0);
        return;
    } else {
        solve_2(id, k - 1);
        vector<ll> t;
        for (ll s : sol[id]) {
            if (!((pow(s, A) ^ B) & mod - 1))
                t.emplace_back(s);
            if (!((pow(s | 1 << k - 1, A) ^ B) & mod - 1))
                t.emplace_back(s | 1 << k - 1);
        }
        swap(sol[id], t);
    }
}

ll BSGS(ll B, ll g, ll mod) {
    unordered_map<ll, ll> map;
    ll L = ceil(sqrt(mod)), t = 1;
    for (ll i = 1; i <= L; ++i) {
        t = t * g % mod;
        map[B * t % mod] = i;
    }
    ll now = 1;
    for (ll i = 1; i <= L; ++i) {
        now = now * t % mod;
        if (map.count(now)) return i * L - map[now];
    }
    assert(0);
}

ll find_primitive_root(ll id) {
    ll phi = mods.phi(id);
    _p.factor_lleger(phi);
    auto check = [&](ll g) {
        for (ll i = 0; i < _p.tot; ++i)
            if (pow(g, phi / _p.prime[i], mods.pk[id]) == 1)
                return 0;
        return 1;
    };
    for (ll g = 2; g < mods.pk[id]; ++g) if (check(g)) return g;
    assert(0);
}

void division(ll id, ll a, ll b, ll mod) { //  ax = b (mod M)
    ll M = mod, g = __gcd(__gcd(a, b), mod);
    a /= g, b /= g, mod /= g;
    if (__gcd(a, mod) > 1) return;
    ll t = b * inverse(a, mod) % mod;
    for (; t < M; t += mod) sol[id].emplace_back(t);
}

void solve_p(ll id, ll B = ::B) {
    ll p = mods.prime[id], e = mods.expo[id], mod = mods.pk[id];
    if (B % mod == 0) {
        ll q = pow(p, (e + A - 1) / A);
        for (ll t = 0; t < mods.pk[id]; t += q)
            sol[id].emplace_back(t);
    } else if (B % p != 0) {
        ll phi = mods.phi(id);
        ll g = find_primitive_root(id), z = BSGS(B, g, mod);
        division(id, A, z, phi);
        for (ll &x : sol[id]) x = pow(g, x, mod);
    } else {
        ll q = 0;
        while (B % p == 0) B /= p, ++q;
        ll pq = pow(p, q);
        if (q % A != 0) return;
        mods.expo[id] -= q, mods.pk[id] /= pq;
        solve_p(id, B);
        mods.expo[id] += q, mods.pk[id] *= pq;
        if (!sol[id].size()) return;

        ll s = pow(p, q - q / A);
        ll t = pow(p, q / A);
        ll u = pow(p, e - q);

        vector<ll> res;
        for (ll y : sol[id]) {
            for (ll i = 0; i < s; ++i)
                res.emplace_back((i * u + y) * t);
        }
        swap(sol[id], res);
    }
}

vector<ll> allans;

void dfs(ll dep, ll ans, ll mod) {
    if (dep == mods.tot) {
        allans.emplace_back(ans);
        return;
    }
    ll p = mods.pk[dep], k = p_inverse(mod % p, dep);
    for (ll a : sol[dep]) {
        ll nxt = (a - ans % p + p) * k % p * mod + ans;
        dfs(dep + 1, nxt, mod * p);
    }
}

int solve() {
    mods.factor_lleger(mod);
    allans.clear();
    for (ll i = mods.tot - 1; ~i; --i) {
        sol[i].clear();
        mods.prime[i] == 2 ? solve_2(i, mods.expo[i]) : solve_p(i);
        if (!sol[i].size()) return 0;
    }
    dfs(0, 0, 1);
    sort(allans.begin(), allans.end());
    cout << allans.size() << '\n';
    for (auto i : allans) cout << i << ' ';
    cout << '\n';
    return 1;
}