cplib-cpp
This documentation is automatically generated by online-judge-tools/verification-helper
number/test/miller-rabin-5e7.test.cpp
- View this file on GitHub
- Last update: 2023-03-10 18:40:39+09:00
- Problem: https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A
Depends on
- Integer factorization （素因数分解）
(number/factorize.hpp)
- Linear sieve （線形篩）
(number/sieve.hpp)
- random/xorshift.hpp
Code
#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A" // DUMMY
#include "../factorize.hpp"
#include "../sieve.hpp"
#include <cassert>
#include <cstdio>
int main() {
int lim = 5e7;
Sieve sieve(lim);
for (int x = 1; x <= lim; x++) {
bool is_prime_1 = (sieve.min_factor[x] == x);
bool is_prime_2 = is_prime(x);
assert(is_prime_1 == is_prime_2);
}
puts("Hello World");
}#line 1 "number/test/miller-rabin-5e7.test.cpp"
#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A" // DUMMY
#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 3 "number/factorize.hpp"
#include <algorithm>
#include <array>
#include <cassert>
#include <numeric>
#include <vector>
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 3 "number/sieve.hpp"
#include <map>
#line 5 "number/sieve.hpp"
// CUT begin
// Linear sieve algorithm for fast prime factorization
// Complexity: O(N) time, O(N) space:
// - MAXN = 10^7: ~44 MB, 80~100 ms (Codeforces / AtCoder GCC, C++17)
// - MAXN = 10^8: ~435 MB, 810~980 ms (Codeforces / AtCoder GCC, C++17)
// Reference:
// [1] D. Gries, J. Misra, "A Linear Sieve Algorithm for Finding Prime Numbers,"
// Communications of the ACM, 21(12), 999-1003, 1978.
// - https://cp-algorithms.com/algebra/prime-sieve-linear.html
// - https://37zigen.com/linear-sieve/
struct Sieve {
std::vector<int> min_factor;
std::vector<int> primes;
Sieve(int MAXN) : min_factor(MAXN + 1) {
for (int d = 2; d <= MAXN; d++) {
if (!min_factor[d]) {
min_factor[d] = d;
primes.emplace_back(d);
}
for (const auto &p : primes) {
if (p > min_factor[d] or d * p > MAXN) break;
min_factor[d * p] = p;
}
}
}
// Prime factorization for 1 <= x <= MAXN^2
// Complexity: O(log x) (x <= MAXN)
// O(MAXN / log MAXN) (MAXN < x <= MAXN^2)
template <class T> std::map<T, int> factorize(T x) const {
std::map<T, int> ret;
assert(x > 0 and
x <= ((long long)min_factor.size() - 1) * ((long long)min_factor.size() - 1));
for (const auto &p : primes) {
if (x < T(min_factor.size())) break;
while (!(x % p)) x /= p, ret[p]++;
}
if (x >= T(min_factor.size())) ret[x]++, x = 1;
while (x > 1) ret[min_factor[x]]++, x /= min_factor[x];
return ret;
}
// Enumerate divisors of 1 <= x <= MAXN^2
// Be careful of highly composite numbers https://oeis.org/A002182/list
// https://gist.github.com/dario2994/fb4713f252ca86c1254d#file-list-txt (n, (# of div. of n)):
// 45360->100, 735134400(<1e9)->1344, 963761198400(<1e12)->6720
template <class T> std::vector<T> divisors(T x) const {
std::vector<T> ret{1};
for (const auto p : factorize(x)) {
int n = ret.size();
for (int i = 0; i < n; i++) {
for (T a = 1, d = 1; d <= p.second; d++) {
a *= p.first;
ret.push_back(ret[i] * a);
}
}
}
return ret; // NOT sorted
}
// Euler phi functions of divisors of given x
// Verified: ABC212 G https://atcoder.jp/contests/abc212/tasks/abc212_g
// Complexity: O(sqrt(x) + d(x))
template <class T> std::map<T, T> euler_of_divisors(T x) const {
assert(x >= 1);
std::map<T, T> ret;
ret[1] = 1;
std::vector<T> divs{1};
for (auto p : factorize(x)) {
int n = ret.size();
for (int i = 0; i < n; i++) {
ret[divs[i] * p.first] = ret[divs[i]] * (p.first - 1);
divs.push_back(divs[i] * p.first);
for (T a = divs[i] * p.first, d = 1; d < p.second; a *= p.first, d++) {
ret[a * p.first] = ret[a] * p.first;
divs.push_back(a * p.first);
}
}
}
return ret;
}
// Moebius function Table, (-1)^{# of different prime factors} for square-free x
// return: [0=>0, 1=>1, 2=>-1, 3=>-1, 4=>0, 5=>-1, 6=>1, 7=>-1, 8=>0, ...] https://oeis.org/A008683
std::vector<int> GenerateMoebiusFunctionTable() const {
std::vector<int> ret(min_factor.size());
for (unsigned i = 1; i < min_factor.size(); i++) {
if (i == 1) {
ret[i] = 1;
} else if ((i / min_factor[i]) % min_factor[i] == 0) {
ret[i] = 0;
} else {
ret[i] = -ret[i / min_factor[i]];
}
}
return ret;
}
// Calculate [0^K, 1^K, ..., nmax^K] in O(nmax)
// Note: **0^0 == 1**
template <class MODINT> std::vector<MODINT> enumerate_kth_pows(long long K, int nmax) const {
assert(nmax < int(min_factor.size()));
assert(K >= 0);
if (K == 0) return std::vector<MODINT>(nmax + 1, 1);
std::vector<MODINT> ret(nmax + 1);
ret[0] = 0, ret[1] = 1;
for (int n = 2; n <= nmax; n++) {
if (min_factor[n] == n) {
ret[n] = MODINT(n).pow(K);
} else {
ret[n] = ret[n / min_factor[n]] * ret[min_factor[n]];
}
}
return ret;
}
};
// Sieve sieve((1 << 20));
#line 5 "number/test/miller-rabin-5e7.test.cpp"
#include <cstdio>
int main() {
int lim = 5e7;
Sieve sieve(lim);
for (int x = 1; x <= lim; x++) {
bool is_prime_1 = (sieve.min_factor[x] == x);
bool is_prime_2 = is_prime(x);
assert(is_prime_1 == is_prime_2);
}
puts("Hello World");
}