結果
| 問題 | No.3333 Consecutive Power Sum (Large) |
| コンテスト | |
| ユーザー |
 Vi24E
|
| 提出日時 | 2026-08-04 02:25:58 |
| 言語 | C++23(gcc16) (gcc 16.1.0 + boost 1.90.0) |
| 結果 |
AC
|
| 実行時間 | 210 ms / 10,000 ms |
| + 433µs | |
| コード長 | 26,459 bytes |
| 記録 | |
| コンパイル時間 | 6,336 ms |
| コンパイル使用メモリ | 403,720 KB |
| 実行使用メモリ | 39,932 KB |
| 最終ジャッジ日時 | 2026-08-04 02:26:13 |
| 合計ジャッジ時間 | 11,580 ms |
|
ジャッジサーバーID （参考情報） |
judge2_0 / judge3_0 |
（要ログイン）
| ファイルパターン | 結果 |
|---|---|
| sample | AC * 3 |
| other | AC * 63 |
ソースコード
#include <bits/stdc++.h>
using namespace std;
using u64 = uint64_t;
using u128 = __uint128_t;
static string to_string_u128(u128 x) {
if (!x) return "0";
char s[64];
int n = 0;
while (x) {
s[n++] = char('0' + x % 10);
x /= 10;
}
reverse(s, s + n);
return string(s, s + n);
}
static u128 parse_u128(const string& s) {
u128 x = 0;
for (char c : s) x = x * 10 + (c - '0');
return x;
}
static u64 isqrt_u128(u128 n) {
u64 x = (u64)sqrt((long double)n);
while ((u128)x * x > n) --x;
while ((u128)(x + 1) * (x + 1) <= n) ++x;
return x;
}
/*
* Montgomery arithmetic modulo an odd integer below 2^127.
* Values are represented modulo R = 2^128.
*/
struct Mont128 {
u128 mod, ninv, r1, r2;
static inline u128 mul_hi(u128 x, u128 y) {
u64 x0 = (u64)x, x1 = (u64)(x >> 64);
u64 y0 = (u64)y, y1 = (u64)(y >> 64);
u128 z11 = (u128)x1 * y1;
u128 z10 = (u128)x1 * y0;
u128 z01 = (u128)x0 * y1;
u128 z00 = (u128)x0 * y0;
u128 mid = (z00 >> 64) + (u64)z10 + (u64)z01;
return z11 + (z10 >> 64) + (z01 >> 64) + (mid >> 64);
}
explicit Mont128(u128 n) : mod(n) {
u128 x = n;
for (int bits = 1; bits < 128; bits <<= 1) {
x *= 2 - x * n;
}
ninv = -x;
r1 = (-n) % n;
r2 = r1;
for (int i = 0; i < 128; ++i) {
r2 <<= 1;
if (r2 >= n) r2 -= n;
}
}
inline u128 mul(u128 x, u128 y) const {
u128 lo = x * y;
u128 m = lo * ninv;
u128 t = mul_hi(x, y) + mul_hi(m, mod) + (lo != 0);
return t >= mod ? t - mod : t;
}
inline u128 init(u128 x) const {
return mul(x % mod, r2);
}
inline u128 norm(u128 x) const {
return mul(x, 1);
}
inline u128 add(u128 x, u128 y) const {
u128 z = x + y;
if (z >= mod) z -= mod;
return z;
}
inline u128 pow_mont(u128 x, u128 e) const {
u128 a = r1;
while (e) {
if (e & 1) a = mul(a, x);
x = mul(x, x);
e >>= 1;
}
return a;
}
inline u128 pow(u128 x, u128 e) const {
return norm(pow_mont(init(x), e));
}
};
static u128 gcd128(u128 a, u128 b) {
while (b) {
u128 r = a % b;
a = b;
b = r;
}
return a;
}
static bool is_prime128(u128 n) {
if (n < 2) return false;
static constexpr uint32_t SMALL[] = {
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37
};
for (uint32_t p : SMALL) {
if (n == p) return true;
if (n % p == 0) return false;
}
Mont128 md(n);
u128 d = n - 1;
int s = 0;
while (!(d & 1)) {
d >>= 1;
++s;
}
auto witness = [&](u128 a) {
if (a % n == 0) return false;
u128 x = md.pow_mont(md.init(a), d);
if (x == md.r1 || x == md.mod - md.r1) return false;
for (int r = 1; r < s; ++r) {
x = md.mul(x, x);
if (x == md.mod - md.r1) return false;
}
return true;
};
// Deterministic for n < 3.317... * 10^24.
static constexpr uint32_t BASES[] = {
2, 3, 5, 7, 11, 13, 17,
19, 23, 29, 31, 37, 41
};
for (uint32_t a : BASES) {
if (witness(a)) return false;
}
return true;
}
static inline u64 mul_mod64(u64 a, u64 b, u64 mod) {
return (u128)a * b % mod;
}
static u64 pow_mod64(u64 a, u64 e, u64 mod) {
u64 r = 1;
while (e) {
if (e & 1) r = mul_mod64(r, a, mod);
a = mul_mod64(a, a, mod);
e >>= 1;
}
return r;
}
/* Brent–Pollard rho for 64-bit recursive cofactors. */
static u64 rho64(u64 n) {
if (n % 2 == 0) return 2;
auto rng = []() -> u64 {
static u64 x =
0x9e3779b97f4a7c15ULL ^
(u64)chrono::high_resolution_clock::now()
.time_since_epoch().count();
x ^= x << 7;
x ^= x >> 9;
return x;
};
while (true) {
u64 y = rng() % (n - 1) + 1;
u64 c = rng() % (n - 1) + 1;
constexpr u64 M = 128;
u64 g = 1, r = 1, q = 1;
u64 x = 0, ys = 0;
auto f = [&](u64 v) {
return (mul_mod64(v, v, n) + c) % n;
};
while (g == 1) {
x = y;
for (u64 i = 0; i < r; ++i) y = f(y);
for (u64 k = 0; k < r && g == 1; k += M) {
ys = y;
u64 lim = min(M, r - k);
for (u64 i = 0; i < lim; ++i) {
y = f(y);
u64 d = x > y ? x - y : y - x;
q = mul_mod64(q, d, n);
}
g = std::gcd(q, n);
}
r <<= 1;
}
if (g == n) {
do {
ys = f(ys);
u64 d = x > ys ? x - ys : ys - x;
g = std::gcd(d, n);
} while (g == 1);
}
if (g != n) return g;
}
}
/* Tonelli–Shanks modulo a small prime. */
static u64 tonelli(uint32_t n, uint32_t p) {
if (p == 2) return n & 1;
if (n == 0) return 0;
if (p % 4 == 3) return pow_mod64(n, (p + 1) / 4, p);
uint32_t q = p - 1;
uint32_t s = 0;
while ((q & 1) == 0) {
q >>= 1;
++s;
}
uint32_t z = 2;
while (pow_mod64(z, (p - 1) / 2, p) != p - 1) ++z;
u64 c = pow_mod64(z, q, p);
u64 x = pow_mod64(n, (q + 1) / 2, p);
u64 t = pow_mod64(n, q, p);
uint32_t m = s;
while (t != 1) {
uint32_t i = 1;
u64 tt = mul_mod64(t, t, p);
while (tt != 1) {
tt = mul_mod64(tt, tt, p);
++i;
}
u64 b = pow_mod64(c, 1ULL << (m - i - 1), p);
x = mul_mod64(x, b, p);
c = mul_mod64(b, b, p);
t = mul_mod64(t, c, p);
m = i;
}
return x;
}
struct FBPrime {
uint32_t p, r1, r2;
uint16_t lg;
};
struct Atom {
u64 x, q;
};
struct Relation {
vector<u64> parity;
array<Atom, 2> a;
uint8_t cnt;
};
struct Partial {
vector<u64> parity;
Atom a;
};
static vector<uint32_t> primes_upto(int bound) {
vector<bool> is_prime(bound + 1, true);
is_prime[0] = is_prime[1] = false;
vector<uint32_t> primes;
for (int i = 2; i <= bound; ++i) {
if (!is_prime[i]) continue;
primes.push_back(i);
if ((int64_t)i * i <= bound) {
for (int j = i * i; j <= bound; j += i) {
is_prime[j] = false;
}
}
}
return primes;
}
/*
* Quadratic sieve specialized for 65–80-bit integers.
*
* Q(x) = (ceil(sqrt(kN)) + x)^2 - kN.
* Within the scanned range Q(x) fits in uint64_t, so almost all sieve and
* trial-division work uses native 64-bit arithmetic.
*/
static u128 quadratic_sieve_factor(
u128 n,
int factor_bound = 7000,
int extra_relations = 48,
u64 max_base = 3500000
) {
constexpr int BLOCK = 1 << 15;
vector<uint32_t> primes = primes_upto(factor_bound);
static constexpr uint32_t MULTIPLIERS[] = {
1, 3, 5, 7, 11, 13, 15, 17, 19, 21,
23, 29, 31, 33, 35, 37, 39, 41, 43, 47
};
uint32_t best_k = 1;
long double best_score = -1e100L;
for (uint32_t k : MULTIPLIERS) {
u128 g = gcd128(n, k);
if (g > 1 && g < n) return g;
long double score = -0.5L * log((long double)k);
for (uint32_t p : primes) {
if (p > 100) break;
u64 a = (u64)((n % p) * (k % p) % p);
if (a == 0) {
score += log((long double)p) / (p - 1);
} else if (
p == 2 ||
pow_mod64(a, (p - 1) / 2, p) == 1
) {
score += 2 * log((long double)p) / (p - 1);
}
}
if (score > best_score) {
best_score = score;
best_k = k;
}
}
u128 kn = n * best_k;
u64 m = isqrt_u128(kn);
if ((u128)m * m < kn) ++m;
if ((u128)m * m == kn) {
u128 g = gcd128(m, n);
if (g > 1 && g < n) return g;
best_k = 1;
kn = n;
m = isqrt_u128(n);
if ((u128)m * m < n) ++m;
}
vector<FBPrime> factor_base;
factor_base.reserve(primes.size() / 2 + 8);
for (uint32_t p : primes) {
if (n % p == 0) return p;
u64 a = (u64)(kn % p);
if (p == 2) {
uint32_t y = a & 1;
uint32_t mm = m & 1;
uint32_t r = (y + 2 - mm) & 1;
factor_base.push_back({p, r, r, 64});
continue;
}
u64 legendre = pow_mod64(a, (p - 1) / 2, p);
if (legendre != 1 && a != 0) continue;
uint32_t y = (uint32_t)tonelli((uint32_t)a, p);
uint32_t mm = m % p;
uint32_t r1 = (y + p - mm) % p;
uint32_t r2 = ((y ? p - y : 0) + p - mm) % p;
factor_base.push_back({
p,
r1,
r2,
(uint16_t)llround(log2((long double)p) * 64)
});
}
const int columns = (int)factor_base.size();
const int parity_words = (columns + 63) >> 6;
vector<Relation> relations;
relations.reserve(columns + extra_relations + 8);
unordered_map<u64, Partial> partials;
partials.reserve((columns + extra_relations) * 8);
vector<uint16_t> sieve(BLOCK);
const u64 large_prime_bound =
(u64)factor_bound * factor_bound;
for (
u64 base = 0;
relations.size() <
(size_t)(columns + extra_relations);
base += BLOCK
) {
fill(sieve.begin(), sieve.end(), 0);
for (const FBPrime& f : factor_base) {
uint32_t p = f.p;
uint32_t base_mod = base % p;
uint32_t start1 = (f.r1 + p - base_mod) % p;
for (uint32_t i = start1; i < BLOCK; i += p) {
sieve[i] += f.lg;
}
if (f.r2 != f.r1) {
uint32_t start2 = (f.r2 + p - base_mod) % p;
for (uint32_t i = start2; i < BLOCK; i += p) {
sieve[i] += f.lg;
}
}
}
for (
uint32_t i = 0;
i < BLOCK &&
relations.size() <
(size_t)(columns + extra_relations);
++i
) {
u64 x = base + i;
u128 qq =
(u128)(m + x) * (m + x) - kn;
if (
qq == 0 ||
qq > numeric_limits<u64>::max()
) {
continue;
}
u64 q = (u64)qq;
int bits = 64 - __builtin_clzll(q);
// Allows repeated small prime factors not represented
// repeatedly in the logarithmic sieve.
if ((int)sieve[i] + 14 * 64 < bits * 64) {
continue;
}
u64 rem = q;
vector<u64> parity(parity_words);
for (int j = 0; j < columns; ++j) {
uint32_t p = factor_base[j].p;
int odd = 0;
while (rem % p == 0) {
rem /= p;
odd ^= 1;
}
if (odd) {
parity[j >> 6] ^=
1ULL << (j & 63);
}
if (rem == 1) break;
}
if (rem == 1) {
Relation relation;
relation.parity = move(parity);
relation.a[0] = {x, q};
relation.cnt = 1;
relations.push_back(move(relation));
} else if (rem <= large_prime_bound) {
/*
* Every prime <= factor_bound that can divide Q(x)
* belongs to the factor base. Therefore a remaining
* value <= B^2 must be prime.
*/
auto it = partials.find(rem);
if (it == partials.end()) {
partials.emplace(
rem,
Partial{move(parity), {x, q}}
);
} else {
Relation relation;
relation.parity = move(parity);
for (int w = 0; w < parity_words; ++w) {
relation.parity[w] ^=
it->second.parity[w];
}
relation.a[0] = it->second.a;
relation.a[1] = {x, q};
relation.cnt = 2;
relations.push_back(move(relation));
partials.erase(it);
}
}
}
if (
base > max_base &&
relations.size() <
(size_t)(columns + extra_relations)
) {
break;
}
}
if (relations.size() <= (size_t)columns) {
return 0;
}
const int relation_count = (int)relations.size();
const int combination_words =
(relation_count + 63) >> 6;
struct BasisRow {
vector<u64> parity;
vector<u64> combination;
bool used = false;
};
vector<BasisRow> basis(columns);
vector<vector<u64>> dependencies;
dependencies.reserve(extra_relations);
for (int i = 0; i < relation_count; ++i) {
vector<u64> parity = relations[i].parity;
vector<u64> combination(combination_words);
combination[i >> 6] |= 1ULL << (i & 63);
bool inserted = false;
for (int column = columns - 1; column >= 0; --column) {
if (
((parity[column >> 6] >>
(column & 63)) & 1ULL) == 0
) {
continue;
}
if (!basis[column].used) {
basis[column].used = true;
basis[column].parity = move(parity);
basis[column].combination =
move(combination);
inserted = true;
break;
}
for (int w = 0; w < parity_words; ++w) {
parity[w] ^= basis[column].parity[w];
}
for (int w = 0; w < combination_words; ++w) {
combination[w] ^=
basis[column].combination[w];
}
}
if (!inserted) {
bool zero = true;
for (u64 x : parity) {
if (x) {
zero = false;
break;
}
}
if (zero) {
dependencies.push_back(move(combination));
}
}
}
Mont128 md(n);
for (const vector<u64>& dependency : dependencies) {
vector<int> exponents(columns);
unordered_map<u64, int> large_exponents;
large_exponents.reserve(16);
u128 X = md.r1;
for (int r = 0; r < relation_count; ++r) {
if (
((dependency[r >> 6] >>
(r & 63)) & 1ULL) == 0
) {
continue;
}
for (int ai = 0; ai < relations[r].cnt; ++ai) {
auto [x, q] = relations[r].a[ai];
X = md.mul(
X,
md.init((u128)(m + x) % n)
);
u64 rem = q;
for (int j = 0; j < columns; ++j) {
uint32_t p = factor_base[j].p;
while (rem % p == 0) {
rem /= p;
++exponents[j];
}
if (rem == 1) break;
}
if (rem > 1) {
++large_exponents[rem];
}
}
}
u128 Y = md.r1;
bool valid = true;
for (int j = 0; j < columns; ++j) {
if (exponents[j] & 1) {
valid = false;
break;
}
if (exponents[j]) {
Y = md.mul(
Y,
md.pow_mont(
md.init(factor_base[j].p),
exponents[j] / 2
)
);
}
}
if (!valid) continue;
for (auto [p, e] : large_exponents) {
if (e & 1) {
valid = false;
break;
}
Y = md.mul(
Y,
md.pow_mont(md.init(p), e / 2)
);
}
if (!valid) continue;
u128 xn = md.norm(X);
u128 yn = md.norm(Y);
u128 difference =
xn > yn ? xn - yn : yn - xn;
u128 divisor = gcd128(difference, n);
if (divisor > 1 && divisor < n) {
return divisor;
}
u128 sum = xn + yn;
if (sum >= n) sum -= n;
divisor = gcd128(sum, n);
if (divisor > 1 && divisor < n) {
return divisor;
}
}
return 0;
}
/* Slow-path fallback for a QS failure. */
static u128 rho128(u128 n) {
if (n % 2 == 0) return 2;
auto rng = []() -> u64 {
static u64 x =
0x243f6a8885a308d3ULL ^
(u64)chrono::high_resolution_clock::now()
.time_since_epoch().count();
x ^= x << 7;
x ^= x >> 9;
return x;
};
Mont128 md(n);
auto random_mod = [&]() {
return (((u128)rng() << 64) | rng()) % n;
};
while (true) {
u128 y = md.init(random_mod());
u128 c = md.init(random_mod() + 1);
u128 x = 0, ys = 0;
constexpr u64 M = 256;
u64 r = 1;
u128 g = 1;
auto f = [&](u128 v) {
return md.add(md.mul(v, v), c);
};
while (g == 1) {
x = y;
for (u64 i = 0; i < r; ++i) y = f(y);
for (u64 k = 0; k < r && g == 1; k += M) {
ys = y;
u128 q = md.r1;
u64 lim = min(M, r - k);
for (u64 i = 0; i < lim; ++i) {
y = f(y);
u128 d = x > y ? x - y : y - x;
q = md.mul(q, d);
}
g = gcd128(q, n);
}
r <<= 1;
}
if (g == n) {
do {
ys = f(ys);
u128 d = x > ys ? x - ys : ys - x;
g = gcd128(d, n);
} while (g == 1);
}
if (g != n) return g;
}
}
static void factor_rec(u128 n, vector<u128>& factors) {
if (n == 1) return;
static constexpr uint32_t SMALL[] = {
2, 3, 5, 7, 11, 13, 17, 19, 23,
29, 31, 37, 41, 43, 47, 53, 59,
61, 67, 71, 73, 79, 83, 89, 97
};
for (uint32_t p : SMALL) {
if (n % p == 0) {
factors.push_back(p);
factor_rec(n / p, factors);
return;
}
}
if (is_prime128(n)) {
factors.push_back(n);
return;
}
u64 square_root = isqrt_u128(n);
if ((u128)square_root * square_root == n) {
factor_rec(square_root, factors);
factor_rec(square_root, factors);
return;
}
u128 divisor;
if (n <= numeric_limits<u64>::max()) {
divisor = rho64((u64)n);
} else {
divisor = quadratic_sieve_factor(
n, 7000, 48, 3500000
);
if (divisor == 0 || divisor == n) {
divisor = quadratic_sieve_factor(
n, 10000, 64, 5000000
);
}
if (divisor == 0 || divisor == n) {
divisor = rho128(n);
}
}
factor_rec(divisor, factors);
factor_rec(n / divisor, factors);
}
static inline u128 pow_cap(
u128 a,
int exponent,
u128 cap
) {
u128 result = 1;
for (int i = 0; i < exponent; ++i) {
if (a && result > cap / a) return cap + 1;
result *= a;
}
return result;
}
/*
* Sum over [l, l+w).
* Only e=1,2,3,4 are used.
*/
static u128 power_sum_small(int e, u128 l, u128 w) {
if (e == 1) {
return l * w + w * (w - 1) / 2;
}
if (e == 2) {
return
l * w * (l + w - 1) +
w * (w - 1) / 2 * (2 * w - 1) / 3;
}
if (e == 3) {
return
w * (2 * l + w - 1) / 2 *
l * (l + w - 1) +
w * w * (w - 1) / 2 *
(2 * l + w - 1) / 2;
}
if (e == 4) {
return
l * w * (l + w - 1) *
(l * (l + w - 1) + w * (w - 1)) +
w * (w - 1) / 2 *
(2 * w - 1) / 3 *
(3 * w * (w - 1) - 1) / 5;
}
abort();
}
using Solution = tuple<int, u128, u128>;
static vector<u128> divisors_limited(
const map<u128, int>& factorization,
u128 limit
) {
vector<pair<u128, int>> factors(
factorization.begin(),
factorization.end()
);
vector<u128> divisors;
auto dfs = [&](auto&& self, int index, u128 value)
-> void
{
if (index == (int)factors.size()) {
divisors.push_back(value);
return;
}
auto [prime, exponent] = factors[index];
for (int e = 0; e <= exponent; ++e) {
self(self, index + 1, value);
if (
e == exponent ||
value > limit / prime
) {
break;
}
value *= prime;
}
};
dfs(dfs, 0, 1);
return divisors;
}
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
string input;
cin >> input;
u128 N = parse_u128(input);
vector<u128> prime_factors;
factor_rec(N, prime_factors);
sort(prime_factors.begin(), prime_factors.end());
map<u128, int> factorization;
for (u128 p : prime_factors) {
++factorization[p];
}
vector<Solution> answers;
/*
* E = 1:
*
* 2N = W(2L+W-1).
* The two factors have opposite parity.
*/
{
map<u128, int> factors_2N = factorization;
++factors_2N[2];
u128 twoN = 2 * N;
u64 limit = isqrt_u128(twoN - 1);
vector<u128> widths =
divisors_limited(factors_2N, limit);
for (u128 width : widths) {
u128 other = twoN / width;
if ((width & 1) == (other & 1)) {
continue;
}
u128 left =
(other - width + 1) / 2;
answers.emplace_back(
1,
left,
left + width - 1
);
}
}
/*
* For every L,W:
*
* W | 6*S(2,L,L+W-1)
* W | 2*S(3,L,L+W-1)
* W | 30*S(4,L,L+W-1)
*
* Therefore, for a solution, W divides D[E]*N.
*/
static constexpr int DENOMINATOR[5] = {
0, 0, 6, 2, 30
};
for (int exponent = 2; exponent <= 4; ++exponent) {
map<u128, int> width_factorization =
factorization;
int d = DENOMINATOR[exponent];
for (int p = 2; p * p <= d; ++p) {
while (d % p == 0) {
++width_factorization[p];
d /= p;
}
}
if (d > 1) ++width_factorization[d];
u64 upper = 1;
while (
power_sum_small(exponent, 1, upper) <= N
) {
upper <<= 1;
}
u64 lower = 0;
while (lower + 1 < upper) {
u64 middle = (lower + upper) / 2;
if (
power_sum_small(
exponent, 1, middle
) <= N
) {
lower = middle;
} else {
upper = middle;
}
}
vector<u128> widths =
divisors_limited(
width_factorization,
lower
);
sort(widths.begin(), widths.end());
u64 maximum_left = 1;
while (
pow_cap(maximum_left, exponent, N) <= N
) {
maximum_left <<= 1;
}
for (u128 width : widths) {
u64 lo = 1;
u64 hi = maximum_left;
while (lo < hi) {
u64 middle = lo + (hi - lo) / 2;
if (
power_sum_small(
exponent,
middle,
width
) >= N
) {
hi = middle;
} else {
lo = middle + 1;
}
}
if (
power_sum_small(
exponent,
lo,
width
) == N
) {
answers.emplace_back(
exponent,
lo,
(u128)lo + width - 1
);
}
}
}
/*
* E >= 5:
* sum of all running times is dominated by N^(1/5).
*/
int maximum_exponent = 0;
for (u128 x = N; x >= 2; x >>= 1) {
++maximum_exponent;
}
for (
int exponent = 5;
exponent <= maximum_exponent;
++exponent
) {
vector<u128> powers(1, 0);
for (u64 i = 1;; ++i) {
u128 value =
pow_cap(i, exponent, N);
if (value > N) break;
powers.push_back(value);
}
size_t right = 1;
u128 sum = 0;
for (
size_t left = 1;
left < powers.size();
++left
) {
if (right < left) {
right = left;
sum = 0;
}
while (
right < powers.size() &&
sum + powers[right] <= N
) {
sum += powers[right];
++right;
}
if (sum == N) {
answers.emplace_back(
exponent,
left,
right - 1
);
}
if (right > left) {
sum -= powers[left];
}
}
}
sort(answers.begin(), answers.end());
cout << answers.size() << '\n';
for (auto [e, l, r] : answers) {
cout
<< e << ' '
<< to_string_u128(l) << ' '
<< to_string_u128(r) << '\n';
}
}