CP-Algorithms Library

This documentation is automatically generated by competitive-verifier/competitive-verifier

View the Project on GitHub cp-algorithms/cp-algorithms-aux

:heavy_check_mark: cp-algo/number_theory/factorize.hpp

Depends on

Required by

Verified with

Code

#ifndef CP_ALGO_MATH_FACTORIZE_HPP
#define CP_ALGO_MATH_FACTORIZE_HPP
#include "primality.hpp"
#include "../util/generator.hpp"
#include "../random/rng.hpp"
namespace cp_algo::math {
    // https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    template<typename _Int>
    auto proper_divisor(_Int m) {
        using Int = std::make_signed_t<_Int>;
        using base = dynamic_modint<Int>;
        return m % 2 == 0 ? 2 : base::with_mod(m, [&]() {
            base t = random::rng();
            auto f = [&](auto x) {
                return x * x + t;
            };
            const base zero = 0, one = 1;
            base x = zero, y = zero;
            base g = one;
            while(g == one) {
                for(int i = 0; i < 64; i++) {
                    x = f(x);
                    y = f(f(y));
                    if(x == y) [[unlikely]] {
                        t = random::rng();
                        x = y = zero;
                    } else {
                        base t = g * (x - y);
                        g = t == zero ? g : t;
                    }
                }
                g = std::gcd(g.getr(), m);
            }
            return g.getr();
        });
    }
    template<typename Int>
    big_generator<Int> factorize(Int m) {
        if(is_prime(m)) {
            co_yield m;
        } else if(m > 1) {
            auto g = proper_divisor(m);
            co_yield std::ranges::elements_of(factorize(g));
            co_yield std::ranges::elements_of(factorize(m / g));
        }
    }
    template<typename Int>
    big_generator<Int> divisors_sqrt(Int m) {
        for(Int i = 1; i * i <= m; i++) {
            if(m % i == 0) {
                co_yield i;
                if(i * i != m) {
                    co_yield m / i;
                }
            }
        }
    }
}
#endif // CP_ALGO_MATH_FACTORIZE_HPP

#line 1 "cp-algo/number_theory/factorize.hpp"


#line 1 "cp-algo/number_theory/primality.hpp"


#line 1 "cp-algo/number_theory/modint.hpp"


#line 1 "cp-algo/math/common.hpp"


#include <functional>
#include <cstdint>
#include <cassert>
#include <bit>
#include <vector>
#include <algorithm>
namespace cp_algo::math {
#ifdef CP_ALGO_MAXN
    const int maxn = CP_ALGO_MAXN;
#else
    const int maxn = 1 << 19;
#endif
    const int magic = 64; // threshold for sizes to run the naive algo

    // Nonnegative 64-bit exponents, with an associative operation and its identity.
    // Windows >1 precompute odd powers only when that saves operations.
    template<int window = 1>
    auto bpow(auto const& x, auto n, auto const& one, auto op) {
        static_assert(window >= 1 && window <= 6);
        if constexpr(window > 1) {
            if(n == 0) {return one;}
            int bits = std::bit_width(uint64_t(n));
            auto low_bit = [&](int high) {
                int low = std::max(0, high - window + 1);
                while(!((n >> low) & 1)) {low++;}
                return low;
            };
            int first = low_bit(bits - 1);
            int cost = (1 << (window - 1)) + first;
            for(int j = first - 1; j >= 0;) {
                if(!((n >> j) & 1)) {j--;}
                else {cost++; j = low_bit(j) - 1;}
            }
            // Do not pay for the table when binary powering uses fewer operations.
            if(cost >= bits + std::popcount(uint64_t(n)) - 2) {return bpow<1>(x, n, one, op);}
            using T = std::decay_t<decltype(x)>;
            std::vector<T> odd;
            odd.reserve(1 << (window - 1));
            odd.push_back(x);
            auto square = op(x, x);
            while(odd.size() < size_t(1 << (window - 1))) {odd.push_back(op(odd.back(), square));}
            auto ans = odd[(n >> first) / 2];
            for(int j = first - 1; j >= 0;) {
                if(!((n >> j) & 1)) {ans = op(ans, ans); j--;}
                else {
                    int low = low_bit(j), length = j - low + 1;
                    auto digit = (n >> low) & ((1u << length) - 1);
                    for(int i = 0; i < length; i++) {ans = op(ans, ans);}
                    ans = op(ans, odd[digit / 2]);
                    j = low - 1;
                }
            }
            return ans;
        } else {
            if (n == 0) {
                return one;
            }
            auto ans = x;
            for(int j = std::bit_width<uint64_t>(n) - 2; ~j; j--) {
                ans = op(ans, ans);
                if((n >> j) & 1) {
                    ans = op(ans, x);
                }
            }
            return ans;
        }
    }
    template<int window = 1>
    auto bpow(auto x, auto n, auto ans) {
        return bpow<window>(x, n, ans, std::multiplies{});
    }
    template<typename T>
    T bpow(T const& x, auto n) {
        return bpow(x, n, T(1));
    }
    inline constexpr auto inv2(auto x) {
        assert(x % 2);
        std::make_unsigned_t<decltype(x)> y = 1;
        while(y * x != 1) {
            y *= 2 - x * y;
        }
        return y;
    }
}

