Compiler Explorer

Source code

#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <string>
#include <random>
#include <cmath>
#include <tuple>
#include <future>
#include <thread>
#include <optional>
#include <set>
#include <map>
#include <unordered_set>
#include <chrono>
#include <type_traits>
#include <atomic>

#define SUPPORT_GMP 1
#define SUPPORT_BOOST 1

#if SUPPORT_GMP
    #include <gmpxx.h>
#endif
#if SUPPORT_BOOST
    #include <boost/multiprecision/cpp_int.hpp>
#endif

#if !SUPPORT_BOOST && (defined(__clang__) || defined(__GNUC__))
    #define SUPPORT_STD_U128 1
#else
    #define SUPPORT_STD_U128 0
#endif

#define ASSERT_MSG(cond, msg) { if (!(cond)) throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg: '" + std::string(msg) + "'."); }
#define ASSERT(cond) ASSERT_MSG(cond, "")
#define DASSERT_MSG(cond, msg) //ASSERT_MSG(cond, msg)
#define DASSERT(cond) DASSERT_MSG(cond, "")
#define OSTR(code) ([&]{ std::ostringstream ss; ss code; return ss.str(); }())
#define COUT(code) { std::lock_guard lock(cout_mux); std::cout code; std::cout << std::flush; }
#define DUMP(var) //{ COUT(<< __LINE__ << ": " << #var << " = (" << (var) << ")" << std::endl); }
#define LN //{ COUT(<< "LN " << __LINE__ << std::endl); }
#define bisizeof(x) (sizeof(x) * 8)

using std::size_t;

using u8  = uint8_t;
using i8  = int8_t;
using u16 = uint16_t;
using i16 = int16_t;
using u32 = uint32_t;
using i32 = int32_t;
using u64 = uint64_t;
using i64 = int64_t;
#if SUPPORT_STD_U128
    using i128 = signed __int128;
    using u128 = unsigned __int128;
#endif
#if SUPPORT_BOOST
    using u128b = boost::multiprecision::uint128_t;
    using i128b = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<128, 128, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    #if !SUPPORT_STD_U128
        using u128 = u128b;
        using i128 = i128b;
    #endif
    using u192 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<192, 192, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i192 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<192, 192, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u256 = boost::multiprecision::uint256_t;
    using i256 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<256, 256, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u384 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<384, 384, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i384 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<384, 384, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u512 = boost::multiprecision::uint512_t;
    using i512 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<512, 512, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u768 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<768, 768, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i768 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<768, 768, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u1024 = boost::multiprecision::uint1024_t;
    using i1024 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<1024, 1024, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u1536 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<1536, 1536, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i1536 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<1536, 1536, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u2048 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<2048, 2048, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i2048 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<2048, 2048, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u3072 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<3072, 3072, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i3072 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<3072, 3072, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
    using u4096 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<4096, 4096, boost::multiprecision::unsigned_magnitude, boost::multiprecision::unchecked, void>>;
    using i4096 = boost::multiprecision::number<boost::multiprecision::cpp_int_backend<4096, 4096, boost::multiprecision::signed_magnitude, boost::multiprecision::unchecked, void>>;
#endif

std::recursive_mutex cout_mux;

template <typename T> struct DWordOf;
template <typename T> using DWordOfT = typename DWordOf<T>::type;

template <> struct DWordOf<u32> : std::type_identity<u64> {};
template <> struct DWordOf<i32> : std::type_identity<i64> {};
#if SUPPORT_STD_U128 || SUPPORT_BOOST
    template <> struct DWordOf<u64> : std::type_identity<u128> {};
    template <> struct DWordOf<i64> : std::type_identity<i128> {};
#endif
#if SUPPORT_BOOST
    template <> struct DWordOf<u128b> : std::type_identity<u256> {};
    template <> struct DWordOf<u192> : std::type_identity<u384> {};
    template <> struct DWordOf<u256> : std::type_identity<u512> {};
    template <> struct DWordOf<u384> : std::type_identity<u768> {};
    template <> struct DWordOf<u512> : std::type_identity<u1024> {};
    template <> struct DWordOf<u768> : std::type_identity<u1536> {};
    template <> struct DWordOf<u1024> : std::type_identity<u2048> {};
    template <> struct DWordOf<u1536> : std::type_identity<u3072> {};
    template <> struct DWordOf<u2048> : std::type_identity<u4096> {};

