Detailed mathematical approach

Problem Summary

The problem revolves around the family of functions

$H(a,b,x)=\frac12-\frac{b\cos(a\pi x)+a\cos(b\pi x)}{2(a+b)}, \qquad 0\le x\le 1.$

For each odd-prime pair \(p<q\) in the interval \([m,n]\), the quantity \(F(p,q,p,2q-p)\) measures the total vertical winding produced when the curve for \(H(p,q,\cdot)\) is compared with the curve for \(H(p,2q-p,\cdot)\). The final target is

$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p),$

and the implementations evaluate the Project Euler instance \(G(500,1000)\), printing the answer to five decimal places.

The key idea is that the winding quantity is not computed by sampling the curves densely. Each curve is compressed to the sequence of its genuine turning values, and the path functional is then reconstructed from overlaps of consecutive monotone height intervals.

Mathematical Approach

The solution is built from three problem-specific facts: the analytic structure of \(H\), a sharp description of its true extrema, and a deterministic walk on the turning-value sequences of two such curves.

Endpoints and derivative factorization

For odd integers \(a\) and \(b\), the curve always starts at height \(0\) and ends at height \(1\):

$H(a,b,0)=0,\qquad H(a,b,1)=1.$

This follows from \(\cos 0=1\) and \(\cos(k\pi)=-1\) when \(k\) is odd. So every curve in the sum runs from the lower-left corner of the height picture to the upper-right one.

Differentiating gives

$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr).$

Using the identity \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\), this becomes

$H'(a,b,x)=\frac{\pi ab}{a+b}\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{(b-a)\pi x}{2}\right).$

Since the prefactor is positive, the sign of the derivative is governed entirely by

$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right).$

Where the possible turning points come from

The zeros of the two trigonometric factors give all critical-point candidates in \((0,1)\):

$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$

and, if \(a\ne b\),

$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$

The implementation collects these candidates, adds the endpoints \(0\) and \(1\), sorts everything, and removes duplicates. At this point we know where \(H'\) can vanish, but we still do not know which of those points are genuine extrema.

Filtering out tangential zeros

A zero of the derivative matters only when the sign of \(H'\) changes across it. If both trigonometric factors vanish at the same position, the derivative can touch zero without reversing direction, and such a point does not create a new piece of the winding path.

The implementation therefore examines each open interval between consecutive candidates, determines the sign of \(H'\) there, and retains only those candidate positions where the sign on the left differs from the sign on the right. This converts the analytic critical-point set into the much smaller set of true turning points.

Turning values are the real state space

Suppose the retained turning positions for \(H(a,b,\cdot)\) are \(x_0<x_1<\dots<x_r\). The solver keeps only the heights

$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$

Each consecutive pair \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) represents one monotone segment of the curve. For the winding computation, the exact \(x\)-coordinate inside such a segment is irrelevant; what matters is only the interval of reachable heights

$I_i^{(a,b)}=\bigl[\min(z_i^{(a,b)},z_{i+1}^{(a,b)}),\ \max(z_i^{(a,b)},z_{i+1}^{(a,b)})\bigr].$

This is the decisive compression step. Once the turning-value sequence is known, the rest of the geometry can be expressed purely in terms of overlaps of these height intervals.

The alternating overlap walk that defines \(F\)

Take two curves, with turning-value sequences \(z^A\) and \(z^B\). Start at the first segment of each sequence, at shared height \(z_{\mathrm{cur}}=0\), and let the first move be upward. If the current active height intervals are \(I_i^A\) and \(I_j^B\), then the next shared height is determined by the overlap of those intervals.

For an upward move, the walk climbs to the top of the overlap:

$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$

For a downward move, it drops to the bottom of the overlap:

$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$

The contribution of that step is always

$|z_{\mathrm{next}}-z_{\mathrm{cur}}|.$

After reaching \(z_{\mathrm{next}}\), whichever segment endpoint has been hit is advanced to the adjacent segment, and the direction flips. Repeating this deterministic rule until both curves have exhausted their final segment yields the total path variation \(F\).

The invariant behind the algorithm is simple but powerful: at every moment, the current shared height lies in the overlap of the two active intervals. That makes the walk unique once the turning-value sequences are fixed.

Worked example: \(H(3,5,x)\)

For \((a,b)=(3,5)\), the derivative factorization produces the interior candidates

$x=\frac14,\ \frac12,\ \frac34.$

However, \(x=\frac12\) is only a tangential zero: the derivative does not change sign there, so this point is discarded. The true turning positions are therefore

$0,\ \frac14,\ \frac34,\ 1,$

with turning values

$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1.$

This example shows exactly why the sign-flip test is needed. The raw zero set of \(H'\) is not yet the right state space; only genuine extrema determine the winding structure.

Worked example: the checkpoint \(F(3,5,3,7)\)

The second curve in that checkpoint has turning-value sequence

$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1.$

The overlap walk therefore begins

$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1,$

alternating between upward and downward moves whenever one of the active segments reaches the top or bottom of the current overlap. Summing all absolute height changes along this walk gives

$F(3,5,3,7)\approx 7.01772,$

which matches the built-in numerical checkpoint exactly.

How the Code Works

Compressing one curve

The implementation evaluates \(H\) in floating-point arithmetic and clamps tiny round-off spillovers back into \([0,1]\). For each parameter pair it builds the candidate set from the two derivative factors, sorts it, removes near-duplicates, samples the derivative sign on the intervals between consecutive candidates, and keeps only the sign-changing positions. Evaluating \(H\) at those positions produces the compact turning-value sequence.

Computing one path functional

Given two turning-value sequences, the implementation stores only the current segment index in each sequence, the current height, the current direction, and the accumulated total variation. Each iteration chooses the next shared height from the overlap formulas above, adds the absolute height change, advances whichever segment endpoint has been reached, and reverses direction. If the walk ever leaves the valid segment range, the result is treated as invalid; otherwise it terminates at height \(1\).

Summing over odd-prime pairs

To evaluate \(G(m,n)\), the implementation generates the odd primes in \([m,n]\) with a sieve and loops over all pairs \(p<q\). Each pair contributes the path variation between \(H(p,q,\cdot)\) and \(H(p,2q-p,\cdot)\). The C++ program performs this directly and can divide the pair list across worker threads. The Java version mirrors the same numerical method in a single runtime. The Python entry point is a thin wrapper that builds and runs the C++ solver, so it inherits exactly the same numerical behavior. The numerical checkpoints \(F(7,17,9,19)\approx 26.79578\) and \(G(3,20)\approx 463.80866\) are used to confirm that the turning-point logic and the overlap walk agree.

Complexity Analysis

Let \(P\) be the number of odd primes in \([m,n]\). Then the outer sum contains \(\binom{P}{2}\) prime pairs. For one curve \(H(a,b,\cdot)\), the number of candidate critical points is \(O(a+b+|b-a|)=O(\max(a,b))\), and the subsequent overlap walk is linear in the combined lengths of the two turning-value sequences. So one pair requires linear work in the scale of its parameters.

If the upper endpoint of the interval is \(N\), a safe worst-case bound is \(O(P^2N)\) time. For the fixed Project Euler input, the parameters lie between 500 and 1000, so the practical cost is dominated by the quadratic sweep over prime pairs rather than by any single path computation. Memory usage is modest: the sieve is linear in the interval limit, and each active worker needs only a few temporary turning-value arrays and one partial sum.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=975
2. Trigonometric identities: Wikipedia - Trigonometric identities
3. Critical point: Wikipedia - Critical point
4. Total variation: Wikipedia - Total variation
5. Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes