Problem 378: Triangle Triples

Project Euler Problem 378 Solution

EulerSolve provides an optimized solution for Project Euler Problem 378, Triangle Triples, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let $$T_n=\frac{n(n+1)}{2},\qquad D_n=d(T_n),$$ where $d(m)$ is the number of positive divisors of $m$. For a given $N$, we must count all triples of indices $$1\le i \lt j \lt k\le N$$ such that $$D_i \gt D_j \gt D_k.$$ In other words, after converting each triangular number $T_n$ into its divisor count $D_n$, the task is to count strictly decreasing subsequences of length three. The implementation targets $N=60{,}000{,}000$ and returns the result modulo $10^{18}$, so both the divisor-count computation and the triple counting must be highly optimized. Mathematical Approach Step 1: Compute $d(T_n)$ from coprime factors Consecutive integers are coprime: $$\gcd(n,n+1)=1.$$ Also, exactly one of $n$ and $n+1$ is even. Therefore we can remove the factor $2$ before applying the divisor function: $$T_n=\begin{cases} \frac{n}{2}(n+1), & n \text{ even},\\[2mm] n\frac{n+1}{2}, & n \text{ odd}. \end{cases}$$ In both cases the two factors are coprime, so the divisor function is multiplicative: $$d(ab)=d(a)d(b)\qquad \text{when } \gcd(a,b)=1.$$ Hence $$D_n=d(T_n)=\begin{cases} d\!\left(\frac{n}{2}\right)d(n+1), & n \text{ even},\\[2mm] d(n)d\!\left(\frac{n+1}{2}\right), & n \text{ odd}. \end{cases}$$ This identity is the core simplification used by every implementation....

Detailed mathematical approach

Problem Summary

Let $T(n)=\frac{n(n+1)}{2}$ be the $n$-th triangle number and let $dT(n)$ be the number of divisors of $T(n)$; for example $T(7)=28$ has the divisors $1,2,4,7,14,28$, so $dT(7)=6$. Let $Tr(n)$ count the triples $(i,j,k)$ with

$$\begin{aligned} &1\le i\lt j\lt k\le n,\\ &dT(i)\gt dT(j)\gt dT(k). \end{aligned}$$

The statement gives $Tr(20)=14$, $Tr(100)=5772$ and $Tr(1000)=11174776$, and asks for $Tr(N)$ with $N=60\,000\,000$, of which only the last $18$ decimal digits are wanted. ^[1]

The whole solution in one sentence: write down $dT(1),\dots,dT(N)$ with one linear sieve, then sweep the list from left to right, keeping in two small tables how many earlier entries and how many earlier decreasing pairs sit above each value, so that every position can report the triples ending at it.

Mathematical Approach

1. The list we have to look at

Only the divisor counts matter, so put $D_n=dT(n)$. The first twelve triangle numbers and their divisor counts are

$$\begin{aligned} T(1),\dots,T(12)&=1,3,6,10,15,21,\\ &\quad\ 28,36,45,55,66,78,\\ D_1,\dots,D_{12}&=1,2,4,4,4,4,\\ &\quad\ 6,9,6,4,8,8, \end{aligned}$$

and continuing to $n=20$ gives

$$D_{13},\dots,D_{20}=4,8,16,8,6,6,8,16,$$

which is the start of OEIS A063440. ^[2] The list is not monotone and has many repetitions; both facts shape everything below. Note that the inequalities in the definition are strict, so two indices with equal divisor counts never appear together in a triple.

2. Splitting $T(n)$ into two coprime factors

Factorising $T(n)$ directly for $6\cdot10^{7}$ values is out of the question, but we never have to. First, consecutive integers are coprime: any common divisor of $n$ and $n+1$ divides their difference $1$, so $\gcd(n,n+1)=1$. Second, exactly one of $n$ and $n+1$ is even, and halving it keeps the two factors coprime. If $n=2k$ then $T(n)=k(2k+1)$ and a common divisor of the two factors divides $2k+1-2k=1$; if $n=2k+1$ then $T(n)=(2k+1)(k+1)$ and it divides $2(k+1)-(2k+1)=1$. The standard advice to halve whichever of the two factors is even before multiplying is exactly this split; what the argument above adds is that the halves are then coprime. ^[3]

3. The divisor count of a product of coprime numbers

