Problem 785: Symmetric Diophantine Equation

Project Euler Problem 785 Solution

EulerSolve provides an optimized solution for Project Euler Problem 785, Symmetric Diophantine Equation, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a bound $N$, define $$S(N)=\sum_{\substack{1\le x\le y\le z\le N\\ \gcd(x,y,z)=1\\ 15(x^2+y^2+z^2)=34(xy+yz+zx)}}(x+y+z).$$ The task is to evaluate $S(10^9)$. The C++, Python, and Java implementations do not search directly in $(x,y,z)$-space. Instead, they generate primitive solutions from a quadratic parametrization, divide out a small class factor, and keep only the parameter pairs that can satisfy the ordering and divisibility conditions. Mathematical Approach The key observation is that the ordered primitive triples on this quadratic surface can be generated from coprime pairs $(u,v)$ lying in a narrow cone. The raw quadratic forms are $$x_0=15u^2-34uv+15v^2,\qquad y_0=105u^2-446uv+473v^2,\qquad z_0=32u^2-176uv+240v^2.$$ After dividing by a class factor $g\in\{1,4,16,19,64,76,304,1216\}$, the counted triple is $$x=\frac{x_0}{g},\qquad y=\frac{y_0}{g},\qquad z=\frac{z_0}{g}.$$ Step 1: Parametrize the Solution Surface The three raw forms are chosen so that substituting them into the symmetric Diophantine equation gives an identity: $$15(x_0^2+y_0^2+z_0^2)=34(x_0y_0+y_0z_0+z_0x_0).$$ Dividing every coordinate by the same positive factor preserves the equation, so the same relation holds for $(x,y,z)$. The primitive condition is enforced by combining $\gcd(u,v)=1$ with the systematic removal of the shared class factor $g$....

Detailed mathematical approach

Problem Summary

For a bound $N$, define

$$S(N)=\sum_{\substack{1\le x\le y\le z\le N\\ \gcd(x,y,z)=1\\ 15(x^2+y^2+z^2)=34(xy+yz+zx)}}(x+y+z).$$

The task is to evaluate $S(10^9)$. The C++, Python, and Java implementations do not search directly in $(x,y,z)$-space. Instead, they generate primitive solutions from a quadratic parametrization, divide out a small class factor, and keep only the parameter pairs that can satisfy the ordering and divisibility conditions.

Mathematical Approach

The key observation is that the ordered primitive triples on this quadratic surface can be generated from coprime pairs $(u,v)$ lying in a narrow cone. The raw quadratic forms are

$$x_0=15u^2-34uv+15v^2,\qquad y_0=105u^2-446uv+473v^2,\qquad z_0=32u^2-176uv+240v^2.$$

After dividing by a class factor $g\in\{1,4,16,19,64,76,304,1216\}$, the counted triple is

$$x=\frac{x_0}{g},\qquad y=\frac{y_0}{g},\qquad z=\frac{z_0}{g}.$$

Step 1: Parametrize the Solution Surface

The three raw forms are chosen so that substituting them into the symmetric Diophantine equation gives an identity:

$$15(x_0^2+y_0^2+z_0^2)=34(x_0y_0+y_0z_0+z_0x_0).$$

Dividing every coordinate by the same positive factor preserves the equation, so the same relation holds for $(x,y,z)$. The primitive condition is enforced by combining $\gcd(u,v)=1$ with the systematic removal of the shared class factor $g$.

Step 2: Restrict to the Branch with $x\le y\le z$

Write $t=u/v$. Then

$$x_0=v^2(15t^2-34t+15)=v^2(5t-3)(3t-5).$$

To obtain positive solutions on the relevant branch we need $t>5/3$. Next,

$$y_0-x_0=2v^2(45t^2-206t+229).$$

The smaller root of $45t^2-206t+229=0$ is

$$h=\frac{103-4\sqrt{19}}{45}.$$

Therefore $x_0\le y_0$ on the ordered branch when

