#1119673 (C++23) No.3332 Consecutive Power Sum (Small)

• 問題: No.3332 Consecutive Power Sum (Small)
• ユーザー: にょぐた
• 提出日時: 2025-09-07 00:29:06
• 言語: C++23 (gcc 13.3.0 + boost 1.87.0)
• 結果: TLE
• 実行時間: -
• コード長: 15,857 bytes
• コンパイル時間: 4,004 ms
• コンパイル使用メモリ: 314,816 KB
• 実行使用メモリ: 19,016 KB
• 最終ジャッジ日時: 2025-11-02 21:09:05
• 合計ジャッジ時間: 10,596 ms
• ジャッジサーバーID （参考情報）: judge2 / judge4

ファイルパターン	結果
sample	TLE * 1 -- * 1
other	-- * 52

コンパイルメッセージ

main.cpp: In function ‘int main()’:
main.cpp:540:33: warning: narrowing conversion of ‘E’ from ‘int’ to ‘__int128 unsigned’ [-Wnarrowing]
  540 |                 ans.push_back({ E, it - cum.begin() + 1, k });
      |                                 ^
main.cpp:540:53: warning: narrowing conversion of ‘(__gnu_cxx::operator-<__int128 unsigned*, std::vector<__int128 unsigned> >(it, cum.std::vector<__int128 unsigned>::begin()) + 1)’ from ‘__gnu_cxx::__normal_iterator<__int128 unsigned*, std::vector<__int128 unsigned> >::difference_type’ {aka ‘long int’} to ‘__int128 unsigned’ [-Wnarrowing]
  540 |                 ans.push_back({ E, it - cum.begin() + 1, k });
      |                                    ~~~~~~~~~~~~~~~~~^~~

ソースコード

#include "bits/stdc++.h"
#include<iostream>
using namespace std;


using ll = long long;
using vll = vector<long long>;
using pll = pair<ll, ll>;
using vpll = vector<pll>;

#define reps(i, a, n) for (ll i = (a); i < (ll)(n); ++i)
#define rep(i, n) reps(i, 0, n)
#define rrep(i, n) reps(i, 1, n + 1)
#define repd(i,n) for(ll i=n-1;i>=0;i--)
#define rrepd(i,n) for(ll i=n;i>=1;i--)
#define repsd(i, a, n) for(ll i=n;i>=a;i--)
#define fore(i,a) for(auto &i:a)

using i16 = std::int16_t;
using u16 = std::uint16_t;
using i32 = std::int32_t;
using u32 = std::uint32_t;
using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = __uint128_t;

using LL = __uint128_t;
istream& operator>>(istream& is, LL& v)
{
    string s;
    is >> s;
    v = 0;
    for (int i = 0; i < (int)s.size(); i++) {
        if (isdigit(s[i])) { v = v * 10 + s[i] - '0'; }
    }
    if (s[0] == '-') { v *= -1; }
    return is;
}
ostream& operator<<(ostream& os, const LL& v)
{
    if (v == 0) { return (os << "0"); }
    LL num = v;
    if (v < 0) {
        os << '-';
        num = -num;
    }
    string s;
    for (; num > 0; num /= 10) { s.push_back((char)(num % 10) + '0'); }
    reverse(s.begin(), s.end());
    return (os << s);
}

LL cbrt128(LL n) {
    LL ng = 0, ok = 1;
    while (ok * ok * ok <= n) ok *= 2;
    while (ok - ng > 1) {
        LL mid = (ok + ng) / 2;
        LL tmp = mid * mid * mid;
        if (tmp >= n) ok = mid;
        else ng = mid;
    }
    return ok;
}

LL sqrt128(LL n) {
    LL ng = 0, ok = 1;
    while (ok * ok <= n) ok *= 2;
    while (ok - ng > 1) {
        LL mid = (ok + ng) / 2;
        LL tmp = mid * mid;
        if (tmp >= n) ok = mid;
        else ng = mid;
    }
    return ok;
}

namespace io
{
    constexpr u128 parse_u128(const std::string &s)
    {
        u128 n = 0;
        for (char c : s) {
            if ('0' <= c && c <= '9') {
                n = n * 10 + (c - '0');
            }
        }
        return n;
    }
    std::istream &operator>>(std::istream &is, u128 &n)
    {
        std::istream::sentry sen(is);
        if (sen) {
            std::string s;
            is >> s;
            n = parse_u128(s);
        }
        return is;
    }
    std::ostream &operator<<(std::ostream &os, u128 n)
    {
        std::ostream::sentry sen(os);
        if (sen) {
            if (n == 0) os << '0';
            std::string ns;
            while (n) {
                int r = n % 10;
                n /= 10;
                ns += char('0' + r);
            }
            std::ranges::reverse(ns);
            os << ns;
        }
        return os;
    }
}  // namespace io

