Problem 374: Maximum Integer Partition Product

Project Euler Problem 374 Solution

EulerSolve provides an optimized solution for Project Euler Problem 374, Maximum Integer Partition Product, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For each positive integer $n$, let $f(n)$ be the largest product obtainable from a partition of $n$ into distinct positive integers, and let $m(n)$ be the number of parts in an optimal partition. Project Euler 374 asks for $$M(N)=\sum_{n=1}^{N} f(n)m(n)$$ at $N=10^{14}$, reduced modulo $$P=982451653.$$ The local C++, Python, and Java solutions all exploit the same fact: optimal partitions fall into simple $\Theta(\sqrt{n})$-sized blocks, so the program never searches over partitions explicitly. Mathematical Approach Step 1: The optimal partition is almost consecutive Write an optimal partition with $k$ parts as $$a_1 \lt a_2 \lt \cdots \lt a_k,\qquad a_1+\cdots+a_k=n.$$ If some adjacent gap is at least $3$, say $a_{i+1}-a_i\ge 3$, then replacing $(a_i,a_{i+1})$ by $(a_i+1,a_{i+1}-1)$ keeps the sum fixed, preserves distinctness, and increases the product because $$ (a_i+1)(a_{i+1}-1)-a_i a_{i+1} = a_{i+1}-a_i-1 \gt 0. $$ Therefore every optimal partition has adjacent differences only $1$ or $2$. For $n\ge 2$, the maximizing partitions are thus “near-consecutive”: a run of consecutive integers with at most one missing value....

Detailed mathematical approach

Problem Summary

An integer partition of $n$ writes $n$ as a sum of positive integers, two sums differing only in the order of their summands counting as the same. ^[2] A partition into distinct parts uses every part at most once, so those of $5$ are $5$, $4+1$ and $3+2$. Let $f(n)$ be the largest product of the parts of such a partition, and $m(n)$ the number of parts of one that attains it. For $n=5$ those products are $5$, $4$ and $6$, so $f(5)=6$ and $m(5)=2$; for $n=10$ the best partition is $10=2+3+5$, so $f(10)=30$, $m(10)=3$ and $f(10)\,m(10)=90$. We must compute

$$S(N)=\sum_{n=1}^{N}f(n)\,m(n),\qquad N=10^{14},$$

modulo $P=982451653$, which the statement calls the fifty millionth prime; it also supplies the checkpoint $S(100)=1683550844462$. ^[1]

The whole solution in one sentence: for every $n\ge2$ the best partition is a run of consecutive integers starting at $2$ or $3$ with at most one value left out, so $f(n)$ is a factorial divided by a small integer, and the $10^{14}$ values of $n$ split into about $1.4\cdot10^{7}$ blocks inside which those factorials add up to one factorial times a short sum of reciprocals.

Mathematical Approach

One structural fact carries everything. We read it off small cases, prove it by four exchange arguments and a count of candidates, then sum block by block.

1. The first nineteen values

Enumerating all partitions into distinct parts gives the maximisers below; in each case exactly one partition attains the largest product. ^[3]

$n=1,2,3,4$: the single part $n$, so $f(n)=n$ and $m(n)=1$;
$n=5,\dots,8$: $2+3$, $2+4$, $3+4$, $3+5$, with $f=6,8,12,15$ and $m(n)=2$;
$n=9,\dots,13$: $2+3+4$, $2+3+5$, $2+4+5$, $3+4+5$, $3+4+6$, with $f=24,30,40,60,72$ and $m(n)=3$;
$n=14,\dots,19$: $2+3+4+5$, $2+3+4+6$, $2+3+5+6$, $2+4+5+6$, $3+4+5+6$, $3+4+5+7$, with $f=120,144,180,240,360,420$ and $m(n)=4$.

The part count is constant on runs of increasing length, and inside a run one value is missing from an interval starting at $2$ or $3$, sliding downwards as $n$ grows.

2. The part 1 is used only for $n=1$

