Problem 154: Exploring Pascal's Pyramid

Project Euler Problem 154 Solution

EulerSolve provides an optimized solution for Project Euler Problem 154, Exploring Pascal's Pyramid, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Problem 154 asks for the number of coefficients in the expansion of $(x+y+z)^{200000}$ that are divisible by $10^{12}$. Every coefficient on the layer $a+b+c=200000$ of Pascal's pyramid has the form $$\binom{200000}{a,b,c}=\frac{200000!}{a!b!c!},\qquad a,b,c\ge 0,\quad a+b+c=200000.$$ Because $10^{12}=2^{12}5^{12}$, a coefficient qualifies exactly when it contains at least twelve factors of 2 and at least twelve factors of 5. The implementations never expand the polynomial; they reduce each candidate to two arithmetic tests on precomputed factorial valuations. Mathematical Approach The problem is really about the prime exponents inside multinomial coefficients. Once those exponents are written in a reusable form, the huge trinomial layer becomes a structured counting problem. Prime valuations are the real divisibility test For any prime $p$, let $v_p(m)$ denote the exponent of $p$ in the prime factorization of $m$. Then for a trinomial coefficient we have $$v_p\!\left(\binom{n}{a,b,c}\right)=v_p(n!)-v_p(a!)-v_p(b!)-v_p(c!),\qquad a+b+c=n.$$ With $n=200000$, divisibility by $10^{12}$ is therefore equivalent to the pair of inequalities $$v_2\!\left(\binom{n}{a,b,c}\right)\ge 12,\qquad v_5\!\left(\binom{n}{a,b,c}\right)\ge 12.$$ No other prime matters....

Detailed mathematical approach

Problem Summary

Balls are stacked into a triangular pyramid in which every ball rests on exactly three balls of the level below. A path starts at the apex and moves down one level at a time, always to one of the three balls directly below, and every ball is labelled with the number of paths that reach it: the apex with $1$, every other ball with the sum of the up to three labels directly above it. The labels form Pascal's pyramid, and those on level $n$ are the coefficients of the expansion of $(x+y+z)^n$. We must count how many coefficients of $(x+y+z)^{200000}$ are multiples of $10^{12}$. ^[1] Every term $x^ay^bz^c$ with $a+b+c=n$ has its own coefficient, so equal values in different places are counted separately. The statement gives no smaller example with its answer.

The whole solution in one sentence: the coefficient of $x^ay^bz^c$ is $n!/(a!\,b!\,c!)$, it is a multiple of $10^{12}$ exactly when it contains each of the prime factors $2$ and $5$ at least twelve times, both counts come from two tables of the exponents of $2$ and $5$ in $m!$, and one pass over the triples $a\le b\le c$, each counted once per rearrangement, adds up the qualifying coefficients.

Mathematical Approach

The coefficients themselves are never computed: the largest ones on level $200000$, at $x^{66666}y^{66667}z^{66667}$ and its rearrangements, have $95419$ decimal digits. Divisibility by $10^{12}$ depends only on the numbers of prime factors $2$ and $5$ of a coefficient, which are small integers.

1. Paths and trinomial coefficients

Give every ball on level $n$ the coordinates $(a,b,c)$, where $a$, $b$ and $c$ count the steps of a path in each of the three downward directions, so that $a+b+c=n$. A path to $(a,b,c)$ is a word of $n$ letters with $a$ letters $x$, $b$ letters $y$ and $c$ letters $z$. If all $n$ letters were different, there would be $n!$ words; making the $a$ copies of $x$ indistinguishable merges them into groups of $a!$, and likewise for $y$ and $z$. So the number of paths is the trinomial coefficient ^[2] ^[3]

$$\binom{n}{a,b,c}=\frac{n!}{a!\,b!\,c!}.$$

