cplib-cpp
This documentation is automatically generated by online-judge-tools/verification-helper
number/test/sieve.stress.test.cpp
- View this file on GitHub
- Last update: 2023-08-05 18:05:47+09:00
- Problem: https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A
Depends on
Code
#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A" // DUMMY
#include "../sieve.hpp"
#include "../modint_runtime.hpp"
#include <algorithm>
#include <cassert>
#include <cstdio>
#include <iostream>
#include <vector>
using namespace std;
struct Case {
int SIEVE_SIZE;
int MAX;
};
int euler_phi(int x) {
int ret = 0;
for (int d = 1; d <= x; d++) ret += (std::__gcd(d, x) == 1);
return ret;
}
void test_divisors(Case testcase) {
const Sieve sieve(testcase.SIEVE_SIZE);
for (int x = 1; x <= testcase.MAX; x++) {
auto divs = sieve.divisors(x);
std::vector<int> is_div(x + 1);
for (auto d : divs) is_div.at(d) = 1;
for (int y = 1; y <= x; y++) assert(is_div.at(y) == (x % y == 0));
}
cerr << "divisors(): passed" << endl;
}
void test_euler_of_divisors(Case testcase) {
const Sieve sieve(testcase.SIEVE_SIZE);
for (int x = 1; x <= testcase.MAX; x++) {
auto div2euler = sieve.euler_of_divisors(x);
for (auto de : div2euler) {
assert(euler_phi(de.first) == de.second);
assert(x % de.first == 0);
}
assert(div2euler.size() == sieve.divisors(x).size());
}
cerr << "euler_of_divisors(): passed" << endl;
}
void test_moebius_table(int HI) {
const Sieve sieve_hi(HI);
const auto answer = sieve_hi.GenerateMoebiusFunctionTable();
assert(int(answer.size()) == HI + 1);
for (int x = 1; x <= HI; x++) {
int ret = 1;
for (auto p : sieve_hi.primes) {
if (x % (p * p) == 0) ret = 0;
if (x % p == 0) ret *= -1;
}
assert(answer[x] == ret);
}
for (int x = 1; x <= HI - 1; x++) {
auto mu = Sieve(x).GenerateMoebiusFunctionTable();
assert(int(mu.size()) == x + 1);
for (int i = 0; i < x + 1; i++) assert(mu[i] == answer[i]);
}
cerr << "GenerateMoebiusFunctionTable(): passed" << endl;
}
void test_enumerate_kth_pows(int hi_sieve_size) {
for (int sieve_size = 1; sieve_size <= hi_sieve_size; sieve_size++) {
const Sieve sieve(sieve_size);
for (int nmax = 1; nmax <= sieve_size; nmax++) {
for (int p = 1; p <= 20; p++) {
using mint = ModIntRuntime;
mint::set_mod(p);
for (int k = 0; k < 10; k++) {
auto ret = sieve.enumerate_kth_pows<mint>(k, nmax);
assert(int(ret.size()) == nmax + 1);
for (int i = 0; i <= nmax; i++) assert(ret[i] == mint(i).pow(k));
}
}
}
}
cerr << "enumerate_kth_pows(): passed" << endl;
}
int main() {
std::vector<Case> cases{{31, 31 * 31}, {100, 1000}, {2, 4}};
for (auto testcase : cases) {
test_divisors(testcase);
test_euler_of_divisors(testcase);
}
test_moebius_table(1000);
test_enumerate_kth_pows(100);
puts("Hello World");
}#line 1 "number/test/sieve.stress.test.cpp"
#define PROBLEM "https://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A" // DUMMY
#line 2 "number/sieve.hpp"
#include <cassert>
#include <map>
#include <vector>
// 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 2 "number/modint_runtime.hpp"
#include <iostream>
#include <set>
#line 5 "number/modint_runtime.hpp"
struct ModIntRuntime {
private:
static int md;
public:
using lint = long long;
static int mod() { return md; }
int val_;
static std::vector<ModIntRuntime> &facs() {
static std::vector<ModIntRuntime> facs_;
return facs_;
}
static int &get_primitive_root() {
static int primitive_root_ = 0;
if (!primitive_root_) {
primitive_root_ = [&]() {
std::set<int> fac;
int v = md - 1;
for (lint i = 2; i * i <= v; i++)
while (v % i == 0) fac.insert(i), v /= i;
if (v > 1) fac.insert(v);
for (int g = 1; g < md; g++) {
bool ok = true;
for (auto i : fac)
if (ModIntRuntime(g).power((md - 1) / i) == 1) {
ok = false;
break;
}
if (ok) return g;
}
return -1;
}();
}
return primitive_root_;
}
static void set_mod(const int &m) {
if (md != m) facs().clear();
md = m;
get_primitive_root() = 0;
}
ModIntRuntime &_setval(lint v) {
val_ = (v >= md ? v - md : v);
return *this;
}
int val() const noexcept { return val_; }
ModIntRuntime() : val_(0) {}
ModIntRuntime(lint v) { _setval(v % md + md); }
explicit operator bool() const { return val_ != 0; }
ModIntRuntime operator+(const ModIntRuntime &x) const {
return ModIntRuntime()._setval((lint)val_ + x.val_);
}
ModIntRuntime operator-(const ModIntRuntime &x) const {
return ModIntRuntime()._setval((lint)val_ - x.val_ + md);
}
ModIntRuntime operator*(const ModIntRuntime &x) const {
return ModIntRuntime()._setval((lint)val_ * x.val_ % md);
}
ModIntRuntime operator/(const ModIntRuntime &x) const {
return ModIntRuntime()._setval((lint)val_ * x.inv().val() % md);
}
ModIntRuntime operator-() const { return ModIntRuntime()._setval(md - val_); }
ModIntRuntime &operator+=(const ModIntRuntime &x) { return *this = *this + x; }
ModIntRuntime &operator-=(const ModIntRuntime &x) { return *this = *this - x; }
ModIntRuntime &operator*=(const ModIntRuntime &x) { return *this = *this * x; }
ModIntRuntime &operator/=(const ModIntRuntime &x) { return *this = *this / x; }
friend ModIntRuntime operator+(lint a, const ModIntRuntime &x) {
return ModIntRuntime()._setval(a % md + x.val_);
}
friend ModIntRuntime operator-(lint a, const ModIntRuntime &x) {
return ModIntRuntime()._setval(a % md - x.val_ + md);
}
friend ModIntRuntime operator*(lint a, const ModIntRuntime &x) {
return ModIntRuntime()._setval(a % md * x.val_ % md);
}
friend ModIntRuntime operator/(lint a, const ModIntRuntime &x) {
return ModIntRuntime()._setval(a % md * x.inv().val() % md);
}
bool operator==(const ModIntRuntime &x) const { return val_ == x.val_; }
bool operator!=(const ModIntRuntime &x) const { return val_ != x.val_; }
bool operator<(const ModIntRuntime &x) const {
return val_ < x.val_;
} // To use std::map<ModIntRuntime, T>
friend std::istream &operator>>(std::istream &is, ModIntRuntime &x) {
lint t;
return is >> t, x = ModIntRuntime(t), is;
}
friend std::ostream &operator<<(std::ostream &os, const ModIntRuntime &x) {
return os << x.val_;
}
lint power(lint n) const {
lint ans = 1, tmp = this->val_;
while (n) {
if (n & 1) ans = ans * tmp % md;
tmp = tmp * tmp % md;
n /= 2;
}
return ans;
}
ModIntRuntime pow(lint n) const { return power(n); }
ModIntRuntime inv() const { return this->pow(md - 2); }
ModIntRuntime fac() const {
int l0 = facs().size();
if (l0 > this->val_) return facs()[this->val_];
facs().resize(this->val_ + 1);
for (int i = l0; i <= this->val_; i++)
facs()[i] = (i == 0 ? ModIntRuntime(1) : facs()[i - 1] * ModIntRuntime(i));
return facs()[this->val_];
}
ModIntRuntime doublefac() const {
lint k = (this->val_ + 1) / 2;
return (this->val_ & 1)
? ModIntRuntime(k * 2).fac() / (ModIntRuntime(2).pow(k) * ModIntRuntime(k).fac())
: ModIntRuntime(k).fac() * ModIntRuntime(2).pow(k);
}
ModIntRuntime nCr(int r) const {
if (r < 0 or this->val_ < r) return ModIntRuntime(0);
return this->fac() / ((*this - r).fac() * ModIntRuntime(r).fac());
}
ModIntRuntime sqrt() const {
if (val_ == 0) return 0;
if (md == 2) return val_;
if (power((md - 1) / 2) != 1) return 0;
ModIntRuntime b = 1;
while (b.power((md - 1) / 2) == 1) b += 1;
int e = 0, m = md - 1;
while (m % 2 == 0) m >>= 1, e++;
ModIntRuntime x = power((m - 1) / 2), y = (*this) * x * x;
x *= (*this);
ModIntRuntime z = b.power(m);
while (y != 1) {
int j = 0;
ModIntRuntime t = y;
while (t != 1) j++, t *= t;
z = z.power(1LL << (e - j - 1));
x *= z, z *= z, y *= z;
e = j;
}
return ModIntRuntime(std::min(x.val_, md - x.val_));
}
};
int ModIntRuntime::md = 1;
#line 4 "number/test/sieve.stress.test.cpp"
#include <algorithm>
#line 6 "number/test/sieve.stress.test.cpp"
#include <cstdio>
#line 9 "number/test/sieve.stress.test.cpp"
using namespace std;
struct Case {
int SIEVE_SIZE;
int MAX;
};
int euler_phi(int x) {
int ret = 0;
for (int d = 1; d <= x; d++) ret += (std::__gcd(d, x) == 1);
return ret;
}
void test_divisors(Case testcase) {
const Sieve sieve(testcase.SIEVE_SIZE);
for (int x = 1; x <= testcase.MAX; x++) {
auto divs = sieve.divisors(x);
std::vector<int> is_div(x + 1);
for (auto d : divs) is_div.at(d) = 1;
for (int y = 1; y <= x; y++) assert(is_div.at(y) == (x % y == 0));
}
cerr << "divisors(): passed" << endl;
}
void test_euler_of_divisors(Case testcase) {
const Sieve sieve(testcase.SIEVE_SIZE);
for (int x = 1; x <= testcase.MAX; x++) {
auto div2euler = sieve.euler_of_divisors(x);
for (auto de : div2euler) {
assert(euler_phi(de.first) == de.second);
assert(x % de.first == 0);
}
assert(div2euler.size() == sieve.divisors(x).size());
}
cerr << "euler_of_divisors(): passed" << endl;
}
void test_moebius_table(int HI) {
const Sieve sieve_hi(HI);
const auto answer = sieve_hi.GenerateMoebiusFunctionTable();
assert(int(answer.size()) == HI + 1);
for (int x = 1; x <= HI; x++) {
int ret = 1;
for (auto p : sieve_hi.primes) {
if (x % (p * p) == 0) ret = 0;
if (x % p == 0) ret *= -1;
}
assert(answer[x] == ret);
}
for (int x = 1; x <= HI - 1; x++) {
auto mu = Sieve(x).GenerateMoebiusFunctionTable();
assert(int(mu.size()) == x + 1);
for (int i = 0; i < x + 1; i++) assert(mu[i] == answer[i]);
}
cerr << "GenerateMoebiusFunctionTable(): passed" << endl;
}
void test_enumerate_kth_pows(int hi_sieve_size) {
for (int sieve_size = 1; sieve_size <= hi_sieve_size; sieve_size++) {
const Sieve sieve(sieve_size);
for (int nmax = 1; nmax <= sieve_size; nmax++) {
for (int p = 1; p <= 20; p++) {
using mint = ModIntRuntime;
mint::set_mod(p);
for (int k = 0; k < 10; k++) {
auto ret = sieve.enumerate_kth_pows<mint>(k, nmax);
assert(int(ret.size()) == nmax + 1);
for (int i = 0; i <= nmax; i++) assert(ret[i] == mint(i).pow(k));
}
}
}
}
cerr << "enumerate_kth_pows(): passed" << endl;
}
int main() {
std::vector<Case> cases{{31, 31 * 31}, {100, 1000}, {2, 4}};
for (auto testcase : cases) {
test_divisors(testcase);
test_euler_of_divisors(testcase);
}
test_moebius_table(1000);
test_enumerate_kth_pows(100);
puts("Hello World");
}