Suppose a partition of some $n\ge2$ into distinct parts contains the part $1$. It has a second part; let $a$ be the largest. Replacing $1$ and $a$ by the single part $a+1$ keeps the sum, keeps all parts distinct because $a+1$ exceeds every remaining part, and raises the product from $a$ to $a+1$. A maximising partition of $n\ge2$ therefore has all parts $\ge2$; the part $1$ only helps when $n=1$.

3. Neighbouring parts differ by at most 2

Let $a\lt b$ be parts with no part strictly between them, and suppose $b-a\ge3$. Then $a+1$ and $b-1$ lie strictly between $a$ and $b$, so neither is a part, and they differ. Replacing $a$ and $b$ by them keeps the sum, while

$$(a+1)(b-1)-ab=b-a-1\ge2\gt0,$$

so the product strictly grows. Hence neighbouring parts differ by $1$ or $2$: the parts fill an interval except for some missing values.

4. At most one value is missing

Suppose two gaps of size $2$ occur: $a$, $a+2$, $c$, $c+2$ are parts while $a+1$ and $c+1$ are not, with $a+1\lt c+1$. Since $a+2$ is a part and $c+1$ is not, $a+2\ne c+1$, so $c\ge a+2$. Now replace $a$ and $c+2$ by $a+1$ and $c+1$: both were missing, they differ, the sum is unchanged, and

$$(a+1)(c+1)-a(c+2)=c-a+1\ge3\gt0.$$

So a maximising partition has at most one gap of size $2$: its parts are all integers of an interval $\{c,c+1,\dots,D\}$ with at most one interior value $h$ removed.

5. The smallest part is 2 or 3

Let $n\ge5$, with maximising parts $\{c,\dots,D\}$ less at most one interior $h$, and assume $c\ge4$. Each case builds a partition of $n$ with a strictly larger product, which is impossible.

$c\ge5$: replace the part $c$ by $2$ and $c-2$. Both are below $c$, hence not parts, and they differ because $c\ne4$; the product is multiplied by $2(c-2)/c\gt1$.
$c=4$ with a missing $h$: then $h\ge5$ and $h+1$ is a part, because $h$ is interior. Replace the parts $4$ and $h+1$ by $2$, $3$ and $h$: the sum is unchanged, the three new values were not parts, and the product is multiplied by $6h/(4h+4)\gt1$.
$c=4$ with nothing missing: the parts are $\{4,\dots,D\}$ with $D\ge5$, since $D=4$ means $n=4$. Now $\{2,3,\dots,D\}\setminus\{5\}$ has the same sum and the product $D!/5$, beating $D!/6$.

So for $n\ge5$ the maximising parts are an interval starting at $2$ or $3$ with at most one interior value removed; the four values $n\le4$ are listed in Section 1.

6. How many parts are possible

By Section 2 a maximising partition of $n\ge2$ has all parts $\ge2$, so the smallest sum of $k$ such parts is

$$T_k=2+3+\cdots+(k+1)=\frac{k(k+3)}{2}.$$

The sequence $T_1=2,T_2=5,T_3=9,T_4=14,\dots$ is strictly increasing, so such a partition has at most $m$ parts, where $m$ is the largest index with $T_m\le n$. Because

$$T_{m+1}-T_m=m+2,$$

the integers $n$ with a given index $m$ form a block of $m+2$ consecutive values, and we write

$$n=T_m+s,\qquad 0\le s\le m+1.$$

Solving $m(m+3)/2\le n$ gives $m=\lfloor(\sqrt{8n+9}-3)/2\rfloor$. Throughout, $m$ without an argument is this block index, $m(n)$ the quantity of the statement; Section 7 shows that they agree.

7. The complete list of candidates

