Problem 404: Crisscross Ellipses

Detailed mathematical approach

Problem Summary

For a positive number \(a\), let \(E_a\) be the ellipse \(x^2+4y^2=4a^2\). It is centred at the origin, with semi-axes \(2a\) along the \(x\)-axis and \(a\) along the \(y\)-axis. ^[2] Let \(E_a^{\prime}\) be the copy of \(E_a\) rotated anticlockwise about the origin by an angle \(\theta\) with \(0^{\circ}\lt\theta\lt90^{\circ}\). The two curves cross in four points. Let \(b\) be the distance from the origin to the two crossings nearest it and \(c\) the distance to the other two. The ordered triplet \((a,b,c)\) is a canonical ellipsoidal triplet when \(a\), \(b\) and \(c\) are all positive integers; the statement offers \((209,247,286)\) as an example. Writing \(C(N)\) for the number of distinct such triplets with \(a\le N\), it gives \(C(10^3)=7\), \(C(10^4)=106\) and \(C(10^6)=11845\), and asks for \(C(10^{17})\). ^[1]

The whole solution in one sentence: the four crossings lie on the two bisectors of the ellipses' axes, which turns the condition on \((a,b,c)\) into the single equation \(4a^2/b^2+4a^2/c^2=5\); its integer solutions are the whole multiples of a short list of basic triplets, so \(C(N)\) is a sum of floor divisions.

Mathematical Approach

1. The ellipse in polar coordinates

Put \(x=r\cos\phi\) and \(y=r\sin\phi\) in \(x^2+4y^2=4a^2\). Using \(\cos^2\phi=1-\sin^2\phi\) the equation collapses to

$$r^2\bigl(1+3\sin^2\phi\bigr)=4a^2 .$$

So the radius in direction \(\phi\) is \(2a/\sqrt{1+3\sin^2\phi}\), which is \(2a\) at \(\phi=0\) and \(a\) at \(\phi=\pi/2\), as it must be. ^[2] Rotating a curve by \(\theta\) adds \(\theta\) to the polar angle of each of its points, so a point \((r,\phi)\) lies on \(E_a^{\prime}\) exactly when \((r,\phi-\theta)\) lies on \(E_a\):

$$r^2\bigl(1+3\sin^2(\phi-\theta)\bigr)=4a^2 .$$

2. Where the two ellipses cross

Each direction \(\phi\) carries exactly one point of \(E_a\) with \(r\gt0\), and that point also lies on \(E_a^{\prime}\) exactly when the two left-hand sides agree, that is when \(\sin^2\phi=\sin^2(\phi-\theta)\). The double-angle and product-to-sum identities ^[3] give

$$\sin^2\phi-\sin^2(\phi-\theta)=\sin(2\phi-\theta)\sin\theta ,$$

and \(\sin\theta\ne0\) because \(0^{\circ}\lt\theta\lt90^{\circ}\). Hence \(\sin(2\phi-\theta)=0\), that is \(2\phi-\theta=k\pi\), so the crossings lie in the four directions

$$\phi=\frac{\theta}{2}+k\cdot\frac{\pi}{2},\qquad k=0,1,2,3 .$$

So there are exactly four crossings, in two opposite pairs. The direction \(\phi=\theta/2\) bisects the two long axes, which point along \(0\) and \(\theta\); \(\phi=\theta/2+\pi/2\) bisects the two short axes.

3. The identity behind the problem

Write \(s=\sin^2(\theta/2)\), so that \(\cos^2(\theta/2)=1-s\) by the Pythagorean identity. ^[3] Substituting the two directions into Section 1 gives the two squared distances

$$\begin{aligned}c^2&=\frac{4a^2}{1+3s},\\[2pt] b^2&=\frac{4a^2}{4-3s} .\end{aligned}$$