Expanding $(x+y+z)^n$ means choosing one of $x$, $y$, $z$ from each of the $n$ factors, so the coefficient of $x^ay^bz^c$ counts the same words; this is the multinomial theorem. The last step of a path to $(a,b,c)$ comes from $(a-1,b,c)$, $(a,b-1,c)$ or $(a,b,c-1)$, whichever of them exist, which is the addition rule of the statement.

Level $2$ holds the rows $1$; $2,2$; $1,2,1$, the six coefficients of $x^2+y^2+z^2+2xy+2xz+2yz$. Level $3$ holds the rows $1$; $3,3$; $3,6,3$; $1,3,3,1$, where $6=3!/(1!\,1!\,1!)$ belongs to $xyz$. ^[4] Level $n$ has $\binom{n+2}{2}=(n+1)(n+2)/2$ places, one for each solution of $a+b+c=n$; on level $200000$ these are $20000300001$ coefficients. ^[2]

2. Divisibility counts prime factors

For a prime $p$ and a positive integer $m$, let $v_p(m)$ be the exponent of $p$ in $m$, the number of times $p$ divides $m$; for example $40=2^3\cdot5$ gives $v_2(40)=3$ and $v_5(40)=1$. Because prime factorisations are unique, exponents add under multiplication, $v_p(rs)=v_p(r)+v_p(s)$, and subtract under exact division. ^[5] For the same reason a number is a multiple of $10^{12}=2^{12}\cdot5^{12}$ exactly when its exponents of $2$ and $5$ are both at least $12$; no other prime matters. For the coefficient $M=\binom{n}{a,b,c}$ we have $M\cdot a!\,b!\,c!=n!$, and hence

$$\begin{aligned}v_p(M)&=v_p(n!)-v_p(a!)\\&\quad-v_p(b!)-v_p(c!).\end{aligned}$$

3. Factors of $p$ in a factorial

Among $1,2,\dots,m$ there are $\lfloor m/p\rfloor$ multiples of $p$, each contributing at least one factor $p$ to $m!$. Of these, $\lfloor m/p^2\rfloor$ are multiples of $p^2$ and contribute a second factor, and so on: a number with exponent $e$ is counted in exactly the first $e$ terms of the sum ^[6]

$$\begin{aligned}v_p(m!)&=\Bigl\lfloor\frac mp\Bigr\rfloor+\Bigl\lfloor\frac m{p^2}\Bigr\rfloor\\&\quad+\Bigl\lfloor\frac m{p^3}\Bigr\rfloor+\cdots,\end{aligned}$$

which is Legendre's formula. For $m=10$ it gives $v_2(10!)=5+2+1=8$ and $v_5(10!)=2$, and indeed $10!=3628800$${}=2^8\cdot3^4\cdot5^2\cdot7$.

Write $m=\sum_i d_ip^i$ in base $p$, with digit sum $s_p(m)=\sum_i d_i$. The quotient $\lfloor m/p^k\rfloor$ is the number with the base-$p$ digits $d_k,d_{k+1},\dots$, so the digit $d_i$ contributes $d_i\,p^{i-k}$ to the $k$-th term for $k=1,\dots,i$, in total $d_i(p^i-1)/(p-1)$ by the geometric sum. Adding over all digits gives the digit form

$$v_p(m!)=\frac{m-s_p(m)}{p-1}.$$

Again for $m=10=1010_2=20_5$: $(10-2)/1=8$ and $(10-2)/4=2$.

4. One coefficient worked out

On level $10$, the coefficient of $x^3y^3z^4$ is $10!/(3!\,3!\,4!)=4200$. With $v_2(3!)=1$, $v_2(4!)=3$ and $v_5(3!)=v_5(4!)=0$, Section 2 gives $v_2(4200)=8-1-1-3=3$ and $v_5(4200)=2-0-0-0=2$, in agreement with $4200=2^3\cdot3\cdot5^2\cdot7$. So $4200$ is a multiple of $100$, but not of $1000$.

