This documentation is automatically generated by online-judge-tools/verification-helper
#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A"
#include <cassert>
#include <iostream>
#include <random>
#include <vector>
#include "math/factorize.hpp"
#include "math/sieve.hpp"
std::vector<std::pair<long long, int>> naive_factorize(long long x) {
std::vector<std::pair<long long, int>> ret;
for (long long p = 2; p * p <= x; p++) {
int e = 0;
while (x % p == 0) x /= p, e++;
if (e) ret.emplace_back(p, e);
}
if (x > 1) ret.emplace_back(x, 1);
return ret;
}
std::vector<long long> naive_divisors(long long x) {
std::vector<long long> ret;
for (long long d = 1; d <= x; d++) {
if (x % d == 0) ret.push_back(d);
}
return ret;
}
int main() {
const int n = 3000;
sieve s(n);
{
std::vector<int> ps;
for (int x = 0; x <= n; x++) {
bool prime = x >= 2;
for (int d = 2; d * d <= x; d++) prime &= x % d != 0;
assert(s.is_prime(x) == prime);
if (prime) ps.push_back(x);
if (x >= 2) assert(s.smallest_prime_factor(x) == naive_factorize(x)[0].first);
if (x >= 1) {
auto f = s.factorize(x);
auto g = naive_factorize(x);
assert(f.size() == g.size());
for (size_t i = 0; i < f.size(); i++) assert(f[i].first == g[i].first && f[i].second == g[i].second);
auto d = s.divisors(x);
auto nd = naive_divisors(x);
assert(d.size() == nd.size());
for (size_t i = 0; i < d.size(); i++) assert(d[i] == nd[i]);
}
}
assert(s.primes() == ps);
}
// sieve(1), sieve(2)
{
sieve s1(1);
assert(s1.primes().empty() && !s1.is_prime(1) && s1.factorize(1).empty() && s1.divisors(1) == std::vector<int>{1});
sieve s2(2);
assert(s2.primes() == std::vector<int>{2} && s2.is_prime(2));
}
// large n: pi(10^7) = 664579
{
sieve big(10000000);
assert(int(big.primes().size()) == 664579);
assert((big.factorize(9999991) == std::vector<std::pair<int, int>>{{9999991, 1}}));
assert(big.divisors(10000000).size() == 8 * 8);
}
// Pollard-rho factorize / divisors (vector form)
{
std::mt19937_64 rng(5);
for (long long x = 1; x <= 2000; x++) {
assert(factorize(x) == naive_factorize(x));
assert(divisors(x) == naive_divisors(x));
}
for (int it = 0; it < 300; it++) {
long long x = 1 + rng() % 1000000000000LL;
assert(factorize(x) == naive_factorize(x));
auto d = divisors(x);
long long cnt = 1;
for (auto [p, e] : factorize(x)) cnt *= e + 1;
assert((long long)d.size() == cnt);
for (size_t i = 0; i < d.size(); i++) {
assert(x % d[i] == 0);
if (i) assert(d[i - 1] < d[i]);
}
}
auto d = divisors(897612484786617600LL); // highly composite, 103680 divisors
assert(d.size() == 103680);
}
std::cout << "Hello World" << std::endl;
return 0;
}#line 1 "test/unit/sieve.test.cpp"
#define PROBLEM "http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=ITP1_1_A"
#include <cassert>
#include <iostream>
#include <random>
#include <vector>
#line 2 "math/factorize.hpp"
#include <algorithm>
#line 4 "math/factorize.hpp"
#include <cstdlib>
#include <numeric>
#include <utility>
#line 8 "math/factorize.hpp"
namespace internal {
long long modmul(long long x, long long y, long long mod) {
return (long long)((__int128_t)x * y % mod);
}
long long modpow(long long a, long long n, long long mod) {
long long ret = 1;
while (n > 0) {
if (n & 1) ret = modmul(ret, a, mod);
a = modmul(a, a, mod);
n >>= 1;
}
return ret;
}
long long rho(long long n) {
long long z = 0;
auto f = [&](long long x) {
long long ret = modmul(x, x, n) + z;
return ret >= n ? ret - n : ret;
};
while (true) {
long long x = ++z;
long long y = f(x);
while (true) {
long long d = std::gcd(std::abs(x - y), n);
if (d == n) break;
if (d > 1) return d;
x = f(x);
y = f(f(y));
}
}
}
} // namespace internal
bool is_prime(long long n) {
assert(n >= 1);
if (n == 1) return false;
if (n % 2 == 0) return n == 2;
long long d = n - 1, s = 0;
while (~d & 1) d >>= 1, s++;
for (long long a : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
if (a >= n) break;
long long x = internal::modpow(a, d, n);
if (x == 1) continue;
bool ok = false;
for (int i = 0; i < s; i++) {
if (x == n - 1) {
ok = true;
break;
}
x = internal::modmul(x, x, n);
}
if (!ok) return false;
}
return true;
}
std::vector<std::pair<long long, int>> factorize(long long n) {
assert(n >= 1);
std::vector<std::pair<long long, int>> ret;
int c = 0;
while (~n & 1) n >>= 1, c++;
if (c > 0) ret.emplace_back(2, c);
std::vector<long long> st = {n};
std::vector<long long> primes;
while (!st.empty()) {
long long p = st.back();
st.pop_back();
if (p == 1) continue;
if (is_prime(p)) {
primes.push_back(p);
continue;
}
long long d = internal::rho(p);
st.push_back(d);
st.push_back(p / d);
}
std::sort(primes.begin(), primes.end());
for (long long p : primes) {
if (!ret.empty() && ret.back().first == p) {
ret.back().second++;
} else {
ret.emplace_back(p, 1);
}
}
return ret;
}
std::vector<long long> divisors(long long n) {
assert(n >= 1);
std::vector<long long> ret = {1};
for (auto [p, e] : factorize(n)) {
int sz = int(ret.size());
long long q = 1;
for (int i = 0; i < e; i++) {
q *= p;
for (int j = 0; j < sz; j++) ret.push_back(ret[j] * q);
}
}
std::sort(ret.begin(), ret.end());
return ret;
}
#line 6 "math/sieve.hpp"
struct sieve {
public:
sieve() : sieve(1) {
}
explicit sieve(int n) : _n(n), lpf(n + 1, 0) {
assert(n >= 1);
for (int i = 2; i <= n; i++) {
if (lpf[i] == 0) {
lpf[i] = i;
_primes.push_back(i);
}
for (int p : _primes) {
if (p > lpf[i] || (long long)p * i > n) break;
lpf[p * i] = p;
}
}
}
bool is_prime(int x) const {
assert(0 <= x && x <= _n);
return x >= 2 && lpf[x] == x;
}
int smallest_prime_factor(int x) const {
assert(2 <= x && x <= _n);
return lpf[x];
}
const std::vector<int>& primes() const {
return _primes;
}
std::vector<std::pair<int, int>> factorize(int x) const {
assert(1 <= x && x <= _n);
std::vector<std::pair<int, int>> ret;
while (x > 1) {
int p = lpf[x];
if (!ret.empty() && ret.back().first == p) {
ret.back().second++;
} else {
ret.emplace_back(p, 1);
}
x /= p;
}
return ret;
}
std::vector<int> divisors(int x) const {
assert(1 <= x && x <= _n);
std::vector<int> ret = {1};
for (auto [p, e] : factorize(x)) {
int sz = int(ret.size());
int q = 1;
for (int i = 0; i < e; i++) {
q *= p;
for (int j = 0; j < sz; j++) ret.push_back(ret[j] * q);
}
}
std::sort(ret.begin(), ret.end());
return ret;
}
private:
int _n;
std::vector<int> lpf, _primes;
};
#line 9 "test/unit/sieve.test.cpp"
std::vector<std::pair<long long, int>> naive_factorize(long long x) {
std::vector<std::pair<long long, int>> ret;
for (long long p = 2; p * p <= x; p++) {
int e = 0;
while (x % p == 0) x /= p, e++;
if (e) ret.emplace_back(p, e);
}
if (x > 1) ret.emplace_back(x, 1);
return ret;
}
std::vector<long long> naive_divisors(long long x) {
std::vector<long long> ret;
for (long long d = 1; d <= x; d++) {
if (x % d == 0) ret.push_back(d);
}
return ret;
}
int main() {
const int n = 3000;
sieve s(n);
{
std::vector<int> ps;
for (int x = 0; x <= n; x++) {
bool prime = x >= 2;
for (int d = 2; d * d <= x; d++) prime &= x % d != 0;
assert(s.is_prime(x) == prime);
if (prime) ps.push_back(x);
if (x >= 2) assert(s.smallest_prime_factor(x) == naive_factorize(x)[0].first);
if (x >= 1) {
auto f = s.factorize(x);
auto g = naive_factorize(x);
assert(f.size() == g.size());
for (size_t i = 0; i < f.size(); i++) assert(f[i].first == g[i].first && f[i].second == g[i].second);
auto d = s.divisors(x);
auto nd = naive_divisors(x);
assert(d.size() == nd.size());
for (size_t i = 0; i < d.size(); i++) assert(d[i] == nd[i]);
}
}
assert(s.primes() == ps);
}
// sieve(1), sieve(2)
{
sieve s1(1);
assert(s1.primes().empty() && !s1.is_prime(1) && s1.factorize(1).empty() && s1.divisors(1) == std::vector<int>{1});
sieve s2(2);
assert(s2.primes() == std::vector<int>{2} && s2.is_prime(2));
}
// large n: pi(10^7) = 664579
{
sieve big(10000000);
assert(int(big.primes().size()) == 664579);
assert((big.factorize(9999991) == std::vector<std::pair<int, int>>{{9999991, 1}}));
assert(big.divisors(10000000).size() == 8 * 8);
}
// Pollard-rho factorize / divisors (vector form)
{
std::mt19937_64 rng(5);
for (long long x = 1; x <= 2000; x++) {
assert(factorize(x) == naive_factorize(x));
assert(divisors(x) == naive_divisors(x));
}
for (int it = 0; it < 300; it++) {
long long x = 1 + rng() % 1000000000000LL;
assert(factorize(x) == naive_factorize(x));
auto d = divisors(x);
long long cnt = 1;
for (auto [p, e] : factorize(x)) cnt *= e + 1;
assert((long long)d.size() == cnt);
for (size_t i = 0; i < d.size(); i++) {
assert(x % d[i] == 0);
if (i) assert(d[i - 1] < d[i]);
}
}
auto d = divisors(897612484786617600LL); // highly composite, 103680 divisors
assert(d.size() == 103680);
}
std::cout << "Hello World" << std::endl;
return 0;
}