template <> struct DWordOf<i128b> : std::type_identity<i256> {};
    template <> struct DWordOf<i192> : std::type_identity<i384> {};
    template <> struct DWordOf<i256> : std::type_identity<i512> {};
    template <> struct DWordOf<i384> : std::type_identity<i768> {};
    template <> struct DWordOf<i512> : std::type_identity<i1024> {};
    template <> struct DWordOf<i768> : std::type_identity<i1536> {};
    template <> struct DWordOf<i1024> : std::type_identity<i2048> {};
    template <> struct DWordOf<i1536> : std::type_identity<i3072> {};
    template <> struct DWordOf<i2048> : std::type_identity<i4096> {};
#endif
#if SUPPORT_GMP
    template <> struct DWordOf<mpz_class> : std::type_identity<mpz_class> {};
#endif

template <typename T> struct SignedOf;
template <typename T> using SignedOfT = typename SignedOf<T>::type;

template <> struct SignedOf<u32> : std::type_identity<i32> {};
template <> struct SignedOf<u64> : std::type_identity<i64> {};
#if SUPPORT_STD_U128
    template <> struct SignedOf<u128> : std::type_identity<i128> {};
#endif
#if SUPPORT_BOOST
    template <> struct SignedOf<u128b> : std::type_identity<i128b> {};
    template <> struct SignedOf<u192> : std::type_identity<i192> {};
    template <> struct SignedOf<u256> : std::type_identity<i256> {};
    template <> struct SignedOf<u384> : std::type_identity<i384> {};
    template <> struct SignedOf<u512> : std::type_identity<i512> {};
    template <> struct SignedOf<u768> : std::type_identity<i768> {};
    template <> struct SignedOf<u1024> : std::type_identity<i1024> {};
    template <> struct SignedOf<u1536> : std::type_identity<i1536> {};
    template <> struct SignedOf<u2048> : std::type_identity<i2048> {};
#endif
#if SUPPORT_GMP
    template <> struct SignedOf<mpz_class> : std::type_identity<mpz_class> {};
#endif

template <typename T> struct IsGmp : std::false_type {};
template <typename T> bool constexpr IsGmpV = IsGmp<T>::value;
#if SUPPORT_GMP
    template <> struct IsGmp<mpz_class> : std::true_type {};
#endif

template <typename T> struct BiSizeS;
template <typename T> size_t constexpr BiSize = BiSizeS<T>::value;

template <> struct BiSizeS<u32> : std::integral_constant<size_t, 32> {};
template <> struct BiSizeS<u64> : std::integral_constant<size_t, 64> {};
#if SUPPORT_STD_U128
    template <> struct BiSizeS<u128> : std::integral_constant<size_t, 128> {};
#endif
#if SUPPORT_BOOST
    template <> struct BiSizeS<u128b> : std::integral_constant<size_t, 128> {};
    template <> struct BiSizeS<u192> : std::integral_constant<size_t, 192> {};
    template <> struct BiSizeS<u256> : std::integral_constant<size_t, 256> {};
    template <> struct BiSizeS<u384> : std::integral_constant<size_t, 384> {};
    template <> struct BiSizeS<u512> : std::integral_constant<size_t, 512> {};
    template <> struct BiSizeS<u768> : std::integral_constant<size_t, 768> {};
    template <> struct BiSizeS<u1024> : std::integral_constant<size_t, 1024> {};
    template <> struct BiSizeS<u1536> : std::integral_constant<size_t, 1536> {};
    template <> struct BiSizeS<u2048> : std::integral_constant<size_t, 2048> {};
#endif
#if SUPPORT_GMP
    template <> struct BiSizeS<mpz_class> : std::integral_constant<size_t, 1ULL << 30> {};
#endif

template <size_t Bits>
using MinType =
    std::conditional_t<Bits <= 32, u32,
    std::conditional_t<Bits <= 64, u64,
#if SUPPORT_BOOST
    std::conditional_t<Bits <= 128, u128b,
    std::conditional_t<Bits <= 192, u192,
    std::conditional_t<Bits <= 256, u256,
    std::conditional_t<Bits <= 384, u384,
    std::conditional_t<Bits <= 512, u512,
    std::conditional_t<Bits <= 768, u768,
    std::conditional_t<Bits <= 1024, u1024,
    std::conditional_t<Bits <= 1536, u1536,
    std::conditional_t<Bits <= 2048, u2048, void>
    >>>>>>>>
#endif
    >>;

template <typename RT, typename T>
auto ToT(T const & x) {
    if constexpr(std::is_same_v<RT, T>)
        return x;
    else
        return RT(x);
}

template <typename T>
T FromU64(u64 x) {
    if constexpr(IsGmpV<T>)
        return (T(u32(x >> 32)) << 32) | u32(x);
    else
        return T(x);
}

template <typename T>
u32 ToU32(T const & x) {
    if constexpr(IsGmpV<T>)
        return x.get_ui();
    else
        return u32(x);
}