If $m=\prod_i p_i^{a_i}$ is the prime factorisation of $m$, a divisor is obtained by choosing each exponent independently in $\{0,1,\dots,a_i\}$, so

$$d(m)=\prod_i (a_i+1).$$

For instance $36=2^2 3^2$ gives $d(36)=3\cdot3=9$. Coprime numbers share no prime, so their exponent lists simply concatenate and $d(ab)=d(a)\,d(b)$ whenever $\gcd(a,b)=1$: the divisor function is multiplicative. ^[4] Combining this with Section 2,

$$\begin{aligned} D_{2k}&=d(k)\,d(2k+1),\\ D_{2k+1}&=d(2k+1)\,d(k+1). \end{aligned}$$

Both checks come out right: $D_7=d(7)\,d(4)=2\cdot3=6$ and $D_8=d(4)\,d(9)=3\cdot3=9$. Every $D_n$ is now a product of two values of $d$ at arguments at most $n+1$, so one table of $d$ up to $N+1$ answers everything.

4. All divisor counts at once with a linear sieve

A linear sieve walks $i=2,3,\dots$, recording for each $i$ its smallest prime factor $\operatorname{lp}(i)$ and the exponent $e(i)$ of that prime. ^[5] An unmarked $i$ is prime, so $\operatorname{lp}(i)=i$, $e(i)=1$, $d(i)=2$. Then, for each known prime $p\le\operatorname{lp}(i)$, it marks $x=ip$ and derives $d(x)$ from $d(i)$ in two cases.

If $p\lt\operatorname{lp}(i)$ then $p\nmid i$, so $x$ adds a brand-new prime with exponent $1$ and $d(x)=2\,d(i)$.
If $p=\operatorname{lp}(i)$, write $i=m\,p^{e}$ with $e=e(i)$ and $p\nmid m$. Then $x=m\,p^{e+1}$, and from $d(i)=d(m)(e+1)$ and $d(x)=d(m)(e+2)$ we get $d(x)=\frac{e+2}{e+1}\,d(i)$. The sieve stops the inner loop here, and that is what marks every composite exactly once: if $q$ is the smallest prime factor of a composite $x$ and $i=x/q$, every prime factor of $i$ is at least $q$, so the loop for $i$ reaches $p=q$ and produces $x$, while no other pair $(i,p)$ can. The pass is therefore linear.

In the second case $e+1$ divides $d(i)$, so $d(i)/(e+1)$ may be computed by integer division before multiplying by $e+2$; no rounding can occur. Both cases on small numbers: $i=9$, $p=2$ give $d(18)=2\,d(9)=6$, and $i=6$, $p=2$ give $d(12)=4\cdot3/2=6$. Three measured facts bound the storage: up to $N+1$ there are $3\,562\,115$ primes, the largest exponent of a smallest prime factor is $25$ (since $2^{25}\lt N\lt 2^{26}$), and the largest divisor count is $d(43\,243\,200)=672$, where $43\,243\,200=2^6 3^3 5^2\cdot7\cdot11\cdot13$. One byte per exponent and $16$ bits per divisor count are ample.

5. What is actually being counted

With $D$ available the question is purely combinatorial: how many strictly decreasing subsequences of length three does $D_1,\dots,D_N$ contain? Testing all triples is hopeless, since

$$\binom{N}{3}=35\,999\,998\,200\,000\,020\,000\,000$$

for $N=6\cdot10^{7}$, and even the $\binom{N}{2}\approx1.8\cdot10^{15}$ pairs are too many to list. The way out is never to look at a triple as a whole.

6. Counting each triple once, at its last index

Fix $k$ and look only at what lies to its left. Define, for every value $v$,

$$C_v^{(k)}=\#\{\,i\lt k\ :\ D_i=v\,\}$$

$$\begin{aligned} P_v^{(k)}=\#\{(i,j):{}&i\lt j\lt k,\\ &D_i\gt D_j=v\} \end{aligned}$$

so $C$ counts earlier entries by value and $P$ counts earlier decreasing pairs by the value of their second member. A triple $(i,j,k)$ is a decreasing pair $(i,j)$ followed by a smaller $D_k$, and $k$ is determined by the triple, so attributing each triple to its largest index counts it exactly once. Hence

$$Tr(N)=\sum_{k=1}^{N}\ \sum_{v\gt D_k}P_v^{(k)}.$$

The tables are maintained by the same idea one level lower. Going from $k$ to $k+1$ adds index $k$ on the left: one new entry at value $D_k$, and one new decreasing pair ending at $k$ for each earlier entry above $D_k$. So

$$\begin{aligned} P_{D_k}&\leftarrow P_{D_k}+\sum_{v\gt D_k}C_v^{(k)},\\ C_{D_k}&\leftarrow C_{D_k}+1, \end{aligned}$$

By induction on $k$ both tables are correct at every step: they are empty at $k=1$, and each step adds precisely what index $k$ contributes. Two remarks. The inner sums run over values strictly above $D_k$, which is where strictness is enforced, so equal divisor counts are skipped. And the running total of $\sum_{v\gt D_k}C_v^{(k)}$ is the number of decreasing pairs, the inversion count of the list: ^[6] this is that classical count with one layer on top, and $m-1$ layers would count decreasing subsequences of length $m$.

7. Worked example: $Tr(20)=14$

Run the sweep on $D_1,\dots,D_{20}$. Most positions contribute nothing, because a triple needs two larger values before it. The positions that do contribute are

$k=10$ with $D_{10}=4$: one triple;
$k=13$ with $D_{13}=4$: three triples;
$k=17$ with $D_{17}=6$: five triples;
$k=18$ with $D_{18}=6$: five triples,

and $1+3+5+5=14$. Take $k=17$ in detail. Just before it, the value table is $C_1=1$, $C_2=1$, $C_4=6$, $C_6=2$, $C_8=4$, $C_9=1$, $C_{16}=1$ (sixteen entries in all), and the pair table is $P_4=8$, $P_6=1$, $P_8=5$, all other $P_v=0$. Since $D_{17}=6$, the values above it are $8$, $9$ and $16$, so position $17$ contributes $P_8+P_9+P_{16}=5+0+0=5$ triples, namely

$$\begin{aligned} &(8,11,17),\ (8,12,17),\\ &(8,14,17),\ (8,16,17),\\ &(15,16,17), \end{aligned}$$

with divisor counts $(9,8,6)$ for the first four and $(16,8,6)$ for the last. The update then adds $C_8+C_9+C_{16}=4+1+1=6$ to $P_6$ and $1$ to $C_6$. After all twenty positions the pair table holds $8+13+7=28$ decreasing pairs and the triple total is $14$, as required.

8. The value axis is tiny

The tables are indexed by divisor counts, not by positions, and that axis is short. Measured over $n\le N$, the largest value is

$$\begin{aligned} D_{59181759}&=d(59181759)\\ &\quad{}\cdot d(29590880)\\ &=120\cdot192=23040, \end{aligned}$$

and the maxima for $N=20,10^2,10^3,\dots,10^7,6\cdot10^{7}$ are $16,40,144,480,1344,3072,7680,23040$: tiny beside $N$, though these are measurements over finite ranges and not an asymptotic law. Moreover only $362$ distinct values of $D_n$ occur at all for $n\le N$. Compressing the axis to those $362$ slots is possible, but direct indexing is already compact: a value $v$ at position $v+1$ needs only $23042$ slots, about $184$ kB per table, which stays in cache. Since every $D_n$ is at least $1$, indexing by $v$ itself would serve just as well; the shift is the code's convention and leaves index $1$ unused.

9. Suffix sums with a Fenwick tree

Each step of Section 6 needs two suffix sums over the value axis and two point updates. A Fenwick tree, or binary indexed tree, does exactly this. Write $b$ for the lowest set bit of an index $t$, computed as t & -t; cell $t$ stores the sum over the block $[t-b+1,t]$. A prefix query subtracts that bit repeatedly, so the prefix through $13$ reads cells $13$, $12$, $8$, covering $[13,13]$, $[9,12]$, $[1,8]$; an update adds it repeatedly, so on an axis of length $16$ a change at $11$ touches cells $11$, $12$, $16$. Either walk takes one step at most per binary digit of the axis length $M$, hence $O(\log M)$ steps. ^[7] It answers prefix queries, so a suffix sum is a difference of two prefixes,

$$\sum_{t=r+1}^{M}a_t=\Big(\sum_{t=1}^{M}a_t\Big)-\sum_{t=1}^{r}a_t,$$

which is why the code asks for the range $[r+1,M]$ with $r=D_k+1$. Two trees are kept, one for $C$ and one for $P$; nothing else about them is ever needed.

