spl

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

View the Project on GitHub Forestedf/spl

:heavy_check_mark: number_theory/primitive_root.hpp

Depends on

Required by

Verified with

Code

#pragma once
#include "factorize.hpp"

// p: prime
unsigned long long find_primitive_root(unsigned long long p) {
    using M = MontgomeryModInt64<20250128>;
    assert(is_prime(p));
    if (p == 2) {
        return 1;
    }
    M::set_mod(p);
    std::vector<unsigned long long> ps = factorize(p - 1);
    ps.erase(std::unique(ps.begin(), ps.end()), ps.end());
    while (true) {
        unsigned long long x = uniform<unsigned long long>(1, p);
        M x_(x), one(1ULL);
        bool ok = true;
        for (unsigned long long q : ps) {
            if (x_.pow((p - 1) / q).x == one.x) {
                ok = false;
                break;
            }
        }
        if (ok) {
            return x;
        }
    }
    return 0;
}
#line 2 "number_theory/factorize.hpp"
#include <algorithm>
#include <vector>
#line 2 "template/random.hpp"
#include <chrono>
#include <random>

#if defined(LOCAL) || defined(FIX_SEED)
std::mt19937_64 mt(123456789);
#else
std::mt19937_64 mt(std::chrono::steady_clock::now().time_since_epoch().count());
#endif

template <typename T>
T uniform(T l, T r) {
    return std::uniform_int_distribution<T>(l, r - 1)(mt);
}
template <typename T>
T uniform(T n) {
    return std::uniform_int_distribution<T>(0, n - 1)(mt);
}
#line 2 "number_theory/primality.hpp"

namespace primality {

using u64 = unsigned long long;
using u128 = __uint128_t;

u64 inv_64(u64 n) {
    u64 r = n;
    for (int i = 0; i < 5; ++i) {
        r *= 2 - n * r;
    }
    return r;
}

// n: odd, < 2^{62}
struct Montgomery64 {
    u64 n, mni, p;
    Montgomery64(u64 n) : n(n), mni(-inv_64(n)), p(-1ULL % n + 1) {}
    u64 mulmr(u64 xr, u64 yr) const {
        u128 z = (u128)xr * yr;
        u64 ret = (z + (u128)((u64)z * mni) * n) >> 64;
        if (ret >= n) {
            ret -= n;
        }
        return ret;
    }
    u64 mr(u64 xr) const {
        u64 ret = (xr + (u128)(xr * mni) * n) >> 64;
        if (ret >= n) {
            ret -= n;
        }
        return ret;
    }
    u64 pow(u64 xr, u64 t) const {
        u64 ret = p;
        while (t) {
            if (t & 1) {
                ret = mulmr(ret, xr);
            }
            xr = mulmr(xr, xr);
            t >>= 1;
        }
        return ret;
    }
};

bool is_prime(u64 n) {
    if (n == 2) {
        return true;
    }
    if (n == 1 || n % 2 == 0) {
        return false;
    }
    u64 s = __builtin_ctzll(n - 1);
    u64 d = (n - 1) >> s;
    u64 base[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
    Montgomery64 mont(n);
    u64 fl = n - mont.p;
    for (u64 b : base) {
        b = mont.mr(b);
        if (!b) {
            continue;
        }
        u64 t = mont.pow(b, d);
        if (t == mont.p) {
            continue;
        }
        u64 i = 0;
        for (; i < s; ++i) {
            if (t == fl) {
                break;
            }
            t = mont.mulmr(t, t);
        }
        if (i == s) {
            return false;
        }
    }
    return true;
}

}  // namespace primality

bool is_prime(unsigned long long n) {
    return primality::is_prime(n);
}
#line 2 "number_theory/montgomery_64.hpp"
#include <cassert>

// mod: odd, < 2^{63}
template <int id>
struct MontgomeryModInt64 {
    using u64 = unsigned long long;
    using u128 = __uint128_t;

    static u64 inv_64(u64 n) {
        u64 r = n;
        for (int i = 0; i < 5; ++i) {
            r *= 2 - n * r;
        }
        return r;
    }

    static u64 mod, neg_inv, sq;
    static void set_mod(u64 m) {
        assert(m % 2 == 1 && m < (1ULL << 63));
        mod = m;
        neg_inv = -inv_64(m);
        sq = -u128(mod) % mod;
    }
    static u64 get_mod() { return mod; }

    static u64 reduce(u128 xr) {
        u64 ret = (xr + u128(u64(xr) * neg_inv) * mod) >> 64;
        if (ret >= mod) {
            ret -= mod;
        }
        return ret;
    }

    using M = MontgomeryModInt64<id>;

