#1179911 (C++23) No.3651 K-th Sum of Divisors

• 問題: No.3651 K-th Sum of Divisors
• ユーザー: kuromath
• 提出日時: 2026-08-28 21:59:02
• 言語: C++23 (gcc 15.2.0 + boost 1.90.0)
• 結果: AC
• 実行時間: 6 ms / 2,000 ms
• コード長: 14,703 bytes
• コンパイル時間: 5,103 ms
• コンパイル使用メモリ: 387,996 KB
• 実行使用メモリ: 6,272 KB
• 最終ジャッジ日時: 2026-08-28 22:00:37
• 合計ジャッジ時間: 6,965 ms
• ジャッジサーバーID （参考情報）: judge1_0 / judge2_0

ファイルパターン	結果
sample	AC * 3
other	AC * 55

ソースコード

#include <bits/stdc++.h>
#if __has_include(<atcoder/all>)
#include <atcoder/all>
using namespace atcoder;
#endif
//#pragma GCC target("avx2")
//#pragma GCC optimize("Ofast")
//#pragma GCC optimize("unroll-loops")
using namespace std;

inline void scan(int &a) { cin >> a; }
inline void scan(unsigned int &a) { cin >> a; }
inline void scan(long long &a) { cin >> a; }
inline void scan(unsigned long long &a) { cin >> a; }
inline void scan(float &a) { cin >> a; }
inline void scan(double &a) { cin >> a; }
inline void scan(long double &a) { cin >> a; }
inline void scan(char &a) { cin >> a; }
inline void scan(string &a) { cin >> a; }
inline void scan(vector<bool> &vec) { for (int i = 0; i < int(vec.size()); i++) { int a; scan(a); vec[i] = a; } }

template<typename T> inline void scan(vector<T> &vec);
template<typename T, size_t size> inline void scan(array<T, size> &vec);
template<typename T> inline void scan(deque<T> &deq);
template<typename T, typename U> inline void scan(pair<T, U> &p);

template<typename T> inline void scan(vector<T> &vec) { for (auto &i : vec) scan(i); }
template<typename T, size_t size> inline void scan(array<T, size> &vec) { for (auto &i : vec) scan(i); }
template<typename T> inline void scan(deque<T> &deq) { for (auto &i : deq) scan(i); }
template<typename T, typename U> inline void scan(pair<T, U> &p) { scan(p.first); scan(p.second); }
template<typename T> inline void scan(T &a) { cin >> a; }
inline void read() { }
template<typename T, typename... U> inline void read(T &a, U &...b) { scan(a); read(b...); }

inline void write(const int &a) { cout << a; }
inline void write(const unsigned int &a) { cout << a; }
inline void write(const long long &a) { cout << a; }
inline void write(const unsigned long long &a) { cout << a; }
inline void write(const float &a) { cout << a; }
inline void write(const double &a) { cout << a; }
inline void write(const long double &a) { cout << a; }
inline void write(const char &a) { cout << a; }
inline void write(const string &a) { cout << a; }

template<typename T> inline void write(const vector<T> &vec);
template<typename T, size_t size> inline void write(const array<T, size> &vec);
template<typename T> inline void write(const deque<T> &deq);
template<typename T, typename U> inline void write(const pair<T, U> &p);

template<typename T> inline void write(const vector<T> &vec) {
    for (int i = 0; i < int(vec.size()); i++) {
        if (i) write(' ');
        write(vec[i]);
    }
}

template<typename T, size_t size> inline void write(const array<T, size> &vec) {
    for (int i = 0; i < int(vec.size()); i++) {
        if (i) write(' ');
        write(vec[i]);
    }
}

template<typename T> inline void write(const deque<T> &deq) {
    for (int i = 0; i < int(deq.size()); i++) {
        if (i) write(' ');
        write(deq[i]);
    }
}

template<typename T, typename U> inline void write(const pair<T, U> &p) {
    write(p.first);
    write(' ');
    write(p.second);
}

template<typename T> inline void write(const T &a) { cout << a; }

inline void print() { cout << '\n'; }
template<typename T> inline void print(const T &a) { write(a); cout << '\n'; }
template<typename T, typename... U> inline void print(const T &a, const U &...b) {
    write(a);
    ((cout << ' ', write(b)), ...);
    write('\n');
}
template<typename T> inline void debug_impl(const char *name, const T &value) {
    cerr << name << " = ";
    write(value);
    cerr << '\n';
}