#line 4 "cp-algo/number_theory/modint.hpp"
#include <iostream>
#line 6 "cp-algo/number_theory/modint.hpp"
namespace cp_algo::math {

    template<typename modint, typename _Int>
    struct modint_base {
        using Int = _Int;
        using UInt = std::make_unsigned_t<Int>;
        static constexpr size_t bits = sizeof(Int) * 8;
        using Int2 = std::conditional_t<bits <= 32, int64_t, __int128_t>;
        using UInt2 = std::conditional_t<bits <= 32, uint64_t, __uint128_t>;
        constexpr static Int mod() {
            return modint::mod();
        }
        constexpr static UInt remod() {
            return modint::remod();
        }
        constexpr static UInt2 modmod() {
            return UInt2(mod()) * mod();
        }
        constexpr modint_base() = default;
        constexpr modint_base(Int2 rr) {
            to_modint().setr(UInt((rr + modmod()) % mod()));
        }
        constexpr modint inv() const {
            return bpow(to_modint(), mod() - 2);
        }
        modint operator-() const {
            modint neg;
            neg.r = std::min(-r, remod() - r);
            return neg;
        }
        modint& operator/=(const modint& t) {
            return to_modint() *= t.inv();
        }
        modint& operator*=(const modint& t) {
            r = UInt(UInt2(r) * t.r % mod());
            return to_modint();
        }
        modint& operator+=(const modint& t) {
            r += t.r;
            r = std::min(r, r - remod());
            return to_modint();
        }
        modint& operator-=(const modint& t) {
            r -= t.r;
            r = std::min(r, r + remod());
            return to_modint();
        }
        modint operator+(const modint& t) const { return modint(to_modint()) += t; }
        modint operator-(const modint& t) const { return modint(to_modint()) -= t; }
        modint operator*(const modint& t) const { return modint(to_modint()) *= t; }
        modint operator/(const modint& t) const { return modint(to_modint()) /= t; }
        // Why <=> doesn't work?..
        auto operator==(const modint& t) const { return to_modint().getr() == t.getr(); }
        auto operator!=(const modint& t) const { return to_modint().getr() != t.getr(); }
        auto operator<=(const modint& t) const { return to_modint().getr() <= t.getr(); }
        auto operator>=(const modint& t) const { return to_modint().getr() >= t.getr(); }
        auto operator<(const modint& t) const { return to_modint().getr() < t.getr(); }
        auto operator>(const modint& t) const { return to_modint().getr() > t.getr(); }
        Int rem() const {
            UInt R = to_modint().getr();
            return R - (R > (UInt)mod() / 2) * mod();
        }
        constexpr void setr(UInt rr) {
            r = rr;
        }
        constexpr UInt getr() const {
            return r;
        }

        // Only use these if you really know what you're doing!
        static uint64_t modmod8() { return uint64_t(8 * modmod()); }
        void add_unsafe(UInt t) { r += t; }
        void pseudonormalize() { r = std::min(r, r - modmod8()); }
        modint const& normalize() {
            if(r >= (UInt)mod()) {
                r %= mod();
            }
            return to_modint();
        }
        void setr_direct(UInt rr) { r = rr; }
        UInt getr_direct() const { return r; }

    protected:
        UInt r;

    private:
        constexpr modint& to_modint() { return static_cast<modint&>(*this); }
        constexpr modint const& to_modint() const { return static_cast<modint const&>(*this); }
    };
    template<typename modint>
    concept modint_type = std::is_base_of_v<modint_base<modint, typename modint::Int>, modint>;
    template<modint_type modint>
    decltype(std::cin)& operator>>(decltype(std::cin)& in, modint& x) {
        typename modint::UInt r;
        auto &res = in >> r;
        x.setr(r);
        return res;
    }
    template<modint_type modint>
    decltype(std::cout)& operator<<(decltype(std::cout)& out, modint const& x) {
        return out << x.getr();
    }

