結果

問題 No.3333 Consecutive Power Sum (Large)
コンテスト
ユーザー にょぐた
提出日時 2025-09-12 01:33:37
言語 C++23
(gcc 13.3.0 + boost 1.87.0)
結果
WA  
実行時間 -
コード長 8,016 bytes
コンパイル時間 3,935 ms
コンパイル使用メモリ 308,872 KB
実行使用メモリ 110,232 KB
最終ジャッジ日時 2025-11-02 21:09:58
合計ジャッジ時間 29,275 ms
ジャッジサーバーID
(参考情報)
judge2 / judge4
このコードへのチャレンジ
(要ログイン)
ファイルパターン 結果
sample AC * 3
other AC * 62 WA * 1
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.cpp:127:17: warning: integer constant is too large for its type
  127 |     if (n < (LL)318665857834031151167461) return miller_test(n, { 2,3,5,7,11,13,17,19,23,29,31,37 });
      |                 ^~~~~~~~~~~~~~~~~~~~~~~~
main.cpp:128:17: warning: integer constant is too large for its type
  128 |     if (n < (LL)3317044064679887385961981) return miller_test(n, { 2,3,5,7,11,13,17,19,23,29,31,37,41 });
      |                 ^~~~~~~~~~~~~~~~~~~~~~~~~
main.cpp: In function ‘int main()’:
main.cpp:298:33: warning: narrowing conversion of ‘E’ from ‘int’ to ‘__int128 unsigned’ [-Wnarrowing]
  298 |                 ans.push_back({ E, it - cum.begin() + 1, k });
      |                                 ^
main.cpp:298: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]
  298 |                 ans.push_back({ E, it - cum.begin() + 1, k });
      |                                    ~~~~~~~~~~~~~~~~~^~~
main.cpp: In function ‘bool is_prime(LL)’:
main.cpp:130:1: warning: control reaches end of non-void function [-Wreturn-type]
  130 | }
      | ^
main.cpp: In function ‘LL find_prime_factor(LL)’:
main.cpp:159:1: warning: control reaches end of non-void function [-Wreturn-type]
  159 | }
      | ^

ソースコード

diff #

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

using LL = __uint128_t;
using vLL = vector<LL>;

#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)


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"); }
    __int128_t 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);
}

struct Montgomery {
    LL n, r2, nInv;

    Montgomery(LL n) {
        LL r2 = ((LL)1 << 124) % n;
        for (int i = 0; i < 132; i++) {
            r2 <<= 1;
            if (r2 >= n) {
                r2 -= n;
            }
        }
        LL nInv = n;
        for (int i = 0; i < 10; ++i) {
            nInv *= 2 - n * nInv;
        }
        this->n = n;
        this->r2 = r2;
        this->nInv = nInv;
    }

    LL high(LL x, LL y) {
        LL xh = x >> 64;
        LL yh = y >> 64;
        LL xl = x & 0xFFFFFFFFFFFFFFFFL;
        LL yl = y & 0xFFFFFFFFFFFFFFFFL;
        return xh * yh + (xh * yl >> 64) + (xl * yh >> 64) + ((((xh * yl & 0xFFFFFFFFFFFFFFFFL) + (xl * yh & 0xFFFFFFFFFFFFFFFFL)) + (xl * yl >> 64)) >> 64);
    }

    LL mulMr(LL x, LL y) {
        return high(x, y) + high(-nInv * x * y, n) + (x * y == 0 ? 0 : 1);
    }

    LL mod(LL x) {
        return x < n ? x : x-n;
    }

    LL mul(LL x, LL y) {
        return mod(mulMr(mulMr(x, r2), y));
    }

    LL pow(LL x, LL y) {
        LL z = mulMr(x, r2);
        LL r = 1;
        while (y > 0) {
            if ((y & 1) == 1) {
                r = mulMr(r, z);
            }
            z = mulMr(z, z);
            y >>= 1;
        }
        return mod(r);
    }
};

bool miller_test(LL n, vLL v) {
    LL s = 0, d = n - 1, y;
    Montgomery mon(n);
    while (d % 2 == 0) s++, d /= 2;
    for(auto a: v){
      LL x = mon.pow(a, d);
      for(int i=0; i<s; i++){
        y = mon.pow(x, 2);
        if(y==1 && x!= 1 && x != n-1) return false;
        x = y;
      }
      if(y!=1) return false;
    }
    return true;
}

