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

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.

Step 2: The block parameter $m$

The baseline near-consecutive partition with $m$ parts is

$$\{2,3,\dots,m+1\},$$

whose sum is

$$T_m = 2+3+\cdots+(m+1)=\frac{m(m+3)}{2}.$$

This is the left endpoint of the block where the optimal partition has exactly $m$ parts. Hence the correct block index for a given $n$ is

$$m=\max\left\{r:\frac{r(r+3)}{2}\le n\right\}.$$

The code computes this with the integer-square-root estimate

$$m\approx \frac{\sqrt{8n+9}-3}{2},$$

and then adjusts by at most a couple of integer steps to remove rounding issues.

Step 3: Closed form inside one block

Write

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

Since

$$T_{m+1}=\frac{(m+1)(m+4)}{2}=T_m+m+2,$$

the full $m$-block is the interval

$$T_m \le n \le T_{m+1}-1=\frac{(m+1)(m+4)}{2}-1.$$

Throughout this whole interval the optimal partition size is constant:

$$m(n)=m.$$

Step 4: The optimal partition shapes

If $0\le s\le m$, the optimal partition is obtained from the consecutive set $\{2,3,\dots,m+2\}$ by deleting exactly one value:

$$\{2,3,\dots,m+2\}\setminus\{m+2-s\}.$$

Its sum is

$$ \left(\sum_{j=2}^{m+2} j\right)-(m+2-s) = \frac{m(m+3)}{2}+s = n, $$

and its product is

$$f(n)=\frac{(m+2)!}{m+2-s}.$$

Multiplying by the part count gives

$$f(n)m(n)=m\frac{(m+2)!}{m+2-s},\qquad 0\le s\le m.$$

If $s=m+1$, the “deleted value” would have to be $1$, so the shape changes to

$$\{3,4,\dots,m+1,m+3\}.$$

Its product is

$$f(n)=\frac{(m+3)!}{2(m+2)},$$

hence

$$f(n)m(n)=m\frac{(m+3)!}{2(m+2)}.$$

This is exactly the piecewise formula implemented by all three solution files. The code keeps the middle case $s=m$ as a separate branch, but mathematically it is simply the denominator-$2$ instance of the first formula.

Worked Example: $n=10$

We have

$$T_3=\frac{3\cdot 6}{2}=9,\qquad T_4=\frac{4\cdot 7}{2}=14,$$

so $m=3$ and $s=10-9=1$. Therefore the optimal partition is

$$\{2,3,4,5\}\setminus\{4\}=\{2,3,5\}.$$

Thus

$$f(10)=2\cdot 3\cdot 5=30,\qquad m(10)=3,\qquad f(10)m(10)=90,$$

which matches the checkpoint in the C++ verifier.

Step 5: Summing one whole block

For a complete block, first sum the $m+1$ terms with $0\le s\le m$:

$$ \sum_{s=0}^{m} m\frac{(m+2)!}{m+2-s} = m(m+2)!\sum_{d=2}^{m+2}\frac{1}{d}, $$

where we reindexed with $d=m+2-s$. The final point of the block contributes

$$m\frac{(m+3)!}{2(m+2)}.$$

So the complete $m$-block contribution is

$$ B_m = m(m+2)!\sum_{d=2}^{m+2}\frac{1}{d} + m\frac{(m+3)!}{2(m+2)}. $$

This is why the implementation only needs a running factorial and a running harmonic-style sum of modular inverses.

Step 6: Modular inverses

Because $P$ is prime and every denominator satisfies $2\le d\le m_{\max}+3 \lt P$, all required inverses exist modulo $P$. The code precomputes them in linear time via

$$\mathrm{inv}[1]=1,\qquad \mathrm{inv}[i]=-\left\lfloor\frac{P}{i}\right\rfloor\mathrm{inv}[P\bmod i]\pmod{P}.$$

This follows from writing $P=qi+r$, so $r\equiv -qi\pmod P$, then multiplying by $r^{-1}i^{-1}$. Once the inverse table is built, every full block is processed in $O(1)$ modular arithmetic.

How the Code Works

The C++, Python, and Java implementations are structurally identical. They first compute $m_{\max}=\max\{m:T_m\le N\}$, then precompute modular inverses up to $m_{\max}+3$. During the main loop they maintain

$$\texttt{fact}=(m+2)!\pmod P,\qquad \texttt{harmonic}=\sum_{d=2}^{m+2}\frac{1}{d}\pmod P.$$

For each $m$, the code knows the block start $T_m$, the block end $T_{m+1}-1$, and how much of that block is still inside the target range. All fully covered blocks use the closed form above; only the last block may be truncated, and the code handles that by summing the needed inverse segment explicitly.

The C++ file also checks the derivation against small exact values:

$$f(5)m(5)=12,\qquad f(10)m(10)=90,\qquad \sum_{n=1}^{100} f(n)m(n)=1683550844462.$$

Complexity Analysis

Since $T_m\sim m^2/2$, the maximal block index satisfies $m_{\max}=\Theta(\sqrt{N})$. Precomputing inverses, advancing the running factorial, and iterating over all blocks therefore costs $O(\sqrt{N})$ time. The inverse table uses $O(\sqrt{N})$ memory. Only the final block can be partial, so there is no hidden extra logarithmic factor.

Footnotes and References

Problem page: https://projecteuler.net/problem=374
Integer partition: Wikipedia — Integer partition
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
Competitive programming reference for the inverse recurrence: cp-algorithms — Modular inverse

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());
    }
}