Fix $k\ge1$ and list every set allowed by Sections 3 to 5 that has exactly $k$ parts. The shape fixes the top end: $k$ parts from $c$ reach $c+k-1$ with nothing removed and $c+k$ with one value removed, so $c=2$ ends at $k+1$ or $k+2$, and $c=3$ at $k+2$ or $k+3$. Two families cover all four cases. The first is $\{2,\dots,k+2\}\setminus\{h\}$ with $2\le h\le k+2$: here $h=k+2$ gives the plain run $\{2,\dots,k+1\}$, $h=2$ the run $\{3,\dots,k+2\}$, and every other $h$ an interval from $2$ with one interior value removed. Its sum is $T_k+(k+2-h)$ and its product $(k+2)!/h$, so it realises each sum $T_k+s'$ with $0\le s'\le k$ exactly once. The second is $\{3,\dots,k+3\}\setminus\{h\}$ with $4\le h\le k+2$, an interval from $3$ with one interior value removed; its sum is $T_k+(2k+3-h)$ and its product $(k+3)!/(2h)$, covering $k+1\le s'\le 2k-1$, again once each. No set of this shape with $k$ parts has any other sum.

Now take $n=T_m+s$ with $n\ge5$ and ask which $k$ occur.

$k\le m-2$ is impossible: the offset would need $s'=n-T_k\ge T_{k+2}-T_k=2k+5$, beyond the largest offset $2k-1$ either family offers.
$k=m-1$ needs $s'=s+m+1\le2m-3$, that is $s\le m-4$, and then the candidate is $\{3,\dots,m+2\}\setminus\{m-s\}$ with the product $(m+2)!/\bigl(2(m-s)\bigr)$.
$k=m$ always gives exactly one candidate: for $0\le s\le m$ it is $\{2,\dots,m+2\}\setminus\{m+2-s\}$ with the product $(m+2)!/(m+2-s)$, and for $s=m+1$ it is $\{3,\dots,m+3\}\setminus\{m+2\}$ with the product $(m+3)!/\bigl(2(m+2)\bigr)$.

When the $k=m-1$ candidate exists, the common factorial $(m+2)!$ cancels from the ratio of the two products, leaving

$$\frac{2(m-s)}{m+2-s},$$

which exceeds $1$ exactly when $s\lt m-2$; the candidate exists only for $s\le m-4$. So the $k=m$ candidate always wins strictly: for every $n\ge5$ the maximising partition is unique, it has $m(n)=m$ parts, and

$$f(n)=\begin{cases}\dfrac{(m+2)!}{m+2-s}, & 0\le s\le m,\\[4pt] \dfrac{(m+3)!}{2(m+2)}, & s=m+1.\end{cases}$$

Uniqueness is what makes the statement's phrase the number of elements well defined. The two lines also give $f(2)=2$, $f(3)=3$, $f(4)=4$, so the formula holds for every $n\ge2$; only $n=1$ stays outside.

8. Worked examples

For $n=10$, $T_3=9\le10\lt T_4=14$ gives $m=3$, hence $s=1\le m$ and the missing value $m+2-s=4$. The partition is $\{2,3,4,5\}\setminus\{4\}=\{2,3,5\}$, with $f(10)=5!/4=30$ and $m(10)=3$ as the statement says.

For $n=13$, still $m=3$ but now $s=4=m+1$, the second line applies: the partition is $\{3,4,5,6\}\setminus\{5\}=\{3,4,6\}$ with $f(13)=6!/(2\cdot5)=72$ and $f(13)\,m(13)=216$, as in Section 1.

For $n=14$ we have $m=4$ and $s=0\le m-4$, the smallest $n$ with two candidates: $\{2,3,4,5\}$ with product $120$ and four parts, and $\{3,5,6\}$ with product $90$ and three parts. Their ratio is $2(m-s)/(m+2-s)=8/6=120/90$, as Section 7 predicts.

9. Summing one whole block

Inside the block of index $m$ the part count is the constant $m$, so the block contributes

$$B_m=m\sum_{s=0}^{m+1}f(T_m+s).$$

For $0\le s\le m$ substitute $d=m+2-s$, running over $2,\dots,m+2$; the value $s=m+1$ is added separately, giving