10. Why $64$ bits suffice, and why only $18$ digits are asked

The pair tree holds exact counts, whose total can never exceed $\binom{N}{2}=1\,799\,999\,970\,000\,000$; at the target it is $790\,912\,595\,361\,069$, so every cell stays below $1.8\cdot10^{15}$. The triple total is a different matter, and one measurement settles it. Lengthening the list can only add triples, so $Tr$ is non-decreasing, and $Tr(10^7)=18\,830\,128\,211\,315\,630\,867$ already exceeds $2^{64}-1=18\,446\,744\,073\,709\,551\,615$; the exact value at $N=6\cdot10^{7}$ therefore cannot be accumulated in $64$ bits at all. (The measured density $Tr(n)/\binom{n}{3}$ climbs from $0.0673$ at $n=10^3$ to $0.1130$ at $n=10^7$.) Asking for the last $18$ digits means asking for the residue modulo $10^{18}$, and since the total is a sum of non-negative exact terms, reducing the running sum is legitimate. A sum below $10^{18}$ plus one term, itself at most $1.8\cdot10^{15}$ (measured maximum $782\,745\,790\,462\,071$), never leaves the unsigned $64$-bit range, ^[8] so one conditional reduction per step suffices.

How the Code Works

The C++ file is the reference implementation. ^[9] build_tau(limit) is Section 4: it allocates lp (32-bit smallest prime factors), cnt (8-bit exponents) and tau (16-bit divisor counts), and in the inner loop writes tau[x] = tau[i] * 2 when p != lp[i] and otherwise tau[x] = tau[i] / (cnt[i] + 1) * (cnt[x] + 1), the factor $\frac{e+2}{e+1}$ with the exact division first. Only tau is returned, so lp and cnt are released before counting starts. main passes kTargetN + 1, because d_triangle needs tau[n + 1].

d_triangle(n, tau) is the parity split of Section 3, and Fenwick is the structure of Section 9: add walks up with idx += idx & -idx, sum_prefix down with idx -= idx & -idx, sum_range subtracting two prefixes.

count_triples_mod(n, tau) carries Sections 6 to 10 in two passes over $1,\dots,n$. The first finds max_d; the second builds count_fw and pair_fw of size max_d + 2 and, for each $i$, reads triples_ending_here from pair_fw above rank = d + 1, adds it to answer with one conditional %= kMod, reads greater_left from count_fw, then updates both trees. Querying first matches the invariant that the tables describe earlier indices only, though it is not forced: both queries are strict suffix sums from rank + 1, which an update at rank leaves untouched, and swapping the steps returns the same totals. Beforehand run_checkpoints verifies d_triangle(7) = 6, $Tr(20)=14$, $Tr(100)=5772$ and $Tr(1000)=11174776$, unless --skip-checkpoints is given.

The Java port holds both trees and the running residue in signed long, which the bound of Section 10 also fits, ^[8] and it carries no checkpoints of its own. Three further differences are visible. It calls buildTau(kTargetN + 2) and loops with i < limit, reaching the same top index $N+1$. It caches every divisor count in dtArr, so the counting pass avoids recomputing them, at the cost of about $240$ MB of extra storage. And its tau is short[] and cnt is byte[], both signed; Java's integer operators do not signal overflow, ^[8] so these would wrap silently, but the measured maxima $672$ and $25$ sit far inside their ranges. And one bound is hard-wired rather than derived: primes is a fixed int[4000000], ample for the $3\,562\,115$ primes below $N+1$ but too small once the limit passes about $6.8\cdot10^{7}$.

The Python file is a bridge, not a port: solve compiles solutionsCpp/Euler378.cpp at -std=c++17 -O2 into a cached binary, rebuilt only when the source is newer, and parse_output returns the printed string. Note that should_skip_cpp_checkpoints searches the source for the literal --skip-checkpoints and finds it, because the C++ file defines that very flag, so the bridge always passes it and the four checkpoints never run on this path. No Python integer arithmetic touches the count, though a native port would need no bound at all, Python integers being unbounded. ^[8] One detail holds for all three: the residue is printed without padding, so a residue below $10^{17}$ would print fewer than $18$ characters, leading zeros and all.

Complexity and Verification

The sieve does $O(N)$ work, since each composite is marked once, and stores $O(N)$ values. The counting pass does a constant number of Fenwick operations per index, each $O(\log M)$ with $M=23042$, so the total is