bool is_prime(LL n) {
    if(n==1) return false;
    if(n<=3) return true;
    if(n%2==0) return false;
    if (n < (LL)2047) return miller_test(n, { 2 });
    if (n < (LL)1373653) return miller_test(n, { 2,3 });
    if (n < (LL)9080191) return miller_test(n, { 31,73 });
    if (n < (LL)25326001) return miller_test(n, { 2,3,5 });
    if (n < (LL)3215031751) return miller_test(n, { 2,3,5,7 });
    if (n < (LL)4759123141) return miller_test(n, { 2,7,61 });
    if (n < (LL)1122004669633) return miller_test(n, { 2,13,23,1662803 });
    if (n < (LL)2152302898747) return miller_test(n, { 2,3,5,7,11 });
    if (n < (LL)3474749660383) return miller_test(n, { 2,3,5,7,11,13 });
    if (n < (LL)341550071728321) return miller_test(n, { 2,3,5,7,11,13,17 });
    if (n < (LL)3825123056546413051) return miller_test(n, { 2,3,5,7,11,13,17,19,23 });
    if (n < (LL)318665857834031151167461) return miller_test(n, { 2,3,5,7,11,13,17,19,23,29,31,37 });
    if (n < (LL)3317044064679887385961981) return miller_test(n, { 2,3,5,7,11,13,17,19,23,29,31,37,41 });
    
}

LL f(LL a, LL c, Montgomery& mon){
  return mon.mod(mon.mul(a,a)+c);
}

LL gcd(LL a, LL b){
  while(a){
    b %=a;
    swap(a,b);
  }
  return b;
}


LL find_prime_factor(LL n){
  Montgomery mon(n);
  for(LL c=0; c<n; c++){
    LL x=0, y=0, g=1;
    while(g==1){
      x = f(x, c, mon);
      y = f(f(y, c, mon), c, mon);
      g = gcd((x>y) ? x-y : y-x, n);
    }
    if(g==n) continue;
    if(is_prime(g)) return g;
    else if(is_prime(n/g)) return n/g;
    else return find_prime_factor(g);
  }
}

map<LL, LL> factorize(LL n) {
    map<LL, LL> mp;
    while (n % 2 == 0) mp[2]++, n/=2;
    if(n==1) return mp;
    while (true) {
        if (is_prime(n)) {
            mp[n]++;
            return mp;
        }
        else {
            LL f = find_prime_factor(n);
            while(n%f==0) mp[f]++, n/=f;
            if(n==1) return mp;
        }
    }
}

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) {
    if (n == 0) return 0;

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

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

void fac2div(vector<pair<LL, LL>>& factors, vector<LL>& divisors, LL val = 1, LL 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<LL> divisor(LL n) {
  auto mp = factorize(2*n);
  vector<pair<LL, LL>> factors(mp.size());
  vector<LL> divisors;
  auto it = mp.begin();
  while(it != mp.end()) factors.push_back({it->first, it->second}), it++;
  fac2div(factors, divisors);
  return divisors;
}

int main()
{
    using namespace std;
    LL n;
    cin >> n;
  vector<array<LL, 3>> ans;

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

    vector<LL> divisors;
    
    // Eが1のとき
    divisors = divisor(2*n);
    // 2*n = d * (2*r-d+1)
    fore(d, divisors) {
        LL r = (2 * n / d + d - 1) / 2;
        if (2 * n != d * (2 * r - d + 1)) continue;
        if (r+1<=d) continue;
        LL l = r - d + 1;
        ans.push_back({ 1,l,r });
    }

    // 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) {
      if(n*12 / d + 1 <d * d) continue;
        LL r = (sqrt128((n*12 / 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;
        if(r+1<=d) continue;
        LL l = r - d + 1;
        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) {
        LL ng = cbrt128(n) + 1, ok = d;
        if (ok > ng) continue;
        while (ng-ok > 1) {
            LL 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;
        LL r = ok;
        LL 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<LL> cum;
        LL sum = 0;

        LL k = 0;
        while (true) {
            LL 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";
    }
}
0