$$\frac{5}{3}<t\le h.$$

On this same interval one also has $y_0\le z_0$, so the search only needs the cone

$$\frac{5}{3}v<u\le hv.$$

Step 3: Extract the Common Class Factor

Introduce the linear forms

$$a=5u-3v,\qquad b=3u-5v.$$

Then the raw coordinates can be rewritten as

$$x_0=ab,$$

$$y_0=\frac{(3a-19b)(a-5b)}{4},$$

$$z_0=\frac{(5a-19b)(a-3b)}{4}.$$

These identities explain why only eight class factors occur. The $2$-power part is

$$g_2=4^{\min(\nu_2(a),\nu_2(b),3)},$$

so $g_2\in\{1,4,16,64\}$. In addition, if $19\mid a$, then all three raw coordinates are divisible by $19$. Hence the full factor is one of

$$g\in\{1,4,16,64\}\times\{1,19\}=\{1,4,16,19,64,76,304,1216\}.$$

Step 4: Turn Divisibility into a Residue Table Modulo $304$

The class factor depends only on

$$a\bmod 16,\qquad b\bmod 16,\qquad a\bmod 19.$$

That is enough because the $2$-power part only needs the powers of $2$ in $a$ and $b$ up to $2^3$, and the $19$-part only asks whether $a\equiv 0\pmod{19}$. Combining modulus $16$ and modulus $19$ gives

$$16\cdot 19=304.$$

So for each class factor $g$ and each residue of $v\bmod 304$, the implementations precompute the admissible residues of $u\bmod 304$. The main search can then skip every parameter pair whose residue class cannot possibly land in the desired primitive branch.

Step 5: Bound the Search by $z\le N$

Still writing $t=u/v$, we have

$$z_0=v^2(32t^2-176t+240).$$

On the interval $\frac{5}{3}<t\le h$, this quadratic factor is minimized at the right endpoint, so

$$z_0\ge c_{\min}v^2,\qquad c_{\min}=32h^2-176h+240.$$

Because the counted triple satisfies $z=z_0/g\le N$, every class obeys

$$v\le \sqrt{\frac{Ng}{c_{\min}}}.$$

This square-root bound is exactly what determines the per-class outer loop limits.

Worked Example

Take $(u,v)=(13,7)$. Then

$$\frac{u}{v}=\frac{13}{7}\approx 1.857,$$

which lies inside $\frac{5}{3}<u/v\le h$. Now

$$a=5u-3v=44,\qquad b=3u-5v=4.$$

Both $a$ and $b$ are divisible by $4$, but not both by $8$, so $g_2=16$. Since $19\nmid 44$, there is no extra factor $19$, hence $g=16$.

The raw forms are

$$x_0=176,\qquad y_0=336,\qquad z_0=1152,$$

and dividing by $16$ gives

$$ (x,y,z)=\left(11,21,72\right). $$

This triple is primitive and ordered, and it satisfies

$$15(11^2+21^2+72^2)=34(11\cdot21+21\cdot72+72\cdot11).$$

Its contribution to the sum is therefore

$$11+21+72=104.$$

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they precompute, for each class factor and each residue of $v\bmod 304$, the admissible residues of $u\bmod 304$. Next they determine class-dependent upper bounds for $v$ from the inequality $z\le N$. For every $v$, they scan only the integers $u$ in the cone $\frac{5}{3}v<u\le hv$ that match one of the allowed residue classes. They then apply three final filters: $\gcd(u,v)=1$, the bound $z_0\le Ng$, and division by the class factor $g$. Every surviving candidate contributes

$$\frac{x_0+y_0+z_0}{g}$$

to the running total. The C++ and Java implementations optionally split the $v$-range across several threads for large $N$; the Python implementation uses the same arithmetic in a single thread.

Complexity Analysis