5. Exponents are carries

Putting the digit form into the formula of Section 2, the terms $n$, $a$, $b$ and $c$ cancel because $a+b+c=n$, and ^[7]

$$\begin{aligned}v_p(M)&=\frac{1}{p-1}\bigl(s_p(a)+s_p(b)\\&\qquad+s_p(c)-s_p(n)\bigr).\end{aligned}$$

The formula has a picture: add $a$, $b$ and $c$ in base $p$ column by column, as at school. In column $i$ the three digits and the carry $k_{i-1}$ from the previous column give the digit $n_i$ of the sum and a new carry $k_i$:

$$\begin{aligned}a_i+b_i+c_i+k_{i-1}&=n_i+p\,k_i,\\k_{-1}&=0.\end{aligned}$$

Beyond the last digit of $n$ nothing is left to carry, so the final carry is $0$. Adding the equations of all columns, every carry appears once on the left and once, multiplied by $p$, on the right, so $s_p(a)+s_p(b)+s_p(c)+K$${}=s_p(n)+pK$, where $K$ is the sum of all carries. Hence $v_p(M)=K$: a coefficient contains $p$ exactly as many times as the carries add up to. For two summands this is Kummer's theorem. A carry is never larger than $2$, because by induction a column sum is at most $3(p-1)+2\lt3p$.

For $3+3+4$ in base $5$, the units column gives $3+3+4=10=0+5\cdot2$, so the digit is $0$ and the carry $2$, and the next column gives $0+0+0+2=2$. The sum is $20_5=10$, and $K=2=v_5(4200)$. In base $2$, with $3=11_2$ and $4=100_2$, the four columns give $1+1+0=2$, $1+1+0+1=3$, $0+0+1+1=2$ and $0+0+0+1=1$, so the digits from the right are $0,1,0,1$, the sum is $1010_2=10$, and the carries $1,1,1$ give $K=3=v_2(4200)$.

6. Why small levels contain no multiple of $10^{12}$

If $n\lt5^L$, then $n$, $a$, $b$ and $c$ have at most $L$ digits in base $5$, only the columns $0,\dots,L-2$ can pass on a carry, and each carry is at most $2$, so $v_5(M)\le2(L-1)$. For $n\lt5^6=15625$ this allows at most $10$ factors $5$. For $15625\le n\lt31250$ the number $n$ has seven digits with leading digit $1$, the top column reads $a_6+b_6+c_6+k_5=1$, so $k_5\le1$ and $v_5(M)\le2\cdot5+1=11$. Hence no level below $31250$ contains a multiple of $5^{12}$, let alone of $10^{12}$.

Level $31250=2000000_5$ is the first that does. There $v_5(M)=12$ requires all six carries to be $2$: the units digits of $a$, $b$ and $c$ must add up to $10$, the digits of each of the next five columns to $8$, and the digits of the top column must be $0$. With digits from $0$ to $4$, the sum $10$ comes from the $6$ rearrangements of $(2,4,4)$ and $(3,3,4)$, and the sum $8$ from the $15$ rearrangements of $(0,4,4)$, $(1,3,4)$, $(2,2,4)$ and $(2,3,3)$. The columns can be chosen independently, and each choice gives one triple $(a,b,c)$, so $6\cdot15^5=4556250$ coefficients on this level contain $5^{12}$; $4017318$ of them also contain $2^{12}$. On level $200000=22400000_5$ there are seven carries with a total of at most $14$, so a multiple of $5^{12}$ may miss at most two units of carry. Accordingly, the test for $5$ rejects more than $97$ percent of the coefficients there, and the test for $2$ only about $1.3$ percent.

7. Two tables and two inequalities

The programs tabulate $V_2(m)=v_2(m!)$ and $V_5(m)=v_5(m!)$ for $0\le m\le n$ with the running sums