template<typename T, typename ... U> inline void debug_impl(const char *names, const T &value, const U &...values) {
    const char *comma = strchr(names, ',');
    cerr.write(names, comma - names);
    cerr << " = ";
    write(value);
    cerr << ", ";
    while (*++comma == ' ');
    debug_impl(comma, values...);
}

using ll = long long;
using ld = long double;
using ull = unsigned long long;
using pii = pair<int, int>;
using pll = pair<ll, ll>;
template<typename T> using pq = priority_queue<T>;
template<typename T> using pq_g = priority_queue<T, vector<T>, greater<T>>;

#define overload3(a, b, c, d, ...) d
#define overload4(a, b, c, d, e, ...) e
#define rep1(n) for (ll i = 0; i < (n); i++)
#define rep2(i, n) for (ll i = 0; i < (n); i++)
#define rep3(i, a, b) for (ll i = (a); i < (b); i++)
#define rep4(i, a, b, c) for (ll i = (a); i < (b); i += (c))
#define rep(...) overload4(__VA_ARGS__, rep4, rep3, rep2, rep1)(__VA_ARGS__)
#define rrep1(n) for (ll i = (n) - 1; i >= 0; i--)
#define rrep2(i, n) for (ll i = (n) - 1; i >= 0; i--)
#define rrep3(i, a, b) for (ll i = (b) - 1; i >= (a); i--)
#define rrep4(i, a, b, c) for (ll i = (a) + ((b) - (a) - 1) / (c) * (c); i >= (a); i -= (c))
#define rrep(...) overload4(__VA_ARGS__, rrep4, rrep3, rrep2, rrep1)(__VA_ARGS__)
#define all1(v) (v).begin(), (v).end()
#define all2(v, n) (v).begin(), (v).begin() + (n)
#define all3(v, a, b) (v).begin() + (a), (v).begin() + (b)
#define all(...) overload3(__VA_ARGS__, all3, all2, all1)(__VA_ARGS__)
#define rall(v) (v).rbegin(), (v).rend()
#define UNIQUE(v) (v).erase(unique((v).begin(), (v).end()), (v).end())
#define NP next_permutation
#define endl '\n'
#define INT(...) int __VA_ARGS__; read(__VA_ARGS__);
#define LL(...) long long __VA_ARGS__; read(__VA_ARGS__);
#define ULL(...) unsigned long long __VA_ARGS__; read(__VA_ARGS__);
#define STR(...) string __VA_ARGS__; read(__VA_ARGS__);
#define CHAR(...) char __VA_ARGS__; read(__VA_ARGS__);
#define LD(...) long double __VA_ARGS__; read(__VA_ARGS__);
#define VEC(type, name, size) vector<type> name(size); read(name);
#define VV(type, name, h, w) vector<vector<type>> name(h, vector<type>(w)); read(name);
#ifdef LOCAL
#define debug(...) debug_impl(#__VA_ARGS__, __VA_ARGS__)
#else
#define debug(...) ((void)0)
#endif

inline void Yes(bool a = true) { cout << (a ? "Yes" : "No") << endl; }
inline void No(bool a = true) { cout << (a ? "No" : "Yes") << endl; }
inline void YES(bool a = true) { cout << (a ? "YES" : "NO") << endl; }
inline void NO(bool a = true) { cout << (a ? "NO" : "YES") << endl; }
inline void Takahashi(bool a = true) { cout << (a ? "Takahashi" : "Aoki") << endl; }
inline void Alice(bool a = true) { cout << (a ? "Alice" : "Bob") << endl; }
inline void Possible(bool a = true) { cout << (a ? "Possible" : "Impossible") << endl; }

constexpr int INF = 1001001001;
constexpr ll LINF = 4004004004004004004LL;
constexpr ld pi = 3.14159265358979323846;
constexpr ld eps = 1e-9;
constexpr int dx[8] = { 1, 0, -1, 0, 1, -1, -1, 1 };
constexpr int dy[8] = { 0, 1, 0, -1, 1, 1, -1, -1 };

