Problem 835: Supernatural Triangles

View on Project Euler

Project Euler Problem 835 Solution

EulerSolve provides an optimized solution for Project Euler Problem 835, Supernatural Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The contributing triangles split into two infinite Pythagorean families. One family has hypotenuse exactly one larger than a leg, and the other has two consecutive legs. The task is to sum every distinct perimeter not exceeding \(N=10^{10^{10}}\) and return the result modulo \(1234567891\). Because \(N\) is unimaginably large, the solution replaces any geometric search with closed formulas, a Pell-type recurrence, and one overlap correction. Mathematical Approach Let \(S(N)\) be the required perimeter sum, and let \(M=1234567891\) be the modulus used at the end. Step 1: First family from Euclid's parametrization For a primitive right triangle we may write $$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,\qquad m>n\ge 1.$$ Impose the condition that the hypotenuse is one larger than one leg: $$c=b+1.$$ Then $$m^2+n^2=2mn+1 \Longrightarrow (m-n)^2=1,$$ so necessarily \(m=n+1\). Substituting gives $$a=2n+1,\qquad b=2n(n+1),\qquad c=2n(n+1)+1.$$ If we set \(t=2n+1\), then \(t\) is odd and the perimeter becomes $$p=t+\frac{t^2-1}{2}+\frac{t^2+1}{2}=t(t+1).$$ So Family I contributes exactly the perimeters \(t(t+1)\) with odd \(t\ge 3\). Step 2: Closed-form sum for Family I The inequality \(t(t+1)\le N\) implies \(t<\sqrt{N}\)....

Detailed mathematical approach

Problem Summary

The contributing triangles split into two infinite Pythagorean families. One family has hypotenuse exactly one larger than a leg, and the other has two consecutive legs. The task is to sum every distinct perimeter not exceeding \(N=10^{10^{10}}\) and return the result modulo \(1234567891\). Because \(N\) is unimaginably large, the solution replaces any geometric search with closed formulas, a Pell-type recurrence, and one overlap correction.

Mathematical Approach

Let \(S(N)\) be the required perimeter sum, and let \(M=1234567891\) be the modulus used at the end.

Step 1: First family from Euclid's parametrization

For a primitive right triangle we may write

$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,\qquad m>n\ge 1.$$

Impose the condition that the hypotenuse is one larger than one leg:

$$c=b+1.$$

Then

$$m^2+n^2=2mn+1 \Longrightarrow (m-n)^2=1,$$

so necessarily \(m=n+1\). Substituting gives

$$a=2n+1,\qquad b=2n(n+1),\qquad c=2n(n+1)+1.$$

If we set \(t=2n+1\), then \(t\) is odd and the perimeter becomes

$$p=t+\frac{t^2-1}{2}+\frac{t^2+1}{2}=t(t+1).$$

So Family I contributes exactly the perimeters \(t(t+1)\) with odd \(t\ge 3\).

Step 2: Closed-form sum for Family I

The inequality \(t(t+1)\le N\) implies \(t<\sqrt{N}\). Here \(N=10^E\) with \(E=10^{10}\), and \(E\) is even, so

$$\sqrt{N}=10^{E/2}.$$

Write

$$s=10^{E/2},\qquad q=\frac{s}{2}.$$

The odd values below \(s\) are \(1,3,\dots,s-1\). Using \(t=2j-1\),

$$\sum_{j=1}^{q}(2j-1)(2j)=\sum_{j=1}^{q}(4j^2-2j)=\frac{q(q+1)(4q-1)}{3}.$$

The term \(t=1\) corresponds to the degenerate triple \((1,0,1)\) with perimeter \(2\), so the valid Family I contribution is

$$S_I(N)=\frac{q(q+1)(4q-1)}{3}-2.$$

Step 3: Second family from consecutive legs

Now look at right triangles whose legs are consecutive, say \(x\) and \(x+1\), with hypotenuse \(c\). Then

$$x^2+(x+1)^2=c^2.$$

Introduce

$$u=2x+1.$$

Since \(u^2=4x^2+4x+1\), the Pythagorean condition becomes the negative Pell equation

$$u^2-2c^2=-1.$$

Its positive solutions are generated by multiplication with \(3+2\sqrt{2}\). Starting from \(7+5\sqrt{2}\), which corresponds to the triangle \((3,4,5)\), we get

$$u_{k+1}=3u_k+4c_k,\qquad c_{k+1}=2u_k+3c_k.$$

The perimeter of such a triangle is

$$p_k=u_k+c_k,$$

so the first values are

$$12,\ 70,\ 408,\ 2378,\dots$$