$$O(N)+O(N\log M)=O(N\log M).$$

Counting the cell visits at the target gives $2\,519\,999\,956$, $42$ per index, $897\,549\,322$ of them in queries and $1\,622\,450\,634$ in updates; that constant is why the sweep, not the sieve, dominates although both are near-linear. Memory is dominated by the sieve: the C++ run peaks at $437$ MB — $436$ MB of arrays ($240$ MB for lp, $120$ MB for tau, $60$ MB for cnt, $16$ MB reserved for the primes) plus process overhead —, and once build_tau returns only tau and the two $184$ kB trees remain; the Java run peaks at $729$ MB because of dtArr. Measured on one machine, the C++ run took $1.7$ s including its checkpoints and the Java run $1.5$ s; such figures depend on the processor and the runtime, unlike the array sizes and operation counts above.

The following was checked for this article. A direct triple loop agrees with the sweep for every $n\le60$, a quadratic pair-then-triple count for every $n\le300$. The parity split of Section 3 matched a plain divisor-count sieve for every $n$ with $T(n)\le80\,002$ and on $400$ random $n$ up to $3\cdot10^{6}$. An implementation that drops the Fenwick trees entirely — plain arrays over the compressed value axis, straight loops for the suffix sums, $128$-bit accumulators — reproduces $Tr(20)=14$, $Tr(100)=5772$, $Tr(1000)=11174776$, $Tr(10^4)=14\,544\,175\,026$, $Tr(10^5)=16\,584\,930\,838\,296$ and $Tr(10^6)=17\,911\,338\,911\,265\,382$ with the matching pair counts. At full size the decreasing-pair total $790\,912\,595\,361\,069$ was recomputed by a square-root decomposition over the value axis instead of a Fenwick tree and came out identical. Finally the published C++ and Java programs print the same $18$ digits, and the C++ passes its four checkpoints.

Footnotes and References

The sources below supply the background facts used above. The coprime split, the counting identity and the accumulator bounds are derived here; the numerical figures are measurements, listed with their checks in the verification section.

Project Euler 378 — Triangle Triples. The official statement defines T(n) = n(n + 1)/2, lets dT(n) count the divisors of T(n) with the example T(7) = 28 and dT(7) = 6, defines Tr(n) as the number of increasing index triples bounded by n whose dT values strictly decrease, and asks for the last 18 digits of Tr(60 000 000). Its stated Tr(20) = 14, Tr(100) = 5772 and Tr(1000) = 11174776 are the checkpoints used throughout.
OEIS A063440 — number of divisors of the n-th triangular number. The entry lists 1, 2, 4, 4, 4, 4, 6, 9, 6, 4, 8, 8, 4, 8, 16, 8, 6, 6, 8, 16, …, the values D₁ to D₂₀ of Sections 1 and 7, and records the parity split of Section 3 as a(2k) = d(k)·d(2k+1), a(2k+1) = d(2k+1)·d(k+1).
Wikipedia — Triangular number. The article derives T(n) = n(n + 1)/2 and recommends halving whichever of n and n + 1 is even before multiplying, but does not say that the two factors are then coprime, which is why Section 2 proves it.
Wikipedia — Divisor function. For a prime factorisation with exponents aᵢ it gives σ₀(n) = ∏ (aᵢ + 1), the formula of Section 3, and states that σₓ is multiplicative, so σₓ(ab) = σₓ(a)σₓ(b) for coprime a and b. It also records that d(n) grows more slowly than any positive power of n, which is why the value axis of Section 8 is far shorter than the index axis; its exact maximum there comes from scanning the computed values.
cp-algorithms — Linear Sieve. The page builds lp[i], the smallest prime factor of i, appending i to the prime list when lp[i] is zero and setting lp[i·p] = p for each prime p up to lp[i], and is linear because every composite has exactly one such representation. It does not discuss divisor counts, so the recurrence of Section 4 is derived here.
Wikipedia — Inversion (discrete mathematics). An inversion is a pair of positions whose entries are out of order, and the inversion number of a sequence of length n lies between 0 and n(n − 1)/2. The lower layer of Section 6 is exactly this count, hence its bound C(N, 2).
cp-algorithms — Fenwick Tree. The page gives the binary indexed tree: point updates and prefix queries in O(log N) time and O(N) memory, a general range as the difference of two prefixes, and the index steps i − (i & −i) and i + (i & −i). It works in prefix terms only, hence Section 9's difference.
cppreference — fixed width integer types documents std::uint16_t, std::uint32_t and std::uint64_t as exact-width unsigned types, so uint64_t holds every value to 18446744073709551615. The Java Language Specification fixes the range of short as −32768 to 32767 in §4.2.1 and states in §4.2.2 that a short operand is widened to int before arithmetic and that the integer operators do not indicate overflow or underflow in any way. Section 10's bound also lies below 2⁶³ − 1, so Java's signed long is safe too. Python — numeric types states that Python integers have unlimited precision.
C++ — Euler378.cpp. The linked immutable revision contains build_tau, d_triangle, Fenwick, count_triples_mod and run_checkpoints, and is the implementation record for this article. The independent comparisons above are extra checks, not tests packaged in that file.