Because \(0^{\circ}\lt\theta/2\lt45^{\circ}\) we have \(0\lt s\lt\tfrac12\), hence \(1+3s\lt4-3s\): the bisector of the long axes carries the far pair and the bisector of the short axes the near pair, so the labelling above is the right one. Adding the two reciprocals now makes the angle disappear completely:

$$\frac{4a^2}{b^2}+\frac{4a^2}{c^2}=(4-3s)+(1+3s)=5 .$$

This one equation is the whole geometric content of the problem. Solved for \(c\) it gives the form the brute force uses, \(c^2=4a^2b^2/(5b^2-4a^2)\).

4. Which triplets can occur

From \(s\in(0,\tfrac12)\) we read off \(1+3s\in(1,\tfrac52)\) and \(4-3s\in(\tfrac52,4)\), hence

$$a\lt b\lt a\sqrt{8/5}\lt c\lt 2a .$$

Conversely, let \(a,b,c\) be positive numbers with \(4a^2/b^2+4a^2/c^2=5\) and \(a\lt b\lt c\). Then \(s:=(4-4a^2/b^2)/3\) satisfies \(s\gt0\) because \(b\gt a\), and \(s\lt\tfrac12\) because \(b\lt c\) forces \(4a^2/b^2\gt\tfrac52\); so a unique \(\theta\in(0^{\circ},90^{\circ})\) produces that triplet. The conditions are therefore equivalent, and

$$\begin{aligned}C(N)=\#\{(a,b,c):\ &a\le N,\ a\lt b\lt c,\\ &4a^2/b^2+4a^2/c^2=5\},\end{aligned}$$

all three entries being positive integers. Dropping \(c\lt2a\) costs nothing: \(b\gt a\) gives \(4a^2/b^2\lt4\), hence \(4a^2/c^2\gt1\).

5. Primitive integer solutions

Set \(X=2a/b\) and \(Y=2a/c\). Section 3 says that \((X,Y)\) is a rational point with \(X,Y\gt0\) on the circle \(X^2+Y^2=5\), and Section 4 that \(1\lt Y\lt X\lt2\). Write it over a common denominator in lowest terms, \(X=U/W\) and \(Y=V/W\) with positive integers and \(\gcd(U,V,W)=1\). Then

$$U^2+V^2=5W^2,\qquad W\lt V\lt U\lt 2W .$$

Such a triple is called primitive. Three facts follow at once.

• Pairwise coprime. If a prime \(p\) divides \(U\) and \(W\) then \(V^2=5W^2-U^2\) gives \(p\mid V\), against \(\gcd(U,V,W)=1\); likewise for \(V\) and \(W\). If \(p\) divides \(U\) and \(V\) then \(p\mid 5W^2\), so either \(p\mid W\) — excluded — or \(p=5\), and then \(25\mid U^2+V^2=5W^2\) forces \(5\mid W\), excluded again.
• \(5\) divides neither \(U\) nor \(V\). If \(5\mid U\) then \(5\mid V^2=5W^2-U^2\), so \(5\mid V\), against the previous point.
• \(W\) is odd and exactly one of \(U,V\) is even. Both even is excluded. Both odd would give \(U^2+V^2\equiv2\pmod 8\), hence \(5W^2\equiv2\) and, multiplying by \(5\) since \(25\equiv1\), \(W^2\equiv2\pmod 8\), which no square is. So \(U^2+V^2\) is odd and \(W\) is odd, and \(UV/2\) is a positive integer.

6. Every triplet is a multiple of a basic one

From \(X=2a/b=U/W\) and \(Y=2a/c=V/W\) we get \(b=2aW/U\) and \(c=2aW/V\). Since \(\gcd(U,W)=1\), integrality of \(b\) forces \(U\mid 2a\); likewise \(V\mid 2a\); and \(\gcd(U,V)=1\) then gives \(UV\mid 2a\). Writing \(2a=kUV\) with \(k\ge1\),

$$a=k\,\frac{UV}{2},\qquad b=kVW,\qquad c=kUW .$$