$$\begin{aligned}V_p(0)&=0,\\V_p(m)&=V_p(m-1)+v_p(m),\end{aligned}$$

which hold because $m!=(m-1)!\cdot m$, and they find $v_p(m)$ by dividing $m$ by $p$ as long as the division is exact. The binary digits of $200000$ add up to $6$ and its base-$5$ digits to $8$, so the digit form gives $V_2(200000)=199994$ and $V_5(200000)=49998$. With the thresholds $T_2=V_2(n)-12=199982$ and $T_5=V_5(n)-12=49986$, Section 2 turns the question for one coefficient into two inequalities between table entries:

$$\begin{aligned}V_2(a)+V_2(b)+V_2(c)&\le T_2,\\V_5(a)+V_5(b)+V_5(c)&\le T_5.\end{aligned}$$

8. Sorting the triple

The coefficient does not change when $a$, $b$ and $c$ are permuted; this is the three-way symmetry of every level of the pyramid. ^[2] So it suffices to test the sorted triples $a\le b\le c$ and to count each passing triple once for every distinct rearrangement: $6$ times if its entries are distinct, $3$ times if exactly two of them are equal, and once if all three are equal. With $c=n-a-b$, the condition $b\le c$ means $b\le\lfloor(n-a)/2\rfloor$, and some $b\ge a$ satisfies it exactly when $a\le(n-a)/2$, that is, when $a\le\lfloor n/3\rfloor$.

The sorted triples are the partitions of $n$ into at most three parts, and there are $\operatorname{round}((n+3)^2/12)$ of them: $14$ for $n=10$ and $3333433334$ for $n=200000$. ^[8] Since $3$ does not divide $200000$, no triple on that level has three equal entries. A triple with two equal entries is the sorted form of $(r,r,200000-2r)$ for exactly one $r$ with $0\le r\le100000$, so there are $100001$ of them, the other $3333333333$ have distinct entries, and indeed $6\cdot3333333333+3\cdot100001$${}=20000300001$, the number of places from Section 1.

9. A complete small example

Which coefficients of $(x+y+z)^{10}$ are multiples of $100=2^2\cdot5^2$? Here $V_2(10)=8$ and $V_5(10)=2$, so with $2$ in place of $12$ the thresholds are $T_2=6$ and $T_5=0$. Since $V_5(m)=0$ exactly for $m\le4$, the condition on $5$ forces $a,b,c\le4$, and the only sorted triples with sum $10$ and entries at most $4$ are $(2,4,4)$ and $(3,3,4)$. With $V_2(2)=V_2(3)=1$ and $V_2(4)=3$, the first triple gives $1+3+3=7\gt6$ and fails, in agreement with $10!/(2!\,4!\,4!)=3150$${}=2\cdot3^2\cdot5^2\cdot7$; the second gives $1+1+3=5\le6$ and passes. The triple $(3,3,4)$ has two equal entries, so exactly $3$ of the $66$ coefficients are multiples of $100$: the number $4200$ at $x^3y^3z^4$, $x^3y^4z^3$ and $x^4y^3z^3$.

For the divisor $10$ ($T_2=7$, $T_5=1$) only $(0,0,10)$, $(0,2,8)$ and $(0,5,5)$ fail, with the coefficients $1$, $45$, which is odd, and $252$, which has no factor $5$, at $3+6+3=12$ places, so $54$ of the $66$ coefficients are multiples of $10$. A direct computation of all $66$ coefficients confirms both counts.

How the Code Works