inline namespace literals
{
    constexpr u128 operator""_u128(const char *s) { return io::parse_u128(s); }
}  // namespace literals

namespace factorization
{
    template <int MOD_BITS>
    struct integer_pack;
    template <>
    struct integer_pack<32> {
        using I0 = i16;
        using U0 = u16;
        using I1 = i32;
        using U1 = u32;
        using I2 = i64;
        using U2 = u64;
        static constexpr int BITS_0 = std::numeric_limits<U0>::digits;
        static constexpr int BITS_1 = std::numeric_limits<U1>::digits;
        static constexpr int BITS_2 = std::numeric_limits<U2>::digits;
    };
    template <>
    struct integer_pack<64> {
        using I0 = i64;
        using U0 = u64;
        using I1 = i64;
        using U1 = u64;
        using I2 = i128;
        using U2 = u128;
        static constexpr int BITS_0 = std::numeric_limits<U0>::digits;
        static constexpr int BITS_1 = std::numeric_limits<U1>::digits;
        static constexpr int BITS_2 = std::numeric_limits<U2>::digits;
    };
    template <>
    struct integer_pack<128> {
        using I0 = i64;
        using U0 = u64;
        using I1 = i128;
        using U1 = u128;
        static constexpr int BITS_0 = std::numeric_limits<U0>::digits;
        static constexpr int BITS_1 = std::numeric_limits<U1>::digits;
    };

    template <int MOD_BITS>
    struct dynamic_mod {
        using pack = integer_pack<MOD_BITS>;
        using U0 = pack::U0;
        using U1 = pack::U1;
        using I1 = pack::I1;
        static constexpr int BITS_0 = pack::BITS_0;
        static constexpr int BITS_1 = pack::BITS_1;

        // N
        U1 mod;
        // R^-1 % N
        U1 R_1;
        // R^1 % N
        U1 R1;
        // R^2 % N
        U1 R2;
        // -(N^-1) % R
        U1 N_;

        // `(x1 * x2) >> BITS_1`
        static U1 multiply_high(U1 x, U1 y)
        {
            U0 hx = x >> BITS_0, lx = x;
            U0 hy = y >> BITS_0, ly = y;
            U1 ans = U1(hx) * hy;
            ans += (U1(hx) * ly) >> BITS_0;
            ans += (U1(lx) * hy) >> BITS_0;
            U1 m = U1(hx * ly) + U1(lx * hy) + ((U1(lx) * ly) >> BITS_0);
            ans += m >> BITS_0;
            return ans;
        }

       private:
        void set_N_()
        {
            U1 n_inv = mod;
            for (int bit = 2; bit < BITS_1; bit <<= 1) {
                n_inv *= 2 - n_inv * mod;
            }
            N_ = -n_inv;
        }
        void set_R1() { R1 = -mod % mod; }
        void set_R2()
        {
            R2 = R1;
            for (int _ = 0; _ < BITS_1; ++_) {
                R2 <<= 1;
                if (R2 >= mod) R2 -= mod;
            }
        }
        void set_R_1() { R_1 = 1 + multiply_high(mod, N_); }

       public:
        dynamic_mod(const U1 &modulo) : mod(modulo)
        {
            assert(~mod >> (BITS_1 - 1));
            assert(mod & 1);
            set_N_();
            set_R1();
            set_R2();
            set_R_1();
        }

        U1 safe_mod(I1 x) const
        {
            x %= I1(mod);
            if (x < 0) x += I1(mod);
            return x;
        }
        // MR(x)
        U1 reduce(const U1 &x) const { return multiply_reduce(x, 1); }
        // MR(x * y)
        U1 multiply_reduce(const U1 &x, const U1 &y) const
        {
            U1 t_ = x * y * N_;
            U1 t = multiply_high(x, y) + multiply_high(t_, mod) + (x * y != 0);
            return t < mod ? t : t - mod;
        }

        U1 from(const U1 &x) const { return multiply_reduce(x % mod, R2); }
        U1 from(U1 &&x) const { return from(x); }
        U1 from_i(const I1 &x) const { return multiply_reduce(safe_mod(x), R2); }
        U1 from_i(I1 &&x) const { return from_i(x); }

        U1 val(const U1 &rx) const { return reduce(rx); }
        U1 pow(U1 rx, U1 e) const
        {
            U1 ans = R1, b = rx;
            while (e) {
                if (e & 1) ans = multiply_reduce(ans, b);
                b = multiply_reduce(b, b);
                e >>= 1;
            }
            return ans;
        }

