Problem 325: Stone Game II
View on Project EulerProject Euler Problem 325 Solution
EulerSolve provides an optimized solution for Project Euler Problem 325, Stone Game II, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We consider ordered piles \((x,y)\) with $$1\le x \lt y \le N.$$ In each move, a player subtracts a positive multiple of the smaller pile from the larger pile, and the result must remain nonnegative. The player who takes the last stone wins. The problem asks for the sum of \(x+y\) over all losing pairs with \(1\le x \lt y\le N\): $$S(N)=\sum (x+y).$$ for the enormous value $$N=10^{16},$$ with the final answer reported modulo $$7^{10}=282475249.$$ Mathematical Approach 1) First isolate the easy winning cases. If \(y\ge 2x\), then the larger pile contains at least two copies of the smaller one, so the current player has several legal Euclid-style subtractions and can force a win immediately. Likewise, if \(x\mid y\), the current player can reduce to a terminal multiple situation. Therefore the only subtle region is $$x \lt y \lt 2x,\qquad x\nmid y.$$ 2) Inside that strip there is only one move. When \(x \lt y \lt 2x\), the only allowed multiple is exactly one copy of the smaller pile. So the game moves deterministically to $$ (x,y)\longrightarrow (y-x,x), $$ after reordering the two piles. This is just one step of Euclid's algorithm. 3) Losing positions form a golden-ratio interval. Write \(y=x+d\) with $$1\le d \lt x.$$ By the previous step, \((x,x+d)\) is losing exactly when \((d,x)\) is winning....
Detailed mathematical approach
Problem Summary
We consider ordered piles \((x,y)\) with
$$1\le x \lt y \le N.$$
In each move, a player subtracts a positive multiple of the smaller pile from the larger pile, and the result must remain nonnegative. The player who takes the last stone wins. The problem asks for the sum of \(x+y\) over all losing pairs with \(1\le x \lt y\le N\):
$$S(N)=\sum (x+y).$$
for the enormous value
$$N=10^{16},$$
with the final answer reported modulo
$$7^{10}=282475249.$$
Mathematical Approach
1) First isolate the easy winning cases.
If \(y\ge 2x\), then the larger pile contains at least two copies of the smaller one, so the current player has several legal Euclid-style subtractions and can force a win immediately. Likewise, if \(x\mid y\), the current player can reduce to a terminal multiple situation. Therefore the only subtle region is
$$x \lt y \lt 2x,\qquad x\nmid y.$$
2) Inside that strip there is only one move.
When \(x \lt y \lt 2x\), the only allowed multiple is exactly one copy of the smaller pile. So the game moves deterministically to
$$ (x,y)\longrightarrow (y-x,x), $$
after reordering the two piles. This is just one step of Euclid's algorithm.
3) Losing positions form a golden-ratio interval.
Write \(y=x+d\) with
$$1\le d \lt x.$$
By the previous step, \((x,x+d)\) is losing exactly when \((d,x)\) is winning. Inductively, the losing positions are characterized by the sharp Beatty boundary
$$x \lt y \le \lfloor \varphi x\rfloor,$$
where
$$\varphi=\frac{1+\sqrt5}{2}.$$
Equivalently, for a fixed \(x\), the admissible differences are
$$1\le d\le a_x,\qquad a_x=\left\lfloor\frac{x}{\varphi}\right\rfloor.$$
The equality of the two forms comes from
$$\varphi=1+\frac{1}{\varphi},\qquad \lfloor \varphi x\rfloor=x+\left\lfloor \frac{x}{\varphi}\right\rfloor.$$
4) Convert the weighted sum to a double sum.
For each fixed \(x\), the losing values are
$$y=x+1,\ x+2,\ \dots,\ x+a_x,$$
unless the global bound \(y\le N\) truncates them. Therefore
$$S(N)=\sum_{x=1}^{N}\sum_{d=1}^{\min(a_x,N-x)}(x+(x+d))$$
or
$$S(N)=\sum_{x=1}^{N}\sum_{d=1}^{\min(a_x,N-x)}(2x+d).$$
5) Prefix-tail split.
Define the cutoff
$$c=\left\lfloor\frac{N+1}{\varphi}\right\rfloor.$$
If \(x\le c\), then the full losing interval still fits below \(N\), because \(x+a_x\le N\). If \(x\gt c\), the interval is cut short by the ceiling \(N\). Hence
$$S(N)=S_{\text{pref}}+S_{\text{tail}}.$$
6) The prefix depends on three Beatty sums.
For \(x\le c\), the inner sum is
$$\sum_{d=1}^{a_x}(2x+d)=2x\,a_x+\frac{a_x(a_x+1)}{2}.$$
Now define
$$G(n)=\sum_{k=1}^{n} a_k,\qquad P(n)=\sum_{k=1}^{n} k\,a_k,\qquad Q(n)=\sum_{k=1}^{n} a_k^2.$$
Then
$$S_{\text{pref}}=\sum_{x=1}^{c}\left(2x\,a_x+\frac{a_x(a_x+1)}{2}\right)=\frac{4P(c)+Q(c)+G(c)}{2}.$$
7) Why \(G,P,Q\) can be computed recursively.
Let
$$m=\left\lfloor\frac{n}{\varphi}\right\rfloor.$$
Classical Beatty partition identities allow us to rewrite the sums up to \(n\) in terms of the same sums up to \(m\). The code uses the exact recurrences
$$G(n)=nm-\frac{m(m+1)}{2}-G(m),$$
$$Q(n)=nm^2-2\sum_{k=1}^{m}k^2-2P(m)+\frac{m(m+1)}{2}+G(m),$$
$$P(n)=\frac{nm(n+1)}{2}-\frac{\left(\sum_{k=1}^{m}k^2+2P(m)+Q(m)+\frac{m(m+1)}{2}+G(m)\right)}{2}.$$
Each recursive step replaces \(n\) by approximately \(n/\varphi\), so the depth is only logarithmic.
8) The tail becomes an ordinary polynomial sum.
For \(x\gt c\), the upper limit is no longer \(a_x\), but \(N-x\). So the inner sum becomes
$$\sum_{d=1}^{N-x}(2x+d)=2x(N-x)+\frac{(N-x)(N-x+1)}{2}.$$
Summing this over \(x=c+1,\dots,N\) reduces the tail to combinations of
$$\sum x,\qquad \sum x^2,$$
which the code evaluates in closed form.
9) Worked example \(N=10\).
The losing pairs are
$$(2,3),(3,4),(4,5),(4,6),(5,6),(5,7),(5,8),(6,7),(6,8),(6,9),$$
$$(7,8),(7,9),(7,10),(8,9),(8,10),(9,10).$$
Their weighted sum is
$$S(10)=211,$$
which matches the checkpoint in the program.
Algorithm
1) Use the golden-ratio characterization of losing positions.
2) Split the total sum at \(c=\lfloor (N+1)/\varphi\rfloor\).
3) Evaluate the prefix via the recursive Beatty sums \(G,P,Q\).
4) Evaluate the tail via polynomial formulas for \(\sum x\) and \(\sum x^2\).
5) Do the same computation modulo \(7^{10}\) for the final output.
Complexity Analysis
Because each Beatty recursion step shrinks \(n\) to roughly \(n/\varphi\), the depth is
$$O(\log N).$$
Each level performs only constant-time big-integer or modular arithmetic, so the whole method is extremely fast even for \(N=10^{16}\).
Checks And Final Result
The code checks
$$S(10)=211,\qquad S(10^4)=230312207313.$$
It also verifies that exact arithmetic and arithmetic modulo \(7^{10}\) agree for \(N=10^6\).
For the target value
$$N=10^{16},$$
the final answer modulo \(7^{10}\) is
$$54672965.$$
Further Reading
- Problem page: https://projecteuler.net/problem=325
- Beatty sequences: https://en.wikipedia.org/wiki/Beatty_sequence
- Euclid-type impartial games: https://en.wikipedia.org/wiki/Euclid's_algorithm#Game_of_Euclid
Problem 325 source code
C++
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <vector>
namespace {
using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;
constexpr u64 kTargetN = 10'000'000'000'000'000ULL;
constexpr u64 kMod = 282'475'249ULL; // 7^10
struct BeattySumsExact {
u128 g = 0; // sum floor(k/phi)
u128 p = 0; // sum k*floor(k/phi)
u128 q = 0; // sum floor(k/phi)^2
};
struct BeattySumsMod {
u64 g = 0;
u64 p = 0;
u64 q = 0;
};
std::string to_string_u128(u128 value) {
if (value == 0) {
return "0";
}
std::string digits;
while (value > 0) {
const int digit = static_cast<int>(value % 10);
digits.push_back(static_cast<char>('0' + digit));
value /= 10;
}
std::reverse(digits.begin(), digits.end());
return digits;
}
u64 isqrt_u128(u128 n) {
if (n == 0) {
return 0;
}
u64 x = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
while ((static_cast<u128>(x) + 1) * (x + 1) <= n) {
++x;
}
while (static_cast<u128>(x) * x > n) {
--x;
}
return x;
}
// floor(n / phi), computed exactly via floor(n*sqrt(5)).
u64 floor_div_phi(u64 n) {
if (n == 0) {
return 0;
}
const u128 nn = static_cast<u128>(n);
const u128 radicand = static_cast<u128>(5) * nn * nn;
const u64 floor_n_sqrt5 = isqrt_u128(radicand);
return (floor_n_sqrt5 - n) / 2;
}
u128 triangular_u128(u64 n) {
return static_cast<u128>(n) * (n + 1) / 2;
}
u128 square_sum_u128(u64 n) {
return static_cast<u128>(n) * (n + 1) * (2 * static_cast<u128>(n) + 1) / 6;
}
u64 mod_add(u64 lhs, u64 rhs, u64 mod) {
lhs += rhs;
if (lhs >= mod) {
lhs -= mod;
}
return lhs;
}
u64 mod_sub(u64 lhs, u64 rhs, u64 mod) {
return (lhs >= rhs) ? (lhs - rhs) : (lhs + mod - rhs);
}
u64 mod_mul(u64 lhs, u64 rhs, u64 mod) {
return static_cast<u64>((static_cast<u128>(lhs) * rhs) % mod);
}
u64 mod_inverse(u64 value, u64 mod) {
i64 t = 0;
i64 new_t = 1;
i64 r = static_cast<i64>(mod);
i64 new_r = static_cast<i64>(value % mod);
while (new_r != 0) {
const i64 quotient = r / new_r;
const i64 next_t = t - quotient * new_t;
t = new_t;
new_t = next_t;
const i64 next_r = r - quotient * new_r;
r = new_r;
new_r = next_r;
}
if (r != 1) {
return 0;
}
if (t < 0) {
t += static_cast<i64>(mod);
}
return static_cast<u64>(t);
}
u64 triangular_mod(u64 n, u64 mod, u64 inv2) {
return mod_mul(mod_mul(n % mod, (n + 1) % mod, mod), inv2, mod);
}
u64 square_sum_mod(u64 n, u64 mod, u64 inv6) {
const u64 a = n % mod;
const u64 b = (n + 1) % mod;
const u64 c = (2 * (n % mod) + 1) % mod;
return mod_mul(mod_mul(mod_mul(a, b, mod), c, mod), inv6, mod);
}
BeattySumsExact beatty_sums_exact(u64 n) {
if (n == 0) {
return BeattySumsExact{};
}
const u64 m = floor_div_phi(n);
const BeattySumsExact child = beatty_sums_exact(m);
const u128 tri_m = triangular_u128(m);
const u128 sq_m = square_sum_u128(m);
const u128 n128 = static_cast<u128>(n);
const u128 m128 = static_cast<u128>(m);
const u128 g = n128 * m128 - tri_m - child.g;
const u128 q = n128 * m128 * m128 - 2 * sq_m - 2 * child.p + tri_m + child.g;
const u128 numerator = sq_m + 2 * child.p + child.q + tri_m + child.g;
const u128 p = n128 * m128 * (n128 + 1) / 2 - numerator / 2;
return BeattySumsExact{g, p, q};
}
BeattySumsMod beatty_sums_mod(u64 n, u64 mod, u64 inv2, u64 inv6) {
if (n == 0) {
return BeattySumsMod{};
}
const u64 m = floor_div_phi(n);
const BeattySumsMod child = beatty_sums_mod(m, mod, inv2, inv6);
const u64 n_mod = n % mod;
const u64 m_mod = m % mod;
const u64 tri_m = triangular_mod(m, mod, inv2);
const u64 sq_m = square_sum_mod(m, mod, inv6);
u64 g = mod_mul(n_mod, m_mod, mod);
g = mod_sub(g, tri_m, mod);
g = mod_sub(g, child.g, mod);
u64 q = mod_mul(n_mod, mod_mul(m_mod, m_mod, mod), mod);
q = mod_sub(q, mod_mul(2 % mod, sq_m, mod), mod);
q = mod_sub(q, mod_mul(2 % mod, child.p, mod), mod);
q = mod_add(q, tri_m, mod);
q = mod_add(q, child.g, mod);
u64 first = mod_mul(mod_mul(n_mod, m_mod, mod), (n + 1) % mod, mod);
first = mod_mul(first, inv2, mod);
u64 numerator = sq_m;
numerator = mod_add(numerator, mod_mul(2 % mod, child.p, mod), mod);
numerator = mod_add(numerator, child.q, mod);
numerator = mod_add(numerator, tri_m, mod);
numerator = mod_add(numerator, child.g, mod);
const u64 second = mod_mul(numerator, inv2, mod);
const u64 p = mod_sub(first, second, mod);
return BeattySumsMod{g, p, q};
}
u128 solve_exact(u64 n) {
const u64 cutoff = floor_div_phi(n + 1);
const BeattySumsExact sums = beatty_sums_exact(cutoff);
const u128 prefix = (4 * sums.p + sums.q + sums.g) / 2;
if (cutoff == n) {
return prefix;
}
const u128 n128 = static_cast<u128>(n);
const u128 count = static_cast<u128>(n - cutoff);
const u128 sum_x = triangular_u128(n) - triangular_u128(cutoff);
const u128 sum_x2 = square_sum_u128(n) - square_sum_u128(cutoff);
const u128 tail_numerator =
count * (n128 * n128 + n128) + (2 * n128 - 1) * sum_x - 3 * sum_x2;
const u128 tail = tail_numerator / 2;
return prefix + tail;
}
u64 solve_mod(u64 n, u64 mod) {
const u64 inv2 = mod_inverse(2, mod);
const u64 inv6 = mod_inverse(6, mod);
const u64 cutoff = floor_div_phi(n + 1);
const BeattySumsMod sums = beatty_sums_mod(cutoff, mod, inv2, inv6);
u64 prefix_numerator = mod_mul(4 % mod, sums.p, mod);
prefix_numerator = mod_add(prefix_numerator, sums.q, mod);
prefix_numerator = mod_add(prefix_numerator, sums.g, mod);
const u64 prefix = mod_mul(prefix_numerator, inv2, mod);
if (cutoff == n) {
return prefix;
}
const u64 n_mod = n % mod;
const u64 count_mod = (n - cutoff) % mod;
const u64 sum_x = mod_sub(triangular_mod(n, mod, inv2),
triangular_mod(cutoff, mod, inv2),
mod);
const u64 sum_x2 = mod_sub(square_sum_mod(n, mod, inv6),
square_sum_mod(cutoff, mod, inv6),
mod);
const u64 n2_plus_n = mod_add(mod_mul(n_mod, n_mod, mod), n_mod, mod);
const u64 coeff = mod_sub(mod_mul(2 % mod, n_mod, mod), 1 % mod, mod);
const u64 term1 = mod_mul(count_mod, n2_plus_n, mod);
const u64 term2 = mod_mul(coeff, sum_x, mod);
const u64 term3 = mod_mul(3 % mod, sum_x2, mod);
u64 tail_numerator = mod_add(term1, term2, mod);
tail_numerator = mod_sub(tail_numerator, term3, mod);
const u64 tail = mod_mul(tail_numerator, inv2, mod);
return mod_add(prefix, tail, mod);
}
bool is_winning_configuration(u64 x, u64 y) {
bool flip = false;
while (true) {
if (y % x == 0 || y / x >= 2) {
return !flip;
}
y -= x; // only possible move when 1 < y/x < 2
if (x > y) {
std::swap(x, y);
}
flip = !flip;
}
}
u128 brute_sum_by_game(u64 n, bool allow_multithreading, unsigned requested_threads = 0) {
if (n < 2) {
return 0;
}
unsigned threads = requested_threads;
if (threads == 0) {
threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 1;
}
}
if (!allow_multithreading || n < 120 || threads <= 1) {
threads = 1;
} else {
threads = std::min<unsigned>(threads, static_cast<unsigned>(n));
}
if (threads == 1) {
u128 total = 0;
for (u64 x = 1; x < n; ++x) {
for (u64 y = x + 1; y <= n; ++y) {
if (!is_winning_configuration(x, y)) {
total += x + y;
}
}
}
return total;
}
std::atomic<u64> next_x{1};
std::vector<u128> partial(threads, 0);
std::vector<std::thread> workers;
workers.reserve(threads);
for (unsigned tid = 0; tid < threads; ++tid) {
workers.emplace_back([&, tid]() {
u128 local = 0;
while (true) {
const u64 x = next_x.fetch_add(1, std::memory_order_relaxed);
if (x >= n) {
break;
}
for (u64 y = x + 1; y <= n; ++y) {
if (!is_winning_configuration(x, y)) {
local += x + y;
}
}
}
partial[tid] = local;
});
}
for (std::thread& worker : workers) {
worker.join();
}
u128 total = 0;
for (const u128 value : partial) {
total += value;
}
return total;
}
bool run_validation_checkpoints(bool allow_multithreading, unsigned requested_threads) {
{
constexpr u64 kBruteN = 220;
const u128 brute = brute_sum_by_game(kBruteN, allow_multithreading, requested_threads);
const u128 fast = solve_exact(kBruteN);
if (brute != fast) {
std::cerr << "Brute-force checkpoint failed for N=" << kBruteN << ": got "
<< to_string_u128(fast) << ", expected " << to_string_u128(brute)
<< "\n";
return false;
}
}
{
constexpr u64 kN = 10;
constexpr u64 kExpected = 211;
const u128 got = solve_exact(kN);
if (got != kExpected) {
std::cerr << "Checkpoint failed for S(10): got " << to_string_u128(got)
<< ", expected " << kExpected << "\n";
return false;
}
}
{
constexpr u64 kN = 10'000;
constexpr u64 kExpected = 230'312'207'313ULL;
const u128 got = solve_exact(kN);
if (got != kExpected) {
std::cerr << "Checkpoint failed for S(10^4): got " << to_string_u128(got)
<< ", expected " << kExpected << "\n";
return false;
}
}
{
constexpr u64 kN = 1'000'000;
const u128 exact = solve_exact(kN);
const u64 exact_mod = static_cast<u64>(exact % kMod);
const u64 fast_mod = solve_mod(kN, kMod);
if (fast_mod != exact_mod) {
std::cerr << "Modulo checkpoint failed for N=" << kN << ": got " << fast_mod
<< ", expected " << exact_mod << "\n";
return false;
}
}
return true;
}
} // namespace
int main() {
unsigned threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 4;
}
if (!run_validation_checkpoints(true, threads)) {
return 1;
}
const u64 answer = solve_mod(kTargetN, kMod);
std::cout << answer << '\n';
return 0;
}
Python
import math
MOD = 282475249
def floor_div_phi(n):
if n == 0:
return 0
radicand = 5 * n * n
floor_n_sqrt5 = math.isqrt(radicand)
return (floor_n_sqrt5 - n) // 2
def triangular_mod(n, mod_val, inv2):
return (n % mod_val) * ((n + 1) % mod_val) % mod_val * inv2 % mod_val
def square_sum_mod(n, mod_val, inv6):
a = n % mod_val
b = (n + 1) % mod_val
c = (2 * (n % mod_val) + 1) % mod_val
return a * b % mod_val * c % mod_val * inv6 % mod_val
def beatty_sums_mod(n, mod_val, inv2, inv6):
if n == 0:
return (0, 0, 0)
m = floor_div_phi(n)
child_g, child_p, child_q = beatty_sums_mod(m, mod_val, inv2, inv6)
n_mod = n % mod_val
m_mod = m % mod_val
tri_m = triangular_mod(m, mod_val, inv2)
sq_m = square_sum_mod(m, mod_val, inv6)
g = (n_mod * m_mod - tri_m - child_g) % mod_val
q = (n_mod * m_mod * m_mod - 2 * sq_m - 2 * child_p + tri_m + child_g) % mod_val
first = n_mod * m_mod % mod_val * ((n + 1) % mod_val) % mod_val * inv2 % mod_val
numerator = (sq_m + 2 * child_p + child_q + tri_m + child_g) % mod_val
second = numerator * inv2 % mod_val
p = (first - second) % mod_val
return ((g + mod_val) % mod_val, (p + mod_val) % mod_val, (q + mod_val) % mod_val)
def solve_mod(n, mod_val):
inv2 = pow(2, -1, mod_val)
inv6 = pow(6, -1, mod_val)
cutoff = floor_div_phi(n + 1)
sums_g, sums_p, sums_q = beatty_sums_mod(cutoff, mod_val, inv2, inv6)
prefix_numerator = (4 * sums_p + sums_q + sums_g) % mod_val
prefix = prefix_numerator * inv2 % mod_val
if cutoff == n:
return prefix
n_mod = n % mod_val
count_mod = (n - cutoff) % mod_val
sum_x = (triangular_mod(n, mod_val, inv2) - triangular_mod(cutoff, mod_val, inv2)) % mod_val
sum_x2 = (square_sum_mod(n, mod_val, inv6) - square_sum_mod(cutoff, mod_val, inv6)) % mod_val
n2_plus_n = (n_mod * n_mod + n_mod) % mod_val
coeff = (2 * n_mod - 1) % mod_val
term1 = count_mod * n2_plus_n % mod_val
term2 = coeff * sum_x % mod_val
term3 = 3 * sum_x2 % mod_val
tail_numerator = (term1 + term2 - term3) % mod_val
tail = tail_numerator * inv2 % mod_val
ans = (prefix + tail) % mod_val
return (ans + mod_val) % mod_val
def solve():
target_n = 10000000000000000
mod_val = 282475249
ans = solve_mod(target_n, mod_val)
return str(ans)
if __name__ == '__main__':
print(solve())
Java
import java.math.BigInteger;
public class Euler325 {
static final long TARGET_N = 10000000000000000L;
static final long MOD = 282475249L;
static long floorDivPhi(long n) {
if (n == 0)
return 0;
BigInteger nn = BigInteger.valueOf(n);
BigInteger radicand = BigInteger.valueOf(5).multiply(nn).multiply(nn);
BigInteger floorNSqrt5 = radicand.sqrt();
return floorNSqrt5.subtract(nn).divide(BigInteger.valueOf(2)).longValue();
}
static long modAdd(long a, long b, long m) {
return (a + b) % m;
}
static long modSub(long a, long b, long m) {
return (a >= b) ? (a - b) : (a + m - b);
}
static long modMul(long a, long b, long m) {
return (a * b) % m;
}
static long modPow(long base, long exp, long m) {
long res = 1;
base %= m;
while (exp > 0) {
if ((exp & 1) != 0)
res = modMul(res, base, m);
base = modMul(base, base, m);
exp >>= 1;
}
return res;
}
static long triangularMod(long n, long m, long inv2) {
return modMul(modMul(n % m, (n + 1) % m, m), inv2, m);
}
static long squareSumMod(long n, long m, long inv6) {
long a = n % m;
long b = (n + 1) % m;
long c = (2 * a + 1) % m;
return modMul(modMul(modMul(a, b, m), c, m), inv6, m);
}
static class BeattySums {
long g, p, q;
BeattySums(long g, long p, long q) {
this.g = g;
this.p = p;
this.q = q;
}
}
static BeattySums beattySumsMod(long n, long mod, long inv2, long inv6) {
if (n == 0)
return new BeattySums(0, 0, 0);
long m = floorDivPhi(n);
BeattySums child = beattySumsMod(m, mod, inv2, inv6);
long nMod = n % mod;
long mMod = m % mod;
long triM = triangularMod(m, mod, inv2);
long sqM = squareSumMod(m, mod, inv6);
long g = modMul(nMod, mMod, mod);
g = modSub(g, triM, mod);
g = modSub(g, child.g, mod);
long q = modMul(nMod, modMul(mMod, mMod, mod), mod);
q = modSub(q, modMul(2, sqM, mod), mod);
q = modSub(q, modMul(2, child.p, mod), mod);
q = modAdd(q, triM, mod);
q = modAdd(q, child.g, mod);
long first = modMul(modMul(nMod, mMod, mod), (n + 1) % mod, mod);
first = modMul(first, inv2, mod);
long numerator = sqM;
numerator = modAdd(numerator, modMul(2, child.p, mod), mod);
numerator = modAdd(numerator, child.q, mod);
numerator = modAdd(numerator, triM, mod);
numerator = modAdd(numerator, child.g, mod);
long second = modMul(numerator, inv2, mod);
long p = modSub(first, second, mod);
return new BeattySums(g, p, q);
}
static long solveMod(long n, long mod) {
BigInteger modBig = BigInteger.valueOf(mod);
long inv2 = BigInteger.valueOf(2).modInverse(modBig).longValue();
long inv6 = BigInteger.valueOf(6).modInverse(modBig).longValue();
long cutoff = floorDivPhi(n + 1);
BeattySums sums = beattySumsMod(cutoff, mod, inv2, inv6);
long prefixNum = modMul(4, sums.p, mod);
prefixNum = modAdd(prefixNum, sums.q, mod);
prefixNum = modAdd(prefixNum, sums.g, mod);
long prefix = modMul(prefixNum, inv2, mod);
if (cutoff == n)
return prefix;
long nMod = n % mod;
long countMod = (n - cutoff) % mod;
long sumX = modSub(triangularMod(n, mod, inv2), triangularMod(cutoff, mod, inv2), mod);
long sumX2 = modSub(squareSumMod(n, mod, inv6), squareSumMod(cutoff, mod, inv6), mod);
long n2PlusN = modAdd(modMul(nMod, nMod, mod), nMod, mod);
long coeff = modSub(modMul(2, nMod, mod), 1, mod);
long term1 = modMul(countMod, n2PlusN, mod);
long term2 = modMul(coeff, sumX, mod);
long term3 = modMul(3, sumX2, mod);
long tailNum = modAdd(term1, term2, mod);
tailNum = modSub(tailNum, term3, mod);
long tail = modMul(tailNum, inv2, mod);
return modAdd(prefix, tail, mod);
}
public static String solve() {
return String.valueOf(solveMod(TARGET_N, MOD));
}
public static void main(String[] args) {
System.out.println(solve());
}
}