Problem 585: Nested Square Roots

View on Project Euler

Project Euler Problem 585 Solution

EulerSolve provides an optimized solution for Project Euler Problem 585, Nested Square Roots, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a given bound \(n\), we count distinct real values of $$\sqrt{x+\sqrt y+\sqrt z},\qquad 0 \lt x \le n,$$ where \(y\) and \(z\) are non-squares and the outer radical can be denested into a finite signed sum of simple square roots. Different triples \((x,y,z)\) may represent the same real number, but that value must be counted only once. The implementations evaluate this count at \(n=5{,}000{,}000\). The program organizes the denestable values into two structural families. The first family comes from denestings with two simple roots. The second family comes from denestings with four simple roots and one minus sign. The final formula counts both families and then removes the four-root cases that collapse back to a value already counted in the two-root family. Mathematical Approach Let \(F(n)\) denote the number of distinct denestable values with integer part \(x\le n\). The implementation builds \(F(n)\) from an arithmetic kernel \(f(s)\) and three counting terms \(C_1(n)\), \(C_2^{\mathrm{ord}}(n)\), and \(C_3(n)\)....

Detailed mathematical approach

Problem Summary

For a given bound \(n\), we count distinct real values of

$$\sqrt{x+\sqrt y+\sqrt z},\qquad 0 \lt x \le n,$$

where \(y\) and \(z\) are non-squares and the outer radical can be denested into a finite signed sum of simple square roots. Different triples \((x,y,z)\) may represent the same real number, but that value must be counted only once. The implementations evaluate this count at \(n=5{,}000{,}000\).

The program organizes the denestable values into two structural families. The first family comes from denestings with two simple roots. The second family comes from denestings with four simple roots and one minus sign. The final formula counts both families and then removes the four-root cases that collapse back to a value already counted in the two-root family.

Mathematical Approach

Let \(F(n)\) denote the number of distinct denestable values with integer part \(x\le n\). The implementation builds \(F(n)\) from an arithmetic kernel \(f(s)\) and three counting terms \(C_1(n)\), \(C_2^{\mathrm{ord}}(n)\), and \(C_3(n)\).

Step 1: Primitive Two-Root Denestings

Start with a value of the form

$$V=\sqrt{k\alpha}+\sqrt{k\beta},\qquad \alpha \gt \beta \ge 1,\qquad \gcd(\alpha,\beta)=1,\qquad k\ge 1.$$

Squaring gives

$$V^2=k(\alpha+\beta)+2k\sqrt{\alpha\beta}.$$

If \(\alpha\beta\) is not a square, then \(V^2\) has exactly one irrational square-root direction, and it fits the required pattern because

$$2k\sqrt{\alpha\beta}=\sqrt{k^2\alpha\beta}+\sqrt{k^2\alpha\beta},$$

with \(k^2\alpha\beta\) still non-square. So each primitive pair \((\alpha,\beta)\) produces one infinite scaling family of admissible values, and the integer part is

$$x=k(\alpha+\beta).$$

For a fixed primitive sum

$$s=\alpha+\beta,$$

the scale \(k\) can be chosen in exactly

$$\left\lfloor\frac{n}{s}\right\rfloor$$

ways.

Step 2: Count the Primitive Kernels \(f(s)\)

For fixed \(s\ge 3\), how many primitive pairs \((\alpha,\beta)\) with \(\alpha+\beta=s\) belong to the first family?

The number of coprime decompositions \(s=\alpha+\beta\) with \(\alpha \gt \beta \ge 1\) is

$$\frac{\varphi(s)}{2},$$

because each reduced residue class modulo \(s\) corresponds to one ordered coprime split, and dividing by \(2\) forgets the order.

Now remove the bad cases where \(\alpha\beta\) is a square. Since \(\gcd(\alpha,\beta)=1\), this happens exactly when

$$\alpha=u^2,\qquad \beta=v^2,\qquad \gcd(u,v)=1.$$

So the excluded pairs are precisely the primitive representations

$$s=u^2+v^2,\qquad u \gt v \ge 1,\qquad \gcd(u,v)=1.$$

If \(R(s)\) denotes the number of those primitive two-square representations, then the kernel used by the implementations is

$$f(s)=\frac{\varphi(s)}{2}-R(s).$$

Therefore the total contribution of the two-root family is

$$C_1(n)=\sum_{s\le n} f(s)\left\lfloor\frac{n}{s}\right\rfloor.$$