        U1 add(const U1 &x, const U1 &y) const
        {
            U1 z = x + y;
            if (z >= mod) z -= mod;
            return z;
        }
        U1 sub(const U1 &x, const U1 &y) const
        {
            U1 z;
            if (__builtin_sub_overflow(x, y, &z)) z += mod;
            return z;
        }
        U1 mul(const U1 &x, const U1 &y) const { return multiply_reduce(x, y); }
        U1 neg(const U1 &x) const { return sub(0, x); }
    };

    bool is_prime(u128 n)
    {
        using U1 = u128;

        if (n < 2) return false;
        for (U1 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
            if (n == p) return true;
            if (n % p == 0) return false;
        }
        if (n < 41 * 41) return true;

        dynamic_mod<128> mont(n);
        const U1 one = mont.R1, neg_one = mont.neg(one);

        auto test_miller_rabin = [&](const std::vector<U1> &bases) {
            // int e = std::countr_zero(n - 1);
            // U1 o = n >> e;
            int e = 0;
            U1 o = n - 1;
            while (~o & 1) {
                o >>= 1;
                e++;
            }
            for (U1 b : bases) {
                U1 x = mont.pow(mont.from(b), o);
                if (x == one) continue;
                for (int ei = 0; ei < e; ++ei) {
                    U1 y = mont.mul(x, x);
                    if (y == one) {
                        if (x == neg_one) break;
                        return false;
                    }
                    x = y;
                    if (ei == e - 1) return false;
                }
            }
            return true;
        };
        if (n < 2047) return test_miller_rabin({2});
        if (n < 9080191) return test_miller_rabin({31, 73});
        if (n < 4759123141) return test_miller_rabin({2, 7, 61});
        if (n < 1122004669633) return test_miller_rabin({2, 13, 23, 1662803});
        if (n < 3770579582154547) return test_miller_rabin({2, 880937, 2570940, 610386380, 4130785767});
        if (n < 18446744073709551616_u128) return test_miller_rabin({2, 325, 9375, 28178, 450775, 9780504, 1795265022});
        if (n < 318665857834031151167461_u128) return test_miller_rabin({2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37});
        if (n < 3317044064679887385961981_u128) return test_miller_rabin({2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41});
        assert(false);
    }

    u128 gcd(u128 a, u128 b)
    {
        while (b) {
            u128 r = a % b;
            std::tie(a, b) = {b, r};
        }
        return a;
    }

    template <int ROUND = 1 << 12>
    std::vector<u128> factorize(u128 n)
    {
        using std::array;
        using std::vector;
        if (n == 1) return {};
        if (is_prime(n)) return {n};

        vector<u128> ans;
        for (u128 p : {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37}) {
            while (n % p == 0) {
                ans.push_back(p);
                n /= p;
            }
        }
        if (n == 1) return ans;
        // here n is odd, n >= 3 (or n > 37)
        auto abs_diff = [](u128 x, u128 y) { return x > y ? x - y : y - x; };

        auto main = [&](auto self, u128 nn) -> vector<u128> {
            if (nn == 1) return {};
            if (is_prime(nn)) return {nn};
            using namespace io;
            dynamic_mod<128> mont(nn);

            auto find_factor = [&]() {
                u128 rc = 0;
                auto f = [&](u128 rx) { return mont.add(mont.mul(rx, rx), rc); };

                while (true) {
                    ++rc;
                    u128 rx = rc, ry = rx;
                    u128 d = 1;

                    // round ごとにみる
                    array<u128, 2> checkpoint{rx, ry};
                    while (d == 1) {
                        u128 combined = 1;
                        for (u128 _ = 0; _ < ROUND; _++) {
                            rx = f(rx);
                            ry = f(f(ry));
                            combined = mont.mul(combined, abs_diff(rx, ry));
                        }
                        d = gcd(combined, nn);
                        if (d == 1) {
                            // この round では見つからず
                            checkpoint = {rx, ry};
                        }
                        if (d != 1 && d != nn) return d;
                    }

                    // 1つずつみる
                    // tie(rx, ry) = checkpoint;
                    rx = checkpoint[0];
                    ry = checkpoint[1];
                    d = 1;
                    while (d == 1) {
                        rx = f(rx);
                        ry = f(f(ry));
                        d = gcd(abs_diff(rx, ry), nn);
                        if (d != 1 && d != nn) return d;
                    }
                    // if (d == nn) continue;
                    // return d;
                }
            };
            u64 d = find_factor();
            vector a1 = self(self, d);
            vector a2 = self(self, nn / d);
            a1.insert(a1.end(), a2.begin(), a2.end());
            return a1;
        };
        vector a = main(main, n);
        ans.insert(ans.end(), a.begin(), a.end());
        return ans;
    }