Conversely every primitive triple and every \(k\ge1\) give a triplet: the entries are integers, they satisfy Section 3's equation because \(U^2+V^2=5W^2\), and \(W\lt V\lt U\) gives \(a\lt b\lt c\); and a triplet determines \((X,Y)\), hence \((U,V,W)\) and \(k\). The correspondence is one-to-one, so with \(a_0=UV/2\):

$$C(N)=\sum_{(U,V,W)}\left\lfloor\frac{N}{a_0}\right\rfloor .$$

All that is left is to list every primitive triple with \(a_0\le N\), once each.

7. All primitive solutions from Gaussian integers

A Gaussian integer is a number \(u+vi\) with \(u,v\) integers and \(i^2=-1\). The whole construction fits on one small example: \((3-2i)^2=5-12i\) and \((2+i)(5-12i)=22-19i\), and its coordinates satisfy \(22^2+19^2=5\cdot13^2\), which Section 6 turns into the statement's triplet \((209,247,286)\).

Why it reaches every primitive triple needs some arithmetic of \(\mathbb{Z}[i]\), the ring of all Gaussian integers. The conjugate of \(z=u+vi\) is \(\bar z=u-vi\) and its norm is \(N(z)=z\bar z=u^2+v^2\); multiplying each factor by its own conjugate shows \(N(zw)=N(z)N(w)\). The units — those \(z\) whose reciprocal is again a Gaussian integer — are \(1,-1,i,-i\), and factorisation into primes is unique up to order and unit factors. The prime \(2\) is a square times a unit, \(2=i(1-i)^2\), and \(1+i=i(1-i)\) is a unit multiple of \(1-i\); a rational prime \(p\equiv1\pmod 4\) splits into two conjugate primes, while one with \(p\equiv3\pmod 4\) stays prime. ^{[4] [5]} The factorisation that matters here is \(5=(2+i)(2-i)\), both factors of norm \(5\). The same square-factor argument proves Euclid's formula for Pythagorean triples, here with the extra factor \(2+i\). ^[6]

For integers \(m,n\) put \(\mu=m-ni\) and \(Z=(2+i)\mu^2\). Then \(N(Z)=5\,N(\mu)^2\), and expanding the product gives \(Z=Y_0+X_0i\) with

$$\begin{aligned}X_0&=m^2-4mn-n^2,\\ Y_0&=2(m^2+mn-n^2),\\ W_0&=m^2+n^2=N(\mu),\end{aligned}$$

so that \(X_0^2+Y_0^2=5W_0^2\) holds identically — the identity is just multiplicativity of the norm. These are exactly the three expressions the programs evaluate.

The converse holds too. Let \((U,V,W)\) be primitive. No prime \(p\equiv3\pmod4\) divides \(W\): such a \(p\) stays prime in \(\mathbb{Z}[i]\), and \(p^2\mid N(U+Vi)\) would force \(p\mid U+Vi\), contradicting \(\gcd(U,V)=1\). Next, \(U^2\equiv-V^2\pmod 5\) and \(5\nmid V\) give \((U/V)^2\equiv4\), so \(U\equiv2\sigma V\pmod 5\) for exactly one sign \(\sigma=\pm1\). Put \(T=U+\sigma Vi\). Both coordinates of \((2-i)T\), namely \(2U+\sigma V\) and \(2\sigma V-U\), are then divisible by \(5\), so \(\zeta=(2-i)T/5\) lies in \(\mathbb{Z}[i]\), has \(N(\zeta)=W^2\), and \(T=(2+i)\zeta\). No rational prime divides \(\zeta\), because such a prime would divide \(T\). For each prime \(p\mid W\), splitting as \(p=\pi\bar\pi\), the exponents of \(\pi\) and \(\bar\pi\) in \(\zeta\) add up to the exponent of \(p\) in \(N(\zeta)=W^2\), which is even, and one of the two is \(0\); so every prime exponent of \(\zeta\) is even and \(\zeta=\varepsilon\mu^2\) for a unit \(\varepsilon\) and some \(\mu\) with \(N(\mu)=W\). A unit multiple of \(T\) only permutes and negates the two coordinates, leaving the pair \(\{U,V\}\) alone, and multiplies \(\varepsilon\) by that unit, so we may take \(\varepsilon=1\) and \(T=(2+i)\mu^2\). Writing \(\mu=m-ni\) we get \(\gcd(m,n)=1\), \(W_0=W\) and \(\{|X_0|,|Y_0|\}=\{U,V\}\). Every primitive triple therefore comes from a coprime pair, with no common factor to remove.