$$\begin{aligned}B_m&=m\,(m+2)!\sum_{d=2}^{m+2}\frac1d\\&\quad+\frac{m\,(m+3)!}{2(m+2)}.\end{aligned}$$

The reciprocal sum is the harmonic number $H_{m+2}$ minus its first term. ^[5] For $m=3$ it equals $77/60$, so the two terms are $3\cdot120\cdot77/60=462$ and $3\cdot720/10=216$, with the total $678$; the direct sum $3(24+30+40+60+72)$ of Section 1 is $678$ as well.

10. Working modulo the fifty millionth prime

These factorials are hopeless to carry exactly: $f(10^{14})$ is $14142136!$ divided by $12392315$, which a sum of logarithms puts at $94981701$ decimal digits, the order Stirling's formula predicts, ^[6] and the grand total is larger still. So everything is reduced modulo the prime $P=982451653$. Division is then multiplication by a modular inverse: the inverse of $d$ is the residue $d^{-1}$ with $d\,d^{-1}\equiv1\pmod P$. Every index here is at most $m+3\le14142137$, far below $P$, so each denominator $d$ of Section 9 is a nonzero residue and invertible. Start from $1^{-1}=1$; for $2\le i\lt P$ write $P=qi+r$ with $q=\lfloor P/i\rfloor$ and $r=P\bmod i$, where $1\le r\lt i$ because the prime $P$ has no divisor $i$ in that range. Then $r\equiv-qi$, and multiplying by $i^{-1}r^{-1}$ gives

$$i^{-1}\equiv-q\,r^{-1}\pmod P.$$

Since $r\lt i$, the inverse on the right is already known, so one loop fills the table. ^[4]

11. The last block

Every block but the last is complete. A cut-off last block keeps its first $\ell$ values with $1\le\ell\le m+1$, whose offsets $s=0,\dots,\ell-1$ become, under $d=m+2-s$, the denominators from $m+3-\ell$ to $m+2$: summing that stretch of inverses in place of the full reciprocal sum truncates the block correctly, and the extra term of $s=m+1$ is added only when the block is complete. For $N=11$, say, the $m=3$ block keeps $n=9,10,11$ and contributes

$$3\cdot120\left(\tfrac15+\tfrac14+\tfrac13\right)=282,$$

which is $3(24+30+40)$. For $N=10^{14}$ the largest index is $m=14142134$, its block starts at $T_m=99999998250179$, and only $1749822$ of the $14142136$ values of that block are at most $N$.

How the Code Works

The C++ program ^[9] opens with two helpers for Section 6. isqrt_u64 takes the floating-point square root as a first guess and corrects it with two loops that compare by division, r + 1 <= n / (r + 1) and r > n / r, so no product can overflow. max_m_for_n evaluates $\lfloor(\sqrt{8n+9}-3)/2\rfloor$ with that exact root and then raises or lowers $m$ until $T_m\le n\lt T_{m+1}$; with an exact root no correction is ever needed. Java's isqrt is the same routine from Math.sqrt, and it is the correction loops that make it exact: even an input below $2^{53}$, hence exactly representable as a double, can have a square root that rounds up past an integer. At the target, $8N+9=800000000000009$ happens to need no correction. Python calls math.isqrt, the floor of the exact square root. ^[8]

value_exact_small returns $f(n)\,m(n)$ exactly, through Section 7 and factorial_u128 in unsigned __int128 ^[7] and serves only run_checkpoints, which verifies $f(5)\,m(5)=12$, $f(10)\,m(10)=90$ and $S(100)=1683550844462$. It branches on $s\le m-1$, $s=m$ and $s=m+1$; the middle branch is the case $d=2$ of the first, so it is redundant but correct. For $n\le100$ the largest factorial it ever evaluates is $14!$, far inside the 128-bit range, and n == 0 and n == 1 return earlier.

