number/primitive_root.hpp

• View this file on GitHub
• Last update: 2025-01-27 19:27:53+09:00
• Include: #include "number/primitive_root.hpp"

Depends on

• hash/splitmix64.hpp
• number/miller_rabin.hpp
• number/pollard_rho.hpp

Verified with

• test/lc_primitive_root.cpp

Code

#include "pollard_rho.hpp"

long long primitive_root(long long p) {
    if (p == 2) return 1;
    if (p == 3) return 2;
    if (p == 5) return 2;

    assert(is_prime(p));

    vector<long long> es;
    long long m = p - 1;
    if ((m & 1) == 0) {
        es.push_back((p - 1) / 2);
    }
    while ((m & 1) == 0) {
        m >>= 1;
    }
    while (m > 1) {
        long long pi = (is_prime(m) ? m : find_prime_factor(m));
        es.push_back((p - 1) / pi);
        while (m % pi == 0) {
            m /= pi;
        }
    }

    long long g = 2;
    for (long long i = 2;; ++i) {
        if (p > 1000) {
            g = splitmix64(i) % p;
            if (g == 0) continue;
        } else {
            g = i;
        }

        bool is = true;
        for (long long e : es) {
            __int128_t po = 1, ge = g;
            while (e) {
                if (e & 1) po = po * ge % p;
                ge = ge * ge % p;
                e >>= 1;
            }
            if (po == 1) {
                is = false;
                break;
            }
        }
        if (is) {
            return g;
        }
    }
}

#line 1 "number/miller_rabin.hpp"
bool miller_rabin(long long n, const vector<long long> &ts) {
    const long long s2 = (n - 1) & (1 - n), d = (n - 1) / s2;
    const auto pow_mod_n = [=](long long b, long long p) {
        long long res = 1;
        while (p > 0) {
            if (p & 1) res = __int128_t(res) * b % n;
            b = __int128_t(b) * b % n;
            p >>= 1;
        }
        return res;
    };
    for (const long long a : ts) {
        if (a >= n) break;
        auto ad = pow_mod_n(a, d);
        if (ad == 1) continue;
        long long t = 1;
        for (; t < s2; t <<= 1) {
            if (ad == n - 1) break;
            ad = __int128_t(ad) * ad % n;
        }
        if (t == s2) return false;
    }
    return true;
}
bool is_prime(long long n) {
    if (n == 1) return false;
    if (n == 2) return true;
    if ((n & 1) == 0) return false;
    if (n < 4759123141) return miller_rabin(n, {2, 7, 61});
    return miller_rabin(n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
}
#line 1 "hash/splitmix64.hpp"
constexpr uint64_t splitmix64(uint64_t x) {
    // http://xorshift.di.unimi.it/splitmix64.c
    x += 0x9e3779b97f4a7c15;
    x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
    x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
    return x ^ (x >> 31);
}
#line 3 "number/pollard_rho.hpp"

// n must not be prime
long long find_prime_factor(long long n) {
    constexpr int m = 200;
    for (long long c_i = 1;; ++c_i) {
        long long c = splitmix64(c_i) % n;
        if (c == 0) continue;

        auto f = [=](long long x) { return (__int128_t(x) * x + c) % n; };
        long long x = 0, y = 0, y_b = 0, q = 1;
        long long g = 1;
        int r = 1, k = 0;
        do {
            x = y;
            while (k < 3 * r / 4) {
                y = f(y);
                ++k;
            }
            while (k < r && g == 1) {
                y_b = y;
                for (int i = 0; i < min(m, r - k); ++i) {
                    y = f(y);
                    q = __int128_t(q) * abs(x - y) % n;
                }
                g = gcd(q, n);
                k += m;
            }
            k = r;
            r <<= 1;
        } while (g == 1);

        if (g == n) {
            g = 1;
            y = y_b;
            do {
                y = f(y);
                g = gcd(abs(x - y), n);
            } while (g == 1);
        }
        if (g == n) continue;

        if (is_prime(g)) return g;
        if (is_prime(n / g)) return n / g;
        return find_prime_factor(min(n / g, g));
    }
    return n;
}

vector<long long> factorize(long long n) {
    vector<long long> ps;
    if (n >= 2) {
        while ((n & 1) == 0) {
            ps.push_back(2);
            n >>= 1;
        }
        if (is_prime(n)) {
            ps.push_back(n);
        } else {
            while (n > 1) {
                long long p = find_prime_factor(n);
                do {
                    ps.push_back(p);
                    n /= p;
                } while (n % p == 0);

                if (is_prime(n)) {
                    ps.push_back(n);
                    break;
                }
            }
        }
    }
    sort(ps.begin(), ps.end());
    return ps;
}
#line 2 "number/primitive_root.hpp"

long long primitive_root(long long p) {
    if (p == 2) return 1;
    if (p == 3) return 2;
    if (p == 5) return 2;

    assert(is_prime(p));

    vector<long long> es;
    long long m = p - 1;
    if ((m & 1) == 0) {
        es.push_back((p - 1) / 2);
    }
    while ((m & 1) == 0) {
        m >>= 1;
    }
    while (m > 1) {
        long long pi = (is_prime(m) ? m : find_prime_factor(m));
        es.push_back((p - 1) / pi);
        while (m % pi == 0) {
            m /= pi;
        }
    }

    long long g = 2;
    for (long long i = 2;; ++i) {
        if (p > 1000) {
            g = splitmix64(i) % p;
            if (g == 0) continue;
        } else {
            g = i;
        }

        bool is = true;
        for (long long e : es) {
            __int128_t po = 1, ge = g;
            while (e) {
                if (e & 1) po = po * ge % p;
                ge = ge * ge % p;
                e >>= 1;
            }
            if (po == 1) {
                is = false;
                break;
            }
        }
        if (is) {
            return g;
        }
    }
}

Back to top page