Problem 378 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using u16 = std::uint16_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;

constexpr int kTargetN = 60000000;
constexpr u64 kMod = 1000000000000000000ULL;  // last 18 digits

struct Fenwick {
    std::vector<u64> bit;

    explicit Fenwick(const int n) : bit(static_cast<std::size_t>(n + 1), 0ULL) {}

    void add(int idx, const u64 value) {
        const int n = static_cast<int>(bit.size()) - 1;
        while (idx <= n) {
            bit[static_cast<std::size_t>(idx)] += value;
            idx += idx & -idx;
        }
    }

    u64 sum_prefix(int idx) const {
        u64 s = 0ULL;
        while (idx > 0) {
            s += bit[static_cast<std::size_t>(idx)];
            idx -= idx & -idx;
        }
        return s;
    }

    u64 sum_range(const int l, const int r) const {
        if (r < l) {
            return 0ULL;
        }
        return sum_prefix(r) - sum_prefix(l - 1);
    }
};

std::vector<u16> build_tau(const int limit) {
    std::vector<u32> lp(static_cast<std::size_t>(limit + 1), 0U);
    std::vector<u16> tau(static_cast<std::size_t>(limit + 1), 0U);
    std::vector<std::uint8_t> cnt(static_cast<std::size_t>(limit + 1), 0U);
    std::vector<int> primes;
    primes.reserve(4000000);

    tau[1] = 1U;
    for (int i = 2; i <= limit; ++i) {
        if (lp[static_cast<std::size_t>(i)] == 0U) {
            lp[static_cast<std::size_t>(i)] = static_cast<u32>(i);
            primes.push_back(i);
            tau[static_cast<std::size_t>(i)] = 2U;
            cnt[static_cast<std::size_t>(i)] = 1U;
        }

        for (const int p : primes) {
            const long long x = 1LL * i * p;
            if (x > limit) {
                break;
            }

            lp[static_cast<std::size_t>(x)] = static_cast<u32>(p);
            if (p == static_cast<int>(lp[static_cast<std::size_t>(i)])) {
                cnt[static_cast<std::size_t>(x)] = static_cast<std::uint8_t>(cnt[static_cast<std::size_t>(i)] + 1U);
                tau[static_cast<std::size_t>(x)] = static_cast<u16>(
                    tau[static_cast<std::size_t>(i)] / (cnt[static_cast<std::size_t>(i)] + 1U) *
                    (cnt[static_cast<std::size_t>(x)] + 1U));
                break;
            } else {
                cnt[static_cast<std::size_t>(x)] = 1U;
                tau[static_cast<std::size_t>(x)] = static_cast<u16>(tau[static_cast<std::size_t>(i)] * 2U);
            }
        }
    }

    return tau;
}

inline int d_triangle(const int n, const std::vector<u16>& tau) {
    if ((n & 1) == 0) {
        return static_cast<int>(tau[static_cast<std::size_t>(n / 2)]) * tau[static_cast<std::size_t>(n + 1)];
    }
    return static_cast<int>(tau[static_cast<std::size_t>(n)]) * tau[static_cast<std::size_t>((n + 1) / 2)];
}