ll modPow(ll a, ll n, ll mod = LINF) { ll ret = 1; ll x = a % mod; while (n) { if (n & 1) ret = ret * x % mod; x = x * x % mod; n >>= 1; } return ret; }
int topbit(int x) { return (x == 0 ? -1 : 31 - __builtin_clz(x)); }
int topbit(ll x) { return (x == 0 ? -1 : 63 - __builtin_clzll(x)); }
int lowbit(int x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int lowbit(ll x) { return (x == 0 ? -1 : __builtin_ctz(x)); }
int popcount(int x) { return __builtin_popcount(x); }
ll popcount(ll x) { return __builtin_popcountll(x); }
ll isqrt(ll x) { ll n = sqrtl(x); while ((n + 1) * (n + 1) <= x) n++; while (n * n > x) n--; return n; }
template<typename T> T floor(T a, T b) { return a / b - (a % b && (a ^ b) < 0); }
template<typename T> T ceil(T a, T b) { return floor(a + b - 1, b); }
template<typename T> T bmod(T a, T b) { return a - b * floor(a, b); }
template<typename T> pair<T, T> divmod(T a, T b) { return { floor(a, b), a - b * floor(a, b) }; }
template<typename T> auto MIN(const vector<T> &a) { return *min_element(all(a)); }
template<typename T> auto MAX(const vector<T> &a) { return *max_element(all(a)); }
template<typename T> auto SUM(const vector<T> &a) { return accumulate(all(a), T(0)); }
template<typename T> bool chmin(T &a, const T &b) { if (a > b) { a = b; return true; } return false; }
template<typename T> bool chmax(T &a, const T &b) { if (a < b) { a = b; return true; } return false; }

void cincout() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    cout << fixed << setprecision(15);
}

struct Montgomery_modint {
    using modint = Montgomery_modint;
    using i64 = int64_t;
    using u64 = uint64_t;
    using u128 = __uint128_t;

    u64 x;
    static inline u64 mod;
    static inline u64 r;
    static inline u64 R2;

    Montgomery_modint() : x(u64(0)) {}
    Montgomery_modint(const int64_t &val) : x(reduce(u128(val % mod + mod) * R2)) {}

    static u64 get_r() {
        u64 ret = mod;
        for (int i = 0; i < 5; i++) ret *= 2 - mod * ret;
        return ret;
    }

    static void set_mod(u64 m) {
        assert(m & 1);
        assert(m < (u64(1) << 62));
        mod = m;
        r = get_r();
        R2 = -u128(mod) % mod;
    }

    static u64 reduce(const u128 &T) {
        return (T + u128(u64(T) * u64(-r)) * mod) >> 64;
    }

    modint &operator+=(const modint &p) {
        if (i64(x += p.x - 2 * mod) < 0) x += 2 * mod;
        return *this;
    }

    modint &operator-=(const modint &p) {
        if (i64(x -= p.x) < 0) x += 2 * mod;
        return *this;
    }

    modint &operator*=(const modint &p) {
        x = reduce(u128(x) * p.x);
        return *this;
    }

    modint &operator/=(const modint &p) {
        *this *= p.inv();
        return *this;
    }

    modint operator-() const { return modint() - *this; }
    modint operator+(const modint &p) const { return modint(*this) += p; }
    modint operator-(const modint &p) const { return modint(*this) -= p; }
    modint operator*(const modint &p) const { return modint(*this) *= p; }
    modint operator/(const modint &p) const { return modint(*this) /= p; }
    bool operator==(const modint &p) const { return (x >= mod ? x - mod : x) == (p.x >= mod ? p.x - mod : p.x); }
    bool operator!=(const modint &p) const { return (x >= mod ? x - mod : x) != (p.x >= mod ? p.x - mod : p.x); }