    template<auto m>
    struct modint: modint_base<modint<m>, decltype(m)> {
        using Base = modint_base<modint<m>, decltype(m)>;
        using Base::Base;
        static constexpr Base::Int mod() { return m; }
        static constexpr Base::UInt remod() { return m; }
        auto getr() const { return Base::r; }
    };

    // Odd moduli up to a quarter of the unsigned word keep Montgomery residues lazily in
    // [0, 2 mod): remod() = 2 mod, and both 4 mod and ab + q * mod fit. Any other modulus keeps
    // fully reduced residues, remod() = mod, which is what the sums of modint_base need then:
    // an even one without the Montgomery form, a wide odd one with a reduction that subtracts
    // high words instead of adding double words that would overflow.
    template<typename Int = int>
    struct dynamic_modint: modint_base<dynamic_modint<Int>, Int> {
        using Base = modint_base<dynamic_modint<Int>, Int>;
        using Base::Base;

        // Out of line, so that the hot path stays as small as it was.
        [[gnu::noinline, gnu::cold]] static Base::UInt m_reduce_reduced(Base::UInt2 ab) {
            if(mod() % 2 == 0) {
                return typename Base::UInt(ab % mod());
            }
            // q * mod has the low word of ab, so the difference of the high words is exact.
            typename Base::UInt q = -(typename Base::UInt(ab) * inverse);
            auto high = typename Base::UInt(ab >> Base::bits);
            auto low = typename Base::UInt(typename Base::UInt2(q) * typename Base::UInt(mod()) >> Base::bits);
            return high >= low ? high - low : high - low + mod();
        }
        static Base::UInt m_reduce(Base::UInt2 ab) {
            if(imod() == 0) [[unlikely]] {
                return m_reduce_reduced(ab);
            }
            typename Base::UInt2 m = typename Base::UInt(ab) * imod();
            return typename Base::UInt((ab + m * mod()) >> Base::bits);
        }
        static Base::UInt m_transform(Base::UInt a) {
            if(mod() % 2 == 0) [[unlikely]] {
                return a;
            } else {
                return m_reduce(a * pw128());
            }
        }
        dynamic_modint& operator*=(const dynamic_modint& t) {
            Base::r = m_reduce(typename Base::UInt2(Base::r) * t.r);
            return *this;
        }
        void setr(Base::UInt rr) {
            Base::r = m_transform(rr);
        }
        Base::UInt getr() const {
            typename Base::UInt res = m_reduce(Base::r);
            return std::min(res, res - mod());
        }
        static Int mod() { return m; }
        static Base::UInt remod() { return rm; }
        static Base::UInt imod() { return im; }
        static Base::UInt2 pw128() { return r2; }
        static void switch_mod(Int nm) {
            m = nm;
            bool lazy = m % 2 && typename Base::UInt(m) <= typename Base::UInt(-1) / 4;
            rm = typename Base::UInt(m) * (lazy ? 2 : 1);
            inverse = m % 2 ? inv2(-m) : 0;
            im = lazy ? inverse : 0;
            r2 = static_cast<Base::UInt>(static_cast<Base::UInt2>(-1) % m + 1);
        }

        // Wrapper for temp switching
        auto static with_mod(Int tmp, auto callback) {
            struct scoped {
                Int prev = mod();
                ~scoped() { switch_mod(prev); }
            } _;
            switch_mod(tmp);
            return callback();
        }

    private:
        static inline thread_local Int m = 1;
        // im: -1 / mod modulo 2^bits for lazy residues and 0 for reduced ones; inverse: the same
        // for every odd mod; rm: the value of remod().
        static inline thread_local Base::UInt im = -1, r2 = 0, inverse = -1, rm = 2;
    };
}