sum_mod is the main computation, following Sections 9 to 11. It starts the answer at $1$ for $n=1$, computes $m_{\max}$, and fills the inverse table inv up to $m_{\max}+3$ by the recurrence of Section 10. It then walks $m=1,\dots,m_{\max}$, maintaining fact as $(m+2)!$ and harmonic as the reciprocal sum through the cursor upto, which only ever increases, so both cost constant time per block. With base set to $m\,(m+2)!$ it computes the block length len from $T_m$, $T_{m+1}-1$ and $N$, and picks one of three branches: a short block multiplies base by the partial inverse sum of Section 11; len == m + 1 multiplies it by the full harmonic; and a complete block adds on top of that the last term $m\,(m+3)!/(2(m+2))$ through inv2 and inv[m + 2]. main accepts --skip-checkpoints and rejects anything else.

The three languages differ only in integer handling. C++ keeps the exact checkpoint arithmetic in unsigned __int128 and widens each modular multiplication to signed __int128 before reducing, which is more conservative than needed: the largest product of two residues is $(P-1)^2\approx9.65\cdot10^{17}$, still inside signed 64-bit range. Java keeps inv as an int[], every residue being below $P\lt2^{31}$, and defers the reduction inside the partial inverse sum, which stays below $(m_{\max}+1)(P-1)\approx1.4\cdot10^{16}$ and so fits a long. Python needs no bounds, its integers being unbounded. Only the C++ file carries the checkpoints.

Complexity and Verification

Since $T_m\sim m^2/2$, the number of blocks is $m_{\max}=\Theta(\sqrt N)$, which is $14142134$ for $N=10^{14}$. The inverse table holds $m_{\max}+3=14142137$ inverses in $m_{\max}+4$ allocated slots, index $0$ staying unused; each block costs a constant number of modular operations, and the truncated block adds at most $m_{\max}$ further additions, so with a modular operation counted as constant time the method runs in $O(\sqrt N)$ time and memory. The table itself occupies $113137104$ bytes of 64-bit entries in C++ and $56568552$ bytes of 32-bit entries in Java, while Python stores references to integer objects. Measured on one desktop machine: C++ $0.65$ s at $115$ MB peak, Java $0.17$ s at $103$ MB, Python $9.7$ s at $581$ MB. All three print the same result.

The following was checked for this article. A brute force over all partitions into distinct parts for $n\le60$ reproduces $f(n)$ and $m(n)$ from Section 7, with a single maximiser for every $n$. A knapsack dynamic program over the parts $1,\dots,n$, which never uses the structure above, agrees for every $n\le2000$, again one maximiser each, and reproduces the checkpoint $S(100)$. The block method with its truncation, modulo $P$, equals a direct term-by-term sum for every $N\le1199$ and for $N=5000$, $12345$, $99999$, $250000$, exercising all three branches. The closed form also agrees for $2\le n\le3000$ with the independent closed form published in the OEIS entry. ^[3] The inverse table satisfies $i\cdot\mathrm{inv}[i]\equiv1\pmod P$, and the block decomposition is complete:

$$1+\sum_{m=1}^{m_{\max}-1}(m+2)+1749822=10^{14}.$$

Footnotes and References

The references below supply the background facts; the structure of the maximising partition, the candidate count, the block formula and the truncation rule are derived explicitly above.