template <typename T>
u64 ToU64(T const & x) {
    if constexpr(IsGmpV<T>)
        return (u64(ToU32<T>(x >> 32)) << 32) | ToU32(x);
    else
        return u64(x);
}

template <typename T = u32>
T GetPrime(size_t i) {
    thread_local std::vector<T> primes = {2, 3};
    while (i >= primes.size())
        for (T p = primes.back() + 2;; p += 2) {
            bool is_prime = true;
            for (auto const & d: primes) {
                if (d * d > p)
                    break;
                if (p % d == 0) {
                    is_prime = false;
                    break;
                }
            }
            if (is_prime) {
                primes.push_back(p);
                break;
            }
        }
    return primes[i];
}

template <typename T = u32>
size_t NextPrimeIdx(T begin) {
    thread_local std::vector<u32> idxs = {0, 0, 1};
    if (begin >= idxs.size())
        for (size_t ip = idxs.back(); begin >= idxs.size(); ++ip) {
            auto const p = GetPrime<T>(ip);
            idxs.resize(p + 1);
            idxs.back() = ip + 1;
        }
    for (size_t i = begin;; ++i)
        if (idxs[i] != 0)
            return idxs[i] - 1;
}

template <typename T>
std::pair<std::vector<T>, bool> FactorTrialDiv(T n, u64 end = u32(-1), u64 begin = 0) {
    // https://en.wikipedia.org/wiki/Trial_division
    ASSERT(end <= u32(-1));
    u64 n64 = 0;
    if constexpr(BiSize<T> <= 64)
        n64 = u64(n);
    else 
        n64 = n < FromU64<T>(u64(-1) >> 2) ? ToU64(n) : (u64(-1) >> 2);
    std::vector<T> fs;
    while ((n & 1) == 0) {
        fs.push_back(2);
        n >>= 1;
    }
    begin = std::max<size_t>(begin, 3);
    bool fully = true;
    for (size_t ip = NextPrimeIdx<u32>(begin);; ++ip) {
        u32 const d = GetPrime<u32>(ip);
        if (u64(d) * d > n64)
            break;
        if (d >= end) {
            fully = false;
            break;
        }
        while (n % d == 0) {
            fs.push_back(d);
            n /= d;
        }
    }
    if (n > 1)
        fs.push_back(n);
    return {std::move(fs), fully};
}

template <typename T>
auto DTMul(T const & a, T const & b) {
    return ToT<DWordOfT<T>>(a) * b;
}

template <typename T>
T MulMod(T const & a, T const & b, T const & c) {
    return ToT<T>((ToT<DWordOfT<T>>(a) * b) % c);
}

template <typename T>
T PowMod(T a, T b, T const & c) {
    T r = 1;
    while (b != 0) {
        if ((b & 1) != 0)
            r = MulMod(r, a, c);
        a = MulMod(a, a, c);
        b >>= 1;
    }
    return r;
}

u64 RandU64() {
    thread_local std::mt19937_64 rng{
        std::random_device{}()
        //123
    };
    return rng();
}

template <typename T>
size_t NumBits(T x) {
    using W = std::conditional_t<IsGmpV<T>, u32, u64>;
    size_t r = 0;
    if constexpr(BiSize<T> > bisizeof(W))
        while (x > W(-1)) {
            r += bisizeof(W);
            x >>= bisizeof(W);
        }
    while (x >= u8(-1)) {
        r += 8;
        x >>= 8;
    }
    while (x > 0) {
        ++r;
        x >>= 1;
    }
    return r;
}

template <typename T>
inline T Shl(T const & x, size_t Bits) {
    if constexpr(IsGmpV<T>)
        return x << Bits;
    else
        return Bits >= BiSize<T> ? ToT<T>(0) : ToT<T>(x << Bits);
}

template <typename T>
T RandRange(T const & begin, T const & end) {
    DASSERT(begin < end);
    using W = std::conditional_t<IsGmpV<T>, u32, u64>;
    T const len = end - begin;
    size_t const nbits = NumBits(len);
    while (true) {
        T r = 0;
        for (size_t i = 0; i < nbits; i += bisizeof(W)) {
            size_t const portion = std::min<size_t>(bisizeof(W), nbits - i);
            if constexpr(IsGmpV<T>)
                r <<= portion;
            else
                r = Shl(r, portion);
            r ^= W(RandU64() & (u64(-1) >> (64 - portion)));
        }
        if (r < len) {
            return begin + r;
        }
    }
}

template <typename T>
T RandPrime(size_t bits) {
    ASSERT(bits >= 3);
    while (true) {
        T const p = T(RandRange<T>(T(1) << (bits - 1), T(1) << bits) | 1);
        if (IsProbablyPrime<T>(p)) {
            return p;
        }
    }
}