The C++ program ^[9] starts in main, which reads --n= (default $200000$) and --skip-checkpoints, runs run_checkpoints unless told to skip it, and prints the count. build_factorial_prime_exponents fills the vectors v2 and v5 with the tables $V_2$ and $V_5$ of Section 7; it strips the factors $2$ of each $i$ with the bit test (x & 1) == 0 and a shift, and the factors $5$ with x % 5. count_coefficients_fast computes need2 $=T_2$ and need5 $=T_5$ and runs the loops of Section 8: $a$ from $0$ to a_max = n / 3, $b$ from $a$ to b_max = (n - a) / 2, and c = n - a - b. A triple that violates one of the two inequalities is skipped with continue. The factor $2$ is tested first, so the test for $5$, which rejects almost every triple (Section 6), runs for nearly all of them; the order does not affect the count. A passing triple adds $1$ if a == b && b == c, $3$ if a == b || b == c, which for a sorted triple means exactly two equal entries, and $6$ otherwise. The values of $a$ are taken in blocks of $64$ (kChunk), whose partial counts are added at the end.

The counts are std::uint64_t, because the level has $20000300001$ places, more than a 32-bit counter can hold, while the table entries fit in int: a sum such as $V_2(a)+V_2(b)+V_2(c)$ never exceeds $V_2(n)=199994$, since the exponent of Section 2 is never negative. count_coefficients_bruteforce loops over all ordered pairs $(a,b)$ with $a+b\le n$ and adds $1$ for every passing triple, without sorting or weights. run_checkpoints requires both functions to agree for $n=40$, $60$, $80$ and $120$, and two differently organised evaluations of the fast count for $n=3000$ to agree.

The Python program mirrors this: count_chunk(a_start, a_end) runs the double loop for one block of values of $a$, and solve adds up the counts of consecutive blocks, $256$ values per block in the default setting. Its integers have no fixed size; it accepts --n and has no checkpoints. The Java program does everything in main with $n=200000$ fixed, with a long counter and no checks.

Complexity and Verification

Building the tables takes $O(n)$ time: the division loops run $V_2(n)+V_5(n)$ times in total, which by the digit form of Section 3 is less than $n+n/4$, and exactly $249992$ times for $n=200000$. The memory is $O(n)$, two arrays of $n+1$ integers. The main loop visits each of the $3333433334$ sorted triples of Section 8 once: about $n^2/12$ steps of constant cost, a sixth of the $20000300001$ coefficients, so the time grows like $n^2$. In tests on a desktop computer the C++ and Java programs finished within seconds and the Python program within a minute. Generating, column by column, only the triples whose base-$5$ carries add up to at least $12$ (Section 6) would avoid testing every sorted triple; the published programs keep the simpler scan.

By Section 6, the levels $40$, $60$, $80$, $120$ and $3000$ of the C++ checkpoints contain no multiple of $10^{12}$, so every count compared there is $0$: the checkpoints confirm that neither the fast nor the brute-force count reports a false coefficient, but never test the weights $3$ and $6$. For this article we therefore ran further checks, which are not part of the published solutions:

For every $0\le n\le120$ and every $k$ from $1$ to $6$, the method of Sections 7 and 8 with $10^k$ in place of $10^{12}$ agrees with a test of the exact coefficients; $342$ of these $726$ counts are positive.
For every triple with $0\le n\le90$ and for $p=2,3,5,7$, the exponent of $p$ in the exact coefficient equals the carry total of Section 5.
On level $31250$ the C++ program with --n=31250, the Python program with --n 31250 and a separate count over all ordered triples, with exponents from digit sums, give $4017318$; the separate count also finds exactly the $4556250$ multiples of $5^{12}$ of Section 6. Level $31249$ gives $0$.
On level $200000$ the three programs print the same number, and so does the separate count over all $20000300001$ ordered triples, which uses neither the running-sum tables nor the weights.

Footnotes and References

The references below supply the background facts used above. The carry form for three summands, the bound that rules out all levels below $31250$, the counts of rearrangements and the operation counts are derived explicitly in this article; the formula for the number of sorted triples is quoted from OEIS A001399.

Project Euler 154 — Exploring Pascal's Pyramid. The official statement labels every ball of a triangular pyramid with the number of paths from the apex, each step going to one of the three balls directly below, notes that level n holds the coefficients of (x + y + z)ⁿ, and asks how many coefficients of (x + y + z)²⁰⁰⁰⁰⁰ are multiples of 10¹².
Wikipedia — Pascal's pyramid. Layer n consists of the coefficients of a trinomial raised to the nth power, each being n! divided by the factorials of the three exponents of its term; every number below the apex is the sum of the three adjacent numbers in the layer above, layer n has as many terms as the (n + 1)th triangular number, and every layer has three-way symmetry.
Wikipedia — Multinomial theorem. The coefficients are the multinomial coefficients n!/(k₁! ⋯ kₘ!), which also count the distinct permutations of a multiset; there are C(n + m − 1, m − 1) of them.
OEIS A046816 — Pascal's tetrahedron, or trinomial coefficients. Its data begin with the levels 0 to 3, 1; 1, 1, 1; 1, 2, 2, 1, 2, 1; 1, 3, 3, 3, 6, 3, 1, 3, 3, 1, and its formula section gives the entries as n!/(i! j! k!).
Wikipedia — p-adic valuation. The p-adic valuation of a nonzero integer is the exponent of the highest power of the prime p that divides it; the valuation of a product is the sum of the valuations, and that of a quotient is their difference.
Wikipedia — Legendre's formula. The exponent of p in n! is the sum of ⌊n/pⁱ⌋ over i ≥ 1, which also equals (n − sₚ(n))/(p − 1), where sₚ(n) is the sum of the base-p digits of n.
Wikipedia — Kummer's theorem. The exponent of p in C(n, m) equals the number of carries when m is added to n − m in base p; for multinomial coefficients the article gives the digit-sum form of Section 5.
OEIS A001399 — partitions of n into at most 3 parts. The entry gives the formula a(n) = round((n + 3)²/12), and its terms begin 1, 1, 2, 3, 4, 5, 7, 8, 10, 12, 14.
C++ — Euler154.cpp. The linked immutable C++ revision contains build_factorial_prime_exponents, count_coefficients_fast, count_coefficients_bruteforce and run_checkpoints. It is the implementation record for this article.

Problem 154 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <string>
#include <thread>
#include <vector>

