結果

問題 No.2318 Phys Bone Maker
コンテスト
ユーザー hliuser1
提出日時 2023-06-09 04:30:14
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
TLE  
実行時間 -
コード長 5,564 bytes
コンパイル時間 3,705 ms
コンパイル使用メモリ 272,160 KB
実行使用メモリ 13,644 KB
最終ジャッジ日時 2024-12-31 11:26:13
合計ジャッジ時間 32,147 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge4
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 40 TLE * 5
権限があれば一括ダウンロードができます

ソースコード

diff # 

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

static constexpr ll Q = 998244353;

template <typename K, typename V>
using Map = unordered_map<K, V>;

template <typename Int>
vector<Int> factorize(Int n, Int i=1) {
    vector<Int> f;
    f.reserve( sqrt(n) );
    for (; i*i < n; ++i)
        if (n % i == 0)
            f.push_back(i);
    i -= (i - (n / i) == 1);
    for (; i >= 1; i--)
        if (n % i == 0)
            f.push_back(n/i);
    return f;
}

struct _big_prime_factorization {
    static uint64_t random_address() { char *p = new char; delete p; return uint64_t(p); }
    inline static auto rng = mt19937_64(chrono::steady_clock::now().time_since_epoch().count() * (random_address() | 1));

    struct barrett_reduction {
        uint64_t mod;
        uint64_t div;
        barrett_reduction(uint64_t m) : mod(m), div(-1LLU / m) {}
        uint64_t operator()(uint64_t a) const {
#ifdef __SIZEOF_INT128__
            uint64_t q = uint64_t(__uint128_t(div) * a >> 64);
            uint64_t r = a - q * mod;
            return uint64_t(r < mod ? r : r - mod);
#endif
            return uint64_t(a % mod);
        }
    };
    static bool miller_rabin(uint64_t n) {
        if (n < 2) return false;
        for (uint64_t p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29})
            if (n % p == 0)
                return n == p;
        auto get_miller_rabin_bases = [&]() -> vector<uint64_t> {
            if (n < 341531) return {9345883071009581737LLU};
            if (n < 1050535501) return {336781006125, 9639812373923155};
            return {4230279247111683200, 14694767155120705706LLU, 16641139526367750375LLU};
        };
        int r = __builtin_ctz(n - 1);
        uint64_t d = (n - 1) >> r;
        barrett_reduction barrett(n);
        auto mod_pow = [&](uint64_t a, uint64_t b) -> uint64_t {
            uint64_t result = 1;
            while (b > 0) {
                if (b & 1) result = barrett(uint64_t(result) * a);
                a = barrett(uint64_t(a) * a);
                b >>= 1;
            }
            return result;
        };
        for (uint64_t a : get_miller_rabin_bases()) {
            if (a % n == 0) continue;
            uint64_t x = mod_pow(uint64_t(a % n), d);
            if (x == 1 || x == n - 1) continue;
            for (int i = 0; i < r - 1 && x != n - 1; i++)
                x = barrett(uint64_t(x) * x);
            if (x != n - 1) return false;
        }
        return true;
    }
    uint64_t binary_gcd(uint64_t a, uint64_t b) {
        if (a == 0 || b == 0) return a + b;
        int common = __builtin_ctzll(a | b);
        b >>= __builtin_ctzll(b);
        do {
            a >>= __builtin_ctzll(a);
            if (a < b) swap(a, b);
            a -= b;
        } while (a != 0);
        return b << common;
    }
    uint64_t pollard_rho(uint64_t n) {
        for (uint64_t p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29})
            if (n % p == 0)
                return p;
        barrett_reduction barrett(n);
        uint64_t increment;
        auto g = [&](uint64_t x) -> uint64_t {
            return barrett(x * x + increment);
        };
        // Choose a jump size much larger than log(n) but much smaller than n^(1/4).
        int jump = (int64_t)(sqrt(log(n) * sqrt(sqrt(n))));
        while (true) {
            increment = uint64_t(rng() % n);
            uint64_t start = uint64_t(rng() % n);
            uint64_t x = start, y = start, p = 1;
            vector<uint64_t> products(jump + 1);
            do {  // https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm#Variants
                products[0] = 1;
                for (int i = 1; i <= jump; i++) {
                    x = g(x);
                    y = g(g(y));
                    products[i] = barrett(uint64_t(products[i - 1]) * (max(x, y) - min(x, y)));
                }
            } while ((p = binary_gcd(products[jump], n)) == 1);
            if (p == n) {
                assert(products[jump] == 0);
                int index = jump;
                while (index > 0 && products[index] == 0) --index;
                p = binary_gcd(products[index], n);
            }
            if (p != 1 && p != n) return p;
        }
    }
    template<typename T>
    Map<T, int> operator()(T n) {
        if (n == 1) return {};
        if (miller_rabin(n)) return {{n, 1}};
        T x = pollard_rho(n);
        Map<T, int> A = this->operator()(x), B = this->operator()(n / x);
        if (A.size() < B.size())
            swap(A, B);
        for (auto [p, cnt] : B)
            A[p] += cnt;
        return A;
    }
} prime_factors;

signed main() {
    ll N;
    cin >> N;
    unordered_map<ll, ll> dp;

    dp[1] = 1;
    auto facts = factorize(N);
    int X = facts.size();

    // dp[x] = sum(dp[y]*Z forall y such that x%y==0) where Z is # of z s.t. lcm(y,z) = x

    vector<Map<ll,int>> prime_facs(X);
    for (int i = 0; i < X; ++i)
        prime_facs[i] = prime_factors(facts[i]);

    auto get = [] (Map<ll, int>& mp, ll key) {
        return mp.contains(key) ? mp[key] : 0;
    };

    for (int i = 1; i < X; ++i) {
        ll x = facts[i];
        auto& pfx = prime_facs[i];
        for (int j = 0; j < i; ++j) {
            ll y = facts[j];
            if (x % y != 0) continue;

            ll ways = dp[y];  // # of z where lcm(y,z) = x
            auto& pfy = prime_facs[j];
            for (auto [p, cnt] : pfx)
                if (get(pfy, p) == cnt)
                    ways = ways * (cnt+1) % Q;

            dp[x] = (dp[x]+ways) % Q;
        }
    }
    cout << dp[N] << '\n';

}
0