Step 4: Linear recurrence for the perimeter sequence

The Pell multiplier has conjugate roots

$$\lambda=3+2\sqrt{2},\qquad \mu=3-2\sqrt{2}.$$

They satisfy \(\lambda+\mu=6\) and \(\lambda\mu=1\), so any sequence built from these Pell solutions satisfies

$$r^2-6r+1=0.$$

In particular, the perimeter sequence obeys

$$p_{k+1}=6p_k-p_{k-1},\qquad p_1=12,\qquad p_2=70.$$

If \(K\) is the largest index with \(p_K\le N\), then Family II contributes

$$S_{II}(N)=\sum_{k=1}^{K}p_k.$$

Step 5: Find the cutoff index without iterating to \(N\)

The recurrence has closed form

$$p_k=\alpha \lambda^k+\beta \mu^k,$$

with \(|\mu|<1\). Therefore \(p_k\) grows like \(\alpha \lambda^k\), and an accurate first estimate is

$$K\approx \frac{\log_{10}N-\log_{10}\alpha}{\log_{10}\lambda}.$$

After that estimate is computed, a few monotone comparisons are enough to adjust \(K\) until \(p_K\le N<p_{K+1}\). This avoids building huge integers such as \(10^{10^{10}}\) explicitly.

Step 6: Combine both families and remove the overlap

The only perimeter counted twice is \(12\), coming from the triangle \((3,4,5)\). In Family I the sides are

$$t,\qquad \frac{t^2-1}{2},\qquad \frac{t^2+1}{2}.$$

To also have consecutive legs we need

$$\frac{t^2-1}{2}=t+1 \Longrightarrow t^2-2t-3=0 \Longrightarrow t=3.$$

Hence the final formula is

$$S(N)=S_I(N)+S_{II}(N)-12 \pmod{1234567891}.$$

Worked Example: \(N=100\)

Family I contributes the perimeters \(12,30,56,90\), whose sum is \(188\).

Family II contributes \(12,70\), whose sum is \(82\).

Subtract the shared perimeter \(12\) once:

$$S(100)=188+82-12=258.$$

This matches the checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations evaluate Family I entirely modulo \(1234567891\). They compute \(10^{E/2}\bmod M\), replace division by \(2\) and \(3\) with modular inverses, and apply the closed formula for \(S_I(N)\) directly.

For Family II, the implementation first estimates the largest admissible index from the dominant root \(3+2\sqrt{2}\), then corrects the index with a few monotone checks. The sum of the recurrence is obtained with fast exponentiation of a \(3\times 3\) transition matrix whose state stores the current perimeter, the previous perimeter, and the running total.

Finally the two partial sums are added modulo \(1234567891\), and the overlap \(12\) is removed exactly once.

Complexity Analysis

Family I is handled in \(O(1)\) modular arithmetic. Family II uses \(O(\log K)\) multiplications of fixed \(3\times 3\) matrices, plus constant-time index correction. Memory usage is \(O(1)\). The running time depends only on the recurrence index, not on the astronomical size of \(N\).

Footnotes and References

  1. Problem page: Project Euler 835
  2. Pythagorean triples: Wikipedia - Pythagorean triple
  3. Pell-type equations: Wikipedia - Pell's equation
  4. Recurrence relations: Wikipedia - Recurrence relation
  5. Matrix exponentiation for linear recurrences: cp-algorithms - Fibonacci numbers and matrix exponentiation

Problem 835 source code

C++

#include <cassert>
#include <cstdint>
#include <iostream>

#include <boost/multiprecision/cpp_dec_float.hpp>

using u64 = std::uint64_t;
using u128 = unsigned __int128;
using Real = boost::multiprecision::cpp_dec_float_100;

static constexpr u64 kMod = 1'234'567'891ULL;

static u64 mod_mul(u64 a, u64 b) {
    return static_cast<u64>((static_cast<u128>(a) * b) % kMod);
}

static u64 mod_pow(u64 a, u64 e, u64 mod) {
    u64 r = 1 % mod;
    a %= mod;
    while (e > 0) {
        if (e & 1ULL) r = static_cast<u64>((static_cast<u128>(r) * a) % mod);
        a = static_cast<u64>((static_cast<u128>(a) * a) % mod);
        e >>= 1ULL;
    }
    return r;
}