The residue precomputation is constant-size, because it only fills tables for $8$ classes and $304$ residues. For a fixed class, the outer variable $v$ runs to $O(\sqrt{N})$, while the raw width of the cone in $u$ is $O(v)$. Thus the unfiltered search region has $O(N)$ area. The residue table removes most candidates before any gcd or bound test is attempted, so the practical constant factor is much smaller than a direct cone scan. Memory usage is $O(1)$ apart from the small residue tables and, in the threaded implementations, one accumulator per worker.

Footnotes and References

Problem page: https://projecteuler.net/problem=785
Diophantine equation: Wikipedia — Diophantine equation
Binary quadratic form: Wikipedia — Binary quadratic form
Modular arithmetic: Wikipedia — Modular arithmetic
Greatest common divisor: Wikipedia — Greatest common divisor

Problem 785 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <pthread.h>
#include <string>
#include <unistd.h>
#include <vector>

using i64 = long long;
using i128 = __int128_t;
using ResidueTable = std::array<std::array<std::vector<int>, 304>, 8>;

static std::string to_string_i128(i128 x) {
    if (x == 0) {
        return "0";
    }
    bool neg = x < 0;
    if (neg) {
        x = -x;
    }
    std::string s;
    while (x > 0) {
        int d = static_cast<int>(x % 10);
        s.push_back(static_cast<char>('0' + d));
        x /= 10;
    }
    if (neg) {
        s.push_back('-');
    }
    std::reverse(s.begin(), s.end());
    return s;
}

static int class_index_from_g(int g) {
    switch (g) {
        case 1: return 0;
        case 4: return 1;
        case 16: return 2;
        case 19: return 3;
        case 64: return 4;
        case 76: return 5;
        case 304: return 6;
        case 1216: return 7;
        default: return -1;
    }
}

static int trailing_zeros_mod16(int x_mod16) {
    x_mod16 &= 15;
    if (x_mod16 == 0) {
        return 4;
    }
    int c = 0;
    while ((x_mod16 & 1) == 0) {
        x_mod16 >>= 1;
        ++c;
    }
    return c;
}

static ResidueTable build_residue_table() {
    ResidueTable table;

    for (int vm = 0; vm < 304; ++vm) {
        for (int um = 0; um < 304; ++um) {
            const int m16 = (5 * um - 3 * vm) & 15;
            const int n16 = (3 * um - 5 * vm) & 15;

            const int tz = std::min(std::min(trailing_zeros_mod16(m16), trailing_zeros_mod16(n16)), 3);
            const int g2 = (tz == 0 ? 1 : (tz == 1 ? 4 : (tz == 2 ? 16 : 64)));

            int m19 = (5 * um - 3 * vm) % 19;
            if (m19 < 0) {
                m19 += 19;
            }
            const bool e19 = (m19 == 0);
            const int g = g2 * (e19 ? 19 : 1);

            const int idx = class_index_from_g(g);
            assert(idx >= 0);
            table[idx][vm].push_back(um);
        }
    }

    return table;
}

static const ResidueTable& get_residue_table() {
    static const ResidueTable table = build_residue_table();
    return table;
}

