Problem 955: Finding Triangles
View on Project EulerProject Euler Problem 955 Solution
EulerSolve provides an optimized solution for Project Euler Problem 955, Finding Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary Let \(T_r=\frac{r(r+1)}2\) be the \(r\)-th triangular number, with \(T_0=0\). The implementations start from \(T_2=3\) and repeatedly look for the next nontrivial identity of the form $$T_n-T_w=T_m,$$ where \(n>m\) and \(w>0\). At each hit, the solver chooses the smallest possible later index \(n\), adds the corresponding offset \(w\) to the running search position, replaces \(m\) by \(n\), and continues. The required output is the accumulated position of the 70th hit, not the triangular value itself. A naive search would test huge numbers of pairs \((n,w)\) for every current \(m\). The local solutions avoid that completely by turning the existence of such a triangle identity into a divisor problem for \(m(m+1)\), so each step becomes: factor two consecutive integers, enumerate divisors, and jump directly to the next valid state. Mathematical Approach Rewriting a triangular difference as a factorization Start from $$T_n-T_w=\frac{n(n+1)-w(w+1)}2=T_m=\frac{m(m+1)}2.$$ Multiplying by 2 and factoring the left-hand side gives $$n(n+1)-w(w+1)=(n-w)(n+w+1)=m(m+1).$$ Define $$P=m(m+1),\qquad u=n-w,\qquad v=n+w+1.$$ Then every nontrivial hit corresponds to a divisor pair $$uv=P.$$ Because \(u+v=2n+1\), the sum \(u+v\) must be odd. Conversely, if \(u\) and \(v\) are divisors of \(P\) with odd sum, then they reconstruct integers \(n\) and \(w\)....
Detailed mathematical approach
Problem Summary
Let \(T_r=\frac{r(r+1)}2\) be the \(r\)-th triangular number, with \(T_0=0\). The implementations start from \(T_2=3\) and repeatedly look for the next nontrivial identity of the form
$$T_n-T_w=T_m,$$
where \(n>m\) and \(w>0\). At each hit, the solver chooses the smallest possible later index \(n\), adds the corresponding offset \(w\) to the running search position, replaces \(m\) by \(n\), and continues. The required output is the accumulated position of the 70th hit, not the triangular value itself.
A naive search would test huge numbers of pairs \((n,w)\) for every current \(m\). The local solutions avoid that completely by turning the existence of such a triangle identity into a divisor problem for \(m(m+1)\), so each step becomes: factor two consecutive integers, enumerate divisors, and jump directly to the next valid state.
Mathematical Approach
Rewriting a triangular difference as a factorization
Start from
$$T_n-T_w=\frac{n(n+1)-w(w+1)}2=T_m=\frac{m(m+1)}2.$$
Multiplying by 2 and factoring the left-hand side gives
$$n(n+1)-w(w+1)=(n-w)(n+w+1)=m(m+1).$$
Define
$$P=m(m+1),\qquad u=n-w,\qquad v=n+w+1.$$
Then every nontrivial hit corresponds to a divisor pair
$$uv=P.$$
Because \(u+v=2n+1\), the sum \(u+v\) must be odd. Conversely, if \(u\) and \(v\) are divisors of \(P\) with odd sum, then they reconstruct integers \(n\) and \(w\).
Recovering the next triangular state
Solving the two linear equations for \(n\) and \(w\) gives
$$n=\frac{u+v-1}{2},\qquad w=\frac{v-u-1}{2}.$$
The condition \(w>0\) is exactly
$$v>u+1.$$
This matters because the divisor pair \((u,v)=(m,m+1)\) always exists and gives \(n=m,\ w=0\), which is only the trivial identity \(T_m-T_0=T_m\). The implementations deliberately reject that case and keep only genuinely later hits.
Why the largest admissible divisor gives the nearest hit
For fixed \(P\) and \(u\le \sqrt{P}\), we have \(v=P/u\), so
$$n(u)=\frac{u+P/u-1}{2}.$$
On the interval \(1\le u\le \sqrt{P}\), increasing \(u\) decreases \(P/u\), and \(n(u)\) becomes smaller. Therefore the smallest later triangle index \(n\) is obtained by taking the largest admissible divisor \(u\) not exceeding \(\sqrt{P}\). That is exactly what the depth-first divisor search is tracking.
The factorization itself is also simplified by the identity
$$\gcd(m,m+1)=1.$$
So instead of factoring the large product \(P\) directly, the implementations factor \(m\) and \(m+1\) separately and merge their prime exponents.
Worked example: from \(T_6\) to the next hit
Take \(m=6\). Then
$$T_6=21,\qquad P=m(m+1)=42.$$
The divisor pairs with \(u\le\sqrt{42}\) are
$$ (1,42),\ (2,21),\ (3,14),\ (6,7). $$
All four pairs have odd sum, but \((6,7)\) gives \(w=0\), so it is the trivial representation and is discarded. Among the remaining pairs, the largest admissible \(u\) is \(u=3\), with \(v=14\). Hence
$$n=\frac{3+14-1}{2}=8,\qquad w=\frac{14-3-1}{2}=5.$$
Indeed,
$$T_8-T_5=36-15=21=T_6.$$
The other admissible pairs would produce later hits such as \(T_{11}-T_9\) and \(T_{21}-T_{20}\), but the algorithm wants the nearest one, so it keeps \(n=8\).
How the Code Works
Factoring \(m\) and \(m+1\)
Each iteration begins with the current triangular index \(m\). The implementations factor \(m\) and \(m+1\) separately using deterministic Miller-Rabin primality checks for 64-bit inputs together with Pollard-Rho for composite splitting. This is much faster than scanning divisors naively, and it exploits the coprimality of consecutive integers.
Enumerating admissible divisor pairs
Once the prime exponents are known, a recursive divisor generator enumerates every divisor \(u\) of \(P=m(m+1)\). Only candidates with \(u\le\sqrt{P}\) are needed, because the partner divisor is \(v=P/u\). For each candidate, the code tests the two number-theoretic conditions already derived:
$$u+v\equiv 1\pmod{2},\qquad v>u+1.$$
The best divisor is simply the largest \(u\) that passes both tests. After that, the code converts \((u,v)\) into the next state \((w,n)\) via the closed formulas above.
Updating the hit chain
The running answer is advanced by the offset \(w\), and the current triangular index is replaced by \(n\). Repeating this step-by-step jump avoids any brute-force scan through failed candidates. The C++, Python, and Java implementations all follow the same mathematics; they differ only in low-level integer handling. They also check a common validation point: the 10th hit occurs at accumulated position \(2964\) and lands on \(T_{1696}=1{,}439{,}056\).
Complexity Analysis
If \(H\) hits are required, the total work is the sum of \(H-1\) transition costs. One transition consists of factoring two consecutive integers and then enumerating the divisors of their product. If \(\tau(P)\) is the number of divisors of \(P=m(m+1)\), the divisor search is \(O(\tau(P))\) because each divisor is generated once.
The harder part is integer factorization. Pollard-Rho is a randomized algorithm, so the cleanest statement is heuristic: in practice the runtime is dominated by a modest number of 64-bit factorizations, and the method is easily fast enough for the requested 70th hit. Memory use stays small, since the program stores only a prime-exponent map, a recursion stack for divisor generation, and a handful of wide integers for \(P\), \(\sqrt{P}\), \(n\), and \(w\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=955
- Triangular number: Wikipedia - Triangular number
- Integer factorization: Wikipedia - Integer factorization
- Pollard's rho algorithm: Wikipedia - Pollard's rho algorithm
- Miller-Rabin primality test: Wikipedia - Miller-Rabin primality test
Problem 955 source code
C++
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <map>
#include <numeric>
#include <random>
#include <vector>
namespace {
using u64 = std::uint64_t;
using u128 = unsigned __int128;
u64 mul_mod_u64(u64 a, u64 b, u64 mod) {
return static_cast<u64>((static_cast<u128>(a) * b) % mod);
}
u64 pow_mod_u64(u64 base, u64 exp, u64 mod) {
u64 result = 1 % mod;
u64 cur = base % mod;
while (exp > 0) {
if ((exp & 1ULL) != 0ULL) {
result = mul_mod_u64(result, cur, mod);
}
cur = mul_mod_u64(cur, cur, mod);
exp >>= 1ULL;
}
return result;
}
bool is_prime_u64(u64 n) {
if (n < 2ULL) {
return false;
}
for (u64 p : {2ULL, 3ULL, 5ULL, 7ULL, 11ULL, 13ULL, 17ULL, 19ULL, 23ULL, 29ULL, 31ULL, 37ULL}) {
if (n == p) {
return true;
}
if (n % p == 0ULL) {
return false;
}
}
u64 d = n - 1ULL;
int s = 0;
while ((d & 1ULL) == 0ULL) {
d >>= 1ULL;
++s;
}
for (u64 a : {2ULL, 325ULL, 9'375ULL, 28'178ULL, 450'775ULL, 9'780'504ULL, 1'795'265'022ULL}) {
if (a % n == 0ULL) {
continue;
}
u64 x = pow_mod_u64(a % n, d, n);
if (x == 1ULL || x == n - 1ULL) {
continue;
}
bool witness = true;
for (int r = 1; r < s; ++r) {
x = mul_mod_u64(x, x, n);
if (x == n - 1ULL) {
witness = false;
break;
}
}
if (witness) {
return false;
}
}
return true;
}
u64 pollard_rho(u64 n, std::mt19937_64& rng) {
if ((n & 1ULL) == 0ULL) {
return 2ULL;
}
if (n % 3ULL == 0ULL) {
return 3ULL;
}
std::uniform_int_distribution<u64> dist(2ULL, n - 2ULL);
while (true) {
const u64 c = dist(rng);
u64 x = dist(rng);
u64 y = x;
u64 d = 1ULL;
auto f = [&](u64 v) {
return (mul_mod_u64(v, v, n) + c) % n;
};
while (d == 1ULL) {
x = f(x);
y = f(f(y));
const u64 diff = x > y ? x - y : y - x;
d = std::gcd(diff, n);
}
if (d != n) {
return d;
}
}
}
void factor_u64(u64 n, std::map<u64, int>& factors, std::mt19937_64& rng) {
if (n == 1ULL) {
return;
}
if (is_prime_u64(n)) {
++factors[n];
return;
}
const u64 d = pollard_rho(n, rng);
factor_u64(d, factors, rng);
factor_u64(n / d, factors, rng);
}
u128 isqrt_u128(u128 x) {
long double approx = std::sqrt(static_cast<long double>(x));
u128 r = static_cast<u128>(approx);
while ((r + 1) <= x / (r + 1)) {
++r;
}
while (r > x / r) {
--r;
}
return r;
}
u128 triangle_u128(u128 m) {
return (m * (m + 1)) / 2;
}
std::pair<u128, u128> next_triangle_state(u128 m, std::mt19937_64& rng) {
const u64 a = static_cast<u64>(m);
const u64 b = static_cast<u64>(m + 1);
std::map<u64, int> fac;
factor_u64(a, fac, rng);
factor_u64(b, fac, rng);
std::vector<std::pair<u64, int>> pf(fac.begin(), fac.end());
const u128 p = m * (m + 1);
const u128 root = isqrt_u128(p);
u128 best_u = 0;
auto dfs = [&](auto&& self, std::size_t idx, u128 cur) -> void {
if (idx == pf.size()) {
if (cur > root) {
return;
}
const u128 v = p / cur;
if (((cur + v) & 1U) == 1U && v > cur + 1 && cur > best_u) {
best_u = cur;
}
return;
}
const auto [prime, exp] = pf[idx];
u128 val = 1;
for (int e = 0; e <= exp; ++e) {
self(self, idx + 1, cur * val);
val *= static_cast<u128>(prime);
}
};
dfs(dfs, 0, 1);
assert(best_u > 0);
const u128 v = p / best_u;
const u128 jump = (v - best_u - 1) / 2;
const u128 next_m = (best_u + v - 1) / 2;
return {jump, next_m};
}
std::pair<u128, u128> solve_kth_triangle_hit(int k) {
std::mt19937_64 rng(0xC0FFEEULL);
int hit = 1;
u128 index = 0;
u128 m = 2;
while (hit < k) {
const auto [jump, next_m] = next_triangle_state(m, rng);
index += jump;
m = next_m;
++hit;
}
return {index, m};
}
std::string to_string_u128(u128 x) {
if (x == 0) {
return "0";
}
std::string s;
while (x > 0) {
const int digit = static_cast<int>(x % 10);
s.push_back(static_cast<char>('0' + digit));
x /= 10;
}
std::reverse(s.begin(), s.end());
return s;
}
void run_validations() {
const auto [idx10, m10] = solve_kth_triangle_hit(10);
assert(idx10 == 2'964);
assert(triangle_u128(m10) == 1'439'056);
}
} // namespace
int main() {
run_validations();
const auto [idx70, _m70] = solve_kth_triangle_hit(70);
std::cout << to_string_u128(idx70) << '\n';
return 0;
}
Python
import math
import random
def mul_mod_u64(a, b, mod):
return (a * b) % mod
def pow_mod_u64(base, exp, mod):
return pow(base, exp, mod)
def is_prime_u64(n):
if n < 2:
return False
for p in [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]:
if n == p:
return True
if n % p == 0:
return False
d = n - 1
s = 0
while (d & 1) == 0:
d >>= 1
s += 1
for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
if a % n == 0:
continue
x = pow_mod_u64(a % n, d, n)
if x == 1 or x == n - 1:
continue
witness = True
for _ in range(1, s):
x = mul_mod_u64(x, x, n)
if x == n - 1:
witness = False
break
if witness:
return False
return True
def pollard_rho(n, rng):
if (n & 1) == 0:
return 2
if n % 3 == 0:
return 3
while True:
c = rng.randint(2, n - 2)
x = rng.randint(2, n - 2)
y = x
d = 1
def f(v):
return (v * v + c) % n
while d == 1:
x = f(x)
y = f(f(y))
diff = x - y if x > y else y - x
d = math.gcd(diff, n)
if d != n:
return d
def factor_u64(n, factors, rng):
if n == 1:
return
if is_prime_u64(n):
factors[n] = factors.get(n, 0) + 1
return
d = pollard_rho(n, rng)
factor_u64(d, factors, rng)
factor_u64(n // d, factors, rng)
def isqrt_u128(x):
if x == 0:
return 0
r = int(math.sqrt(x))
while (r + 1) * (r + 1) <= x:
r += 1
while r * r > x:
r -= 1
return r
def triangle_u128(m):
return (m * (m + 1)) // 2
def next_triangle_state(m, rng):
a = m
b = m + 1
fac = {}
factor_u64(a, fac, rng)
factor_u64(b, fac, rng)
pf = list(fac.items())
p = m * (m + 1)
root = isqrt_u128(p)
best_u = [0]
def dfs(idx, cur):
if idx == len(pf):
if cur > root:
return
v = p // cur
if ((cur + v) & 1) == 1 and v > cur + 1 and cur > best_u[0]:
best_u[0] = cur
return
prime, exp = pf[idx]
val = 1
for _ in range(exp + 1):
dfs(idx + 1, cur * val)
val *= prime
dfs(0, 1)
v = p // best_u[0]
jump = (v - best_u[0] - 1) // 2
next_m = (best_u[0] + v - 1) // 2
return jump, next_m
def solve_kth_triangle_hit(k):
rng = random.Random(0xC0FFEE)
hit = 1
index = 0
m = 2
while hit < k:
jump, next_m = next_triangle_state(m, rng)
index += jump
m = next_m
hit += 1
return index, m
def solve():
idx70, _ = solve_kth_triangle_hit(70)
return str(idx70)
if __name__ == "__main__":
idx10, m10 = solve_kth_triangle_hit(10)
assert idx10 == 2964
assert triangle_u128(m10) == 1439056
print(solve())
Java
import java.math.BigInteger;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Random;
public class Euler955 {
static long mulModU64(long a, long b, long mod) {
return BigInteger.valueOf(a).multiply(BigInteger.valueOf(b)).mod(BigInteger.valueOf(mod)).longValue();
}
static long powModU64(long base, long exp, long mod) {
long result = 1 % mod;
long cur = base % mod;
while (exp > 0) {
if ((exp & 1L) != 0L) {
result = mulModU64(result, cur, mod);
}
cur = mulModU64(cur, cur, mod);
exp >>= 1L;
}
return result;
}
static boolean isPrimeU64(long n) {
if (n < 2L) {
return false;
}
long[] smallPrimes = { 2L, 3L, 5L, 7L, 11L, 13L, 17L, 19L, 23L, 29L, 31L, 37L };
for (long p : smallPrimes) {
if (n == p)
return true;
if (n % p == 0L)
return false;
}
long d = n - 1L;
int s = 0;
while ((d & 1L) == 0L) {
d >>= 1L;
++s;
}
long[] aValues = { 2L, 325L, 9375L, 28178L, 450775L, 9780504L, 1795265022L };
for (long a : aValues) {
if (a % n == 0L)
continue;
long x = powModU64(a % n, d, n);
if (x == 1L || x == n - 1L)
continue;
boolean witness = true;
for (int r = 1; r < s; ++r) {
x = mulModU64(x, x, n);
if (x == n - 1L) {
witness = false;
break;
}
}
if (witness)
return false;
}
return true;
}
static long gcd(long a, long b) {
while (b != 0) {
long t = b;
b = a % b;
a = t;
}
return a;
}
static long nextLong(Random rng, long n) {
long bits, val;
do {
bits = (rng.nextLong() << 1) >>> 1;
val = bits % n;
} while (bits - val + (n - 1) < 0L);
return val;
}
static long pollardRho(long n, Random rng) {
if ((n & 1L) == 0L)
return 2L;
if (n % 3L == 0L)
return 3L;
while (true) {
long c = nextLong(rng, n - 3) + 2;
long x = nextLong(rng, n - 3) + 2;
long y = x;
long d = 1L;
while (d == 1L) {
x = (mulModU64(x, x, n) + c) % n;
y = (mulModU64(y, y, n) + c) % n;
y = (mulModU64(y, y, n) + c) % n;
long diff = x > y ? x - y : y - x;
d = gcd(diff, n);
}
if (d != n)
return d;
}
}
static void factorU64(long n, Map<Long, Integer> factors, Random rng) {
if (n == 1L)
return;
if (isPrimeU64(n)) {
factors.put(n, factors.getOrDefault(n, 0) + 1);
return;
}
long d = pollardRho(n, rng);
factorU64(d, factors, rng);
factorU64(n / d, factors, rng);
}
static BigInteger isqrtU128(BigInteger x) {
if (x.equals(BigInteger.ZERO))
return BigInteger.ZERO;
BigInteger s = new BigInteger(String.format("%.0f", Math.sqrt(x.doubleValue())));
while (s.add(BigInteger.ONE).pow(2).compareTo(x) <= 0) {
s = s.add(BigInteger.ONE);
}
while (s.pow(2).compareTo(x) > 0) {
s = s.subtract(BigInteger.ONE);
}
return s;
}
static BigInteger triangleU128(BigInteger m) {
return m.multiply(m.add(BigInteger.ONE)).divide(BigInteger.valueOf(2));
}
static class Pair {
BigInteger first;
BigInteger second;
Pair(BigInteger first, BigInteger second) {
this.first = first;
this.second = second;
}
}
static class Factor {
long prime;
int exp;
Factor(long prime, int exp) {
this.prime = prime;
this.exp = exp;
}
}
static Pair nextTriangleState(BigInteger m, Random rng) {
long a = m.longValue();
long b = m.add(BigInteger.ONE).longValue();
Map<Long, Integer> fac = new HashMap<>();
factorU64(a, fac, rng);
factorU64(b, fac, rng);
List<Factor> pf = new ArrayList<>();
for (Map.Entry<Long, Integer> entry : fac.entrySet()) {
pf.add(new Factor(entry.getKey(), entry.getValue()));
}
BigInteger p = m.multiply(m.add(BigInteger.ONE));
BigInteger root = isqrtU128(p);
BigInteger[] bestU = { BigInteger.ZERO };
dfs(pf, 0, BigInteger.ONE, p, root, bestU);
BigInteger v = p.divide(bestU[0]);
BigInteger jump = v.subtract(bestU[0]).subtract(BigInteger.ONE).divide(BigInteger.valueOf(2));
BigInteger nextM = bestU[0].add(v).subtract(BigInteger.ONE).divide(BigInteger.valueOf(2));
return new Pair(jump, nextM);
}
static void dfs(List<Factor> pf, int idx, BigInteger cur, BigInteger p, BigInteger root, BigInteger[] bestU) {
if (idx == pf.size()) {
if (cur.compareTo(root) > 0)
return;
BigInteger v = p.divide(cur);
if (cur.add(v).testBit(0) && v.compareTo(cur.add(BigInteger.ONE)) > 0 && cur.compareTo(bestU[0]) > 0) {
bestU[0] = cur;
}
return;
}
Factor factor = pf.get(idx);
BigInteger val = BigInteger.ONE;
for (int e = 0; e <= factor.exp; ++e) {
dfs(pf, idx + 1, cur.multiply(val), p, root, bestU);
val = val.multiply(BigInteger.valueOf(factor.prime));
}
}
static Pair solveKthTriangleHit(int k) {
Random rng = new Random(0xC0FFEE);
int hit = 1;
BigInteger index = BigInteger.ZERO;
BigInteger m = BigInteger.valueOf(2);
while (hit < k) {
Pair next = nextTriangleState(m, rng);
index = index.add(next.first);
m = next.second;
++hit;
}
return new Pair(index, m);
}
public static String solve() {
return solveKthTriangleHit(70).first.toString();
}
public static void main(String[] args) {
Pair res10 = solveKthTriangleHit(10);
if (!res10.first.equals(BigInteger.valueOf(2964))
|| !triangleU128(res10.second).equals(BigInteger.valueOf(1439056))) {
System.out.println("Validation failed");
return;
}
System.out.println(solve());
}
}