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},\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. For example, $T_7=28=4\cdot 7$, so

$$D_7=d(4)d(7)=3\cdot 2=6,$$

which matches the explicit checkpoint in the C++ solver.

Step 2: Build all divisor counts with a linear sieve

Because the formula above only needs $d(x)$ for integers up to $N+1$, the program first precomputes the divisor-count table once. The arrays in build_tau store:

$lp[i]$: the smallest prime dividing $i$,

$cnt[i]$: the exponent of $lp[i]$ inside $i$,

$d(i)$: the divisor count itself.

$$i=m p^e,\qquad p\nmid m,$$

then

$$d(i)=d(m)(e+1).$$

When the sieve extends $i$ to $x=i p$, there are two cases:

$$d(x)=\begin{cases} d(i)\cdot 2, & p\ne lp[i],\\[2mm] d(i)\cdot \dfrac{e+2}{e+1}, & p=lp[i]. \end{cases}$$

The second formula is exactly what the source code writes as

$$d(x)=\frac{d(i)}{cnt[i]+1}\,(cnt[x]+1).$$

This produces all values up to $N+1$ in linear total work, after which each $D_n$ can be evaluated in constant time from the parity split above.

Step 3: Recast the problem as decreasing subsequences

Now consider the sequence

$$D_1,D_2,\dots,D_N.$$

We want the number of triples $i \lt j \lt k$ with strict decrease. This is a dynamic-programming problem over values rather than over explicit triples.

For each processed position, maintain:

$C(r)$: how many earlier indices have divisor-count rank $r$,

$P(r)$: how many decreasing pairs $(i,j)$ seen so far end with rank $r$.

The code does not need coordinate compression. It first scans all $D_n$ to find

$$M=\max_{1\le n\le N} D_n,$$

then uses the direct 1-based rank

$$r=D_n+1.$$

That is why the Fenwick trees are sized as max_d + 2 in both C++ and Java.

Step 4: Fenwick-tree recurrence

Suppose the current position is $k$ and its rank is $r$. Every earlier value larger than $D_k$ creates a new decreasing pair ending at $k$, so

$$\text{greater\_left}(k)=\sum_{t=r+1}^{M+1} C(t).$$

Similarly, every earlier decreasing pair whose second value is still larger than $D_k$ can be extended to a triple ending at $k$, so

$$\text{triples\_ending\_here}(k)=\sum_{t=r+1}^{M+1} P(t).$$

After computing those two quantities, the update rules are

$$P(r)\leftarrow P(r)+\text{greater\_left}(k),$$

$$C(r)\leftarrow C(r)+1.$$

A Fenwick tree gives each suffix sum and each point update in $O(\log M)$. Two such trees are enough: one for single elements and one for decreasing pairs.

Worked Example: $N=20$

The first twenty values are

$$\left(D_1,\dots,D_{20}\right)=\left(1,2,4,4,4,4,6,9,6,4,8,8,4,8,16,8,6,6,8,16\right).$$

Because the inequalities are strict, equal values never contribute. For instance, $(8,9,10)$ is valid because

$$D_8=9,\qquad D_9=6,\qquad D_{10}=4,$$

and $(15,16,17)$ is valid because

$$16 \gt 8 \gt 6.$$

The total number of valid triples up to $20$ is

$$14,$$

which is exactly the checkpoint count_triples_mod(20) == 14 in the C++ source.

How the Code Works

The C++ implementation is the reference solver. It builds tau up to kTargetN + 1, verifies the checkpoints d_triangle(7)=6, Tr(20)=14, Tr(100)=5772, and Tr(1000)=11174776, then performs two passes over $1,\dots,N$: the first finds max_d, and the second performs the Fenwick-based dynamic programming.

The Java version mirrors the same mathematics, but stores every $D_n$ in an array dtArr during the first pass so that the second pass does not recompute them. The Python file is intentionally thin: it compiles and runs the C++ solver, then parses the produced answer. So the mathematical logic is shared across all three languages, with C++ providing the canonical implementation details.

Complexity Analysis

Let $N$ be the target limit and let $M=\max_{1\le n\le N} D_n$. The divisor-count sieve runs in $O(N)$ time and uses $O(N)$ memory for the sieve arrays. The triple-counting phase performs two Fenwick queries and two Fenwick updates per index, so it costs $O(N\log M)$ time and $O(M)$ memory for the trees. Overall, the algorithm is dominated by

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

time, with $O(N)+O(M)$ memory. In practice $M$ is far smaller than $N$, which is why this approach is feasible for $N=60{,}000{,}000$.

Footnotes and References

Problem page: https://projecteuler.net/problem=378
Triangular numbers: Wikipedia — Triangular number
Divisor function: Wikipedia — Divisor function
Fenwick tree: cp-algorithms — Fenwick Tree
Linear sieve: cp-algorithms — Linear Sieve

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