This documentation is automatically generated by online-judge-tools/verification-helper
View the Project on GitHub hitonanode/cplib-cpp
#define PROBLEM "https://judge.yosupo.jp/problem/primitive_root" #include "../primitive_root.hpp" #include <iostream> using namespace std; int main() { cin.tie(nullptr), ios::sync_with_stdio(false); int Q; cin >> Q; while (Q--) { long long p; cin >> p; cout << find_smallest_primitive_root(p) << '\n'; } }
#line 1 "number/test/primitive_root_1e18.test.cpp" #define PROBLEM "https://judge.yosupo.jp/problem/primitive_root" #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; } #line 3 "number/test/primitive_root_1e18.test.cpp" #include <iostream> using namespace std; int main() { cin.tie(nullptr), ios::sync_with_stdio(false); int Q; cin >> Q; while (Q--) { long long p; cin >> p; cout << find_smallest_primitive_root(p) << '\n'; } }