Project Euler 374 — Maximum Integer Partition Product. The official statement defines f(n) as the maximum product of the parts of a partition of n into distinct parts and m(n) as the number of elements of such a partition, gives f(5) = 6, m(5) = 2, f(10) = 30, m(10) = 3 and the checkpoint 1683550844462 for the sum of f(n)·m(n) over 1 ≤ n ≤ 100, and asks for that sum over 1 ≤ n ≤ 10^14 modulo 982451653, which it calls the 50 millionth prime. Its primality was checked here.
Wikipedia — Partition (number theory). It defines a partition of n as a way of writing n as a sum of positive integers, states that two sums differing only in the order of their summands are the same partition, and defines partitions with distinct parts, in which no part repeats.
OEIS A034893 — maximum product of a partition into distinct parts. The entry tabulates the maximal products 1, 1, 2, 3, 4, 6, 8, 12, 15, 24 from offset 0, the values f(n) of Section 1, comments that the maximiser lies among the partitions without a part 1 that have the most parts, and gives a closed form in k = ⌊(√(8n+9) − 1)/2⌋, checked here against Section 7 for all 2 ≤ n ≤ 3000.
cp-algorithms — Modular inverse. The page derives the linear-time table of inverses from m = ka + r, giving r ≡ −ka and hence a⁻¹ ≡ −k·r⁻¹, and notes that it needs a prime modulus, since otherwise r need not be invertible.
Wikipedia — Harmonic number. The n-th harmonic number is the sum of the reciprocals of the first n natural numbers, which for n = m + 2, less its first term, is the reciprocal sum of Section 9.
Wikipedia — Stirling's approximation. It gives n! ∼ √(2πn)·(n/e)ⁿ and ln n! = n ln n − n + O(ln n), which fixes the order of the factorials in Section 10.
GCC manual — 128-bit integers. The types __int128 and unsigned __int128 are 128 bits wide and exist on targets with an integer mode that wide; the C++ file uses the unsigned form for the exact checkpoints and the signed form to widen modular multiplications.
Python documentation — integer mathematical functions. isqrt(n) returns the integer square root of a nonnegative n, the floor of its exact square root, equivalently the greatest a with a² ≤ n; the published Python file reaches it through the classic math.isqrt alias.
C++ — Euler374.cpp. The linked immutable revision contains isqrt_u64, max_m_for_n, factorial_u128, value_exact_small, sum_mod and run_checkpoints as described above, with the --skip-checkpoints switch; it is the implementation record for this article.

Problem 374 source code

C++

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>
#include <cmath>
#include <functional>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr i64 kMod = 982451653LL;
constexpr u64 kLimit = 100000000000000ULL;

u64 isqrt_u64(const u64 n) {
    u64 r = static_cast<u64>(std::sqrt(static_cast<long double>(n)));
    while ((r + 1ULL) <= n / (r + 1ULL)) {
        ++r;
    }
    while (r > n / r) {
        --r;
    }
    return r;
}

u64 max_m_for_n(const u64 n) {
    // Largest m with m(m+3)/2 <= n.
    u64 m = (isqrt_u64(8ULL * n + 9ULL) - 3ULL) / 2ULL;
    while ((m + 1ULL) * (m + 4ULL) / 2ULL <= n) {
        ++m;
    }
    while (m > 0ULL && m * (m + 3ULL) / 2ULL > n) {
        --m;
    }
    return m;
}

u128 factorial_u128(const u64 n) {
    u128 f = 1U;
    for (u64 i = 2; i <= n; ++i) {
        f *= static_cast<u128>(i);
    }
    return f;
}

u128 value_exact_small(const u64 n) {
    if (n == 0ULL) {
        return 0U;
    }
    if (n == 1ULL) {
        return 1U;
    }

    const u64 m = max_m_for_n(n);
    const u64 t = m * (m + 3ULL) / 2ULL;
    const u64 s = n - t;

    const u128 fact_m2 = factorial_u128(m + 2ULL);
    if (s <= m - 1ULL) {
        return static_cast<u128>(m) * (fact_m2 / static_cast<u128>(m + 2ULL - s));
    }
    if (s == m) {
        return static_cast<u128>(m) * (fact_m2 / 2U);
    }
    return static_cast<u128>(m) * (factorial_u128(m + 3ULL) / (2U * static_cast<u128>(m + 2ULL)));
}