static u64 sum_family_I_mod(u64 half_exp) {
    // N = 10^(2*half_exp), s = 10^half_exp, q = s/2
    const u64 s_mod = mod_pow(10, half_exp, kMod);
    const u64 inv2 = (kMod + 1) / 2;
    const u64 inv3 = mod_pow(3, kMod - 2, kMod);
    const u64 q = mod_mul(s_mod, inv2);

    // Sum over odd t = 1..s-1: t(t+1) = q(q+1)(4q-1)/3
    u64 part = q;
    part = mod_mul(part, (q + 1) % kMod);
    u64 four_q_minus_one = (mod_mul(4 % kMod, q) + kMod - 1) % kMod;
    part = mod_mul(part, four_q_minus_one);
    part = mod_mul(part, inv3);

    // remove t=1 term (degenerate triangle)
    part = (part + kMod - 2) % kMod;
    return part;
}

struct Mat3 {
    u64 a[3][3]{};
};

static Mat3 mat_mul(const Mat3& A, const Mat3& B) {
    Mat3 C{};
    for (int i = 0; i < 3; ++i) {
        for (int k = 0; k < 3; ++k) {
            if (A.a[i][k] == 0) continue;
            for (int j = 0; j < 3; ++j) {
                if (B.a[k][j] == 0) continue;
                C.a[i][j] = (C.a[i][j] + static_cast<u64>((static_cast<u128>(A.a[i][k]) * B.a[k][j]) % kMod)) % kMod;
            }
        }
    }
    return C;
}

static Mat3 mat_pow(Mat3 base, u64 exp) {
    Mat3 res{};
    for (int i = 0; i < 3; ++i) res.a[i][i] = 1;
    while (exp > 0) {
        if (exp & 1ULL) res = mat_mul(res, base);
        base = mat_mul(base, base);
        exp >>= 1ULL;
    }
    return res;
}

static u64 sum_family_II_mod(u64 K) {
    if (K == 0) return 0;
    if (K == 1) return 12;

    // p_{k+1} = 6 p_k - p_{k-1}
    // s_{k+1} = s_k + p_{k+1}
    Mat3 A{};
    A.a[0][0] = 6;
    A.a[0][1] = kMod - 1;
    A.a[0][2] = 0;

    A.a[1][0] = 1;
    A.a[1][1] = 0;
    A.a[1][2] = 0;

    A.a[2][0] = 6;
    A.a[2][1] = kMod - 1;
    A.a[2][2] = 1;

    Mat3 P = mat_pow(A, K - 2);

    // state at k=2: [p2, p1, s2] = [70,12,82]
    const u64 v[3] = {70, 12, 82};
    u64 out[3] = {0, 0, 0};

    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 3; ++j) {
            out[i] = (out[i] + static_cast<u64>((static_cast<u128>(P.a[i][j]) * v[j]) % kMod)) % kMod;
        }
    }

    return out[2];
}

static u64 K_max_for_power10_N(u64 exp10) {
    // p_k = alpha * lambda^k + beta * mu^k, with lambda = 3 + 2*sqrt(2)
    const Real sqrt2 = sqrt(Real(2));
    const Real lambda = Real(3) + Real(2) * sqrt2;
    const Real mu = Real(3) - Real(2) * sqrt2;
    const Real p1 = 12;
    const Real p2 = 70;

    const Real alpha = (p2 - p1 * mu) / (lambda * (lambda - mu));

    const Real log10_lambda = log(lambda) / log(Real(10));
    const Real log10_alpha = log(alpha) / log(Real(10));

    const Real E = Real(exp10);
    Real kest = (E - log10_alpha) / log10_lambda;
    u64 K = kest.convert_to<u64>();

    auto log10_pk = [&](u64 k) -> Real {
        return Real(k) * log10_lambda + log10_alpha;
    };

    while (log10_pk(K + 1) <= E) ++K;
    while (K > 0 && log10_pk(K) > E) --K;

    return K;
}

static u64 S_bruteforce(u64 N) {
    u64 sum = 0;

    // Family I: (t, (t^2-1)/2, (t^2+1)/2), t odd >=3, perimeter t(t+1)
    for (u64 t = 3;; t += 2) {
        u128 p = static_cast<u128>(t) * (t + 1);
        if (p > N) break;
        sum += static_cast<u64>(p);
    }

    // Family II: consecutive legs sequence perimeters
    u64 u = 7, c = 5;  // first valid gives perimeter 12
    while (true) {
        u64 p = u + c;
        if (p > N) break;
        sum += p;
        u64 nu = 3 * u + 4 * c;
        u64 nc = 2 * u + 3 * c;
        u = nu;
        c = nc;
    }

    // overlap at 3-4-5 only
    sum -= 12;
    return sum;
}

int main() {
    assert(S_bruteforce(100) == 258);
    assert(S_bruteforce(10000) == 172004);

    const u64 exp10 = 10'000'000'000ULL;
    const u64 K = K_max_for_power10_N(exp10);

    const u64 sumI = sum_family_I_mod(exp10 / 2);
    const u64 sumII = sum_family_II_mod(K);

    u64 ans = (sumI + sumII) % kMod;
    ans = (ans + kMod - 12) % kMod;

    std::cout << ans << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 1234567891; exp10 = 10000000000

    def mod_pow(a, e, m=MOD):
        r = 1; a %= m
        while e > 0:
            if e & 1: r = r*a%m
            a = a*a%m; e >>= 1
        return r

    # Family I
    half = exp10 // 2
    s = mod_pow(10, half); inv2 = (MOD+1)//2; inv3 = mod_pow(3, MOD-2)
    q = s * inv2 % MOD
    part = q * ((q+1)%MOD) % MOD * ((4*q%MOD + MOD - 1)%MOD) % MOD * inv3 % MOD
    part = (part + MOD - 2) % MOD
    sumI = part

    # Family II: p_k = 6*p_{k-1} - p_{k-2}, matrix exponentiation
    def mat_mul(A, B):
        C = [[0]*3 for _ in range(3)]
        for i in range(3):
            for k in range(3):
                if A[i][k] == 0: continue
                for j in range(3): C[i][j] = (C[i][j] + A[i][k]*B[k][j]) % MOD
        return C
    def mat_pow(base, exp):
        res = [[1 if i==j else 0 for j in range(3)] for i in range(3)]
        while exp > 0:
            if exp & 1: res = mat_mul(res, base)
            base = mat_mul(base, base); exp >>= 1
        return res

    # Find K_max: p_k ~ alpha * lambda^k
    sqrt2 = math.sqrt(2); lam = 3 + 2*sqrt2
    p1, p2 = 12, 70
    alpha = (p2 - p1*(3-2*sqrt2)) / (lam*(lam-(3-2*sqrt2)))
    log10_lam = math.log10(lam); log10_alpha = math.log10(alpha)
    K = int((exp10 - log10_alpha) / log10_lam)
    while (K+1)*log10_lam + log10_alpha <= exp10: K += 1
    while K > 0 and K*log10_lam + log10_alpha > exp10: K -= 1

    if K <= 0: sumII = 0
    elif K == 1: sumII = 12
    else:
        A = [[6, MOD-1, 0], [1, 0, 0], [6, MOD-1, 1]]
        P = mat_pow(A, K-2); v = [70, 12, 82]
        sumII = sum(P[2][j]*v[j] for j in range(3)) % MOD

    ans = (sumI + sumII + MOD - 12) % MOD
    return str(ans)

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

Java

import java.math.BigDecimal;
import java.math.MathContext;
import java.math.RoundingMode;

public class Euler835 {

    static final long kMod = 1234567891L;

    static long modMul(long a, long b) {
        return (a * b) % kMod;
    }

    static long modPow(long base, long exp) {
        long r = 1 % kMod;
        base %= kMod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                r = modMul(r, base);
            base = modMul(base, base);
            exp >>= 1;
        }
        return r;
    }

    static long sumFamilyIMod(long halfExp) {
        long sMod = modPow(10, halfExp);
        long inv2 = (kMod + 1) / 2;
        long inv3 = modPow(3, kMod - 2);
        long q = modMul(sMod, inv2);

        long part = q;
        part = modMul(part, (q + 1) % kMod);
        long fourQMinusOne = (modMul(4, q) + kMod - 1) % kMod;
        part = modMul(part, fourQMinusOne);
        part = modMul(part, inv3);

        part = (part + kMod - 2) % kMod;
        return part;
    }

    static long[][] matMul(long[][] A, long[][] B) {
        long[][] C = new long[3][3];
        for (int i = 0; i < 3; ++i) {
            for (int k = 0; k < 3; ++k) {
                if (A[i][k] == 0)
                    continue;
                for (int j = 0; j < 3; ++j) {
                    if (B[k][j] == 0)
                        continue;
                    C[i][j] = (C[i][j] + modMul(A[i][k], B[k][j])) % kMod;
                }
            }
        }
        return C;
    }

    static long[][] matPow(long[][] base, long exp) {
        long[][] res = new long[3][3];
        for (int i = 0; i < 3; ++i)
            res[i][i] = 1;
        while (exp > 0) {
            if ((exp & 1) == 1)
                res = matMul(res, base);
            base = matMul(base, base);
            exp >>= 1;
        }
        return res;
    }

    static long sumFamilyIIMod(long K) {
        if (K == 0)
            return 0;
        if (K == 1)
            return 12;

        long[][] A = new long[3][3];
        A[0][0] = 6;
        A[0][1] = kMod - 1;
        A[0][2] = 0;
        A[1][0] = 1;
        A[1][1] = 0;
        A[1][2] = 0;
        A[2][0] = 6;
        A[2][1] = kMod - 1;
        A[2][2] = 1;

        long[][] P = matPow(A, K - 2);
        long[] v = { 70, 12, 82 };
        long[] out = { 0, 0, 0 };

        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 3; ++j) {
                out[i] = (out[i] + modMul(P[i][j], v[j])) % kMod;
            }
        }

        return out[2];
    }

    static BigDecimal sqrt(BigDecimal A, final int SCALE) {
        BigDecimal x0 = new BigDecimal("0");
        BigDecimal x1 = new BigDecimal(Math.sqrt(A.doubleValue()));
        while (!x0.equals(x1)) {
            x0 = x1;
            x1 = A.divide(x0, SCALE, RoundingMode.HALF_UP).add(x0)
                    .divide(new BigDecimal("2"), SCALE, RoundingMode.HALF_UP);
        }
        return x1;
    }

    // Since we only need logarithms we can compute log10 directly or just use
    // double for the integer part.
    // E = 10^10 so log10 calculations need enough precision to not make an
    // off-by-one error.
    // ~100 decimal places is enough. Let's use BigInteger/BigDecimal for high
    // precision root finding.
    // Note: K \approx 10^10 / log10(lambda)
    // Actually, we can use a known approximation or BigDecimal.
    // Java doesn't have Math.log10(BigDecimal). We can use a Taylor series or
    // Newton's method, or simply compute log(x) by using high precision double
    // since 10^10 isn't that large. Wait, double has 53 bits (15 digits). We need
    // more than 10 digits of log10(lambda) to get K exact. We need at least 15
    // digits. Double is borderline.
    // We can use BigDecimal for high precision constants.
    // lambda = 3 + 2 * sqrt(2) =
    // 5.828427124746190097603377448419396157139343750753896...
    // log10(lambda) = 0.765551370675727192892973719114704381335...
    // alpha = (70 - 12(3-2sqrt2)) / (8sqrt2) = (34 + 24sqrt2) / (8sqrt2) = 3 +
    // (17/4)sqrt2 = 3 + 4.25 * sqrt(2) = 9.01040764008565...
    // log10(alpha) = 0.954743209598287515...

    static long kMaxForPower10N(long exp10) {
        MathContext mc = new MathContext(100, RoundingMode.HALF_UP);
        BigDecimal sqrt2 = sqrt(new BigDecimal("2"), 100);
        BigDecimal lambda = new BigDecimal("3").add(new BigDecimal("2").multiply(sqrt2));
        BigDecimal mu = new BigDecimal("3").subtract(new BigDecimal("2").multiply(sqrt2));
        BigDecimal p1 = new BigDecimal("12");
        BigDecimal p2 = new BigDecimal("70");

        BigDecimal alpha = p2.subtract(p1.multiply(mu)).divide(lambda.multiply(lambda.subtract(mu)), mc);

        // We need to compute log10(lambda) and log10(alpha) to high precision.
        // We can cheat by using precalculated high precision values.
        BigDecimal log10Lambda = new BigDecimal("0.7655513706757271928929737191147043813350284560759086");
        BigDecimal log10Alpha = new BigDecimal("0.9547432095982875151520625345719335967272288304245084");

        BigDecimal E = new BigDecimal(exp10);
        BigDecimal kest = E.subtract(log10Alpha).divide(log10Lambda, mc);
        long K = kest.longValue();

        // evaluate log10_pk(K)
        BigDecimal log10pkK = new BigDecimal(K).multiply(log10Lambda).add(log10Alpha);
        BigDecimal log10pkKp1 = new BigDecimal(K + 1).multiply(log10Lambda).add(log10Alpha);

        while (log10pkKp1.compareTo(E) <= 0) {
            K++;
            log10pkKp1 = new BigDecimal(K + 1).multiply(log10Lambda).add(log10Alpha);
        }
        while (K > 0) {
            log10pkK = new BigDecimal(K).multiply(log10Lambda).add(log10Alpha);
            if (log10pkK.compareTo(E) <= 0)
                break;
            K--;
        }

        return K;
    }

    public static String solve() {
        long exp10 = 10000000000L;
        long K = kMaxForPower10N(exp10);

        long sumI = sumFamilyIMod(exp10 / 2);
        long sumII = sumFamilyIIMod(K);

        long ans = (sumI + sumII) % kMod;
        ans = (ans + kMod - 12) % kMod;

        return Long.toString(ans);
    }

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