double Time() {
    static auto const gtb = std::chrono::high_resolution_clock::now();
    return std::chrono::duration_cast<std::chrono::duration<double>>(
        std::chrono::high_resolution_clock::now() - gtb).count();
}

template <typename T>
bool FermatPRP(T const & n, size_t ntrials = 32) {
    // https://en.wikipedia.org/wiki/Fermat_primality_test
    for (size_t i = 0; i < ntrials; ++i) {
        if (PowMod<T>(RandRange<T>(2, n - 1), n - 1, n) != 1)
            return false;
    }
    return true;
}

template <typename T>
bool IsProbablyPrime(T const & n) {
    if (n <= 16)
        return n == 2 || n == 3 || n == 5 || n == 7 || n == 11 || n == 13;
    if (n < (1 << 20))
        return FactorTrialDiv<u32>(ToU32(n)).first.size() == 1;
    if (FactorTrialDiv<T>(n, 1 << 5).first.size() >= 2)
        return false;
    return FermatPRP<T>(n);
}

template <typename T>
std::vector<T> GenPrimes(u64 const & end) {
    // https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
    std::vector<bool> composite(end / 2);
    for (size_t p = 3; p * p < end; p += 2) {
        if (composite[p / 2])
            continue;
        for (size_t i = p * p; i < end; i += p * 2)
            composite[i / 2] = true;
    }
    std::vector<T> primes = {2};
    for (size_t i = 3; i < end; i += 2)
        if (!composite[i / 2])
            primes.push_back(T(i));
    primes.shrink_to_fit();
    return primes;
}

template <typename T>
T GCD(T const & a0, T const & b0) {
    // https://en.wikipedia.org/wiki/Euclidean_algorithm
    DASSERT(a0 >= 0 && b0 >= 0);
    T a = a0, b = b0;
    while (b != 0)
        std::tie(a, b) = std::make_tuple(b, ToT<T>(a % b));
    DASSERT_MSG(a > 0 && a0 % a == 0 && b0 % a == 0, "a0 " + IntToStr(a0) + " b0 " + IntToStr(b0) +
        " a " + IntToStr(a) + " a0 % a " + IntToStr<T>(a0 % a) + " b0 % b " + IntToStr<T>(b0 % b));
    return a;
}

template <typename T, typename ST = SignedOfT<T>>
std::tuple<T, ST, ST> EGCD(T const & a, T const & b) {
    // https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
    using DST = DWordOfT<ST>;
    DASSERT(a >= 0 && b >= 0);
    T ro = a, r = b, q = 0, tmp0 = 0;
    ST so = 1, s = 0, tmp1 = 0;
    while (r != 0) {
        q = ro / r;
        tmp0 = r;
        r = ToT<T>(ro - q * r);
        ro = tmp0;
        tmp1 = s;
        s = ToT<ST>(so - ToT<ST>(q) * s);
        so = tmp1;
    }
    ST const to = ToT<ST>((ToT<DST>(ro) - ToT<DST>(a) * so) / ToT<ST>(b));
    DASSERT(ro == ToT<DST>(a) * so + ToT<DST>(b) * to);
    return std::make_tuple(ro, so, to);
}

template <typename T>
class ModInvExc : std::runtime_error {
public:
    ModInvExc() = delete;
    ModInvExc(T const & _g)
        : std::runtime_error("ModInvExc"), g(_g) {}
    T g = 0;
};

template <typename T>
T ModInv(T a, T const & n) {
    // https://en.wikipedia.org/wiki/Modular_multiplicative_inverse
    DASSERT(a >= 0);
    a %= n;
    auto [g, s, t] = EGCD<T>(a, n);
    using ST = std::decay_t<decltype(s)>;
    if (g != 1)
        throw ModInvExc<T>(g);
    return ToT<T>(ToT<ST>(n) + s) % n;
}

template <typename T>
std::pair<T, T> EcAdd(T const & N, T const & A, T const & B,
        T const & X0, T const & Y0, T const & X1, T const & Y1) {
    // https://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
    if (X0 == X1 && Y0 == Y1) {
        T const l = MulMod<T>((MulMod<T>(X0, X0, N) * 3 + A), ModInv<T>(Y0 * 2, N), N);
        T const x = (MulMod<T>(l, l, N) + N * 2 - X0 * 2) % N;
        T const y = (MulMod<T>(l, X0 + N - x, N) + N - Y0) % N;
        return std::make_pair(x, y);
    } else {
        T const l = MulMod<T>(Y1 + N - Y0, ModInv<T>(X1 + N - X0, N), N);
        T const x = (MulMod<T>(l, l, N) + N * 2 - X0 - X1) % N;
        T const y = (MulMod<T>(l, X0 + N - x, N) + N - Y0) % N;
        return std::make_pair(x, y);
    }
}

