cplib-cpp

This documentation is automatically generated by online-judge-tools/verification-helper

View the Project on GitHub hitonanode/cplib-cpp

:heavy_check_mark: Primitive root modulo $n$ （原始根の発見）
(number/primitive_root.hpp)

与えられた整数 $n$ に対して，乗法群 $\mathbb{Z}^{\times}_{n}$ に原始根が存在する場合に最小の原始根を計算する．$n \le 10^{18}$ を想定．

乗法群の原始根とは

$\mod n$ 上の整数のうち，乗法の逆元が存在するもののみをとって作る群を乗法群 $\mathbb{Z}^{\times}_{n}$ と呼ぶ．

乗法群 $\mathbb{Z}^{\times}_{n}$ の元 $r$ について，$r^0, r^1, r^2, \dots$ の集合が $\mathbb{Z}^{\times}_{n}$ と（集合として）一致するとき，$r$ を乗法群の原始根と呼ぶ．

$\mathbb{Z}^\times_{n}$ に原始根が存在するような $n$ は限られている．具体的には，$n = 2, 4, p^k, 2p^k$ のみである（$p$ は奇素数，$k$ は正の整数）．

最小の原始根について

本コードの

    for (long long g = 1; g < n; g++) {

の部分について，g をインクリメンタルに調べることで最小の原始根を見つけている．$n$ が素数でかつ $10^{16}$ 以下の範囲では， $n = 6525032504501281$ のとき最小の原始根が $417$ となってこれが最大という（参考文献 [2]）．

上記のようなケースに対して，$g$ について乱択すれば（「最小の」原始根を諦める代わりに）実用上の高速化が可能と思われる [3]．

問題例

参考文献・リンク

Depends on

Verified with

Code

#pragma once
#include "bare_mod_algebra.hpp"
#include "factorize.hpp"

// Find smallest primitive root for given number n （最小の原始根探索）
// n must be 2 / 4 / p^k / 2p^k (p: odd prime, k > 0)
// (https://en.wikipedia.org/wiki/Primitive_root_modulo_n)
//
// Complexity: maybe O(sqrt(n)), if n is
// prime Algorithm: http://kirika-comp.hatenablog.com/entry/2018/03/12/210446 Verified:
// - https://yukicoder.me/submissions/405938
// - https://judge.yosupo.jp/problem/primitive_root
// - SRM 840 Div.1 900 https://community.topcoder.com/stat?c=problem_statement&pm=17877
// Sample:
//  - 998244353 ( = (119 << 23) + 1 ) -> 3
//  - 163577857 ( = (39 << 22) + 1 ) -> 23
//  - 2 -> 1
//  - 1 -> -1
long long find_smallest_primitive_root(long long n) {
    std::vector<long long> fac;
    const long long phi = FactorizeLonglong.euler_phi(n);
    for (long long q : FactorizeLonglong(phi)) {
        if (fac.empty() or fac.back() != q) fac.push_back(q);
    }

    for (long long g = 1; g < n; g++) {
        if (std::__gcd(n, g) != 1) continue;
        if (pow_mod<long long, __int128>(g, phi, n) != 1) return -1;
        bool ok = true;
        for (auto pp : fac) {
            if (pow_mod<long long, __int128>(g, phi / pp, n) == 1) {
                ok = false;
                break;
            }
        }
        if (ok) return g;
    }
    return -1;
}
#line 2 "number/bare_mod_algebra.hpp"
#include <algorithm>
#include <cassert>
#include <tuple>
#include <utility>
#include <vector>

// CUT begin
// Solve ax+by=gcd(a, b)
template <class Int> Int extgcd(Int a, Int b, Int &x, Int &y) {
    Int d = a;
    if (b != 0) {
        d = extgcd(b, a % b, y, x), y -= (a / b) * x;
    } else {
        x = 1, y = 0;
    }
    return d;
}
// Calculate a^(-1) (MOD m) s if gcd(a, m) == 1
// Calculate x s.t. ax == gcd(a, m) MOD m
template <class Int> Int mod_inverse(Int a, Int m) {
    Int x, y;
    extgcd<Int>(a, m, x, y);
    x %= m;
    return x + (x < 0) * m;
}

// Require: 1 <= b
// return: (g, x) s.t. g = gcd(a, b), xa = g MOD b, 0 <= x < b/g
template <class Int> /* constexpr */ std::pair<Int, Int> inv_gcd(Int a, Int b) {
    a %= b;
    if (a < 0) a += b;
    if (a == 0) return {b, 0};
    Int s = b, t = a, m0 = 0, m1 = 1;
    while (t) {
        Int u = s / t;
        s -= t * u, m0 -= m1 * u;
        auto tmp = s;
        s = t, t = tmp, tmp = m0, m0 = m1, m1 = tmp;
    }
    if (m0 < 0) m0 += b / s;
    return {s, m0};
}

template <class Int>
/* constexpr */ std::pair<Int, Int> crt(const std::vector<Int> &r, const std::vector<Int> &m) {
    assert(r.size() == m.size());
    int n = int(r.size());
    // Contracts: 0 <= r0 < m0
    Int r0 = 0, m0 = 1;
    for (int i = 0; i < n; i++) {
        assert(1 <= m[i]);
        Int r1 = r[i] % m[i], m1 = m[i];
        if (r1 < 0) r1 += m1;
        if (m0 < m1) {
            std::swap(r0, r1);
            std::swap(m0, m1);
        }
        if (m0 % m1 == 0) {
            if (r0 % m1 != r1) return {0, 0};
            continue;
        }
        Int g, im;
        std::tie(g, im) = inv_gcd<Int>(m0, m1);

        Int u1 = m1 / g;
        if ((r1 - r0) % g) return {0, 0};

        Int x = (r1 - r0) / g % u1 * im % u1;
        r0 += x * m0;
        m0 *= u1;
        if (r0 < 0) r0 += m0;
    }
    return {r0, m0};
}

// 蟻本 P.262
// 中国剰余定理を利用して，色々な素数で割った余りから元の値を復元
// 連立線形合同式 A * x = B mod M の解
// Requirement: M[i] > 0
// Output: x = first MOD second (if solution exists), (0, 0) (otherwise)
template <class Int>
std::pair<Int, Int>
linear_congruence(const std::vector<Int> &A, const std::vector<Int> &B, const std::vector<Int> &M) {
    Int r = 0, m = 1;
    assert(A.size() == M.size());
    assert(B.size() == M.size());
    for (int i = 0; i < (int)A.size(); i++) {
        assert(M[i] > 0);
        const Int ai = A[i] % M[i];
        Int a = ai * m, b = B[i] - ai * r, d = std::__gcd(M[i], a);
        if (b % d != 0) {
            return std::make_pair(0, 0); // 解なし
        }
        Int t = b / d * mod_inverse<Int>(a / d, M[i] / d) % (M[i] / d);
        r += m * t;
        m *= M[i] / d;
    }
    return std::make_pair((r < 0 ? r + m : r), m);
}

template <class Int = int, class Long = long long> Int pow_mod(Int x, long long n, Int md) {
    static_assert(sizeof(Int) * 2 <= sizeof(Long), "Watch out for overflow");
    if (md == 1) return 0;
    Int ans = 1;
    while (n > 0) {
        if (n & 1) ans = (Long)ans * x % md;
        x = (Long)x * x % md;
        n >>= 1;
    }
    return ans;
}
#line 2 "random/xorshift.hpp"
#include <cstdint>

// CUT begin
uint32_t rand_int() // XorShift random integer generator
{
    static uint32_t x = 123456789, y = 362436069, z = 521288629, w = 88675123;
    uint32_t t = x ^ (x << 11);
    x = y;
    y = z;
    z = w;
    return w = (w ^ (w >> 19)) ^ (t ^ (t >> 8));
}
double rand_double() { return (double)rand_int() / UINT32_MAX; }
#line 4 "number/factorize.hpp"
#include <array>
#line 6 "number/factorize.hpp"
#include <numeric>
#line 8 "number/factorize.hpp"

namespace SPRP {
// http://miller-rabin.appspot.com/
const std::vector<std::vector<__int128>> bases{
    {126401071349994536},                              // < 291831
    {336781006125, 9639812373923155},                  // < 1050535501 (1e9)
    {2, 2570940, 211991001, 3749873356},               // < 47636622961201 (4e13)
    {2, 325, 9375, 28178, 450775, 9780504, 1795265022} // <= 2^64
};
inline int get_id(long long n) {
    if (n < 291831) {
        return 0;
    } else if (n < 1050535501) {
        return 1;
    } else if (n < 47636622961201)
        return 2;
    else { return 3; }
}
} // namespace SPRP

// Miller-Rabin primality test
// https://ja.wikipedia.org/wiki/%E3%83%9F%E3%83%A9%E3%83%BC%E2%80%93%E3%83%A9%E3%83%93%E3%83%B3%E7%B4%A0%E6%95%B0%E5%88%A4%E5%AE%9A%E6%B3%95
// Complexity: O(lg n) per query
struct {
    long long modpow(__int128 x, __int128 n, long long mod) noexcept {
        __int128 ret = 1;
        for (x %= mod; n; x = x * x % mod, n >>= 1) ret = (n & 1) ? ret * x % mod : ret;
        return ret;
    }
    bool operator()(long long n) noexcept {
        if (n < 2) return false;
        if (n % 2 == 0) return n == 2;
        int s = __builtin_ctzll(n - 1);

        for (__int128 a : SPRP::bases[SPRP::get_id(n)]) {
            if (a % n == 0) continue;
            a = modpow(a, (n - 1) >> s, n);
            bool may_composite = true;
            if (a == 1) continue;
            for (int r = s; r--; a = a * a % n) {
                if (a == n - 1) may_composite = false;
            }
            if (may_composite) return false;
        }
        return true;
    }
} is_prime;

struct {
    // Pollard's rho algorithm: find factor greater than 1
    long long find_factor(long long n) {
        assert(n > 1);
        if (n % 2 == 0) return 2;
        if (is_prime(n)) return n;
        long long c = 1;
        auto f = [&](__int128 x) -> long long { return (x * x + c) % n; };

        for (int t = 1;; t++) {
            for (c = 0; c == 0 or c + 2 == n;) c = rand_int() % n;
            long long x0 = t, m = std::max(n >> 3, 1LL), x, ys, y = x0, r = 1, g, q = 1;
            do {
                x = y;
                for (int i = r; i--;) y = f(y);
                long long k = 0;
                do {
                    ys = y;
                    for (int i = std::min(m, r - k); i--;)
                        y = f(y), q = __int128(q) * std::abs(x - y) % n;
                    g = std::__gcd<long long>(q, n);
                    k += m;
                } while (k < r and g <= 1);
                r <<= 1;
            } while (g <= 1);
            if (g == n) {
                do {
                    ys = f(ys);
                    g = std::__gcd(std::abs(x - ys), n);
                } while (g <= 1);
            }
            if (g != n) return g;
        }
    }

    std::vector<long long> operator()(long long n) {
        std::vector<long long> ret;
        while (n > 1) {
            long long f = find_factor(n);
            if (f < n) {
                auto tmp = operator()(f);
                ret.insert(ret.end(), tmp.begin(), tmp.end());
            } else
                ret.push_back(n);
            n /= f;
        }
        std::sort(ret.begin(), ret.end());
        return ret;
    }
    long long euler_phi(long long n) {
        long long ret = 1, last = -1;
        for (auto p : this->operator()(n)) ret *= p - (last != p), last = p;
        return ret;
    }
} FactorizeLonglong;
#line 4 "number/primitive_root.hpp"

// Find smallest primitive root for given number n （最小の原始根探索）
// n must be 2 / 4 / p^k / 2p^k (p: odd prime, k > 0)
// (https://en.wikipedia.org/wiki/Primitive_root_modulo_n)
//
// Complexity: maybe O(sqrt(n)), if n is
// prime Algorithm: http://kirika-comp.hatenablog.com/entry/2018/03/12/210446 Verified:
// - https://yukicoder.me/submissions/405938
// - https://judge.yosupo.jp/problem/primitive_root
// - SRM 840 Div.1 900 https://community.topcoder.com/stat?c=problem_statement&pm=17877
// Sample:
//  - 998244353 ( = (119 << 23) + 1 ) -> 3
//  - 163577857 ( = (39 << 22) + 1 ) -> 23
//  - 2 -> 1
//  - 1 -> -1
long long find_smallest_primitive_root(long long n) {
    std::vector<long long> fac;
    const long long phi = FactorizeLonglong.euler_phi(n);
    for (long long q : FactorizeLonglong(phi)) {
        if (fac.empty() or fac.back() != q) fac.push_back(q);
    }

    for (long long g = 1; g < n; g++) {
        if (std::__gcd(n, g) != 1) continue;
        if (pow_mod<long long, __int128>(g, phi, n) != 1) return -1;
        bool ok = true;
        for (auto pp : fac) {
            if (pow_mod<long long, __int128>(g, phi / pp, n) == 1) {
                ok = false;
                break;
            }
        }
        if (ok) return g;
    }
    return -1;
}
Back to top page