    modint pow(u64 n) const {
        modint ret(1), mul(*this);
        while (n > 0) {
            if (n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }

    modint inv() const {
        assert(x != 0);
        return pow(mod - 2);
    }

    u64 val() const {
        u64 ret = reduce(u128(x));
        return ret >= mod ? ret - mod : ret;
    }

    static u64 get_mod() {
        return mod;
    }

    friend istream &operator>>(istream &is, modint &a) {
        int64_t x;
        is >> x;
        a = modint(x);
        return is;
    }

    friend ostream &operator<<(ostream &os, const modint &a) {
        return os << a.val();
    }
};

namespace fast_factorize {
    using mint = Montgomery_modint;
    using u64 = uint64_t;
    mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

    bool MillerRabin(u64 n, const vector<u64> &A) {
        mint::set_mod(n);
        u64 s = 0, d = n - 1;
        while (~d & 1) {
            s++;
            d >>= 1;
        }

        for (u64 a : A) {
            if (n <= a) break;
            mint x = mint(a).pow(d);
            if (x != 1) {
                int t = 0;
                while (t < s) {
                    if (x == n - 1) break;
                    x *= x;
                    t++;
                }
                if (t == s) return false;
            }
        }
        return true;
    }

    bool primetest(u64 n) {
        if (n <= 1) return false;
        if (n == 2) return true;
        if (~n & 1) return false;
        if (n < 4759123141LL) return MillerRabin(n, {2, 7, 61});
        return MillerRabin(n, {2, 325, 9375, 28178, 450775, 9780504, 1795265022});
    }

    u64 PollardRho(u64 n) {
        if (~n & 1) return 2;
        if (primetest(n)) return n;
        if (mint::get_mod() != n) mint::set_mod(n);
        mint c, unit = 1;
        auto f = [&](mint x) { return x * x + c; };
        auto random = [&](u64 l, u64 r) {
            return uniform_int_distribution<u64>(l, r)(rng);
        };

        while (1) {
            mint x, y, z, q = unit;
            constexpr int m = 128;
            u64 g = 1;
            c = random(1, n - 1);
            y = random(1, n - 1);
            for (int r = 1; g == 1; r <<= 1) {
                x = y;
                for (int i = 0; i < r; i++) y = f(y);
                for (int k = 0; k < r && g == 1; k += m) {
                    z = y;
                    for (int i = 0; i < min(m, r - k); i++) {
                        y = f(y);
                        q *= x - y;
                    }
                    g = gcd(q.val(), n);
                }
            }
            if (g != n) return g;
            do {
                z = f(z);
                g = gcd((x - z).val(), n);
            } while (g == 1);
        }
    }

    vector<u64> inner_factorize(u64 n) {
        if (n <= 1) return {};
        u64 p = PollardRho(n);
        if (p == n) return {p};
        auto l = inner_factorize(p);
        auto r = inner_factorize(n / p);
        l.append_range(r);
        return l;
    }

    vector<u64> factorize(u64 n) {
        auto F = inner_factorize(n);
        sort(F.begin(), F.end());
        return F;
    }

    vector<u64> divisors(u64 n) {
        if (n == 0) return {};
        auto pf = factorize(n);
        vector<pair<u64, int>> F;
        F.reserve(int(pf.size()));
        for (u64 &p : pf) {
            if (F.empty() || F.back().first != p) {
                F.emplace_back(p, 1);
            } else {
                F.back().second++;
            }
        }

        vector<u64> d = {1};
        for (auto &[p, e] : F) {
            int sz = d.size();
            u64 pw = 1;
            for (int i = 0; i < e; i++) {
                pw *= p;
                for (int j = 0; j < sz; j++) {
                    d.push_back(d[j] * pw);
                }
            }
        }
        sort(d.begin(), d.end());
        return d;
    }
};

using fast_factorize::primetest;
using fast_factorize::factorize;
using fast_factorize::divisors;

constexpr int P = 100003;

int f(ll x) {
    int s = 0;
    for (auto d : divisors(x)) {
        s += d;
    }
    return s % P;
}

void solve() {
    LL(N, K);
    K--;
    if (K == 0) {
        print(N);
        return;
    }

    modint::set_mod(100003);
    vector<int> dp(P), ndp(P);
    for (int i = 1; i < P; i++) {
        for (int j = i; j < P; j += i) {
            dp[j] += i;
            if (dp[j] >= P) dp[j] -= P;
        }
    }

    int ans = f(N);
    K--;
    for (int i = 0; i <= 59; i++) {
        ndp.assign(P, 0);
        if (K >> i & 1) {
            ans = dp[ans];
        }
        for (int j = 1; j < P; j++) ndp[j] = dp[dp[j]];
        swap(dp, ndp);
    }
    print(ans);
}

int main() {
    cincout();
    int t = 1;
    // cin >> t;
    while (t--) solve();
    return 0;
}