template <typename T, typename KT>
std::pair<T, T> EcMul(T const & N, T const & A, T const & B, T X, T Y, KT k) {
    // https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
    DASSERT(k >= 1);
    if (k == 1)
        return {X, Y};
    --k;
    T BX = X, BY = Y;
    while (k != 0) {
        if ((k & 1) != 0) {
            std::tie(X, Y) = EcAdd<T>(N, A, B, X, Y, BX, BY);
        }
        std::tie(BX, BY) = EcAdd<T>(N, A, B, BX, BY, BX, BY);
        k >>= 1;
    }
    return std::make_pair(X, Y);
}

template <typename T>
void SplitFactors(std::set<T> & fs) {
    fs.erase(1);
    while (true) {
        std::vector<T> fsv(fs.begin(), fs.end());
        std::set<T> fsn;
        bool splited = false;
        size_t i = 0, j = 0;
        for (i = 0; i < fsv.size(); ++i) {
            for (j = i + 1; j < fsv.size(); ++j) {
                T const g = GCD<T>(fsv[i], fsv[j]);
                if (g <= 1)
                    continue;
                fsn.insert(g);
                fsn.insert(fsv[i] / g);
                fsn.insert(fsv[j] / g);
                splited = true;
                break;
            }
            if (splited)
                break;
        }
        if (!splited)
            break;
        fsn.erase(1);
        fsn.insert(fsv.begin() +       0, fsv.begin() +          i);
        fsn.insert(fsv.begin() + (i + 1), fsv.begin() +          j);
        fsn.insert(fsv.begin() + (j + 1), fsv.begin() + fsv.size());
        fs = std::move(fsn);
    }
}

template <typename T>
std::string IntToStr(T x) {
    bool const neg = x < 0;
    if constexpr(std::is_same_v<T, SignedOfT<T>>)
        if (neg)
            x = -x;
    std::string res;
    u32 constexpr mod = 1'000'000'000;
    while (x != 0) {
        res = std::to_string(ToU32<T>(x % mod)) + res;
        x /= mod;
    }
    return (neg ? "-" : "") + res;
}

template <typename ItT>
std::string IntsToStr(ItT begin, ItT end) {
    std::ostringstream ss;
    ss << "[";
    for (auto it = begin; it != end; ++it) {
        ss << IntToStr(*it);
        if (std::next(it) != end)
            ss << ", ";
    }
    ss << "]";
    return ss.str();
}

template <typename T, typename NT>
T NumPower(NT const & p, T const & bound2) {
    T mp = p, t = 0;
    while (true) {
        t = mp * p;
        if (t >= bound2)
            return mp;
        mp = t;
    }
}

template <typename T>
constexpr T ParseNum(char const * s) {
    size_t len = 0;
    for (len = 0; s[len]; ++len);
    T r = 0;
    for (size_t i = 0; i < len; ++i) {
        r *= 10;
        r += s[i] - '0';
    }
    return r;
}

static double ecm_succ_avg = 0, ecm_succ_avg_cnt = 0, ecm_stats_reported = 0,
    tdiv_time_avg = 0, tdiv_time_avg_cnt = 0, tdiv_fully = 0, tdiv_fully_tot = 0;
static std::string ecm_stats, tdiv_stats, total_stats;

