【Gym103261D】FFT Algorithm (Carmichael function + Pollard Rho)

题目链接

题目大意

给定 $m,k(m\le 10^{18},15\le k\le 23)$ ，求一个数 $w$ 满足：

• $w^{2^k}\equiv 1(\mathrm{mod}m)$
• $w^{2^{k-1}}\not\equiv 1(\mathrm{mod}m)$

思路

Carmichael 函数 λ ( n ) = max ⁡ { o ( a ) ∣ a ⊥ n } \lambda(n)=\max\{o(a)\mid a\perp n\} λ(n)=max{o(a)∣a⊥n}

其中

• $\lambda (2)=1$
• $\lambda (4)=2$
• $\lambda (2^k)=2^{k-2}(k\ge 3)$
• $\lambda (p^k)=p^{k-1}(p-1)(p\ge 3)$
• $\lambda (uv)=\text{lcm}(\lambda (u),\lambda (v))(u\perp v)$

原根存在的充要条件是 $\lambda (n)=\phi (n)$ ，即 $n=1,2,4,p^k,2p^k(p\ge 3)$

答案存在当且仅当 $2^k|\lambda (m)$ 。

为了求 $\lambda (m)$ ，需要对 $m$ 进行质因数分解。可以采用 Pollard Rho 算法：

1. 利用 Miller-Rabin 判断本身是否是质数
2. 如果不是质数，随机一个序列 $x,f(x),f(f(x)),...$ ，其中 $f(x)=x^2+c$ （这样可以保证 $x_1-x_2|f(x_1)-f(x_2)$），序列最后会形成一个 $\rho$ 形。我们可以搞两个变量 $x_1,x_2$ ，初始都设成 $x$ ，然后每次迭代时令 $x_1=f(x_1),x_2=f(f(x_2))$ ，依次检验即可。如果遇到套圈的情况，就换个种子重新生成。（实际操作时，一般让 $x=2,c=1$ ，127步为一组判断。如果套圈就令 $c++$ ）期望复杂度 $O(n^{\frac{1}{4}})$

求出来后，随机一些 $x$ ，令 $w=x^{\frac{\lambda (m)}{2^k}}$ ，多试一些就能出答案。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
ll prime[12] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
ll mul(ll x, ll y, ll mod) {  // can also use __int128 on Linux
    ll ret = 0;
    while (y) {
        if (y & 1) ret = (ret + x) % mod;
        x = x * 2 % mod;
        y >>= 1;
    }
    return ret;
}
ll pw(ll x, ll y, ll mod) {
    ll ret = 1;
    while (y) {
        if (y & 1) ret = mul(ret, x, mod);
        x = mul(x, x, mod);
        y >>= 1;
    }
    return ret;
}
bool isprime(ll x) {
    if (x == 1) return false;
    for (int i = 0; i < 12; ++i) {
        if (x == prime[i]) return true;  // special(itself)
        ll t = x - 1;
        if (pw(prime[i], t, x) != 1) return false;  // failure
        while (t % 2 == 0) {
            t >>= 1;
            ll now = pw(prime[i], t, x);
            if (now == x - 1) break;
            if (now != 1) return false;  // failure
        }
    }
    return true;
}
ll gcd(ll x, ll y) { return y ? gcd(y, x % y) : x; }
ll lcm(ll x, ll y) { return x / gcd(x, y) * y; }
ll pollard_rho(ll x) {
    if (x % 2 == 0) return 2;  // special(unable to generate 2)
    ll x1 = 2, x2 = 2, c = 1;
    while (1) {
        ll y1 = x1, y2 = x2, y = 1;
        for (int i = 1; i <= 127; ++i) {
            y1 = (mul(y1, y1, x) + c) % x;
            y2 = (mul(y2, y2, x) + c) % x;
            y2 = (mul(y2, y2, x) + c) % x;
            y = mul(y, y2 - y1 + x,
                    x);  // to generate possible (times of factors)
        }
        y = gcd(y, x);
        if (y == x) {  // rejudge one by one
            for (int i = 1; i <= 127; ++i) {
                x1 = (mul(x1, x1, x) + c) % x;
                x2 = (mul(x2, x2, x) + c) % x;
                x2 = (mul(x2, x2, x) + c) % x;
                y = gcd(x2 - x1 + x, x);
                if (y == x)
                    ++c;
                else if (y > 1)
                    return y;
            }
        } else if (y > 1)
            return y;
        x1 = y1;
        x2 = y2;
    }
}
ll sta[1005], top = 0;
void decompose(ll x) {
    if (x == 1) return;
    if (isprime(x)) {
        sta[++top] = x;
        return;
    }
    ll t = pollard_rho(x);
    decompose(t);
    decompose(x / t);
}
ll m;
int k;
ll solve() {
    int cnt = 0;
    ll x = m;
    while (x % 2 == 0) {
        x >>= 1;
        ++cnt;
    }
    ll lambda;
    if (cnt <= 1)
        lambda = 1;
    else if (cnt == 2)
        lambda = 2;
    else
        lambda = (1ll << cnt - 2);
    top = 0;
    decompose(x);
    sort(sta + 1, sta + top + 1);
    for (int i = 1, last; i <= top; i = last + 1) {
        ll tmp = sta[i] - 1;
        last = i;
        while (last < top && sta[last + 1] == sta[i]) {
            ++last;
            tmp *= sta[i];
        }
        lambda = lcm(lambda, tmp);
    }
    // printf("%lld\n", lambda);
    if (lambda & ((1ll << k) - 1))
        return -1;
    else {
        for (ll x = 1; x < m; x++) {
            if (gcd(x, m) > 1) continue;
            ll w = pw(x, lambda >> k, m);
            if (pw(w, 1ll << k - 1, m) != 1) return w;
        }
    }
    return -2;
}
int main() {
    scanf("%lld%d", &m, &k);
    printf("%lld", solve());

    return 0;
}