namespace {

using u64 = std::uint64_t;

struct Options {
    int n = 200000;
    int threads = static_cast<int>(std::thread::hardware_concurrency());
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    if (options.threads <= 0) {
        options.threads = 1;
    }
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--n=", options.n) ||
            parse_int_after_prefix(arg, "--threads=", options.threads)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.n >= 1 && options.threads >= 1;
}

void build_factorial_prime_exponents(const int n, std::vector<int>& v2, std::vector<int>& v5) {
    v2.assign(static_cast<std::size_t>(n + 1), 0);
    v5.assign(static_cast<std::size_t>(n + 1), 0);
    for (int i = 1; i <= n; ++i) {
        int x = i;
        int add2 = 0;
        while ((x & 1) == 0) {
            x >>= 1;
            ++add2;
        }

        x = i;
        int add5 = 0;
        while (x % 5 == 0) {
            x /= 5;
            ++add5;
        }

        v2[static_cast<std::size_t>(i)] = v2[static_cast<std::size_t>(i - 1)] + add2;
        v5[static_cast<std::size_t>(i)] = v5[static_cast<std::size_t>(i - 1)] + add5;
    }
}

u64 count_coefficients_fast(const int n, int threads) {
    std::vector<int> v2;
    std::vector<int> v5;
    build_factorial_prime_exponents(n, v2, v5);

    const int need2 = v2[static_cast<std::size_t>(n)] - 12;
    const int need5 = v5[static_cast<std::size_t>(n)] - 12;

    const int a_max = n / 3;
    threads = std::max(1, std::min(threads, a_max + 1));

    std::atomic<int> next_a{0};
    constexpr int kChunk = 64;

    std::vector<u64> partial(static_cast<std::size_t>(threads), 0);
    std::vector<std::thread> pool;
    pool.reserve(static_cast<std::size_t>(threads));

    for (int t = 0; t < threads; ++t) {
        pool.emplace_back([&, t]() {
            u64 local = 0;
            while (true) {
                const int start = next_a.fetch_add(kChunk, std::memory_order_relaxed);
                if (start > a_max) {
                    break;
                }
                const int end = std::min(a_max, start + kChunk - 1);
                for (int a = start; a <= end; ++a) {
                    const int a2 = v2[static_cast<std::size_t>(a)];
                    const int a5 = v5[static_cast<std::size_t>(a)];
                    const int b_max = (n - a) / 2;
                    for (int b = a; b <= b_max; ++b) {
                        const int c = n - a - b;
                        if (a2 + v2[static_cast<std::size_t>(b)] + v2[static_cast<std::size_t>(c)] > need2) {
                            continue;
                        }
                        if (a5 + v5[static_cast<std::size_t>(b)] + v5[static_cast<std::size_t>(c)] > need5) {
                            continue;
                        }

                        if (a == b && b == c) {
                            local += 1;
                        } else if (a == b || b == c) {
                            local += 3;
                        } else {
                            local += 6;
                        }
                    }
                }
            }
            partial[static_cast<std::size_t>(t)] = local;
        });
    }

    for (auto& th : pool) {
        th.join();
    }

    u64 answer = 0;
    for (u64 v : partial) {
        answer += v;
    }
    return answer;
}

u64 count_coefficients_bruteforce(const int n) {
    std::vector<int> v2;
    std::vector<int> v5;
    build_factorial_prime_exponents(n, v2, v5);

    const int need2 = v2[static_cast<std::size_t>(n)] - 12;
    const int need5 = v5[static_cast<std::size_t>(n)] - 12;

    u64 answer = 0;
    for (int a = 0; a <= n; ++a) {
        const int a2 = v2[static_cast<std::size_t>(a)];
        const int a5 = v5[static_cast<std::size_t>(a)];
        for (int b = 0; b <= n - a; ++b) {
            const int c = n - a - b;
            if (a2 + v2[static_cast<std::size_t>(b)] + v2[static_cast<std::size_t>(c)] > need2) {
                continue;
            }
            if (a5 + v5[static_cast<std::size_t>(b)] + v5[static_cast<std::size_t>(c)] > need5) {
                continue;
            }
            ++answer;
        }
    }
    return answer;
}

bool run_checkpoints() {
    for (int n : {40, 60, 80, 120}) {
        const u64 slow = count_coefficients_bruteforce(n);
        const u64 fast = count_coefficients_fast(n, 1);
        if (slow != fast) {
            std::cerr << "Checkpoint failed at n=" << n << ": brute=" << slow
                      << ", fast=" << fast << '\n';
            return false;
        }
    }

    const int thread_check_n = 3000;
    const u64 single = count_coefficients_fast(thread_check_n, 1);
    unsigned hw = std::thread::hardware_concurrency();
    if (hw == 0) {
        hw = 2;
    }
    const int multi_threads = static_cast<int>(std::min<unsigned>(hw, 8));
    const u64 multi = count_coefficients_fast(thread_check_n, multi_threads);
    if (single != multi) {
        std::cerr << "Thread consistency failed at n=" << thread_check_n
                  << ": single=" << single << ", multi=" << multi << '\n';
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << count_coefficients_fast(options.n, options.threads) << '\n';
    return 0;
}

Python

from __future__ import annotations

import argparse
import multiprocessing as mp
import os
from concurrent.futures import ProcessPoolExecutor


G_V2: list[int] = []
G_V5: list[int] = []
G_N = 0
G_NEED2 = 0
G_NEED5 = 0


def build_factorial_prime_exponents(n: int) -> tuple[list[int], list[int]]:
    v2 = [0] * (n + 1)
    v5 = [0] * (n + 1)
    for i in range(1, n + 1):
        x = i
        add2 = 0
        while (x & 1) == 0:
            x >>= 1
            add2 += 1

        x = i
        add5 = 0
        while x % 5 == 0:
            x //= 5
            add5 += 1

        v2[i] = v2[i - 1] + add2
        v5[i] = v5[i - 1] + add5
    return v2, v5


def init_worker(v2: list[int], v5: list[int], n: int, need2: int, need5: int) -> None:
    global G_V2, G_V5, G_N, G_NEED2, G_NEED5
    G_V2 = v2
    G_V5 = v5
    G_N = n
    G_NEED2 = need2
    G_NEED5 = need5


def count_chunk(a_start: int, a_end: int) -> int:
    v2 = G_V2
    v5 = G_V5
    n = G_N
    need2 = G_NEED2
    need5 = G_NEED5

    local = 0
    for a in range(a_start, a_end):
        a2 = v2[a]
        a5 = v5[a]
        b_max = (n - a) // 2
        for b in range(a, b_max + 1):
            c = n - a - b
            if a2 + v2[b] + v2[c] > need2:
                continue
            if a5 + v5[b] + v5[c] > need5:
                continue
            if a == b == c:
                local += 1
            elif a == b or b == c:
                local += 3
            else:
                local += 6
    return local


def solve(n: int = 200000, threads: int | None = None) -> int:
    v2, v5 = build_factorial_prime_exponents(n)
    need2 = v2[n] - 12
    need5 = v5[n] - 12
    a_max = n // 3

    if threads is None:
        threads = min(8, os.cpu_count() or 1)
    threads = max(1, threads)

    if threads == 1:
        init_worker(v2, v5, n, need2, need5)
        return count_chunk(0, a_max + 1)

    ctx = mp.get_context("fork") if "fork" in mp.get_all_start_methods() else mp.get_context()
    chunk_size = 256
    tasks = []
    for start in range(0, a_max + 1, chunk_size):
        end = min(a_max + 1, start + chunk_size)
        tasks.append((start, end))

    total = 0
    with ProcessPoolExecutor(
        max_workers=threads,
        mp_context=ctx,
        initializer=init_worker,
        initargs=(v2, v5, n, need2, need5),
    ) as pool:
        futures = [pool.submit(count_chunk, start, end) for start, end in tasks]
        for fut in futures:
            total += fut.result()
    return total


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--n", type=int, default=200000)
    parser.add_argument("--threads", type=int, default=min(8, os.cpu_count() or 1))
    args = parser.parse_args()
    print(solve(args.n, args.threads))


if __name__ == "__main__":
    main()

Java

public class Euler154 {
    public static void main(String[] args) {
        int n = 200000;
        int[] v2 = new int[n + 1], v5 = new int[n + 1];
        for (int i = 1; i <= n; i++) {
            int x = i, a2 = 0;
            while (x % 2 == 0) {
                x /= 2;
                a2++;
            }
            x = i;
            int a5 = 0;
            while (x % 5 == 0) {
                x /= 5;
                a5++;
            }
            v2[i] = v2[i - 1] + a2;
            v5[i] = v5[i - 1] + a5;
        }
        int need2 = v2[n] - 12, need5 = v5[n] - 12;
        long count = 0;
        int aMax = n / 3;
        for (int a = 0; a <= aMax; a++) {
            int a2v = v2[a], a5v = v5[a];
            int bMax = (n - a) / 2;
            for (int b = a; b <= bMax; b++) {
                int c = n - a - b;
                if (a2v + v2[b] + v2[c] > need2)
                    continue;
                if (a5v + v5[b] + v5[c] > need5)
                    continue;
                if (a == b && b == c)
                    count += 1;
                else if (a == b || b == c)
                    count += 3;
                else
                    count += 6;
            }
        }
        System.out.println(count);
    }
}