8. The window \(m/2\lt n\lt m\)

Which pairs land in the arc \(W\lt V\lt U\lt2W\)? Suppose \(m\gt n\gt m/2\gt0\). Then \(X_0\lt0\), \(Y_0\gt0\), and four factorisations settle it:

• \(|X_0|-W_0=2m(2n-m)\gt0\);
• \(2W_0-|X_0|=(3m-n)(m-n)\gt0\);
• \(Y_0-W_0=(m+3n)(m-n)\gt0\);
• \(2W_0-Y_0=2n(2n-m)\gt0\).

So both \(|X_0|\) and \(Y_0\) lie strictly between \(W_0\) and \(2W_0\). Dividing all three by a common factor leaves these ratios alone, so after reduction and sorting the triple satisfies \(W\lt V\lt U\lt2W\). Sorting only fixes the names: the product \(UV\), hence \(a_0\), does not care.

The same computation, read as a statement about \(r=n/m\) alone and allowing any positive \(m,n\), gives the converse: the two ratios \(|X_0|/W_0\) and \(|Y_0|/W_0\) both lie in \((1,2)\) exactly when \(r\in(\tfrac12,1)\) or \(r\gt3\), the factors governing that second branch being \((r-1)(r-3)\) for \(|X_0|\) and \((r-3)(r+1)\) for \(|Y_0|\). Pairs with \(mn=0\) are out, since they give \(|X_0|=W_0\).

9. The common factor and the duplicates

Let \(g=\gcd(|X_0|,Y_0,W_0)\). Then \(g\in\{1,2,5,10\}\). For an odd prime \(p\mid g\): \(p\) divides \(W_0+X_0=2m(m-2n)\) and \(2W_0-Y_0=2n(2n-m)\). If \(p\nmid(m-2n)\) then \(p\mid m\) and \(p\mid n\), against \(\gcd(m,n)=1\); so \(m\equiv2n\pmod p\), whence \(W_0\equiv5n^2\equiv0\) and, as \(p\nmid n\), \(p=5\). And \(25\nmid g\): otherwise \(25\mid Z\), which forces \((2+i)(2-i)^2\mid\mu^2\), hence \(5\mid\mu\) and \(5\mid\gcd(m,n)\). For the factor \(2\): if \(m\) and \(n\) have different parity then \(W_0\) is odd; if both are odd then \(W_0\equiv2\pmod 8\), so \(4\nmid W_0\), while \(X_0\) and \(Y_0\) are both even. So the \(2\)-part of \(g\) is \(1\) or \(2\) and its \(5\)-part \(1\) or \(5\), leaving exactly the four values claimed.

A pair with \(5\mid g\) is a duplicate and is skipped. In \(\mathbb{Z}[i]\), \(5\mid Z\) means \((2-i)\mid\mu\), so \(\mu=(2-i)\nu\) and \(Z/5=(2-i)\nu^2\), which is the conjugate of \((2+i)\bar\nu^2\); conjugation leaves \(\{|X_0|,|Y_0|\}\) unchanged, so the reduced triple is one that a smaller pair already produces. The smallest example in the window is \((m,n)=(7,6)\), giving \((|X_0|,Y_0,W_0)=(155,110,85)\) and \(g=5\); its reduced triple \((31,22,17)\) comes from \((5,3)\). And \(g=10\) behaves the same way: \((9,7)\) gives \((220,190,130)\), again the triple of \((3,2)\).

