library

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

View the Project on GitHub gyouzasushi/library

:heavy_check_mark: test/unit/sieve.test.cpp

Depends on

Code

#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;
}
Back to top page