    u64 x;
    MontgomeryModInt64() : x(0) {}
    MontgomeryModInt64(u64 _x) : x(reduce(u128(_x) * sq)) {}

    u64 val() const { return reduce(u128(x)); }

    M &operator+=(M rhs) {
        if ((x += rhs.x) >= mod) {
            x -= mod;
        }
        return *this;
    }
    M &operator-=(M rhs) {
        if ((x -= rhs.x) >= mod) {
            x += mod;
        }
        return *this;
    }
    M &operator*=(M rhs) {
        x = reduce(u128(x) * rhs.x);
        return *this;
    }
    M operator+(M rhs) const { return M(*this) += rhs; }
    M operator-(M rhs) const { return M(*this) -= rhs; }
    M operator*(M rhs) const { return M(*this) *= rhs; }

    M pow(u64 t) const {
        M ret(1);
        M self = *this;
        while (t) {
            if (t & 1) {
                ret *= self;
            }
            self *= self;
            t >>= 1;
        }
        return ret;
    }
    M inv() const {
        assert(x);
        return this->pow(mod - 2);
    }

    M &operator/=(M rhs) {
        *this /= rhs.inv();
        return *this;
    }
    M operator/(M rhs) const { return M(*this) /= rhs; }
};

template <int id> unsigned long long MontgomeryModInt64<id>::mod = 1;
template <int id> unsigned long long MontgomeryModInt64<id>::neg_inv = 1;
template <int id> unsigned long long MontgomeryModInt64<id>::sq = 1;
#line 7 "number_theory/factorize.hpp"

namespace factorize_impl {

unsigned long long bgcd(unsigned long long x, unsigned long long y) {
    if (x == 0) {
        return y;
    }
    if (y == 0) {
        return x;
    }
    int n = __builtin_ctzll(x);
    int m = __builtin_ctzll(y);
    x >>= n;
    y >>= m;
    while (x != y) {
        if (x > y) {
            x = (x - y) >> __builtin_ctzll(x - y);
        } else {
            y = (y - x) >> __builtin_ctzll(y - x);
        }
    }
    return x << (n < m ? n : m);
}

template <typename T>
unsigned long long rho(unsigned long long n, unsigned long long c) {
    T cc(c);
    auto f = [cc](T x) -> T {
        return x * x + cc;
    };
    T y(2);
    T x = y;
    T z = y;
    T p(1);
    unsigned long long g = 1;
    constexpr int M = 128;
    for (int r = 1; g == 1; r *= 2) {
        x = y;
        for (int i = 0; i < r && g == 1; i += M) {
            z = y;
            for (int j = 0; j < r - i && j < M; ++j) {
                y = f(y);
                p *= y - x;
            }
            g = bgcd(p.val(), n);
        }
    }
    if (g == n) {
        do {
            z = f(z);
            g = bgcd((z - x).val(), n);
        } while (g == 1);
    }
    return g;
}

unsigned long long find_factor(unsigned long long n) {
    using M = MontgomeryModInt64<20250127>;
    M::set_mod(n);
    while (true) {
        unsigned long long c = uniform(n);
        unsigned long long g = rho<M>(n, c);
        if (g != n) {
            return g;
        }
    }
    return 0;
}

void factor_inner(unsigned long long n, std::vector<unsigned long long> &ps) {
    if (is_prime(n)) {
        ps.push_back(n);
        return;
    }
    if (n % 2 == 0) {
        ps.push_back(2);
        factor_inner(n / 2, ps);
        return;
    }
    unsigned long long m = find_factor(n);
    factor_inner(m, ps);
    factor_inner(n / m, ps);
}

}

std::vector<unsigned long long> factorize(unsigned long long n) {
    if (n <= 1) {
        return std::vector<unsigned long long>();
    }
    std::vector<unsigned long long> ps;
    factorize_impl::factor_inner(n, ps);
    std::sort(ps.begin(), ps.end());
    return ps;
}
#line 3 "number_theory/primitive_root.hpp"

// p: prime
unsigned long long find_primitive_root(unsigned long long p) {
    using M = MontgomeryModInt64<20250128>;
    assert(is_prime(p));
    if (p == 2) {
        return 1;
    }
    M::set_mod(p);
    std::vector<unsigned long long> ps = factorize(p - 1);
    ps.erase(std::unique(ps.begin(), ps.end()), ps.end());
    while (true) {
        unsigned long long x = uniform<unsigned long long>(1, p);
        M x_(x), one(1ULL);
        bool ok = true;
        for (unsigned long long q : ps) {
            if (x_.pow((p - 1) / q).x == one.x) {
                ok = false;
                break;
            }
        }
        if (ok) {
            return x;
        }
    }
    return 0;
}
Back to top page