#line 6 "cp-algo/number_theory/primality.hpp"
namespace cp_algo::math {
    // https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
    template<typename _Int>
    bool is_prime(_Int m) {
        using Int = std::make_signed_t<_Int>;
        using UInt = std::make_unsigned_t<Int>;
        if(m == 1 || m % 2 == 0) {
            return m == 2;
        }
        // m - 1 = 2^s * d
        int s = std::countr_zero(UInt(m - 1));
        auto d = (m - 1) >> s;
        using base = dynamic_modint<Int>;
        auto test = [&](base x) {
            x = bpow(x, d);
            if(std::abs(x.rem()) <= 1) {
                return true;
            }
            for(int i = 1; i < s && x != -1; i++) {
                x *= x;
            }
            return x == -1;
        };
        return base::with_mod(m, [&]() {
#ifdef CP_ALGO_NUMBER_THEORY_PRIMALITY_BASES_HPP
            uint16_t base2 = 7, base3 = 61;
            if (m != uint32_t(m)) {
                base2 = base_table1[uint32_t(m * 0xAD625B89) >> 18];
                base3 = base_table2[base2 >> 13];
            }
            return test(2) && test(base2) && test(base3);
#else
            return std::ranges::all_of(std::array{2, 325, 9375, 28178, 450775, 9780504, 1795265022}, test);
#endif
        });
    }
}

#line 1 "cp-algo/util/generator.hpp"



#line 1 "cp-algo/util/big_alloc.hpp"



#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#line 10 "cp-algo/util/big_alloc.hpp"
#include <string>
#include <cstddef>
#line 14 "cp-algo/util/big_alloc.hpp"
#include <forward_list>

// Single macro to detect POSIX platforms (Linux, Unix, macOS)
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#define CP_ALGO_USE_MMAP 1
#include <sys/mman.h>
#else
#define CP_ALGO_USE_MMAP 0
#endif

namespace cp_algo {
    template<typename T, size_t Align = 32>
    class big_alloc {
        static_assert(Align >= alignof(void*), "Align must be at least pointer-size");
        static_assert(std::has_single_bit(Align), "Align must be a power of two");

    public:
        using value_type = T;
        template<class U>
        struct rebind {
            using other = big_alloc<U, Align>;
        };
        constexpr bool operator==(const big_alloc&) const = default;

        big_alloc() noexcept = default;
        template<typename U, std::size_t A>
        big_alloc(const big_alloc<U, A>&) noexcept {}

        [[nodiscard]] T* allocate(std::size_t n) {
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
            if(padded >= MEGABYTE) {
                void* raw = mmap(nullptr, padded, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0);
                madvise(raw, padded, MADV_HUGEPAGE);
                return static_cast<T*>(raw);
            }
#endif
            return static_cast<T*>(::operator new(padded, std::align_val_t(align)));
        }

        void deallocate(T* p, std::size_t n) noexcept {
            if(!p) {
                return;
            }
            std::size_t padded = round_up(n * sizeof(T));
            std::size_t align = std::max<std::size_t>(alignof(T), Align);
#if CP_ALGO_USE_MMAP
            if(padded >= MEGABYTE) {
                munmap(p, padded);
                return;
            }
#endif
            ::operator delete(p, padded, std::align_val_t(align));
        }

    private:
        static constexpr std::size_t MEGABYTE = 1 << 20;
        static constexpr std::size_t round_up(std::size_t x) noexcept {
            return (x + Align - 1) / Align * Align;
        }
    };

    template<typename T>
    using big_vector = std::vector<T, big_alloc<T>>;
    template<typename T>
    using big_basic_string = std::basic_string<T, std::char_traits<T>, big_alloc<T>>;
    template<typename T>
    using big_deque = std::deque<T, big_alloc<T>>;
    template<typename T>
    using big_stack = std::stack<T, big_deque<T>>;
    template<typename T>
    using big_queue = std::queue<T, big_deque<T>>;
    template<typename T>
    using big_priority_queue = std::priority_queue<T, big_vector<T>>;
    template<typename T>
    using big_forward_list = std::forward_list<T, big_alloc<T>>;
    using big_string = big_basic_string<char>;

    template<typename Key, typename Value, typename Compare = std::less<Key>>
    using big_map = std::map<Key, Value, Compare, big_alloc<std::pair<const Key, Value>>>;
    template<typename T, typename Compare = std::less<T>>
    using big_multiset = std::multiset<T, Compare, big_alloc<T>>;
    template<typename T, typename Compare = std::less<T>>
    using big_set = std::set<T, Compare, big_alloc<T>>;
}


#line 5 "cp-algo/util/generator.hpp"
#include <generator>

namespace cp_algo {
    template<typename Ref, typename V = void>
    using big_generator = std::generator<Ref, V, big_alloc<std::byte>>;
}

namespace std::ranges {
    template<typename Ref, typename V>
    elements_of(cp_algo::big_generator<Ref, V>&&) -> elements_of<cp_algo::big_generator<Ref, V>&&, cp_algo::big_alloc<std::byte>>;
}