i64 sum_mod(const u64 n) {
    if (n == 0ULL) {
        return 0;
    }

    i64 answer = 1;  // n = 1
    const u64 m_max = max_m_for_n(n);
    const u64 inv_limit = m_max + 3ULL;

    std::vector<i64> inv(static_cast<std::size_t>(inv_limit + 1ULL), 0);
    inv[1] = 1;
    for (u64 i = 2ULL; i <= inv_limit; ++i) {
        inv[static_cast<std::size_t>(i)] =
            (kMod - static_cast<i64>((static_cast<__int128>(kMod / static_cast<i64>(i)) *
                                      inv[static_cast<std::size_t>(kMod % static_cast<i64>(i))]) %
                                     kMod)) %
            kMod;
    }

    const i64 inv2 = (kMod + 1LL) / 2LL;

    i64 fact = 1;          // running factorial up to `upto`
    i64 harmonic = 0;      // sum_{d=2..upto} inv[d] mod kMod
    u64 upto = 1ULL;

    for (u64 m = 1ULL; m <= m_max; ++m) {
        const u64 target = m + 2ULL;
        while (upto < target) {
            ++upto;
            fact = static_cast<i64>((static_cast<__int128>(fact) * static_cast<i64>(upto)) % kMod);
            if (upto >= 2ULL) {
                harmonic += inv[static_cast<std::size_t>(upto)];
                if (harmonic >= kMod) {
                    harmonic -= kMod;
                }
            }
        }

        const u64 start = m * (m + 3ULL) / 2ULL;
        const u64 end_full = (m + 1ULL) * (m + 4ULL) / 2ULL - 1ULL;
        const u64 len = std::min(n, end_full) - start + 1ULL;

        const i64 mm = static_cast<i64>(m % static_cast<u64>(kMod));
        const i64 base = static_cast<i64>((static_cast<__int128>(mm) * fact) % kMod);

        i64 add = 0;
        if (len <= m) {
            const u64 d1 = m + 3ULL - len;
            const u64 d2 = m + 2ULL;
            i64 partial_h = 0;
            for (u64 d = d1; d <= d2; ++d) {
                partial_h += inv[static_cast<std::size_t>(d)];
                if (partial_h >= kMod) {
                    partial_h -= kMod;
                }
            }
            add = static_cast<i64>((static_cast<__int128>(base) * partial_h) % kMod);
        } else if (len == m + 1ULL) {
            add = static_cast<i64>((static_cast<__int128>(base) * harmonic) % kMod);
        } else {  // len == m + 2
            add = static_cast<i64>((static_cast<__int128>(base) * harmonic) % kMod);
            const i64 fact_m3 = static_cast<i64>((static_cast<__int128>(fact) * static_cast<i64>(m + 3ULL)) % kMod);
            i64 extra = static_cast<i64>((static_cast<__int128>(mm) * fact_m3) % kMod);
            extra = static_cast<i64>((static_cast<__int128>(extra) * inv2) % kMod);
            extra = static_cast<i64>((static_cast<__int128>(extra) * inv[static_cast<std::size_t>(m + 2ULL)]) % kMod);
            add += extra;
            if (add >= kMod) {
                add -= kMod;
            }
        }

        answer += add;
        if (answer >= kMod) {
            answer -= kMod;
        }
    }

    return answer;
}