u64 count_triples_mod(const int n, const std::vector<u16>& tau) {
    int max_d = 0;
    for (int i = 1; i <= n; ++i) {
        const int d = d_triangle(i, tau);
        if (d > max_d) {
            max_d = d;
        }
    }

    const int max_rank = max_d + 2;
    Fenwick count_fw(max_rank);
    Fenwick pair_fw(max_rank);

    u64 answer = 0ULL;
    for (int i = 1; i <= n; ++i) {
        const int d = d_triangle(i, tau);
        const int rank = d + 1;

        const u64 triples_ending_here = pair_fw.sum_range(rank + 1, max_rank);
        answer += triples_ending_here;
        if (answer >= kMod) {
            answer %= kMod;
        }

        const u64 greater_left = count_fw.sum_range(rank + 1, max_rank);
        pair_fw.add(rank, greater_left);
        count_fw.add(rank, 1ULL);
    }

    return answer % kMod;
}

bool run_checkpoints(const std::vector<u16>& tau) {
    if (d_triangle(7, tau) != 6) {
        std::cerr << "Checkpoint failed: dT(7)\n";
        return false;
    }
    if (count_triples_mod(20, tau) != 14ULL) {
        std::cerr << "Checkpoint failed: Tr(20)\n";
        return false;
    }
    if (count_triples_mod(100, tau) != 5772ULL) {
        std::cerr << "Checkpoint failed: Tr(100)\n";
        return false;
    }
    if (count_triples_mod(1000, tau) != 11174776ULL) {
        std::cerr << "Checkpoint failed: Tr(1000)\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;
        }
    }

    const std::vector<u16> tau = build_tau(kTargetN + 1);

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

    const u64 answer = count_triples_mod(kTargetN, tau);
    std::cout << answer << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


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

Java

public class Euler378 {

    static final int kTargetN = 60000000;
    static final long kMod = 1000000000000000000L;

    static class Fenwick {
        long[] bit;

        Fenwick(int n) {
            bit = new long[n + 1];
        }

        void add(int idx, long val) {
            while (idx < bit.length) {
                bit[idx] += val;
                idx += idx & -idx;
            }
        }

        long sumPrefix(int idx) {
            long s = 0;
            while (idx > 0) {
                s += bit[idx];
                idx -= idx & -idx;
            }
            return s;
        }

        long sumRange(int l, int r) {
            if (r < l)
                return 0;
            return sumPrefix(r) - sumPrefix(l - 1);
        }
    }

    static short[] buildTau(int limit) {
        int[] lp = new int[limit];
        short[] tau = new short[limit];
        byte[] cnt = new byte[limit];
        int[] primes = new int[4000000];
        int primeCount = 0;

        tau[1] = 1;
        for (int i = 2; i < limit; i++) {
            if (lp[i] == 0) {
                lp[i] = i;
                primes[primeCount++] = i;
                tau[i] = 2;
                cnt[i] = 1;
            }

            int lpi = lp[i];
            short tauI = tau[i];
            byte cntI = cnt[i];

            for (int j = 0; j < primeCount; j++) {
                int p = primes[j];
                long x = (long) i * p;
                if (x >= limit)
                    break;

                lp[(int) x] = p;
                if (p == lpi) {
                    cnt[(int) x] = (byte) (cntI + 1);
                    tau[(int) x] = (short) (tauI / (cntI + 1) * (cntI + 2));
                    break;
                } else {
                    cnt[(int) x] = 1;
                    tau[(int) x] = (short) (tauI * 2);
                }
            }
        }
        return tau;
    }

    static int dTriangle(int n, short[] tau) {
        if ((n & 1) == 0) {
            return tau[n / 2] * tau[n + 1];
        }
        return tau[n] * tau[(n + 1) / 2];
    }

    static String solve() {
        short[] tau = buildTau(kTargetN + 2);

        int maxD = 0;
        int[] dtArr = new int[kTargetN + 1];
        for (int i = 1; i <= kTargetN; i++) {
            int d = dTriangle(i, tau);
            dtArr[i] = d;
            if (d > maxD)
                maxD = d;
        }

        int maxRank = maxD + 2;
        Fenwick countFw = new Fenwick(maxRank);
        Fenwick pairFw = new Fenwick(maxRank);

        long answer = 0;
        for (int i = 1; i <= kTargetN; i++) {
            int rank = dtArr[i] + 1;

            long triplesEndingHere = pairFw.sumRange(rank + 1, maxRank);
            answer = (answer + triplesEndingHere) % kMod;

            long greaterLeft = countFw.sumRange(rank + 1, maxRank);
            pairFw.add(rank, greaterLeft);
            countFw.add(rank, 1);
        }

        return Long.toString(answer);
    }

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