    std::map<u128, int> to_freq(const std::vector<u128> &primes)
    {
        std::map<u128, int> freq;
        for (u128 p : primes) {
            freq[p]++;
        }
        return freq;
    }

    // [d for d divides n if pred(d)]
    // pred(d) == true  <==>  small enough
    template <typename T>
    requires requires(T pred, u128 l) {
        { pred(l) } -> std::same_as<bool>;
    }
    std::vector<u128> list_divisors(const std::map<u128, int> &pe, T pred)
    {
        if (!pred(1)) return {};
        std::vector<u128> ans{1};
        for (auto [p, e] : pe) {
            for (int i = 0, s = ans.size(); i < s; ++i) {
                u128 d = ans[i];
                for (int ei = 1; ei <= e; ++ei) {
                    d *= p;
                    if (!pred(d)) break;
                    ans.push_back(d);
                }
            }
        }
        return ans;
    }
}  // namespace factorization


u128 pow_int(u128 x, u128 p) {
    u128 res = 1;
    for (int i = 0; i < p; i++) res *= x;
    return res;
}

void fac2div(vector<pair<u128, u128>>& factors, vector<u128>& divisors, u128 val = 1, u128 idx = 0) {
    if (idx == factors.size()) divisors.push_back(val);
    else {
        fac2div(factors, divisors, val, idx + 1);
        rep(i, factors[idx].second) {
            val *= factors[idx].first;
            fac2div(factors, divisors, val, idx + 1);
        }
    }
}

// 約数列挙
vector<u128> divisor(u128 n) {
  auto mp = factorization::to_freq(factorization::factorize(2 * n));
  vector<pair<u128, u128>> factors(mp.size());
  vector<u128> divisors;
  auto it = mp.begin();
  while(it != mp.end()) factors.push_back({it->first, it->second}), it++;
  fac2div(factors, divisors);
  return divisors;
}
/*
vll divisor(ll n) {
    vll ret;
    for (ll i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            ret.push_back(i);
            if (i * i != n) ret.push_back(n / i);
        }
    }
    sort(begin(ret), end(ret));
    return ret;
}
*/



int main() {

    u128 n;
    cin >> n;
    vector<array<u128, 3>> ans;

    // d = r-l+1 によって場合分け

    vector<u128> divisors;

    // Eが2のとき
    divisors = divisor(6 * n);
    // 6*n = d * (2*d*d-6*d*r-3*d+6*r*r+6*r+1)
    // r = sqrt((12*n/d-d*d+1)/3)+d-1/2
    fore(d, divisors) {
        u128 r = (sqrt128((12 * n / d - d * d + 1) / 3) + d - 1) / 2;
        if (6 * n != d * (2 * d * d - 6 * d * r - 3 * d + 6 * r * r + 6 * r + 1)) continue;
        u128 l = r - d + 1;
        if(l<=0) continue;
        ans.push_back({ 2,l,r });
    }

    // Eが3のとき
    divisors = divisor(4 * n);
    // 4*n = d * (2*r-d+1)*(d*d-2*d*r-d+2*r*r+2*r)
    fore(d, divisors) {
        u128 ng = cbrt128(n) + 1, ok = d;
        if (ok > ng) continue;
        while (ng-ok > 1) {
            u128 mid = ok + ng >> 1;
            if (4 * n >= mid * mid * (mid + 1) * (mid + 1) - (mid - d) * (mid - d) * (mid - d + 1) * (mid - d + 1)) ok = mid;
            else ng = mid;
        }
        if (4 * n != ok * ok * (ok + 1) * (ok + 1) - (ok - d) * (ok - d) * (ok - d + 1) * (ok - d + 1)) continue;
        u128 r = ok;
        u128 l = r - d + 1;
        ans.push_back({ 3,l,r });
    }

    // Eが4以上のとき
    int E = 4;
    while (true) {
        if (pow_int(2, E) > n) break;
        vector<u128> cum;
        u128 sum = 0;

        u128 k = 0;
        while (true) {
            u128 tmp = pow_int(k, E);
            if (tmp > n) break;
            sum += tmp;
            cum.push_back(sum);
            auto it = lower_bound(cum.begin(), cum.end(), sum - n);
            if (it != cum.end() && *it == sum - n) {
                ans.push_back({ E, it - cum.begin() + 1, k });
            }

            k++;
        }


        E++;
    }

    

    sort(ans.begin(), ans.end());

    cout << ans.size() << "\n";
    for (auto e : ans) {
        cout << e[0] << " " << e[1] << " " << e[2] << "\n";
    }

}