#line 1 "cp-algo/random/rng.hpp"


#include <chrono>
#include <random>
namespace cp_algo::random {
    std::mt19937_64 gen(
        std::chrono::steady_clock::now().time_since_epoch().count()
    );
    uint64_t rng() {
        return gen();
    }
}

#line 6 "cp-algo/number_theory/factorize.hpp"
namespace cp_algo::math {
    // https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
    template<typename _Int>
    auto proper_divisor(_Int m) {
        using Int = std::make_signed_t<_Int>;
        using base = dynamic_modint<Int>;
        return m % 2 == 0 ? 2 : base::with_mod(m, [&]() {
            base t = random::rng();
            auto f = [&](auto x) {
                return x * x + t;
            };
            const base zero = 0, one = 1;
            base x = zero, y = zero;
            base g = one;
            while(g == one) {
                for(int i = 0; i < 64; i++) {
                    x = f(x);
                    y = f(f(y));
                    if(x == y) [[unlikely]] {
                        t = random::rng();
                        x = y = zero;
                    } else {
                        base t = g * (x - y);
                        g = t == zero ? g : t;
                    }
                }
                g = std::gcd(g.getr(), m);
            }
            return g.getr();
        });
    }
    template<typename Int>
    big_generator<Int> factorize(Int m) {
        if(is_prime(m)) {
            co_yield m;
        } else if(m > 1) {
            auto g = proper_divisor(m);
            co_yield std::ranges::elements_of(factorize(g));
            co_yield std::ranges::elements_of(factorize(m / g));
        }
    }
    template<typename Int>
    big_generator<Int> divisors_sqrt(Int m) {
        for(Int i = 1; i * i <= m; i++) {
            if(m % i == 0) {
                co_yield i;
                if(i * i != m) {
                    co_yield m / i;
                }
            }
        }
    }
}


#ifndef CP_ALGO_MATH_FACTORIZE_HPP
#define CP_ALGO_MATH_FACTORIZE_HPP
#include "primality.hpp"
#include "../util/generator.hpp"
#include "../random/rng.hpp"
namespace cp_algo::math{ template<typename _Int> auto proper_divisor(_Int m){ using Int=std::make_signed_t<_Int>; using base=dynamic_modint<Int>; return m%2==0?2:base::with_mod(m,[&](){ base t=random::rng(); auto f= [&](auto x){ return x*x+t; }; const base zero=0,one=1; base x=zero,y=zero; base g=one; while(g==one){ for(int i=0;i<64;i++){ x=f(x); y=f(f(y)); if(x==y)[[unlikely]]{ t=random::rng(); x=y=zero; }else{ base t=g*(x-y); g=t==zero?g:t; } } g=std::gcd(g.getr(),m); } return g.getr(); }); } template<typename Int> big_generator<Int>factorize(Int m){ if(is_prime(m)){ co_yield m; }else if(m>1){ auto g=proper_divisor(m); co_yield std::ranges::elements_of(factorize(g)); co_yield std::ranges::elements_of(factorize(m/g)); } } template<typename Int> big_generator<Int>divisors_sqrt(Int m){ for(Int i=1;i*i<=m;i++){ if(m%i==0){ co_yield i; if(i*i!=m){ co_yield m/i; } } } } }
#endif