template <typename T>
std::pair<std::vector<T>, bool> FactorECM(T const & N0, bool verbose = false, size_t bound_inc = 1 << 9, size_t tdiv_limit = 1 << 14, double max_time = 1ULL << 30) {
    // https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization
    // https://stackoverflow.com/a/75765566/941531
    auto total_time = Time();
    
    size_t const tdiv_endA = tdiv_limit, tdiv_endB = tdiv_endA * 4;
    size_t const nthreads = std::thread::hardware_concurrency();
    size_t const ncurves_per_bound = (std::max<size_t>(1 << 3, nthreads * 3) + nthreads - 1) / nthreads * nthreads;
    
    if (verbose) {
        COUT(<< "Initial N to factor with ECM: " << IntToStr(N0) << " (" << NumBits(N0) << "-bit)" << std::endl
             << "Number of threads " << nthreads << "." << std::endl);
    }
    
    size_t const start_bound = bound_inc;
    
    std::map<T, bool> is_prime;
    std::unordered_set<u64> tdiv_done;
    auto IsPrime = [&](T const & x){
        auto it = is_prime.find(x);
        if (it != is_prime.end())
            return it->second;
        return is_prime.insert(std::pair<T, bool>{x, IsProbablyPrime<T>(x)}).first->second;
    };
    
    std::set<T> factors;
    bool needs_main = true, timed_out = false;
    
    {
        auto tdiv_time = Time();
        auto const [fs, fully] = FactorTrialDiv<T>(N0, tdiv_endA);
        tdiv_time = Time() - tdiv_time;
        if (fs.size() == 1) {
            ASSERT(fs.front() == N0);
            factors.insert(N0);
        } else {
            ASSERT(fs.size() >= 1);
            if (verbose) {
                COUT(<< "Factors TrialDivA: " << IntsToStr(fs.begin(), fs.begin() + (fully ? fs.size() : fs.size() - 1)) << std::endl);
            }
            factors.insert(fs.begin(), fs.end());
            SplitFactors(factors);
        }
        tdiv_fully_tot += 1;
        if (fully) {
            for (size_t i = 0; i < fs.size(); ++i)
                is_prime[fs[i]] = true;
            needs_main = false;
        } else
            for (size_t i = 0; i + 1 < fs.size(); ++i)
                is_prime[fs[i]] = true;
        {
            tdiv_time_avg += tdiv_time;
            tdiv_time_avg_cnt += 1;
            if (std::llround(tdiv_time_avg_cnt) % 50 == 0)
                tdiv_stats = OSTR(<< std::fixed << std::setprecision(2) << std::setw(6) << 100.0 * tdiv_fully / tdiv_fully_tot) +
                    "% " + OSTR(<< std::setw(4) << std::llround(1'000'000 * tdiv_time_avg / tdiv_time_avg_cnt)) +"mcs";
        }
    }
    
    static T const limit64 = FromU64<T>(u64(-10));
    bool used_ecm = false;
    
    for (size_t imain = 0;; ++imain) {
        if (!needs_main)
            break;
        
        T N = 0;
        for (auto const & f: factors)
            if (!IsPrime(f)) {
                N = f;
                break;
            }
        if (N == 0)
            break;
        
        if (verbose) {
            COUT(<< "Factoring " << IntToStr(N) << std::endl);
        }
        
        if (N < limit64) {
            u64 const N64 = ToU64(N);
            if (!tdiv_done.count(N64)) {
                auto const [fs, fully] = FactorTrialDiv<u64>(N64, tdiv_endB, tdiv_endA);
                tdiv_done.insert(N64);
                if (fs.size() >= 2) {
                    ASSERT(fs.size() >= 1);
                    if (verbose) {
                        COUT(<< "Factors TrialDivB: " << IntsToStr(fs.begin(), fs.begin() + (fully ? fs.size() : fs.size() - 1)) << std::endl);
                    }
                    for (auto const & f: fs)
                        factors.insert(FromU64<T>(f));
                    SplitFactors(factors);
                    continue;
                }
            }
        }
        
        size_t bound = start_bound;
        
        for (size_t icycle = 0, icurve = 0;; ++icycle, bound += bound_inc, icurve += ncurves_per_bound) {
            if (Time() - total_time >= max_time) {
                timed_out = true;
                break;
            }
            
            if (verbose) {
                COUT(<< "Curves [" << std::setw(4) << icurve << ", " << std::setw(4) << (icurve + ncurves_per_bound) << "), bound 2^"
                    << std::fixed << std::setprecision(3) << std::setw(7) << std::log2(bound));
            }

auto ecm_time = Time();
            used_ecm = true;
            auto const primes = GenPrimes<u32>(bound);
            
            static std::atomic<double> avg_ecm_time = 0, avg_ecm_time_cnt = 0;
            std::vector<std::future<std::optional<T>>> asyncs;
            for (size_t ithr = 0; ithr < ncurves_per_bound; ++ithr)
                asyncs.push_back(std::async([&, ithr]()->std::optional<T>{
                    auto ecm_time = Time();
                    T X = 0, Y = 0, A = 0, B = 0;
                    while (true) {
                        std::tie(X, Y, A) = std::make_tuple(RandRange<T>(0, N), RandRange<T>(0, N), RandRange<T>(0, N));
                        B = T((DTMul<T>(Y, Y) + DTMul<T>(N, N) * 4 - DTMul<T>(X, MulMod<T>(X, X, N)) - DTMul<T>(A, X)) % N);
                        if ((DTMul<T>(A, MulMod<T>(A, A, N)) * 4 + DTMul<T>(N, N) * 27 - DTMul<T>(B, B) * 27) % N == 0)
                            continue;
                        break;
                    }
                    
                    size_t const bound2 = size_t(std::sqrt(bound) + 1.01);
                    
                    for (size_t iprime = 0; iprime < primes.size(); ++iprime) {
                        auto const p = primes[iprime];
                        auto const k = NumPower<u64>(p, bound2);
                        try {
                            std::tie(X, Y) = EcMul(N, A, B, X, Y, k);
                        } catch (ModInvExc<T> const & ex) {
                            auto const g = ex.g;
                            ASSERT_MSG(g > 1, "g " + IntToStr<T>(g));
                            if (g == N)
                                break;
                            return g;
                        }
                    }
                    
                    ecm_time = Time() - ecm_time;
                    avg_ecm_time += ecm_time;
                    avg_ecm_time_cnt += 1;
                    
                    return std::nullopt;
                }));
            std::set<T> fsn;
            bool any_factored = false, any_unfactored = false;
            for (auto & e: asyncs) {
                auto const r = e.get();
                if (r) {
                    ecm_succ_avg += 1;
                    fsn.insert(*r);
                }
                ecm_succ_avg_cnt += 1;
            }
            ecm_time = Time() - ecm_time;
            if (verbose) {
                COUT(<< ", " << std::fixed << std::setprecision(3) << ecm_time << " sec" << std::endl);
            }
            if (avg_ecm_time_cnt - ecm_stats_reported >= 50) {
                ecm_stats = OSTR(<< std::fixed << std::setprecision(1) << std::setw(4) << 100.0 * ecm_succ_avg / ecm_succ_avg_cnt) + "% " +
                    OSTR(<< std::setw(4) << std::llround(1'000'000 * avg_ecm_time / avg_ecm_time_cnt)) + "mcs";
                ecm_stats_reported = avg_ecm_time_cnt;
            }
            if (!fsn.empty()) {
                SplitFactors(fsn);
                if (verbose) {
                    COUT(<< "Factors ECM: " << IntsToStr(fsn.begin(), fsn.end()) << std::endl);
                }
                factors.insert(fsn.begin(), fsn.end());
                SplitFactors(factors);
                break;
            }
        }
        if (timed_out)
            break;
    }
    
    std::vector<T> all_factors;
    
    {
        auto N = N0;
        for (auto const & f: factors) {
            ASSERT(N % f == 0);
            do {
                all_factors.push_back(f);
                N /= f;
            } while (N % f == 0);
        }
        ASSERT(N == 1);
    }
    
    bool fs_all_primes = true;
    for (auto const & f: all_factors)
        if (!IsPrime(f)) {
            fs_all_primes = false;
            break;
        }
    if (fs_all_primes && !used_ecm)
        tdiv_fully += 1;
    if (!timed_out) {
        ASSERT_MSG(fs_all_primes, "N0 " + IntToStr(N0) + " factors " + IntsToStr(all_factors.begin(), all_factors.end()));
    }
    
    total_time = Time() - total_time;
    
    if (!timed_out) {
        static double tot_time_avg = 0, tot_time_avg_cnt = 0;
        tot_time_avg += total_time;
        tot_time_avg_cnt += 1;
        total_stats = OSTR(<< std::setw(4) << std::llround(1'000'000 * tot_time_avg / tot_time_avg_cnt)) + "mcs";
    }
    
    if (verbose && !timed_out) {
        COUT(<< "Fully factored. N " << IntToStr(N0) << ": " << IntsToStr(all_factors.begin(), all_factors.end())
            << std::endl << "Time " << std::fixed << std::setprecision(3) << total_time << " sec." << std::endl);
    }
    
    return {all_factors, !timed_out};
}