10. Nothing is missed and nothing is counted twice

Start from the representation of Section 7, where \(N(\mu)=W\) and nothing has to be divided out. A unit multiple of \(\mu\) only changes the sign of \(Z\), and each unit orbit of \(\mu=m-ni\) has one representative with \(m,n\gt0\), so take both positive; \(mn=0\) is out, since it gives \(|X_0|=W_0\). Since \(W=m^2+n^2\) is odd, \(m\) and \(n\) have opposite parity. By Section 8 the ratio \(n/m\) now lies in the window or above \(3\). In the second case multiply \(\mu\) by \(1+i\), which sends the pair to \((m+n,\,n-m)\) and \(Z\) to \(2iZ\). Both new entries are positive and odd, and a common divisor would divide \(2m\) and \(2n\), so they are coprime; and \(n\gt3m\) gives

$$\frac12\lt\frac{n-m}{m+n}\lt1 .$$

The denominator has doubled to \(2W\), so the common factor is now exactly \(2\) and reduction returns the same primitive triple. Every primitive triple of the arc therefore has a retained representation in the window.

It has only one. For a retained window pair put \(\nu=\mu\) when \(g=1\) and \(\nu=\mu/(1+i)\) when \(g=2\) — both \(m,n\) are then odd, so \(1+i\) divides \(\mu\) — so that \(N(\nu)=W\) and \(T=(2+i)\nu^2\) has \(\{U,V\}\) as its absolute coordinates. Suppose two retained pairs give the same primitive triple. Their values \(T_1,T_2\) then differ by a unit, or by a unit and conjugation. Each \(T_j\) is divisible by \(2+i\) but not by \(2-i\), since both together would make \(5\) divide \(U\) and \(V\), which Section 5 forbids; conjugation exchanges the two factors, so the conjugate case cannot occur. Hence \(\nu_2^2\) is a unit times \(\nu_1^2\), and unique factorisation makes \(\nu_2\) a unit times \(\nu_1\).

If the two pairs have the same \(g\), their \(\mu\) also differ by a unit, and as each orbit has one representative with positive entries the pairs coincide. If the factors are \(1\) and \(2\), then \(\mu_2\) is a unit times \((1+i)\mu_1\), whose positive representative is \((m_1-n_1,\,m_1+n_1)\) with ratio \((m_1+n_1)/(m_1-n_1)\gt3\), outside the window. The mixed case is therefore impossible, and the enumeration over the window, with the \(5\mid g\) pairs skipped, visits every primitive triple of the arc exactly once.

11. Where the search stops

Only pairs with \(a_0\le N\) contribute. With \(r=n/m\in(\tfrac12,1)\) the two numerators are \(|X_0|=m^2(r^2+4r-1)\) and \(Y_0=2m^2(1+r-r^2)\), and since \(g\le2\) after the skip,

$$\begin{aligned}a_0=\frac{|X_0|\,Y_0}{2g^2}&\ge\frac{|X_0|\,Y_0}{8}\\ &=\frac{m^4}{4}\,h(r)\end{aligned}$$