#line 1 "cp-algo/number_theory/factorize.hpp"
#line 1 "cp-algo/number_theory/primality.hpp"
#line 1 "cp-algo/number_theory/modint.hpp"
#line 1 "cp-algo/math/common.hpp"
#include <functional>
#include <cstdint>
#include <cassert>
#include <bit>
#include <vector>
#include <algorithm>
namespace cp_algo::math{
#ifdef CP_ALGO_MAXN
const int maxn=CP_ALGO_MAXN;
#else
const int maxn=1<<19;
#endif
const int magic=64; template<int window=1> auto bpow(auto const&x,auto n,auto const&one,auto op){ static_assert(window>=1&&window<=6); if constexpr(window>1){ if(n==0){return one;} int bits=std::bit_width(uint64_t(n)); auto low_bit= [&](int high){ int low=std::max(0,high-window+1); while(!((n>>low)&1)){low++;} return low; }; int first=low_bit(bits-1); int cost=(1<<(window-1))+first; for(int j=first-1;j>=0;){ if(!((n>>j)&1)){j--;} else{cost++;j=low_bit(j)-1;} } if(cost>=bits+std::popcount(uint64_t(n))-2){return bpow<1>(x,n,one,op);} using T=std::decay_t<decltype(x)>; std::vector<T>odd; odd.reserve(1<<(window-1)); odd.push_back(x); auto square=op(x,x); while(odd.size()<size_t(1<<(window-1))){odd.push_back(op(odd.back(),square));} auto ans=odd[(n>>first)/2]; for(int j=first-1;j>=0;){ if(!((n>>j)&1)){ans=op(ans,ans);j--;} else{ int low=low_bit(j),length=j-low+1; auto digit=(n>>low)&((1u<<length)-1); for(int i=0;i<length;i++){ans=op(ans,ans);} ans=op(ans,odd[digit/2]); j=low-1; } } return ans; }else{ if(n==0){ return one; } auto ans=x; for(int j=std::bit_width<uint64_t>(n)-2;~j;j--){ ans=op(ans,ans); if((n>>j)&1){ ans=op(ans,x); } } return ans; } } template<int window=1> auto bpow(auto x,auto n,auto ans){ return bpow<window>(x,n,ans,std::multiplies{}); } template<typename T> T bpow(T const&x,auto n){ return bpow(x,n,T(1)); } inline constexpr auto inv2(auto x){ assert(x%2); std::make_unsigned_t<decltype(x)>y=1; while(y*x!=1){ y*=2-x*y; } return y; } }
#line 4 "cp-algo/number_theory/modint.hpp"
#include <iostream>
#line 6 "cp-algo/number_theory/modint.hpp"
namespace cp_algo::math{ template<typename modint,typename _Int> struct modint_base{ using Int=_Int; using UInt=std::make_unsigned_t<Int>; static constexpr size_t bits=sizeof(Int)*8; using Int2=std::conditional_t<bits<=32,int64_t,__int128_t>; using UInt2=std::conditional_t<bits<=32,uint64_t,__uint128_t>; constexpr static Int mod(){ return modint::mod(); } constexpr static UInt remod(){ return modint::remod(); } constexpr static UInt2 modmod(){ return UInt2(mod())*mod(); } constexpr modint_base()=default; constexpr modint_base(Int2 rr){ to_modint().setr(UInt((rr+modmod())%mod())); } constexpr modint inv()const{ return bpow(to_modint(),mod()-2); } modint operator-()const{ modint neg; neg.r=std::min(-r,remod()-r); return neg; } modint&operator/=(const modint&t){ return to_modint()*=t.inv(); } modint&operator*=(const modint&t){ r=UInt(UInt2(r)*t.r%mod()); return to_modint(); } modint&operator+=(const modint&t){ r+=t.r; r=std::min(r,r-remod()); return to_modint(); } modint&operator-=(const modint&t){ r-=t.r; r=std::min(r,r+remod()); return to_modint(); } modint operator+(const modint&t)const{return modint(to_modint())+=t;} modint operator-(const modint&t)const{return modint(to_modint())-=t;} modint operator*(const modint&t)const{return modint(to_modint())*=t;} modint operator/(const modint&t)const{return modint(to_modint())/=t;} auto operator==(const modint&t)const{return to_modint().getr()==t.getr();} auto operator!=(const modint&t)const{return to_modint().getr()!=t.getr();} auto operator<=(const modint&t)const{return to_modint().getr()<=t.getr();} auto operator>=(const modint&t)const{return to_modint().getr()>=t.getr();} auto operator<(const modint&t)const{return to_modint().getr()<t.getr();} auto operator>(const modint&t)const{return to_modint().getr()>t.getr();} Int rem()const{ UInt R=to_modint().getr(); return R-(R>(UInt)mod()/2)*mod(); } constexpr void setr(UInt rr){ r=rr; } constexpr UInt getr()const{ return r; } static uint64_t modmod8(){return uint64_t(8*modmod());} void add_unsafe(UInt t){r+=t;} void pseudonormalize(){r=std::min(r,r-modmod8());} modint const&normalize(){ if(r>=(UInt)mod()){ r%=mod(); } return to_modint(); } void setr_direct(UInt rr){r=rr;} UInt getr_direct()const{return r;} protected: UInt r; private: constexpr modint&to_modint(){return static_cast<modint&>(*this);} constexpr modint const&to_modint()const{return static_cast<modint const&>(*this);} }; template<typename modint> concept modint_type=std::is_base_of_v<modint_base<modint,typename modint::Int>,modint>; template<modint_type modint> decltype(std::cin)&operator>>(decltype(std::cin)&in,modint&x){ typename modint::UInt r; auto&res=in>>r; x.setr(r); return res; } template<modint_type modint> decltype(std::cout)&operator<<(decltype(std::cout)&out,modint const&x){ return out<<x.getr(); } template<auto m> struct modint:modint_base<modint<m>,decltype(m)>{ using Base=modint_base<modint<m>,decltype(m)>; using Base::Base; static constexpr Base::Int mod(){return m;} static constexpr Base::UInt remod(){return m;} auto getr()const{return Base::r;} }; template<typename Int=int> struct dynamic_modint:modint_base<dynamic_modint<Int>,Int>{ using Base=modint_base<dynamic_modint<Int>,Int>; using Base::Base; [[gnu::noinline,gnu::cold]]static Base::UInt m_reduce_reduced(Base::UInt2 ab){ if(mod()%2==0){ return typename Base::UInt(ab%mod()); } typename Base::UInt q= -(typename Base::UInt(ab)*inverse); auto high=typename Base::UInt(ab>>Base::bits); auto low=typename Base::UInt(typename Base::UInt2(q)*typename Base::UInt(mod())>>Base::bits); return high>=low?high-low:high-low+mod(); } static Base::UInt m_reduce(Base::UInt2 ab){ if(imod()==0)[[unlikely]]{ return m_reduce_reduced(ab); } typename Base::UInt2 m=typename Base::UInt(ab)*imod(); return typename Base::UInt((ab+m*mod())>>Base::bits); } static Base::UInt m_transform(Base::UInt a){ if(mod()%2==0)[[unlikely]]{ return a; }else{ return m_reduce(a*pw128()); } } dynamic_modint&operator*=(const dynamic_modint&t){ Base::r=m_reduce(typename Base::UInt2(Base::r)*t.r); return*this; } void setr(Base::UInt rr){ Base::r=m_transform(rr); } Base::UInt getr()const{ typename Base::UInt res=m_reduce(Base::r); return std::min(res,res-mod()); } static Int mod(){return m;} static Base::UInt remod(){return rm;} static Base::UInt imod(){return im;} static Base::UInt2 pw128(){return r2;} static void switch_mod(Int nm){ m=nm; bool lazy=m%2&&typename Base::UInt(m)<=typename Base::UInt(-1)/4; rm=typename Base::UInt(m)*(lazy?2:1); inverse=m%2?inv2(-m):0; im=lazy?inverse:0; r2=static_cast<Base::UInt>(static_cast<Base::UInt2>(-1)%m+1); } auto static with_mod(Int tmp,auto callback){ struct scoped{ Int prev=mod(); ~scoped(){switch_mod(prev);} }_; switch_mod(tmp); return callback(); } private: static inline thread_local Int m=1; static inline thread_local Base::UInt im= -1,r2=0,inverse= -1,rm=2; }; }
#line 6 "cp-algo/number_theory/primality.hpp"
namespace cp_algo::math{ template<typename _Int> bool is_prime(_Int m){ using Int=std::make_signed_t<_Int>; using UInt=std::make_unsigned_t<Int>; if(m==1||m%2==0){ return m==2; } int s=std::countr_zero(UInt(m-1)); auto d=(m-1)>>s; using base=dynamic_modint<Int>; auto test= [&](base x){ x=bpow(x,d); if(std::abs(x.rem())<=1){ return true; } for(int i=1;i<s&&x!= -1;i++){ x*=x; } return x== -1; }; return base::with_mod(m,[&](){
#ifdef CP_ALGO_NUMBER_THEORY_PRIMALITY_BASES_HPP
uint16_t base2=7,base3=61; if(m!=uint32_t(m)){ base2=base_table1[uint32_t(m*0xAD625B89)>>18]; base3=base_table2[base2>>13]; } return test(2)&&test(base2)&&test(base3);
#else
return std::ranges::all_of(std::array{2,325,9375,28178,450775,9780504,1795265022},test);
#endif
}); } }
#line 1 "cp-algo/util/generator.hpp"
#line 1 "cp-algo/util/big_alloc.hpp"
#include <set>
#include <map>
#include <deque>
#include <stack>
#include <queue>
#line 10 "cp-algo/util/big_alloc.hpp"
#include <string>
#include <cstddef>
#line 14 "cp-algo/util/big_alloc.hpp"
#include <forward_list>
#if defined(__linux__) || defined(__unix__) || (defined(__APPLE__) && defined(__MACH__))
#define CP_ALGO_USE_MMAP 1
#include <sys/mman.h>
#else
#define CP_ALGO_USE_MMAP 0
#endif
namespace cp_algo{ template<typename T,size_t Align=32> class big_alloc{ static_assert(Align>=alignof(void*),"Align must be at least pointer-size"); static_assert(std::has_single_bit(Align),"Align must be a power of two"); public: using value_type=T; template<class U> struct rebind{ using other=big_alloc<U,Align>; }; constexpr bool operator==(const big_alloc&)const=default; big_alloc()noexcept=default; template<typename U,std::size_t A> big_alloc(const big_alloc<U,A>&)noexcept{} [[nodiscard]]T*allocate(std::size_t n){ std::size_t padded=round_up(n*sizeof(T)); std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){ void*raw=mmap(nullptr,padded,PROT_READ|PROT_WRITE,MAP_PRIVATE|MAP_ANONYMOUS,-1,0); madvise(raw,padded,MADV_HUGEPAGE); return static_cast<T*>(raw); }
#endif
return static_cast<T*>(::operator new(padded,std::align_val_t(align))); } void deallocate(T*p,std::size_t n)noexcept{ if(!p){ return; } std::size_t padded=round_up(n*sizeof(T)); std::size_t align=std::max<std::size_t>(alignof(T),Align);
#if CP_ALGO_USE_MMAP
if(padded>=MEGABYTE){ munmap(p,padded); return; }
#endif
::operator delete(p,padded,std::align_val_t(align)); } private: static constexpr std::size_t MEGABYTE=1<<20; static constexpr std::size_t round_up(std::size_t x)noexcept{ return(x+Align-1)/Align*Align; } }; template<typename T> using big_vector=std::vector<T,big_alloc<T>>; template<typename T> using big_basic_string=std::basic_string<T,std::char_traits<T>,big_alloc<T>>; template<typename T> using big_deque=std::deque<T,big_alloc<T>>; template<typename T> using big_stack=std::stack<T,big_deque<T>>; template<typename T> using big_queue=std::queue<T,big_deque<T>>; template<typename T> using big_priority_queue=std::priority_queue<T,big_vector<T>>; template<typename T> using big_forward_list=std::forward_list<T,big_alloc<T>>; using big_string=big_basic_string<char>; template<typename Key,typename Value,typename Compare=std::less<Key>> using big_map=std::map<Key,Value,Compare,big_alloc<std::pair<const Key,Value>>>; template<typename T,typename Compare=std::less<T>> using big_multiset=std::multiset<T,Compare,big_alloc<T>>; template<typename T,typename Compare=std::less<T>> using big_set=std::set<T,Compare,big_alloc<T>>; }
#line 5 "cp-algo/util/generator.hpp"
#include <generator>
namespace cp_algo{ template<typename Ref,typename V=void> using big_generator=std::generator<Ref,V,big_alloc<std::byte>>; } namespace std::ranges{ template<typename Ref,typename V> elements_of(cp_algo::big_generator<Ref,V>&&)->elements_of<cp_algo::big_generator<Ref,V>&&,cp_algo::big_alloc<std::byte>>; }
#line 1 "cp-algo/random/rng.hpp"
#include <chrono>
#include <random>
namespace cp_algo::random{ std::mt19937_64 gen( std::chrono::steady_clock::now().time_since_epoch().count() ); uint64_t rng(){ return gen(); } }
#line 6 "cp-algo/number_theory/factorize.hpp"
namespace cp_algo::math{ template<typename _Int> auto proper_divisor(_Int m){ using Int=std::make_signed_t<_Int>; using base=dynamic_modint<Int>; return m%2==0?2:base::with_mod(m,[&](){ base t=random::rng(); auto f= [&](auto x){ return x*x+t; }; const base zero=0,one=1; base x=zero,y=zero; base g=one; while(g==one){ for(int i=0;i<64;i++){ x=f(x); y=f(f(y)); if(x==y)[[unlikely]]{ t=random::rng(); x=y=zero; }else{ base t=g*(x-y); g=t==zero?g:t; } } g=std::gcd(g.getr(),m); } return g.getr(); }); } template<typename Int> big_generator<Int>factorize(Int m){ if(is_prime(m)){ co_yield m; }else if(m>1){ auto g=proper_divisor(m); co_yield std::ranges::elements_of(factorize(g)); co_yield std::ranges::elements_of(factorize(m/g)); } } template<typename Int> big_generator<Int>divisors_sqrt(Int m){ for(Int i=1;i*i<=m;i++){ if(m%i==0){ co_yield i; if(i*i!=m){ co_yield m/i; } } } } }
Back to top page