bool run_checkpoints() {
    if (value_exact_small(5ULL) != 12U) {
        std::cerr << "Checkpoint failed: f(5)*m(5)\n";
        return false;
    }
    if (value_exact_small(10ULL) != 90U) {
        std::cerr << "Checkpoint failed: f(10)*m(10)\n";
        return false;
    }

    u128 sum_100 = 0U;
    for (u64 n = 1ULL; n <= 100ULL; ++n) {
        sum_100 += value_exact_small(n);
    }
    if (sum_100 != 1683550844462ULL) {
        std::cerr << "Checkpoint failed: sum up to 100\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << sum_mod(kLimit) << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 982451653
    LIMIT = 100_000_000_000_000

    def isqrt(n):
        return math.isqrt(n)

    def max_m_for_n(n):
        m = (isqrt(8*n + 9) - 3) // 2
        while (m+1)*(m+4)//2 <= n:
            m += 1
        while m > 0 and m*(m+3)//2 > n:
            m -= 1
        return m

    m_max = max_m_for_n(LIMIT)
    inv_limit = m_max + 3

    inv = [0] * (inv_limit + 1)
    inv[1] = 1
    for i in range(2, inv_limit + 1):
        inv[i] = (MOD - MOD // i * inv[MOD % i] % MOD) % MOD

    inv2 = (MOD + 1) // 2
    answer = 1  # n=1
    fact = 1
    harmonic = 0
    upto = 1

    for m in range(1, m_max + 1):
        target = m + 2
        while upto < target:
            upto += 1
            fact = fact * upto % MOD
            if upto >= 2:
                harmonic = (harmonic + inv[upto]) % MOD

        start = m * (m + 3) // 2
        end_full = (m + 1) * (m + 4) // 2 - 1
        ln = min(LIMIT, end_full) - start + 1

        mm = m % MOD
        base = mm * fact % MOD

        if ln <= m:
            d1 = m + 3 - ln
            d2 = m + 2
            partial_h = 0
            for d in range(d1, d2 + 1):
                partial_h = (partial_h + inv[d]) % MOD
            add = base * partial_h % MOD
        elif ln == m + 1:
            add = base * harmonic % MOD
        else:  # ln == m + 2
            add = base * harmonic % MOD
            fact_m3 = fact * (m + 3) % MOD
            extra = mm * fact_m3 % MOD * inv2 % MOD * inv[m + 2] % MOD
            add = (add + extra) % MOD

        answer = (answer + add) % MOD

    return str(answer)

if __name__ == '__main__':
    print(solve())

Java

public class Euler374 {
    static final long kMod = 982451653L;

    static long isqrt(long n) {
        long approx = (long) Math.sqrt(n);
        while ((approx + 1) <= n / (approx + 1))
            approx++;
        while (approx > n / approx)
            approx--;
        return approx;
    }

    static long maxMForN(long n) {
        long m = (isqrt(8 * n + 9) - 3) / 2;
        while ((m + 1) * (m + 4) / 2 <= n)
            m++;
        while (m > 0 && m * (m + 3) / 2 > n)
            m--;
        return m;
    }

    static String solve() {
        long n = 100000000000000L;
        if (n == 0)
            return "0";

        long answer = 1;

        long mMax = maxMForN(n);
        int invLimit = (int) (mMax + 3);

        int[] inv = new int[invLimit + 1];
        inv[1] = 1;
        for (int i = 2; i <= invLimit; i++) {
            inv[i] = (int) ((kMod - (kMod / i) * inv[(int) (kMod % i)] % kMod) % kMod);
        }

        long inv2 = (kMod + 1) / 2;

        long fact = 1;
        long harmonic = 0;
        int upto = 1;

        for (long m = 1; m <= mMax; m++) {
            long target = m + 2;
            while (upto < target) {
                upto++;
                fact = (fact * upto) % kMod;
                if (upto >= 2) {
                    harmonic = (harmonic + inv[upto]) % kMod;
                }
            }

            long start = m * (m + 3) / 2;
            long endFull = (m + 1) * (m + 4) / 2 - 1;
            long len = Math.min(n, endFull) - start + 1;

            long mm = m % kMod;
            long base = (mm * fact) % kMod;

            long add = 0;
            if (len <= m) {
                int d1 = (int) (m + 3 - len);
                int d2 = (int) (m + 2);
                long partialH = 0;
                for (int d = d1; d <= d2; d++) {
                    partialH += inv[d];
                }
                partialH %= kMod;
                add = (base * partialH) % kMod;
            } else if (len == m + 1) {
                add = (base * harmonic) % kMod;
            } else {
                add = (base * harmonic) % kMod;
                long factM3 = (fact * (m + 3)) % kMod;
                long extra = (mm * factM3) % kMod;
                extra = (extra * inv2) % kMod;
                extra = (extra * inv[(int) (m + 2)]) % kMod;
                add = (add + extra) % kMod;
            }

            answer = (answer + add) % kMod;
        }

        return Long.toString(answer);
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}