template <typename T>
T SumOfProperDivisors(T const & Na, std::vector<T> const & fsA) {
    std::map<T, size_t> fsm;
    for (auto const & f: fsA)
        ++fsm[f];
    std::vector<std::pair<T, size_t>> fsmv(fsm.begin(), fsm.end());
    T Nb = 0;
    for (size_t i = 0;; ++i) {
        T prod = 1;
        auto j = i;
        for (auto const & [f, c]: fsmv) {
            for (size_t k = 0, kend = j % (c + 1); k < kend; ++k)
                prod *= f;
            j /= c + 1;
        }
        if (j != 0)
            break;
        if (prod >= Na)
            break;
        Nb += prod;
    }
    return Nb;
}

double AmicableChance(size_t bits) {
    // https://sech.me/ap/
    double const
        ref_digits = 20, ref_cnt = 2249600, ref_coeff = 2249600.0 / 1046355.0,
        num = std::exp2(double(bits) - 0.5),
        cnt = ref_cnt / std::pow(ref_coeff, ref_digits - std::log10(num)),
        prob = cnt / std::exp2(bits);
    return prob;
}

template <typename T>
void FindAmicable(size_t bits_amicable_min, size_t bits_amicable_max) {
    // https://en.wikipedia.org/wiki/Amicable_numbers
    ASSERT(bits_amicable_min <= bits_amicable_max);
    ASSERT(bits_amicable_max + 5 <= BiSize<T>);
    double const amicable_prob = AmicableChance(bits_amicable_max);
    ecm_succ_avg = 0; ecm_succ_avg_cnt = 0; ecm_stats_reported = 0;
    tdiv_time_avg = 0; tdiv_time_avg_cnt = 0; tdiv_fully = 0; tdiv_fully_tot = 0;
    size_t dm = 0, ds = 0, dl = 0, dots_len = 0, dl_before = 0, min_prec = 1000;
    double report_time = -1000, factor_time_avg = 0, factor_time_avg_cnt = 0;
    auto Dots = [&]{
        auto const tdelta = Time() - report_time;
        size_t const prec = std::max<int>(2, -std::llround(std::floor(std::log10(
            std::max(1.0e-20, amicable_prob * 4.0 * (ptrdiff_t(dl) - ptrdiff_t(dl_before))))))) - 2;
        min_prec = std::min(min_prec, prec);
        std::ostringstream ss;
        bool constexpr show_tdiv_ecm_stats = 1;
        ss  << "  ." << std::left << std::setw(8) << dm << " _" << std::setw(8) << ds << " *" << std::setw(8) << dl
            << ", " << std::right << std::setw(6) << std::llround(factor_time_avg_cnt <= 0.01 ? 0 : factor_time_avg_cnt / factor_time_avg) << " nums/sec";
        if constexpr(show_tdiv_ecm_stats)
            ss << ", tdiv " << tdiv_stats << ", ecm " << ecm_stats;
        ss << ", total " << total_stats << ", done " << std::fixed << std::setprecision(min_prec)
            << std::right << std::setw(3) << 100.0 * amicable_prob * dl << "%    ";
        dl_before = dl;
        return ss.str();
    };
    auto ShowDots = [&]{
        if (Time() - report_time <= 1.0)
            return;
        COUT(<< std::string(dots_len, '\x08'));
        dots_len = 0;
        auto const dots = Dots();
        dots_len = dots.size();
        COUT(<< dots);
        report_time = Time();
    };
    std::set<T> done;
    while (true) {
        T Na = RandRange<T>(T(1) << (bits_amicable_min - 1), T(1) << bits_amicable_max);
        auto tb = Time();
        auto const [fsA, fullyA] = FactorECM<T>(Na, false, 1 << 5, bits_amicable_max <= 32 ? (1 << 8) : (1 << 14), 0.01);
        tb = Time() - tb;
        if (fullyA) { factor_time_avg += tb; factor_time_avg_cnt += 1; }
        ++dm;
        ShowDots();
        if (!fullyA || fsA.size() <= 1)
            continue;
        T Nb = SumOfProperDivisors<T>(Na, fsA);
        tb = Time();
        auto const [fsB, fullyB] = FactorECM<T>(Nb, false, 1 << 5, bits_amicable_max <= 32 ? (1 << 8) : (1 << 14), 0.05);
        tb = Time() - tb;
        if (fullyB) { factor_time_avg += tb; factor_time_avg_cnt += 1; }
        ++ds;
        ShowDots();
        if (!fullyB || fsB.size() <= 1)
            continue;
        ++dl;
        ShowDots();
        T const Nc = SumOfProperDivisors<T>(Nb, fsB);
        if (Nc == Na) {
            if (Na > Nb)
                std::swap(Na, Nb);
            if (!done.count(Na)) {
                COUT(<< "found (" << IntToStr(Na) << ", " << IntToStr(Nb) << ")" << std::endl);
                done.insert(Na);
                dm = 0; ds = 0; dl = 0; dots_len = 0; dl_before = 0; report_time = Time();
                return; // For GodBolt only
            }
        }
    }
}

void Test() {
    {
        size_t constexpr bits_factor = 80;
        using T = MinType<bits_factor + 5>;
        //T const N = T(RandPrime<T>((bits_factor + 1) / 2) * RandPrime<T>(bits_factor / 2));
        T const N = RandRange<T>(T(1) << (bits_factor - 1), T(1) << bits_factor);
        auto const fs = FactorECM<T>(N, true);
    }
    {
        size_t constexpr bits_amicable_min = 7, bits_amicable_max = 23;
        using T = MinType<bits_amicable_max + 5>;
        FindAmicable<T>(bits_amicable_min, bits_amicable_max);
    }
}

int main() {
    try {
        Test();
        return 0;
    } catch (std::exception const & ex) {
        COUT(<< "Exception: " << ex.what() << std::endl);
        return -1;
    }
}