Step 3: A Genuine Four-Root Family

The second family combines two primitive kernels. Take two primitive pairs \((\alpha,\beta)\) and \((\gamma,\delta)\), each counted by \(f\), and form

$$W=\sqrt{k\alpha\gamma}+\sqrt{k\alpha\delta}+\sqrt{k\beta\gamma}-\sqrt{k\beta\delta}.$$

Expanding \(W^2\) gives a crucial cancellation:

$$W^2=k(\alpha+\beta)(\gamma+\delta)+2k(\alpha-\beta)\sqrt{\gamma\delta}+2k(\gamma-\delta)\sqrt{\alpha\beta}.$$

The mixed term involving \(\sqrt{\alpha\beta\gamma\delta}\) disappears because of the sign pattern. This leaves exactly two square-root directions, so \(W\) is another denesting of the required form.

If we write

$$s=\alpha+\beta,\qquad t=\gamma+\delta,$$

then the integer part of \(W^2\) is

$$x=kst.$$

Hence each ordered pair of primitive kernels \((s,t)\) contributes

$$\left\lfloor\frac{n}{st}\right\rfloor$$

scaled values, and the ordered count is

$$C_2^{\mathrm{ord}}(n)=\sum_{s\le n}\sum_{t\le n} f(s)f(t)\left\lfloor\frac{n}{st}\right\rfloor.$$

The superscript “ord” matters: swapping the two kernels does not change the final value, so this sum counts the generic four-root construction twice.

Step 4: Detect the Four-Root Values that Collapse to Two Roots

Not every ordered four-root construction is genuinely new. Some of them collapse to only two independent squarefree radicals and are therefore already included in \(C_1(n)\).

Write the two primitive pairs as

$$\alpha=\sigma_1 a^2,\qquad \beta=\sigma_2 b^2,\qquad \gamma=\tau_1 c^2,\qquad \delta=\tau_2 d^2,$$

where \(\sigma_1,\sigma_2,\tau_1,\tau_2\) are squarefree. The four radicals inside \(W\) then have squarefree parts

$$\sigma_1\tau_1,\qquad \sigma_1\tau_2,\qquad \sigma_2\tau_1,\qquad \sigma_2\tau_2.$$

A collision with the two-root family occurs exactly when these four squarefree products merge into only two distinct squarefree directions. The implementation parameterizes those collisions by squarefree, pairwise-coprime factors \(\rho_1,\rho_2,\rho_3,\rho_4\) such that

$$\sigma_1=\rho_1\rho_2,\qquad \sigma_2=\rho_3\rho_4,\qquad \tau_1=\rho_1\rho_3,\qquad \tau_2=\rho_2\rho_4.$$

Then

$$\sqrt{\sigma_1\tau_1}=\rho_1\sqrt{\rho_2\rho_3},\qquad \sqrt{\sigma_2\tau_2}=\rho_4\sqrt{\rho_2\rho_3},$$

and likewise

$$\sqrt{\sigma_1\tau_2}=\rho_2\sqrt{\rho_1\rho_4},\qquad \sqrt{\sigma_2\tau_1}=\rho_3\sqrt{\rho_1\rho_4}.$$

So the supposed four-root expression is really a two-root value. Let \(C_3(n)\) be the number of ordered constructions of exactly this collapsing kind. These must be removed from the ordered four-root count before the final symmetry factor is applied.

Step 5: Final Formula

After subtracting the collapsing ordered constructions and then forgetting the order of the two primitive kernels, the implementations compute

$$\boxed{F(n)=C_1(n)+\frac{C_2^{\mathrm{ord}}(n)-C_3(n)}{2}.}$$

This is the exact decomposition used by the C++, Python, and Java implementations.

Worked Example

A simple two-root example is

$$V=\sqrt5+\sqrt2.$$

Then

$$V^2=7+2\sqrt{10}=7+\sqrt{10}+\sqrt{10}.$$

Here the primitive kernel is \(s=5+2=7\), and because \(5\cdot2=10\) is not a square, this value belongs to the first family. Every scale \(k\) with \(7k\le n\) gives another valid value from the same kernel.

A genuine four-root example uses \((\alpha,\beta)=(2,1)\) and \((\gamma,\delta)=(3,1)\):

$$W=\sqrt6+\sqrt2+\sqrt3-1.$$

Its square is

$$W^2=(2+1)(3+1)+2(2-1)\sqrt3+2(3-1)\sqrt2=12+2\sqrt3+4\sqrt2.$$