static i128 accumulate_stride(i64 N, const std::array<int, 8>& vmax, int tid, int step) {
    constexpr std::array<int, 8> G = {1, 4, 16, 19, 64, 76, 304, 1216};
    const auto& residue_table = get_residue_table();

    const long double lo = 5.0L / 3.0L;
    const long double hi = (103.0L - 4.0L * std::sqrt(19.0L)) / 45.0L;

    i128 ans = 0;

    for (int ci = 0; ci < 8; ++ci) {
        const int g = G[ci];
        for (int v = tid + 1; v <= vmax[ci]; v += step) {
            const i64 umin = static_cast<i64>(std::floor(lo * static_cast<long double>(v))) + 1;
            const i64 umax = static_cast<i64>(std::floor(hi * static_cast<long double>(v)));
            if (umin > umax) {
                continue;
            }

            const auto& residues = residue_table[ci][v % 304];
            for (int r : residues) {
                i64 first = umin;
                const int delta = (r - static_cast<int>(first % 304) + 304) % 304;
                first += delta;

                for (i64 u = first; u <= umax; u += 304) {
                    if (std::gcd(u, static_cast<i64>(v)) != 1) {
                        continue;
                    }

                    const i64 uu = u;
                    const i64 vv = v;

                    const i64 z0 = 32LL * uu * uu - 176LL * uu * vv + 240LL * vv * vv;
                    if (z0 > N * static_cast<i64>(g)) {
                        continue;
                    }

                    const i64 x0 = 15LL * uu * uu - 34LL * uu * vv + 15LL * vv * vv;
                    const i64 y0 = 105LL * uu * uu - 446LL * uu * vv + 473LL * vv * vv;

                    const i64 x = x0 / g;
                    const i64 y = y0 / g;
                    const i64 z = z0 / g;

                    ans += static_cast<i128>(x) + static_cast<i128>(y) + static_cast<i128>(z);
                }
            }
        }
    }

    return ans;
}

struct ThreadTask785 {
    i64 N = 0;
    std::array<int, 8> vmax{};
    int tid = 0;
    int step = 1;
    i128 partial = 0;
};

static void* thread_worker_785(void* arg) {
    auto* t = static_cast<ThreadTask785*>(arg);
    t->partial = accumulate_stride(t->N, t->vmax, t->tid, t->step);
    return nullptr;
}

static i128 solve_fast(i64 N) {
    constexpr std::array<int, 8> G = {1, 4, 16, 19, 64, 76, 304, 1216};
    const long double hi = (103.0L - 4.0L * std::sqrt(19.0L)) / 45.0L;
    const long double cmin = 32.0L * hi * hi - 176.0L * hi + 240.0L;

    std::array<int, 8> vmax{};
    for (int i = 0; i < 8; ++i) {
        const long double lim = std::sqrt((static_cast<long double>(N) * static_cast<long double>(G[i])) / cmin);
        vmax[i] = static_cast<int>(std::floor(lim)) + 3;
    }

    long cpu_count = ::sysconf(_SC_NPROCESSORS_ONLN);
    int thread_count = (cpu_count > 1 ? static_cast<int>(cpu_count) : 1);
    if (thread_count > 16) {
        thread_count = 16;
    }
    if (N < 1'000'000LL || thread_count <= 1) {
        return accumulate_stride(N, vmax, 0, 1);
    }

    std::vector<pthread_t> threads((std::size_t)thread_count);
    std::vector<ThreadTask785> tasks((std::size_t)thread_count);

    for (int t = 0; t < thread_count; ++t) {
        tasks[(std::size_t)t].N = N;
        tasks[(std::size_t)t].vmax = vmax;
        tasks[(std::size_t)t].tid = t;
        tasks[(std::size_t)t].step = thread_count;
        const int rc = ::pthread_create(&threads[(std::size_t)t], nullptr, thread_worker_785, &tasks[(std::size_t)t]);
        assert(rc == 0);
    }

    i128 ans = 0;
    for (int t = 0; t < thread_count; ++t) {
        const int rc = ::pthread_join(threads[(std::size_t)t], nullptr);
        assert(rc == 0);
        ans += tasks[(std::size_t)t].partial;
    }
    return ans;
}