with \(h(r)=(r^2+4r-1)(1+r-r^2)\). Expanding, \(h(r)=-r^4-3r^3+6r^2+3r-1\), so \(h''(r)=-6(2r-1)(r+2)\), which is negative for \(r\gt\tfrac12\). Hence \(h'\) decreases on \([\tfrac12,1]\), and \(h'(1)=2\gt0\), so \(h'\gt0\) there and \(h\) increases from \(h(\tfrac12)=\tfrac{25}{16}\) to \(h(1)=4\). Therefore \(a_0\ge\frac{25}{64}m^4\), and a pair can contribute only when

$$m^4\le\frac{64}{25}\,N .$$

In particular \(m\le(4N)^{1/4}\) is a safe cut-off. The programs use \(\lfloor(4N)^{1/4}\rfloor+2\), two values of slack above it. At \(N=10^{17}\) the fourth root is \(25148.67\), so the cut-off is \(m\le25148\) and the loop reaches \(m=25150\); the largest \(m\) that actually contributes is \(22493\).

12. Worked example: \(C(1000)=7\)

Here \((4N)^{1/4}\approx7.95\), so the programs look at \(m\le9\) and only a handful of pairs matter. Three of them survive.

• \((m,n)=(3,2)\): \((X_0,Y_0,W_0)=(-19,22,13)\), \(g=1\), so \((U,V,W)=(22,19,13)\), \(a_0=22\cdot19/2=209\), \(b_0=19\cdot13=247\) and \(c_0=22\cdot13=286\) — the statement's example, at \(\theta\approx75.96^{\circ}\). It contributes \(\lfloor1000/209\rfloor=4\).
• \((m,n)=(5,3)\): \((-44,62,34)\), \(g=2\), so \((U,V,W)=(31,22,17)\) and \((a_0,b_0,c_0)=(341,374,527)\), with \(\theta\approx56.62^{\circ}\). It contributes \(\lfloor1000/341\rfloor=2\).
• \((m,n)=(4,3)\): \((-41,38,25)\), \(g=1\), so \((U,V,W)=(41,38,25)\) and \((a_0,b_0,c_0)=(779,950,1025)\), with \(\theta\approx82.74^{\circ}\). It contributes \(\lfloor1000/779\rfloor=1\).

The seven triplets are \((209,247,286)\) with its multiples by \(2,3,4\), then \((341,374,527)\) with its double, and \((779,950,1025)\). A direct search over \(a\le1000\) finds the same seven and no others, so \(C(1000)=7\). The smallest basic triplet already has \(a_0=209\), which is why \(C(100)=0\); the next one after the three above comes from \((7,5)\) with \(a_0=1711\). As a check of Section 3, \(4\cdot209^2=174724\) while \(247^2=61009\) and \(286^2=81796\), and indeed \(174724/61009+174724/81796=5\); also \(247\lt209\sqrt{8/5}\lt286\).

How the Code Works

The C++ program reads its options first. parse_arguments accepts --n= to replace the default limit \(10^{17}\) and --skip-checkpoints to suppress the self-checks; parse_u64_after_prefix parses the digits itself and raises std::overflow_error instead of wrapping.

fourth_root_bound implements Section 11. It starts from a floating-point estimate of \((4N)^{1/4}\), walks down while \(m^4\) exceeds \(4N\) and up while \((m+1)^4\) does not — so the result is exact whatever the floating-point error was — and adds \(2\) as slack. The fourth powers and \(4N\) itself are formed in u128, the file's alias for __uint128_t, ^[8] so that an overshooting estimate cannot overflow the comparison; \(4N\) alone already needs \(59\) bits.

count_triplets_fast is the enumeration of Sections 7 to 11. For each \(m\) from \(2\) to the bound it runs \(n\) from \(\lfloor m/2\rfloor+1\) to \(m-1\) — exactly the window \(m/2\lt n\lt m\) — skips \(n\) with \(\gcd(m,n)\ne1\) using std::gcd, ^[7] and forms x_raw, y_abs and w from the formulas of Section 7. Because \(X_0\) is negative throughout the window it is computed in a signed 128-bit type and passed through abs_i128; since \(n\lt m\) we have \(0\lt Y_0\lt4m^2\) and \(0\lt W_0\lt2m^2\), which at \(m\le25150\) stay well inside 64 bits. Two more calls give \(g\); pairs with g % 5 == 0 are dropped as the duplicates of Section 9; the numerators are divided by \(g\) and swapped into \(u\ge v\); and a_primitive \(=uv/2\) is formed in __uint128_t. If it exceeds the limit the family is empty, otherwise limit / a_primitive joins a 128-bit running sum, and the total is returned as a 64-bit unsigned integer. That it fits needs no computation: the window holds fewer than \(m/2\) pairs for each \(m\ge3\), none at \(m=2\), and \(a_0\ge25m^4/64\), so

$$C(N)\lt\frac{32N}{25}\sum_{m\ge3}\frac1{m^3}\lt\frac{4N}{25},$$

by \(\sum_{m\ge3}m^{-3}\lt\int_2^{\infty}x^{-3}dx=\tfrac18\) — far inside both the unsigned 64-bit C++ result and the signed Java sum.

count_triplets_bruteforce is the independent check. For each \(a\) it bounds \(b\) by \(b^2\le(8a^2-1)/5\), which is exactly \(b\lt a\sqrt{8/5}\) from Section 4, computes \(c^2=4a^2b^2/(5b^2-4a^2)\) when the division is exact, and accepts \(b\lt c\lt2a\) when \(c^2\) is a perfect square. isqrt_u128 corrects a floating-point square-root estimate in both directions, so it is exact. run_checkpoints then verifies \(C(100)=0\) and \(C(10^4)=106\) and compares the fast method with the brute force at \(N=5000\), throwing on any mismatch.

The Python version does the same arithmetic in one function. Its integers have unlimited precision, ^[9] so no size analysis is needed: fourth_root_bound, the math.gcd filter, the g % 5 skip and total += limit // a_prim read almost like the mathematics. It has a fixed limit and no self-checks.

The Java version follows the C++ one with a hand-written Euclidean gcd and keeps everything in long, whose range is \(-2^{63}\) to \(2^{63}-1\). ^[9] This is sound but not generous: the largest product \(u\cdot v\) at \(N=10^{17}\) is \(3200614060901447098\), about a third of Long.MAX_VALUE, since \(UV\le10\,m^4\) is of order \(40N\); the fourth powers also fit, as \(25150^4\lt4.1\cdot10^{17}\). The test p * p > 0 in fourthRootBound rejects a wrapped product that turned non-positive, but it is not a general overflow test: a wrapped long product can stay positive, as \(70000^4\) does. Raising the limit needs fresh bounds for \(4N\), the fourth powers, \(uv\) and the sum. Like Python, the Java version runs no checks.

Complexity and Verification

The work is one pass over the window. With \(M=\lfloor(4N)^{1/4}\rfloor+2\) the window holds

$$\begin{aligned}P(M)&=\sum_{m=2}^{M}\Bigl(m-1-\Bigl\lfloor\frac m2\Bigr\rfloor\Bigr)\\ &=\Bigl\lfloor\frac{(M-1)^2}{4}\Bigr\rfloor=\tfrac12\sqrt N+O(N^{1/4})\end{aligned}$$

pairs, the closed form following by grouping even and odd \(m\). Counting arithmetic on a fixed number of machine words as unit cost, the method takes \(O(\sqrt N\log N)\) time — at most three gcd evaluations per pair, each with a logarithmic step count — and keeps a constant number of integers, no list of triplets. In Python those integers grow to \(O(\log N)\) bits, so the bound counts arithmetic steps, not bit operations. For \(N=10^{17}\) the exact figures are \(M=25150\), \(P(M)=158\,118\,050\) pairs visited, \(96\,133\,502\) of them coprime, \(80\,112\,017\) surviving the \(5\mid g\) skip and \(31\,881\,671\) contributing a family with \(a_0\le N\); the coprime fraction matches \(6/\pi^2\) to three digits. On one desktop machine a pass takes about \(12\) seconds for the C++ program at -O2 and about \(36\) seconds for the Python one on CPython 3.

Beyond the packaged checks, the following were verified for this article. A direct search over all triplets agrees with the fast method for every \(N\le5000\), where \(C(3000)=27\) and \(C(5000)=50\); the C++, Python and Java programs, run at reduced limits, all return \(C(10^3)=7\), \(C(10^4)=106\), \(C(10^5)=1153\) and \(C(10^6)=11845\), reproducing the three statement checkpoints. The crossing directions and radii of Section 2 were checked numerically against both ellipse equations for several \((a,\theta)\), and the angle recovered from \((209,247,286)\) reproduces \(b\) and \(c\). Independently of the parametrisation, a direct search for primitive solutions in the arc \(W\lt V\lt U\lt2W\) finds \(327\) with \(W\le3000\), and the window enumeration produces exactly that set, each triple once, confirming Section 10; tracing each back to its coprime pair, the ratio \(n/m\) always fell in \((\tfrac12,1)\) or above \(3\), as Section 8 claims. Over the window pairs with \(m\lt600\) the factor \(g\) took only the values \(1,2,5,10\), with frequencies \(30360\), \(15188\), \(6075\), \(3046\) — close to \(10:5:2:1\). Among the pairs kept after the \(5\mid g\) skip the smallest \(a_0/m^4\) was \(0.3919336\), at \((m,n)=(597,299)\), just above \(25/64=0.390625\) as Section 11 requires; counting the skipped pairs too it drops to about \(0.0159\), so the bound really does need that skip.

Footnotes and References

The references below supply the background facts used above. The polar description of the crossings, the identity of Section 3, the reduction to primitive triples, the completeness of the window and the fourth-root bound are derived here.

1. Project Euler 404 — Crisscross Ellipses. The official statement defines E_a by x² + 4y² = 4a² and E_a′ as its anticlockwise rotation by θ degrees with 0° < θ < 90°, calls (a, b, c) canonical when a and the two intersection distances b and c are all positive integers, and gives C(10³) = 7, C(10⁴) = 106 and C(10⁶) = 11845 before asking for C(10¹⁷).
2. Wikipedia — Ellipse. It gives x²/a² + y²/b² = 1 with a and b the semi-major and semi-minor axes, and the polar form about the centre, r(θ) = ab/√((b cos θ)² + (a sin θ)²); Section 1 is that formula for semi-axes 2a and a.
3. Wikipedia — List of trigonometric identities. It states the Pythagorean identity sin²θ + cos²θ = 1, which makes the angle cancel in Section 3, and the product-to-sum identities, which with the double-angle formula give the squared-sine difference of Section 2.
4. Wikipedia — Gaussian integer. It states that ℤ[i] is a Euclidean domain and hence a unique factorisation domain, that its units are 1, −1, i, −i, that the norm of a + bi is a² + b², that 2 is the unique ramified prime with 2 = i(1 − i)², and that a rational prime splits into two conjugate Gaussian primes when it is 1 mod 4 and is inert when it is 3 mod 4.
5. Wikipedia — Fermat's theorem on sums of two squares. An odd prime p is a sum of two squares if and only if p ≡ 1 (mod 4). Hence the denominator W of Section 5 has only prime factors 1 mod 4 and is a norm from ℤ[i]: 13, 17 and 25 occur in the worked example, 3, 7 and 11 never.
6. Wikipedia — Pythagorean triple. It gives Euclid's formula m² − n², 2mn, m² + n² for primitive triples when gcd(m, n) = 1 with exactly one of m, n even, and derives it in ℤ[i] by writing a + bi = ε(m + ni)² with ε a unit. Section 7 is that argument with the extra factor 2 + i.
7. cppreference — std::gcd. Declared in <numeric> since C++17, it returns the greatest common divisor of the absolute values of its two integer arguments, and zero only when both are zero.
8. GCC manual — 128-bit integers. __int128 and unsigned __int128 are GCC extensions supported on targets with an integer mode wide enough to hold 128 bits; the C++ version uses them for the bound's fourth powers, the negative X₀, the product uv and the running total.
9. Java — Long documents Long.MAX_VALUE as 2⁶³ − 1 and Long.MIN_VALUE as −2⁶³, the range the Java port's products and sums must fit. Python — numeric types states that integers have unlimited precision, so Python needs no such analysis.