So the integer part is \(12=3\cdot4\), which shows why the second family is indexed by products of primitive kernel sums.

How the Code Works

The implementations first build a totient sieve up to \(n\). From that array they initialize

$$f(s)=\frac{\varphi(s)}{2}$$

for \(s\ge 3\), then subtract one for every primitive representation \(s=u^2+v^2\) with \(u \gt v\) and \(\gcd(u,v)=1\). That produces the kernel sequence \(f(s)\).

Next they form prefix sums of \(f\). The first family \(C_1(n)\) is then a direct divisor-style sum. The ordered second-family term \(C_2^{\mathrm{ord}}(n)\) is not evaluated by a quadratic double loop. Instead, the code groups equal quotients of \(\left\lfloor n/i\right\rfloor\) into harmonic intervals and combines them with prefix sums, which turns the convolution into a much smaller collection of block sums.

The overlap term \(C_3(n)\) is handled separately. The code precomputes squarefree flags, coprimality lookup tables, and squares up to \(\lfloor\sqrt n\rfloor\). It then enumerates only the parameter combinations that can force a four-root value to collapse into two squarefree directions, using many early inequality breaks to prune the search. The C++ implementation parallelizes this expensive overlap phase; the Python version mirrors the same mathematics more directly, and the Java version reuses the same optimized core computation.

Complexity Analysis

Building the totient sieve and the kernel \(f\) takes \(O(n\log\log n)\) time and \(O(n)\) memory. The first family sum \(C_1(n)\) is \(O(n)\). The ordered convolution \(C_2^{\mathrm{ord}}(n)\) is accelerated by harmonic blocking and prefix sums, so the code avoids the naive \(O(n^2)\) double loop and works with a subquadratic number of quotient blocks instead.