static i128 solve_bruteforce(int N) {
    i128 ans = 0;
    for (int x = 1; x <= N; ++x) {
        for (int y = x; y <= N; ++y) {
            const i64 A = 15;
            const i64 B = -34LL * (x + y);
            const i64 C = 15LL * (1LL * x * x + 1LL * y * y) - 34LL * x * y;
            const i64 D = B * B - 4LL * A * C;
            const i64 t = static_cast<i64>(std::sqrt(static_cast<long double>(D)));
            if (t * t != D) {
                continue;
            }
            for (i64 sgn : {-1LL, 1LL}) {
                const i64 num = -B + sgn * t;
                if (num % (2LL * A) != 0) {
                    continue;
                }
                const i64 z = num / (2LL * A);
                if (z < y || z > N) {
                    continue;
                }
                const i64 lhs = 15LL * (1LL * x * x + 1LL * y * y + z * z);
                const i64 rhs = 34LL * (1LL * x * y + 1LL * y * z + 1LL * z * x);
                if (lhs != rhs) {
                    continue;
                }
                if (std::gcd(std::gcd(x, y), static_cast<int>(z)) != 1) {
                    continue;
                }
                ans += static_cast<i128>(x + y + z);
            }
        }
    }
    return ans;
}

int main() {
    assert(to_string_i128(solve_fast(100)) == "184");
    assert(solve_fast(200) == solve_bruteforce(200));

    const i128 ans = solve_fast(1'000'000'000LL);
    std::cout << to_string_i128(ans) << '\n';
    return 0;
}

Python

import math

def solve():
    N = 1000000000
    G = [1, 4, 16, 19, 64, 76, 304, 1216]
    lo_r = 5/3
    hi_r = (103 - 4*math.sqrt(19)) / 45

    def cidx(g):
        return {1:0, 4:1, 16:2, 19:3, 64:4, 76:5, 304:6, 1216:7}[g]

    def tz16(x):
        x &= 15
        if x == 0: return 4
        c = 0
        while (x & 1) == 0: x >>= 1; c += 1
        return c

    # Build residue table
    table = [[[] for _ in range(304)] for _ in range(8)]
    for vm in range(304):
        for um in range(304):
            m16 = (5*um - 3*vm) & 15; n16 = (3*um - 5*vm) & 15
            tz = min(tz16(m16), tz16(n16), 3)
            g2 = [1, 4, 16, 64][tz]
            m19 = (5*um - 3*vm) % 19
            e19 = m19 == 0
            g = g2 * (19 if e19 else 1)
            table[cidx(g)][vm].append(um)

    cmin = 32*hi_r*hi_r - 176*hi_r + 240
    ans = 0
    for ci in range(8):
        g = G[ci]; vmax = int(math.sqrt(N * g / cmin)) + 3
        for v in range(1, vmax+1):
            umin = int(math.floor(lo_r * v)) + 1
            umax = int(math.floor(hi_r * v))
            if umin > umax: continue
            for r in table[ci][v % 304]:
                delta = (r - umin % 304 + 304) % 304
                first = umin + delta
                for u in range(first, umax+1, 304):
                    if math.gcd(u, v) != 1: continue
                    z0 = 32*u*u - 176*u*v + 240*v*v
                    if z0 > N * g: continue
                    x0 = 15*u*u - 34*u*v + 15*v*v
                    y0 = 105*u*u - 446*u*v + 473*v*v
                    ans += x0//g + y0//g + z0//g

    return str(ans)

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

Java

import java.util.ArrayList;

public class Euler785 {

    static int classIndexFromG(int g) {
        switch (g) {
            case 1:
                return 0;
            case 4:
                return 1;
            case 16:
                return 2;
            case 19:
                return 3;
            case 64:
                return 4;
            case 76:
                return 5;
            case 304:
                return 6;
            case 1216:
                return 7;
            default:
                return -1;
        }
    }

    static int trailingZerosMod16(int xMod16) {
        xMod16 &= 15;
        if (xMod16 == 0)
            return 4;
        int c = 0;
        while ((xMod16 & 1) == 0) {
            xMod16 >>= 1;
            c++;
        }
        return c;
    }

    @SuppressWarnings("unchecked")
    static ArrayList<Integer>[][] buildResidueTable() {
        ArrayList<Integer>[][] table = new ArrayList[8][304];
        for (int i = 0; i < 8; i++) {
            for (int j = 0; j < 304; j++) {
                table[i][j] = new ArrayList<>();
            }
        }

        for (int vm = 0; vm < 304; ++vm) {
            for (int um = 0; um < 304; ++um) {
                int m16 = (5 * um - 3 * vm) & 15;
                int n16 = (3 * um - 5 * vm) & 15;

                int tz = Math.min(trailingZerosMod16(m16), trailingZerosMod16(n16));
                tz = Math.min(tz, 3);
                int g2 = (tz == 0 ? 1 : (tz == 1 ? 4 : (tz == 2 ? 16 : 64)));

                int m19 = (5 * um - 3 * vm) % 19;
                if (m19 < 0)
                    m19 += 19;
                boolean e19 = (m19 == 0);
                int g = g2 * (e19 ? 19 : 1);

                int idx = classIndexFromG(g);
                table[idx][vm].add(um);
            }
        }
        return table;
    }

    static final ArrayList<Integer>[][] residueTable = buildResidueTable();

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = a % b;
            a = b;
            b = t;
        }
        return a;
    }

    static long accumulateStride(long N, int[] vmax, int tid, int step) {
        int[] G = { 1, 4, 16, 19, 64, 76, 304, 1216 };
        double lo = 5.0 / 3.0;
        double hi = (103.0 - 4.0 * Math.sqrt(19.0)) / 45.0;

        long ans = 0; // The total answer fits in a 64 bit integer for N=10^9

        for (int ci = 0; ci < 8; ++ci) {
            int g = G[ci];
            for (int v = tid + 1; v <= vmax[ci]; v += step) {
                long umin = (long) Math.floor(lo * v) + 1;
                long umax = (long) Math.floor(hi * v);
                if (umin > umax)
                    continue;

                ArrayList<Integer> residues = residueTable[ci][v % 304];
                for (int r : residues) {
                    long first = umin;
                    int delta = (r - (int) (first % 304) + 304) % 304;
                    first += delta;

                    for (long u = first; u <= umax; u += 304) {
                        if (gcd(u, v) != 1)
                            continue;

                        long uu = u;
                        long vv = v;

                        long z0 = 32L * uu * uu - 176L * uu * vv + 240L * vv * vv;
                        if (z0 > N * (long) g)
                            continue;

                        long x0 = 15L * uu * uu - 34L * uu * vv + 15L * vv * vv;
                        long y0 = 105L * uu * uu - 446L * uu * vv + 473L * vv * vv;

                        long x = x0 / g;
                        long y = y0 / g;
                        long z = z0 / g;

                        ans += x + y + z;
                    }
                }
            }
        }

        return ans;
    }

    static long solveFast(long N) {
        int[] G = { 1, 4, 16, 19, 64, 76, 304, 1216 };
        double hi = (103.0 - 4.0 * Math.sqrt(19.0)) / 45.0;
        double cmin = 32.0 * hi * hi - 176.0 * hi + 240.0;

        int[] vmax = new int[8];
        for (int i = 0; i < 8; ++i) {
            double lim = Math.sqrt((N * (double) G[i]) / cmin);
            vmax[i] = (int) Math.floor(lim) + 3;
        }

        int threadCount = Runtime.getRuntime().availableProcessors();
        if (threadCount > 16)
            threadCount = 16;
        if (threadCount <= 1 || N < 1000000) {
            return accumulateStride(N, vmax, 0, 1);
        }

        Thread[] threads = new Thread[threadCount];
        long[] partials = new long[threadCount];

        for (int t = 0; t < threadCount; t++) {
            final int tid = t;
            final int step = threadCount;
            threads[t] = new Thread(() -> {
                partials[tid] = accumulateStride(N, vmax, tid, step);
            });
            threads[t].start();
        }

        long ans = 0;
        for (int t = 0; t < threadCount; t++) {
            try {
                threads[t].join();
                ans += partials[t];
            } catch (InterruptedException e) {
                e.printStackTrace();
            }
        }
        return ans;
    }

    public static String solve() {
        return Long.toString(solveFast(1000000000L));
    }

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