Problem 621: Expressing an Integer as the Sum of Triangular Numbers
View on Project EulerProject Euler Problem 621 Solution
EulerSolve provides an optimized solution for Project Euler Problem 621, Expressing an Integer as the Sum of Triangular Numbers, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let $T_k=\frac{k(k+1)}{2}\qquad (k\ge 0)$ be the \(k\)-th triangular number. The function \(G(n)\) counts ordered triples \((a,b,c)\) of nonnegative integers such that $T_a+T_b+T_c=n.$ The task is to evaluate \(G(17526\cdot 10^9)\). Directly enumerating triangular triples is hopeless at that scale, so the solution converts the problem into counting representations by three squares and then uses a class-number formula. Mathematical Approach The implementation turns the additive problem for triangular numbers into an arithmetic problem for quadratic forms. Step 1: Transform Triangular Numbers into Odd Squares The identity $8T_k+1=4k(k+1)+1=(2k+1)^2$ shows that every triangular number corresponds to an odd square. Therefore $n=T_a+T_b+T_c$ is equivalent to $8n+3=(2a+1)^2+(2b+1)^2+(2c+1)^2.$ If we write $N=8n+3,$ then \(G(n)\) is the number of ordered representations of \(N\) as a sum of three positive odd squares. Step 2: Relate \(G(n)\) to the Classical Three-Squares Counting Function Define $r_3(N)=\#\left\{(x,y,z)\in\mathbb{Z}^3:x^2+y^2+z^2=N\right\},$ where order and signs are both counted....
Detailed mathematical approach
Problem Summary
Let
$T_k=\frac{k(k+1)}{2}\qquad (k\ge 0)$
be the \(k\)-th triangular number. The function \(G(n)\) counts ordered triples \((a,b,c)\) of nonnegative integers such that
$T_a+T_b+T_c=n.$
The task is to evaluate \(G(17526\cdot 10^9)\). Directly enumerating triangular triples is hopeless at that scale, so the solution converts the problem into counting representations by three squares and then uses a class-number formula.
Mathematical Approach
The implementation turns the additive problem for triangular numbers into an arithmetic problem for quadratic forms.
Step 1: Transform Triangular Numbers into Odd Squares
The identity
$8T_k+1=4k(k+1)+1=(2k+1)^2$
shows that every triangular number corresponds to an odd square. Therefore
$n=T_a+T_b+T_c$
is equivalent to
$8n+3=(2a+1)^2+(2b+1)^2+(2c+1)^2.$
If we write
$N=8n+3,$
then \(G(n)\) is the number of ordered representations of \(N\) as a sum of three positive odd squares.
Step 2: Relate \(G(n)\) to the Classical Three-Squares Counting Function
Define
$r_3(N)=\#\left\{(x,y,z)\in\mathbb{Z}^3:x^2+y^2+z^2=N\right\},$
where order and signs are both counted. Since \(N\equiv 3 \pmod{8}\), each square in any representation must be congruent to \(1 \pmod{8}\), because the only way to obtain \(3\) from three quadratic residues modulo \(8\) is
$1+1+1\equiv 3 \pmod{8}.$
So every valid \(x,y,z\) is odd, and none of them can be zero. Each positive odd triple corresponds to exactly \(2^3=8\) signed triples in \(r_3(N)\). Hence
$G(n)=\frac{r_3(8n+3)}{8}.$
Step 3: Use the Arithmetic Formula for \(r_3(N)\)
Factor \(N\) as
$N=f^2m,$
where \(m\) is squarefree. Because \(N\equiv 3\pmod{8}\), the squarefree part \(m\) is odd and satisfies \(m\equiv 3\pmod{4}\). Set
$D=-m.$
Then \(D\) is a negative fundamental discriminant, and the implementation evaluates
$r_3(N)=\frac{12\cdot 2h(D)\left(1-\left(\frac{D}{2}\right)\right)}{w(D)}\sum_{d\mid f}\mu(d)\left(\frac{D}{d}\right)\sigma\left(\frac{f}{d}\right).$
Here \(h(D)\) is the class number of discriminant \(D\), \(\mu\) is the Möbius function, \(\sigma\) is the sum-of-divisors function, \(\left(\frac{D}{d}\right)\) is the Jacobi or Kronecker character, and \(w(D)\) is the number of units in the quadratic order:
$w(D)=\begin{cases} 6,& D=-3,\\ 4,& D=-4,\\ 2,& \text{otherwise.} \end{cases}$
The sum is multiplicative in the square part \(f\), which is why the factorization of \(N\) is the first major step of the program.
Step 4: Compute \(h(D)\) by Counting Reduced Binary Quadratic Forms
For a negative fundamental discriminant \(D\), the class number can be obtained by counting reduced positive definite binary quadratic forms
$ax^2+bxy+cy^2,\qquad b^2-4ac=D.$
The reduced conditions used by the implementation are
$|b|\le a\le c,$
$b\ge 0\quad\text{when}\quad |b|=a\ \text{or}\ a=c.$
Since \(D<0\), one only needs to scan
$1\le a\le \sqrt{\frac{|D|}{3}}.$
For each odd \(a\), the congruence condition coming from the discriminant is
$b^2\equiv D \pmod{4a}.$
The implementation factors \(a\), finds square roots of \(D\) modulo each prime power, lifts them to higher powers, combines the local roots with the Chinese remainder theorem, and then reconstructs the admissible values of \(b\) and \(c\). Every valid reduced form contributes exactly one class.
Step 5: Worked Example \(n=9\)
For \(n=9\), we have
$N=8\cdot 9+3=75.$
The positive odd-square representations of \(75\) are
$75=1^2+5^2+7^2=5^2+5^2+5^2.$
The first pattern has \(3!=6\) orderings, and each ordering has \(8\) sign choices in \(r_3(75)\). The second pattern contributes \(1\cdot 8\). Therefore
$r_3(75)=6\cdot 8+1\cdot 8=56,$
so
$G(9)=\frac{56}{8}=7.$
In triangular-number language, those seven ordered representations are the six permutations of \(0+3+6\) and the single representation \(3+3+3\).
Step 6: Final Evaluation Strategy
For the target input, the program forms
$N=8\cdot (17526\cdot 10^9)+3=140208000000003,$
factors \(N\), extracts \(f\) and \(m\), computes \(h(-m)\), evaluates the divisor sum in the formula for \(r_3(N)\), and finally divides by \(8\) to obtain \(G(n)\). This avoids any search over triangular numbers themselves.
How the Code Works
The C++, Python, and Java implementations follow the same pipeline. First they factor \(N=8n+3\) with a 64-bit primality test and Pollard-Rho splitting. From that factorization they separate the squarefree part \(m\) from the square part \(f^2\), which gives the discriminant \(D=-m\).
Next they compute the class number \(h(D)\). To do that efficiently, the implementation enumerates candidate values of \(a\), factors each \(a\) with a small-prime sieve, solves the modular square-root conditions prime-power by prime-power, merges the local solutions with the Chinese remainder theorem, and counts the reduced forms that satisfy the discriminant equation.
Once \(h(D)\) is known, the code evaluates the divisor sum
$\sum_{d\mid f}\mu(d)\left(\frac{D}{d}\right)\sigma\left(\frac{f}{d}\right)$
by iterating over the squarefree divisors encoded by the prime factors of \(f\). That produces \(r_3(N)\), and the final answer is \(r_3(N)/8\). The implementations also include checkpoint values such as \(G(9)=7\), \(G(1000)=78\), and \(G(10^6)=2106\).
Complexity Analysis
The brute-force search space for three triangular numbers grows far too quickly, so the arithmetic method is essential. Integer factorization is handled by Miller-Rabin and Pollard-Rho, which are very fast in practice for a single 64-bit input. After factorization, the dominant deterministic work is the class-number computation.
If
$A=\left\lfloor\sqrt{\frac{|D|}{3}}\right\rfloor,$
then the sieve used inside the class-number computation requires \(O(A)\) memory and \(O(A\log\log A)\) preprocessing time. The subsequent enumeration scans odd \(a\le A\); each step performs a small factorization of \(a\), a few modular square-root lifts, and a small CRT merge. In practice this is easily manageable for a single target value and is incomparably faster than enumerating triangular triples directly.
Footnotes and References
- Problem page: https://projecteuler.net/problem=621
- Triangular numbers: Wikipedia — Triangular number
- Three-square theorem: Wikipedia — Legendre's three-square theorem
- Binary quadratic forms: Wikipedia — Binary quadratic form
- Class numbers: Wikipedia — Class number
- Jacobi symbol: Wikipedia — Jacobi symbol