The practical bottleneck is the overlap correction \(C_3(n)\). Its exact worst-case count is awkward to state cleanly because it depends on several nested squarefree and coprimality constraints, but the implementation reduces it sharply through precomputed tables, squarefree filtering, and early stopping inequalities. Overall memory usage remains \(O(n)\), dominated by the arithmetic arrays up to \(n\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=585
  2. Nested radical: Wikipedia - Nested radical
  3. Euler's totient function: Wikipedia - Euler's totient function
  4. Coprime integers: Wikipedia - Coprime integers
  5. Sum of two squares: Wikipedia - Sum of two squares
  6. Square-free integer: Wikipedia - Square-free integer

Problem 585 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <numeric>
#include <thread>
#include <vector>

using int64 = long long;

static std::vector<int> totient_sieve(int n) {
    std::vector<int> phi(n + 1);
    std::iota(phi.begin(), phi.end(), 0);
    for (int i = 2; i <= n; ++i) {
        if (phi[i] == i) {
            for (int j = i; j <= n; j += i) {
                phi[j] -= phi[j] / i;
            }
        }
    }
    return phi;
}

// f(s) = phi(s)/2 - R(s), where R(s) counts primitive i^2 + j^2 representations.
static std::vector<int> build_f(int n) {
    std::vector<int> phi = totient_sieve(n);
    std::vector<int> f(n + 1, 0);
    for (int s = 3; s <= n; ++s) {
        f[s] = phi[s] / 2;
    }

    int lim = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    for (int i = 2; i <= lim; ++i) {
        int i2 = i * i;
        for (int j = 1; j < i; ++j) {
            if (std::gcd(i, j) != 1) continue;
            int s = i2 + j * j;
            if (s > n) break;
            f[s] -= 1;
        }
    }
    return f;
}

static std::vector<int64> prefix_sum(const std::vector<int>& f) {
    std::vector<int64> pref(f.size(), 0);
    for (size_t i = 1; i < f.size(); ++i) {
        pref[i] = pref[i - 1] + f[i];
    }
    return pref;
}

static int64 partial_sum_floor(int64 m, const std::vector<int64>& pref) {
    int64 ret = 0;
    for (int64 i = 1; i <= m;) {
        int64 q = m / i;
        int64 j = m / q + 1;
        ret += q * (pref[j - 1] - pref[i - 1]);
        i = j;
    }
    return ret;
}

static int64 count_type1(int n, const std::vector<int>& f) {
    int64 cnt = 0;
    for (int s = 3; s <= n; ++s) {
        if (f[s] == 0) continue;
        cnt += static_cast<int64>(f[s]) * (n / s);
    }
    return cnt;
}

static int64 count_type2_ordered(int n, const std::vector<int64>& pref) {
    int64 cnt = 0;
    for (int64 i = 1; i <= n;) {
        int64 q = n / i;
        int64 j = n / q + 1;
        int64 sum_f = pref[j - 1] - pref[i - 1];
        cnt += sum_f * partial_sum_floor(q, pref);
        i = j;
    }
    return cnt;
}

static std::vector<char> squarefree_flags(int n) {
    std::vector<char> ok(n + 1, 1);
    ok[0] = 0;
    for (int p = 2; static_cast<int64>(p) * p <= n; ++p) {
        int64 p2 = static_cast<int64>(p) * p;
        for (int64 j = p2; j <= n; j += p2) {
            ok[static_cast<size_t>(j)] = 0;
        }
    }
    return ok;
}

static std::vector<unsigned char> pairwise_coprime_small(int lim) {
    const int stride = lim + 1;
    std::vector<unsigned char> table(static_cast<size_t>(stride) * stride, 0U);
    for (int a = 1; a <= lim; ++a) {
        for (int b = a; b <= lim; ++b) {
            const unsigned char v = (std::gcd(a, b) == 1) ? 1U : 0U;
            table[static_cast<size_t>(a) * stride + b] = v;
            table[static_cast<size_t>(b) * stride + a] = v;
        }
    }
    return table;
}

static std::vector<unsigned char> pairwise_coprime_values(const std::vector<int>& values) {
    const int m = static_cast<int>(values.size());
    std::vector<unsigned char> table(static_cast<size_t>(m) * m, 0U);
    for (int i = 0; i < m; ++i) {
        for (int j = i; j < m; ++j) {
            const unsigned char v = (std::gcd(values[static_cast<size_t>(i)],
                                              values[static_cast<size_t>(j)]) == 1) ? 1U : 0U;
            table[static_cast<size_t>(i) * m + j] = v;
            table[static_cast<size_t>(j) * m + i] = v;
        }
    }
    return table;
}

static std::vector<unsigned char> coprime_w_to_values(const std::vector<int>& values, int lim) {
    const int m = static_cast<int>(values.size());
    const int stride = lim + 1;
    std::vector<unsigned char> table(static_cast<size_t>(m) * stride, 0U);
    for (int i = 0; i < m; ++i) {
        const int v = values[static_cast<size_t>(i)];
        for (int w = 1; w <= lim; ++w) {
            table[static_cast<size_t>(i) * stride + w] = (std::gcd(w, v) == 1) ? 1U : 0U;
        }
    }
    return table;
}

static int64 count_overlap(int n) {
    int lim = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    auto is_sf = squarefree_flags(lim);

    std::vector<int> sf;
    sf.reserve(lim);
    for (int i = 1; i <= lim; ++i) {
        if (is_sf[i]) sf.push_back(i);
    }

    std::vector<int64> squares(lim + 1, 0);
    for (int i = 0; i <= lim; ++i) squares[i] = static_cast<int64>(i) * i;

    const int sf_count = static_cast<int>(sf.size());
    const int w_stride = lim + 1;
    const auto sf_pair_coprime = pairwise_coprime_values(sf);
    const auto sf_w_coprime = coprime_w_to_values(sf, lim);
    const auto w_pair_coprime = pairwise_coprime_small(lim);

    auto worker = [&](int start, int end) {
        int64 acc = 0;
        for (int p_idx = start; p_idx < end; ++p_idx) {
            const int p = sf[static_cast<size_t>(p_idx)];
            if (static_cast<int64>(p + 1) * (p + 1) > n) break;
            const unsigned char* cp =
                &sf_w_coprime[static_cast<size_t>(p_idx) * w_stride];
            for (int q_idx = 0; q_idx < sf_count; ++q_idx) {
                const int q = sf[static_cast<size_t>(q_idx)];
                if (sf_pair_coprime[static_cast<size_t>(p_idx) * sf_count + q_idx] == 0U) continue;
                int64 pq = static_cast<int64>(p) * static_cast<int64>(q);
                if ((pq + 1) * (p + q) > n) break;
                const unsigned char* cq =
                    &sf_w_coprime[static_cast<size_t>(q_idx) * w_stride];

                for (int r_idx = 0; r_idx < sf_count; ++r_idx) {
                    const int r = sf[static_cast<size_t>(r_idx)];
                    if (sf_pair_coprime[static_cast<size_t>(p_idx) * sf_count + r_idx] == 0U) continue;
                    if (sf_pair_coprime[static_cast<size_t>(q_idx) * sf_count + r_idx] == 0U) continue;
                    int64 pr = static_cast<int64>(p) * static_cast<int64>(r);
                    if ((pq + r) * (pr + q) > n) break;
                    const unsigned char* cr =
                        &sf_w_coprime[static_cast<size_t>(r_idx) * w_stride];

                    for (int s_idx = 0; s_idx < sf_count; ++s_idx) {
                        const int s = sf[static_cast<size_t>(s_idx)];
                        if (sf_pair_coprime[static_cast<size_t>(p_idx) * sf_count + s_idx] == 0U) continue;
                        if (sf_pair_coprime[static_cast<size_t>(q_idx) * sf_count + s_idx] == 0U) continue;
                        if (sf_pair_coprime[static_cast<size_t>(r_idx) * sf_count + s_idx] == 0U) continue;
                        int64 rs = static_cast<int64>(r) * static_cast<int64>(s);
                        int64 qs = static_cast<int64>(q) * static_cast<int64>(s);
                        if ((pq + rs) * (pr + qs) > n) break;
                        const unsigned char* cs =
                            &sf_w_coprime[static_cast<size_t>(s_idx) * w_stride];

                        int64 u = pq;
                        int64 v = rs;
                        int64 a = pr;
                        int64 b = qs;
                        if (u == v || a == b) continue;

                        int64 max_x = n / (a + b);
                        for (int w1 = 1; w1 <= lim; ++w1) {
                            int64 uw1 = u * squares[w1];
                            if (uw1 + v > max_x) break;
                            if (cr[w1] == 0U || cs[w1] == 0U) continue;
                            for (int w2 = 1; w2 <= lim; ++w2) {
                                int64 vw2 = v * squares[w2];
                                int64 x = uw1 + vw2;
                                if (x > max_x) break;
                                if (uw1 <= vw2) continue;
                                if (cp[w2] == 0U || cq[w2] == 0U) continue;
                                if (w_pair_coprime[static_cast<size_t>(w1) * w_stride + w2] == 0U) continue;

                                int64 m = n / x;
                                for (int w3 = 1; w3 <= lim; ++w3) {
                                    int64 aw3 = a * squares[w3];
                                    if (aw3 + b > m) break;
                                    if (cq[w3] == 0U || cs[w3] == 0U) continue;
                                    for (int w4 = 1; w4 <= lim; ++w4) {
                                        int64 bw4 = b * squares[w4];
                                        int64 y = aw3 + bw4;
                                        if (y > m) break;
                                        if (aw3 <= bw4) continue;
                                        if (cp[w4] == 0U || cr[w4] == 0U) continue;
                                        if (w_pair_coprime[static_cast<size_t>(w3) * w_stride + w4] == 0U) continue;
                                        acc += m / y;
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        return acc;
    };

    unsigned threads = std::max(1u, std::thread::hardware_concurrency());
    if (threads == 1 || n < 200000) {
        return worker(0, static_cast<int>(sf.size()));
    }

    std::vector<int64> partial(threads, 0);
    std::vector<std::thread> pool;
    pool.reserve(threads);
    for (unsigned t = 0; t < threads; ++t) {
        int start = static_cast<int>((static_cast<int64>(sf.size()) * t) / threads);
        int end = static_cast<int>((static_cast<int64>(sf.size()) * (t + 1)) / threads);
        pool.emplace_back([&, t, start, end]() {
            partial[t] = worker(start, end);
        });
    }
    for (auto& th : pool) th.join();

    int64 total = 0;
    for (int64 v : partial) total += v;
    return total;
}

static int64 compute_f_value(int n) {
    auto f = build_f(n);
    auto pref = prefix_sum(f);

    int64 cnt2 = count_type1(n, f);
    int64 cnt1 = count_type2_ordered(n, pref);
    int64 cnt3 = count_overlap(n);

    return cnt2 + (cnt1 - cnt3) / 2;
}

int main(int argc, char** argv) {
    int n = 5000000;
    if (argc > 1) n = std::atoi(argv[1]);

    std::vector<std::pair<int, int64>> checks = {
        {10, 17}, {15, 46}, {20, 86}, {30, 213}, {100, 2918}, {5000, 11134074}
    };
    for (const auto& c : checks) {
        int64 got = compute_f_value(c.first);
        if (got != c.second) {
            std::cerr << "Validation failed at n=" << c.first << ": got=" << got
                      << " expected=" << c.second << "\n";
            return 1;
        }
    }

    std::cout << compute_f_value(n) << "\n";
    return 0;
}

Python

import math

def solve():
    N = 5000000

    # Totient sieve
    phi = list(range(N+1))
    for i in range(2, N+1):
        if phi[i]==i:
            for j in range(i, N+1, i): phi[j]-=phi[j]//i

    # f(s) = phi(s)/2 - R(s) where R counts primitive i^2+j^2 reps
    f = [0]*(N+1)
    for s in range(3, N+1): f[s] = phi[s]//2
    lim = int(math.isqrt(N))
    for i in range(2, lim+1):
        i2 = i*i
        for j in range(1, i):
            if math.gcd(i,j)!=1: continue
            s = i2+j*j
            if s>N: break
            f[s] -= 1

    # Prefix sum
    pref = [0]*(N+1)
    for i in range(1, N+1): pref[i] = pref[i-1]+f[i]

    def partial_sum_floor(m):
        ret=0; i=1
        while i<=m:
            q=m//i; j=m//q+1
            ret+=q*(pref[j-1]-pref[i-1]); i=j
        return ret

    # Type 1: sum f(s)*floor(n/s)
    cnt1 = 0
    for s in range(3, N+1):
        if f[s]==0: continue
        cnt1 += f[s]*(N//s)

    # Type 2 ordered
    cnt2 = 0; i=1
    while i<=N:
        q=N//i; j=N//q+1
        sf = pref[j-1]-pref[i-1]
        cnt2 += sf*partial_sum_floor(q); i=j

    # Overlap (simplified - 4-tuple coprime enumeration)
    def sq_flags(n):
        ok=[True]*(n+1); ok[0]=False
        for p in range(2, int(math.isqrt(n))+1):
            p2=p*p
            for j in range(p2,n+1,p2): ok[j]=False
        return ok

    lim2 = int(math.isqrt(N))
    isf = sq_flags(lim2)
    sf = [i for i in range(1,lim2+1) if isf[i]]
    sqs = [i*i for i in range(lim2+1)]

    cnt3 = 0
    for pi in range(len(sf)):
        p=sf[pi]
        if (p+1)*(p+1)>N: break
        for qi in range(len(sf)):
            q=sf[qi]; pq=p*q
            if (pq+1)*(p+q)>N: break
            if math.gcd(p,q)!=1: continue
            for ri in range(len(sf)):
                r=sf[ri]; pr=p*r
                if (pq+r)*(pr+q)>N: break
                if math.gcd(p,r)!=1 or math.gcd(q,r)!=1: continue
                for si in range(len(sf)):
                    s=sf[si]; rs=r*s; qs=q*s
                    if (pq+rs)*(pr+qs)>N: break
                    if math.gcd(p,s)!=1 or math.gcd(q,s)!=1 or math.gcd(r,s)!=1: continue
                    u=pq; v=rs; a2=pr; b2=qs
                    if u==v or a2==b2: continue
                    mx=N//(a2+b2)
                    for w1 in range(1,lim2+1):
                        uw1=u*sqs[w1]
                        if uw1+v>mx: break
                        if math.gcd(r,w1)!=1 or math.gcd(s,w1)!=1: continue
                        for w2 in range(1,lim2+1):
                            vw2=v*sqs[w2]; x2=uw1+vw2
                            if x2>mx: break
                            if uw1<=vw2: continue
                            if math.gcd(p,w2)!=1 or math.gcd(q,w2)!=1: continue
                            if math.gcd(w1,w2)!=1: continue
                            m2=N//x2
                            for w3 in range(1,lim2+1):
                                aw3=a2*sqs[w3]
                                if aw3+b2>m2: break
                                if math.gcd(q,w3)!=1 or math.gcd(s,w3)!=1: continue
                                for w4 in range(1,lim2+1):
                                    bw4=b2*sqs[w4]; y2=aw3+bw4
                                    if y2>m2: break
                                    if aw3<=bw4: continue
                                    if math.gcd(p,w4)!=1 or math.gcd(r,w4)!=1: continue
                                    if math.gcd(w3,w4)!=1: continue
                                    cnt3+=m2//y2

    return str(cnt1 + (cnt2-cnt3)//2)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler585 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler585.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler585_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler585 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler585 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler585 C++ bridge produced empty output.");
        }
        return parsed;
    }

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