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.
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.
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).$$
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.
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.
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.
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.
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.
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.
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.
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\).
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.
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.
Im Mittelpunkt der Aufgabe steht die Funktionenfamilie
$$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.$$
Für jedes Paar ungerader Primzahlen \(p<q\) im Intervall \([m,n]\) misst \(F(p,q,p,2q-p)\), wie stark sich der zugehörige Höhenpfad windet, wenn man die Kurve \(H(p,q,\cdot)\) mit \(H(p,2q-p,\cdot)\) vergleicht. Gesucht ist also
$$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p).$$
Die Implementierungen berechnen den Project-Euler-Fall \(G(500,1000)\) und geben das Ergebnis mit fünf Nachkommastellen aus.
Entscheidend ist, dass die Pfadgröße nicht durch dichtes Abtasten der Kurven bestimmt wird. Jede Kurve wird auf die Folge ihrer echten Wendewerte komprimiert; daraus lässt sich der gesamte Windungspfad allein über Überlappungen benachbarter monotoner Höhenintervalle rekonstruieren.
Die Lösung ruht auf drei problemspezifischen Bausteinen: der analytischen Struktur von \(H\), einer präzisen Beschreibung der echten Extremstellen und einem deterministischen Lauf auf den Wendewertfolgen zweier solcher Kurven.
Für ungerade ganze Zahlen \(a\) und \(b\) beginnt die Kurve stets bei Höhe \(0\) und endet bei Höhe \(1\):
$$H(a,b,0)=0,\qquad H(a,b,1)=1.$$
Der Grund ist \(\cos 0=1\) und \(\cos(k\pi)=-1\) für ungerades \(k\). Jede in der Summe vorkommende Kurve verbindet also die unteren und oberen Randhöhen.
Die Ableitung lautet
$$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr).$$
Mit \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\) folgt
$$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).$$
Damit wird das Vorzeichen von \(H'\) vollständig durch
$$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right)$$
bestimmt.
Aus den Nullstellen der beiden Faktoren erhält man alle Kandidaten für kritische Punkte in \((0,1)\):
$$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$$
und, falls \(a\ne b\), zusätzlich
$$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$$
Die Implementierung sammelt diese Punkte, ergänzt \(0\) und \(1\), sortiert die Liste und entfernt Dubletten. Damit kennt man alle Stellen, an denen \(H'\) verschwinden kann, aber noch nicht die echten Extremstellen.
Eine Nullstelle der Ableitung ist nur dann relevant, wenn das Vorzeichen von \(H'\) dort tatsächlich wechselt. Verschwinden beide trigonometrischen Faktoren gleichzeitig, kann die Ableitung Null werden, ohne dass die Kurve ihre Bewegungsrichtung ändert. Solche Stellen erzeugen keinen neuen Abschnitt des Windungspfads.
Deshalb untersucht die Implementierung jedes offene Intervall zwischen zwei benachbarten Kandidaten, bestimmt dort das Vorzeichen von \(H'\) und behält nur diejenigen Kandidaten, an denen links und rechts unterschiedliche Vorzeichen vorliegen. So wird aus der analytischen Kandidatenmenge die deutlich kleinere Menge der echten Wendestellen.
Seien \(x_0<x_1<\dots<x_r\) die behaltenen Wendestellen von \(H(a,b,\cdot)\). Gespeichert werden nur die zugehörigen Höhen
$$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$$
Jedes aufeinanderfolgende Paar \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) beschreibt einen monotonen Kurvenabschnitt. Für die Windungsgröße ist die genaue \(x\)-Position innerhalb dieses Abschnitts unwichtig; entscheidend ist nur das erreichbare Höhenintervall
$$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].$$
Hier liegt die eigentliche Kompression des Problems: Kennt man die Wendewertfolge, lässt sich die restliche Geometrie allein mit Überlappungen dieser Höhenintervalle ausdrücken.
Betrachte zwei Kurven mit Wendewertfolgen \(z^A\) und \(z^B\). Man startet auf dem ersten Segment jeder Folge bei gemeinsamer Höhe \(z_{\mathrm{cur}}=0\); die erste Bewegung geht nach oben. Sind die aktuell aktiven Höhenintervalle \(I_i^A\) und \(I_j^B\), dann wird die nächste gemeinsame Höhe nur aus deren Überlappung bestimmt.
Bei einer Aufwärtsbewegung steigt der Pfad bis zum oberen Rand der Überlappung:
$$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$$
Bei einer Abwärtsbewegung fällt er bis zum unteren Rand:
$$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$$
Der Beitrag dieses Schritts ist stets
$$|z_{\mathrm{next}}-z_{\mathrm{cur}}|.$$
Danach wird jeder Abschnitt, dessen passender Rand gerade erreicht wurde, auf den benachbarten Abschnitt weitergeschaltet, und die Bewegungsrichtung wechselt. Wiederholt man diese deterministische Regel bis beide Kurven ihren letzten Abschnitt verlassen haben, erhält man die gesamte Pfadvariation \(F\).
Die tragende Invariante lautet: Die aktuelle gemeinsame Höhe liegt immer in der Überlappung der beiden aktiven Intervalle. Dadurch ist der Lauf eindeutig, sobald die beiden Wendewertfolgen feststehen.
Für \((a,b)=(3,5)\) liefert die Faktorisierung der Ableitung die inneren Kandidaten
$$x=\frac14,\ \frac12,\ \frac34.$$
Der Punkt \(x=\frac12\) ist jedoch nur eine tangentiale Nullstelle; dort wechselt das Vorzeichen der Ableitung nicht. Die echten Wendestellen sind daher
$$0,\ \frac14,\ \frac34,\ 1,$$
mit den Wendewerten
$$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1.$$
Genau dieses Beispiel zeigt, warum der Vorzeichenwechsel-Test unverzichtbar ist: Die rohe Nullstellenmenge von \(H'\) ist noch nicht der richtige Zustandsraum.
Die zweite Kurve in diesem Prüfbeispiel besitzt die Wendewertfolge
$$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1.$$
Der Überlappungslauf beginnt also mit
$$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1,$$
wobei sich Aufwärts- und Abwärtsschritte immer dann abwechseln, wenn einer der aktiven Abschnitte den oberen oder unteren Rand der aktuellen Überlappung erreicht. Die Summe aller absoluten Höhenänderungen ergibt
$$F(3,5,3,7)\approx 7.01772,$$
genau wie in der numerischen Selbstprüfung der Implementierungen.
Die Implementierung wertet \(H\) mit Gleitkommaarithmetik aus und klemmt winzige Rundungsabweichungen wieder in das Intervall \([0,1]\). Für jedes Parameterpaar wird die Kandidatenmenge aus den beiden Ableitungsfaktoren aufgebaut, sortiert, mit einer engen Toleranz dedupliziert, auf Vorzeichenwechsel von \(H'\) getestet und anschließend nur an den behaltenen Stellen ausgewertet. Das Ergebnis ist die kompakte Wendewertfolge.
Für zwei Wendewertfolgen speichert die Implementierung nur den aktuellen Segmentindex beider Folgen, die aktuelle Höhe, die Bewegungsrichtung und die bisher aufsummierte Variation. In jedem Schritt wird die nächste gemeinsame Höhe mit den Überlappungsformeln bestimmt, die absolute Höhenänderung addiert, jeder getroffene Segmentrand vorgerückt und die Richtung umgedreht. Verlässt der Lauf jemals den gültigen Segmentbereich, wird das Ergebnis als ungültig behandelt; sonst endet er bei Höhe \(1\).
Zur Berechnung von \(G(m,n)\) erzeugt die Implementierung die ungeraden Primzahlen in \([m,n]\) mit einem Sieb und durchläuft dann alle Paare \(p<q\). Jedes Paar liefert die Pfadvariation zwischen \(H(p,q,\cdot)\) und \(H(p,2q-p,\cdot)\). Das C++-Programm kann diese Paarliste auf mehrere Worker-Threads verteilen. Die Java-Version bildet dieselbe Numerik in einer Laufzeitumgebung nach. Der Python-Einstieg ist nur eine dünne Hülle, die den C++-Löser baut und startet; deshalb stimmen alle drei Zugänge numerisch überein. Die Prüfwerte \(F(7,17,9,19)\approx 26.79578\) und \(G(3,20)\approx 463.80866\) sichern ab, dass Wendepunktlogik und Überlappungslauf korrekt zusammenspielen.
Sei \(P\) die Anzahl ungerader Primzahlen in \([m,n]\). Dann enthält die äußere Summe \(\binom{P}{2}\) Primzahlpaare. Für eine einzelne Kurve \(H(a,b,\cdot)\) ist die Zahl der kritischen Kandidaten \(O(a+b+|b-a|)=O(\max(a,b))\), und der anschließende Überlappungslauf ist linear in den kombinierten Längen der beiden Wendewertfolgen. Ein einzelnes Paar kostet also lineare Arbeit in der Größenordnung seiner Parameter.
Ist \(N\) die obere Intervallgrenze, erhält man damit als sichere Schranke \(O(P^2N)\) Zeit. Im festen Project-Euler-Fall liegen die Parameter zwischen 500 und 1000; praktisch dominiert also der quadratische Durchlauf über alle Primzahlpaare. Der Speicherbedarf bleibt klein: Das Sieb ist linear in der Intervallgrenze, und jeder Worker benötigt nur einige temporäre Wendewertarrays sowie eine partielle Summe.
Bu problemin merkezinde şu fonksiyon ailesi vardır:
$$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.$$
\([m,n]\) aralığındaki her \(p<q\) tek asal çifti için \(F(p,q,p,2q-p)\), \(H(p,q,\cdot)\) eğrisi ile \(H(p,2q-p,\cdot)\) eğrisi karşılaştırıldığında ortaya çıkan dikey sarılma miktarını ölçer. Son hedef
$$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p)$$
toplamıdır ve uygulamalar Project Euler için istenen \(G(500,1000)\) değerini beş ondalık basamakla üretir.
Kilit nokta şu: bu büyüklük eğrileri çok sık örnekleyerek hesaplanmıyor. Her eğri, gerçek dönüm değerlerinin dizisine indirgeniyor; ardından yol fonksiyoneli, ardışık monoton yükseklik aralıklarının kesişimleri üzerinden yeniden kuruluyor.
Çözüm üç problem-özel gözleme dayanır: \(H\)'nin analitik yapısı, gerçek ekstrem noktalarının nasıl ayıklandığı ve iki dönüm-değeri dizisi üzerinde çalışan deterministik yürüyüş.
\(a\) ve \(b\) tek tam sayılar olduğunda eğri her zaman yükseklik \(0\)'dan başlayıp \(1\)'de biter:
$$H(a,b,0)=0,\qquad H(a,b,1)=1.$$
Bunun nedeni \(\cos 0=1\) ve tek \(k\) için \(\cos(k\pi)=-1\) olmasıdır. Dolayısıyla toplamda yer alan her eğri yükseklik resminde alt uçtan üst uca gider.
Türev
$$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr)$$
şeklindedir. \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\) özdeşliğiyle
$$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)$$
elde edilir. Bu yüzden türevin işareti tamamen
$$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right)$$
çarpanı tarafından belirlenir.
İki trigonometrik çarpanın kökleri, \((0,1)\) içindeki bütün kritik nokta adaylarını verir:
$$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$$
ve \(a\ne b\) ise ayrıca
$$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$$
Uygulama bu adayları toplar, uç noktalar \(0\) ve \(1\)'i ekler, listeyi sıralar ve tekrarları siler. Böylece \(H'\)'nin sıfır olabileceği bütün yerler bulunur; fakat bunların hepsi gerçek ekstremum değildir.
Türevin bir kökü ancak işaret değiştiriyorsa önemlidir. İki trigonometrik çarpan aynı yerde sıfır olursa türev sıfıra değip yön değiştirmeyebilir; böyle bir nokta sarılan yolun yeni bir bölümünü oluşturmaz.
Bu yüzden uygulama, art arda gelen adaylar arasındaki her açık aralıkta \(H'\)'nin işaretini kontrol eder ve yalnızca sol ve sağ işaretlerin farklı olduğu adayları tutar. Böylece analitik kritik nokta kümesi, gerçek dönüm noktalarının çok daha küçük kümesine indirgenir.
\(H(a,b,\cdot)\) için korunan dönüm noktaları \(x_0<x_1<\dots<x_r\) olsun. Saklanan şey yalnızca bunlardaki yüksekliklerdir:
$$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$$
Ardışık her \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) çifti eğrinin bir monoton parçasını temsil eder. Sarılma hesabı için bu parçanın içindeki tam \(x\)-konumu gerekli değildir; önemli olan yalnızca erişilebilir yükseklik aralığıdır:
$$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].$$
Problemi çözümleyen sıkıştırma tam burada yapılır. Dönüm-değeri dizisi bilindiğinde, geri kalan geometri yalnızca bu yükseklik aralıklarının kesişimleriyle ifade edilebilir.
İki eğri alalım; dönüm-değeri dizileri \(z^A\) ve \(z^B\) olsun. Her iki dizinin ilk parçasında, ortak yükseklik \(z_{\mathrm{cur}}=0\) noktasından başlanır ve ilk hareket yukarı yönlü kabul edilir. Aktif yükseklik aralıkları \(I_i^A\) ve \(I_j^B\) ise, bir sonraki ortak yükseklik yalnızca bu iki aralığın kesişiminden çıkar.
Yukarı giderken yürüyüş kesişimin üst ucuna çıkar:
$$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$$
Aşağı giderken ise alt uca iner:
$$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$$
Bu adımın katkısı her zaman
$$|z_{\mathrm{next}}-z_{\mathrm{cur}}|$$
olur. \(z_{\mathrm{next}}\)'e ulaşıldığında hangi aktif parçanın ilgili ucu vurulduysa o parça komşu segmente ilerletilir ve yön ters çevrilir. Bu kural iki eğri de son parçasını bitirene kadar uygulanınca toplam yol varyasyonu \(F\) elde edilir.
Algoritmanın temel invarianti şudur: mevcut ortak yükseklik her zaman iki aktif aralığın kesişiminde kalır. Bu sayede dönüm-değeri dizileri sabitlenince yürüyüş de tek anlamlı hale gelir.
\((a,b)=(3,5)\) için türevin çarpanlara ayrılması şu iç adayları verir:
$$x=\frac14,\ \frac12,\ \frac34.$$
Fakat \(x=\frac12\) yalnızca teğetsel bir sıfırdır; türev burada işaret değiştirmez. Bu nedenle gerçek dönüm noktaları
$$0,\ \frac14,\ \frac34,\ 1$$
olur ve karşılık gelen dönüm değerleri
$$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1$$
şeklindedir. Bu örnek, işaret-değişimi filtresinin neden vazgeçilmez olduğunu doğrudan gösterir.
Bu denetimdeki ikinci eğrinin dönüm-değeri dizisi
$$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1$$
olur. Kesişim yürüyüşü de şu şekilde başlar:
$$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1.$$
Aktif segmentlerden biri geçerli kesişimin üst ya da alt ucuna her ulaştığında yön değişir. Tüm mutlak yükseklik farkları toplandığında
$$F(3,5,3,7)\approx 7.01772$$
elde edilir; bu değer uygulamalardaki sayısal kontrolle birebir aynıdır.
Uygulama \(H\)'yi kayan noktalı aritmetik ile hesaplar ve çok küçük yuvarlama taşmalarını yeniden \([0,1]\) aralığına sıkıştırır. Her parametre çifti için iki türev çarpanından aday küme oluşturulur, sıralanır, çok yakın tekrarlar temizlenir, ardışık adaylar arasındaki aralıklarda türev işareti örneklenir ve yalnızca işaret değiştiren noktalar tutulur. Bu noktalarda \(H\)'nin hesaplanması kompakt dönüm-değeri dizisini verir.
İki dönüm-değeri dizisi verildiğinde uygulama yalnızca her iki dizideki aktif segment indekslerini, mevcut yüksekliği, yönü ve biriken toplam varyasyonu tutar. Her adımda üstteki kesişim ya da alttaki kesişim formülüyle yeni ortak yükseklik bulunur, mutlak fark eklenir, ilgili segment uçlarına çarpan taraflar ilerletilir ve yön ters çevrilir. Yürüyüş geçerli segment aralığının dışına çıkarsa sonuç geçersiz kabul edilir; aksi halde yükseklik \(1\)'de biter.
\(G(m,n)\)'yi hesaplamak için uygulama önce \([m,n]\) içindeki tek asalları bir elek ile üretir, sonra bütün \(p<q\) çiftlerini dolaşır. Her çift, \(H(p,q,\cdot)\) ile \(H(p,2q-p,\cdot)\) arasındaki yol varyasyonunu katkı olarak ekler. C++ programı bunu doğrudan yapar ve isterse çift listesini iş parçacıklarına bölebilir. Java sürümü aynı sayısal yöntemi tek bir çalışma zamanında tekrarlar. Python girişi ise C++ çözücüyü derleyip çalıştıran ince bir kabuktur; bu yüzden üç giriş de aynı sayısal davranışı paylaşır. \(F(7,17,9,19)\approx 26.79578\) ve \(G(3,20)\approx 463.80866\) denetimleri, dönüm noktası mantığı ile kesişim yürüyüşünün uyumlu olduğunu doğrular.
\(P\), \([m,n]\) aralığındaki tek asal sayısı olsun. Dış toplamda \(\binom{P}{2}\) asal çifti vardır. Tek bir \(H(a,b,\cdot)\) eğrisi için kritik aday sayısı \(O(a+b+|b-a|)=O(\max(a,b))\) mertebesindedir; ardından gelen kesişim yürüyüşü de iki dönüm-değeri dizisinin toplam uzunluğunda doğrusaldır. Dolayısıyla bir çiftin maliyeti, parametre ölçeğinde doğrusaldır.
Aralığın üst sınırı \(N\) ise güvenli üst sınır \(O(P^2N)\) zamandır. Sabit Project Euler girdisinde parametreler 500 ile 1000 arasındadır; pratikte baskın maliyet, tek tek yol hesabından çok asal çiftleri üzerindeki karesel taramadır. Bellek kullanımı küçüktür: elek aralığın üst sınırında lineerdir, her çalışan ise yalnızca birkaç geçici dönüm-değeri dizisi ve bir kısmi toplam taşır.
La tarea gira en torno a la familia de funciones
$$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.$$
Para cada par de primos impares \(p<q\) dentro del intervalo \([m,n]\), la cantidad \(F(p,q,p,2q-p)\) mide la variación vertical acumulada cuando se compara la curva \(H(p,q,\cdot)\) con la curva \(H(p,2q-p,\cdot)\). El objetivo final es
$$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p),$$
y las implementaciones evalúan el caso de Project Euler \(G(500,1000)\), mostrando la respuesta con cinco decimales.
La idea clave es que la cantidad de “winding” no se obtiene muestreando las curvas en una malla fina. Cada curva se reduce a la sucesión de sus verdaderos valores de giro, y a partir de esa compresión el funcional de camino se reconstruye usando las intersecciones de los intervalos de altura de los tramos monótonos consecutivos.
La solución se apoya en tres hechos específicos del problema: la forma analítica de \(H\), la caracterización exacta de sus extremos reales y un recorrido determinista sobre dos secuencias de valores de giro.
Si \(a\) y \(b\) son impares, la curva siempre empieza en altura \(0\) y termina en altura \(1\):
$$H(a,b,0)=0,\qquad H(a,b,1)=1.$$
Esto se deduce de \(\cos 0=1\) y de \(\cos(k\pi)=-1\) cuando \(k\) es impar. En otras palabras, cada curva del problema conecta el nivel inferior con el superior.
Derivando obtenemos
$$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr).$$
La identidad \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\) la convierte en
$$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).$$
Por tanto, el signo de la derivada queda controlado por
$$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right).$$
Los ceros de los dos factores trigonométricos producen todos los candidatos críticos dentro de \((0,1)\):
$$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$$
y, si \(a\ne b\), además
$$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$$
La implementación reúne esos candidatos, añade \(0\) y \(1\), los ordena y elimina duplicados. En ese momento conocemos todos los puntos donde \(H'\) puede anularse, pero todavía no sabemos cuáles son extremos verdaderos.
Un cero de la derivada solo importa si el signo de \(H'\) cambia a través de él. Si los dos factores trigonométricos se anulan en la misma posición, la derivada puede tocar cero sin cambiar de dirección, y ese punto no genera un nuevo tramo para el camino sinuoso.
Por eso la implementación examina cada intervalo abierto entre candidatos consecutivos, calcula allí el signo de \(H'\) y conserva únicamente los candidatos donde el signo de la izquierda difiere del signo de la derecha. Así se pasa del conjunto analítico de candidatos al conjunto mucho más pequeño de giros reales.
Sean \(x_0<x_1<\dots<x_r\) las posiciones retenidas para \(H(a,b,\cdot)\). El solver guarda solo las alturas
$$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$$
Cada par consecutivo \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) representa un tramo monótono de la curva. Para el cálculo de la trayectoria ya no hace falta la coordenada \(x\) exacta dentro de ese tramo; solo importa el intervalo de alturas alcanzables
$$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].$$
Aquí está la compresión decisiva del problema: una vez conocida la sucesión de valores de giro, el resto de la geometría se puede describir únicamente mediante intersecciones de esos intervalos de altura.
Tomemos dos curvas con secuencias de valores de giro \(z^A\) y \(z^B\). Se empieza en el primer tramo de cada una, con altura compartida \(z_{\mathrm{cur}}=0\), y el primer movimiento es ascendente. Si los intervalos activos son \(I_i^A\) e \(I_j^B\), la siguiente altura compartida viene dada únicamente por su intersección.
Si el movimiento actual es ascendente, el camino sube hasta la parte superior de la intersección:
$$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$$
Si el movimiento es descendente, baja hasta la parte inferior:
$$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$$
La contribución del paso es siempre
$$|z_{\mathrm{next}}-z_{\mathrm{cur}}|.$$
Después de alcanzar \(z_{\mathrm{next}}\), todo tramo cuyo extremo relevante haya sido tocado avanza al tramo vecino, y la dirección se invierte. Repitiendo esta regla determinista hasta que ambas curvas agoten su último tramo se obtiene \(F\).
La invariante esencial es que la altura actual siempre permanece dentro de la intersección de los dos intervalos activos. Por eso, fijadas las dos secuencias de giro, el recorrido ya no tiene ambigüedad.
Para \((a,b)=(3,5)\), la factorización de la derivada da como candidatos interiores
$$x=\frac14,\ \frac12,\ \frac34.$$
Sin embargo, \(x=\frac12\) es solo un cero tangencial: la derivada no cambia de signo allí. Por tanto, las posiciones reales de giro son
$$0,\ \frac14,\ \frac34,\ 1,$$
con valores de giro
$$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1.$$
Este ejemplo muestra con claridad por qué el filtro por cambio de signo es indispensable: el conjunto bruto de ceros de \(H'\) todavía no es el estado correcto para el problema.
La segunda curva de ese control produce la secuencia
$$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1.$$
El recorrido por intersecciones empieza entonces como
$$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1,$$
alternando pasos hacia arriba y hacia abajo cada vez que uno de los tramos activos alcanza la cota superior o inferior de la intersección actual. La suma de todas las variaciones absolutas de altura es
$$F(3,5,3,7)\approx 7.01772,$$
exactamente el valor que usan las implementaciones como verificación numérica.
La implementación evalúa \(H\) en coma flotante y recorta pequeños desbordamientos numéricos de vuelta a \([0,1]\). Para cada par de parámetros construye el conjunto de candidatos a partir de los dos factores de la derivada, lo ordena, elimina casi duplicados, muestrea el signo de la derivada entre candidatos consecutivos y conserva solo los puntos donde sí hay cambio de signo. Evaluar \(H\) en esos puntos retenidos produce la secuencia compacta de valores de giro.
Dadas dos secuencias de valores de giro, la implementación solo necesita el índice del tramo activo en cada secuencia, la altura actual, la dirección actual y la variación acumulada. En cada iteración calcula la siguiente altura compartida con las fórmulas de intersección, añade la diferencia absoluta de altura, avanza los extremos de tramo que han sido alcanzados y cambia la dirección. Si el recorrido llegara a salir del rango válido de segmentos, el resultado se considera inválido; en caso contrario termina en la altura \(1\).
Para evaluar \(G(m,n)\), la implementación genera con una criba todos los primos impares de \([m,n]\) y recorre todos los pares \(p<q\). Cada par aporta la variación de trayectoria entre \(H(p,q,\cdot)\) y \(H(p,2q-p,\cdot)\). El programa en C++ hace esto de forma directa y puede repartir la lista de pares entre varios hilos de trabajo. La versión Java reproduce el mismo método numérico en una sola ejecución. La entrada de Python es un envoltorio mínimo que compila y ejecuta el solver en C++, por lo que hereda exactamente el mismo comportamiento numérico. Los controles \(F(7,17,9,19)\approx 26.79578\) y \(G(3,20)\approx 463.80866\) comprueban que la extracción de giros y el recorrido por intersecciones están sincronizados.
Sea \(P\) el número de primos impares en \([m,n]\). Entonces la suma exterior contiene \(\binom{P}{2}\) pares. Para una curva \(H(a,b,\cdot)\), el número de candidatos críticos es \(O(a+b+|b-a|)=O(\max(a,b))\), y el recorrido posterior es lineal en la longitud combinada de las dos secuencias de giro. Así, el costo de un solo par es lineal en la escala de sus parámetros.
Si \(N\) es el extremo superior del intervalo, una cota segura es \(O(P^2N)\) tiempo. En la instancia fija de Project Euler, los parámetros están entre 500 y 1000, de modo que en la práctica domina el barrido cuadrático sobre pares de primos. El uso de memoria es reducido: la criba es lineal en el límite del intervalo y cada trabajador solo mantiene unas pocas secuencias temporales de giro y una suma parcial.
这道题围绕如下函数族展开：
$$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.$$
对于区间 \([m,n]\) 中每一对奇素数 \(p<q\)，量 \(F(p,q,p,2q-p)\) 描述了把 \(H(p,q,\cdot)\) 与 \(H(p,2q-p,\cdot)\) 放在一起比较时，公共高度路径所产生的总“绕行”量。最终要求的是
$$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p),$$
而实现代码计算的是 Project Euler 指定的目标 \(G(500,1000)\)，并以五位小数输出结果。
真正重要的点在于：这个量不是靠对曲线做高密度采样得到的。每条曲线都会先被压缩成“真实转折值”序列，然后再根据相邻单调高度区间之间的重叠关系，把整条 winding path 重新构造出来。
整个解法依赖三个与本题直接相关的对象：\(H\) 的解析结构、如何从导数零点中筛出真正的极值点，以及如何在两条转折值序列之间做一个确定性的交替行走。
当 \(a,b\) 都是奇整数时，曲线总是从高度 \(0\) 出发，在高度 \(1\) 结束：
$$H(a,b,0)=0,\qquad H(a,b,1)=1.$$
原因很直接：\(\cos 0=1\)，而奇数倍的 \(\pi\) 满足 \(\cos(k\pi)=-1\)。因此本题中出现的每一条曲线都把高度区间 \([0,1]\) 的两个端点连接起来。
对 \(H\) 求导可得
$$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr).$$
再利用恒等式 \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\)，就得到
$$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).$$
前面的系数始终为正，所以导数符号完全由下面这个乘积决定：
$$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right).$$
上面两个三角因子的零点给出了 \((0,1)\) 内全部临界点候选：
$$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$$
如果 \(a\ne b\)，还要再加上
$$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$$
实现先把这些点收集起来，再加入端点 \(0\) 和 \(1\)，排序并去重。做到这一步时，我们已经知道导数可能为零的所有位置，但还不知道哪些位置真的是曲线方向发生改变的地方。
导数为零并不自动意味着它是有用的转折点。只有当 \(H'\) 的符号在该点两侧发生变化时，这个点才会改变路径结构。如果两个三角因子在同一点同时为零，导数可能只是轻轻碰到零，却没有真正改变上升或下降方向。
因此实现会查看相邻候选点之间每个开区间内的导数符号，只保留那些左右符号不同的候选点。这样，解析上得到的“大候选集”就被压缩成了真正影响 winding path 的“小转折集”。
设 \(H(a,b,\cdot)\) 保留下来的转折位置是 \(x_0<x_1<\dots<x_r\)。程序并不继续保留全部函数图像，而只保留这些位置上的高度
$$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$$
相邻的两个值 \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) 就表示了曲线的一个单调片段。对于后续的路径计算，这个片段内部的精确 \(x\) 坐标已经不重要了；真正重要的是它能覆盖的高度区间
$$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].$$
这一步就是本题的核心压缩：一旦转折值序列确定，剩余的几何结构就可以完全用这些高度区间之间的重叠关系来描述。
现在考虑两条曲线，它们的转折值序列分别记作 \(z^A\) 与 \(z^B\)。从两条序列的第一个片段出发，初始共享高度取 \(z_{\mathrm{cur}}=0\)，并规定第一步向上。若当前活动区间是 \(I_i^A\) 和 \(I_j^B\)，那么下一次公共高度完全由这两个区间的交叠部分决定。
如果当前是向上走，那么下一高度是交叠部分的上端点：
$$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$$
如果当前是向下走，那么下一高度是交叠部分的下端点：
$$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$$
这一小步对函数值的贡献始终是
$$|z_{\mathrm{next}}-z_{\mathrm{cur}}|.$$
到达 \(z_{\mathrm{next}}\) 之后，哪个片段刚好碰到了对应的端点，就把哪个片段推进到相邻片段，然后把方向翻转。不断重复，直到两条曲线都走完最后一个片段，就得到了总路径变差 \(F\)。
这个算法的关键不变量是：当前共享高度始终落在两个活动区间的重叠部分里。因此，只要两条转折值序列固定，整条行走路径就是唯一的，不存在分叉选择。
当 \((a,b)=(3,5)\) 时，由导数分解得到的内部候选点是
$$x=\frac14,\ \frac12,\ \frac34.$$
但是 \(x=\frac12\) 只是一个“切到零”的点，不会引起导数符号变化，因此必须删掉。真正保留下来的转折位置是
$$0,\ \frac14,\ \frac34,\ 1,$$
对应的转折值为
$$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1.$$
这个例子精确说明了为什么不能把导数零点原封不动地拿来用。真正决定路径结构的是“发生方向变化的零点”，而不是“所有零点”。
该校验中的第二条曲线对应的转折值序列是
$$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1.$$
于是交替行走一开始会经过
$$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1.$$
每当某个活动片段撞到当前交叠区间的上界或下界，方向就翻转一次。把整个过程中所有高度变化的绝对值累加起来，就得到
$$F(3,5,3,7)\approx 7.01772,$$
这与实现中的数值校验完全一致。
实现使用浮点数计算 \(H\)，并把非常小的舍入越界重新钳回 \([0,1]\)。对每组参数，它先根据导数的两个因子生成候选点集合，排序、去重，再在相邻候选点之间采样导数符号，只保留真正发生符号翻转的位置。最后在这些保留位置上求 \(H\) 的值，得到紧凑的转折值序列。
有了两条转折值序列之后，程序只维护四类状态：两边当前所在的片段编号、当前公共高度、当前方向，以及累计的总变差。每一步都利用上面的交叠公式求出下一公共高度，加入绝对高度差，推进刚刚触到端点的片段，并翻转方向。如果行走过程超出了合法片段范围，就把这次计算视为失败；否则它会在高度 \(1\) 处正常结束。
为了计算 \(G(m,n)\)，实现先用筛法生成区间 \([m,n]\) 内的全部奇素数，再枚举所有 \(p<q\) 的组合。每一对都要比较 \(H(p,q,\cdot)\) 与 \(H(p,2q-p,\cdot)\) 的路径变差。C++ 程序直接执行这套流程，并且可以把素数对列表切分给多个工作线程。Java 版本在自己的运行时里复现相同的数值步骤。Python 入口则是一个很薄的包装层，它负责构建并运行 C++ 求解器，因此三种入口共享同一套最终数值行为。实现里还用 \(F(7,17,9,19)\approx 26.79578\) 和 \(G(3,20)\approx 463.80866\) 作为校验点，确保转折点提取与路径行走逻辑没有偏移。
设 \([m,n]\) 中奇素数的个数为 \(P\)。那么外层求和共有 \(\binom{P}{2}\) 对素数。对于单条曲线 \(H(a,b,\cdot)\)，临界点候选数是 \(O(a+b+|b-a|)=O(\max(a,b))\)，之后的交替行走则与两条转折值序列的总长度成线性关系。因此，单个素数对的代价在线性数量级内。
若区间上界记为 \(N\)，则一个安全的总时间上界是 \(O(P^2N)\)。在固定的 Project Euler 输入里，参数都落在 500 到 1000 之间，所以实际运行时主要由素数对的二次枚举主导，而不是由某一次单独的路径计算主导。空间使用很节制：筛法需要关于区间上界的线性布尔数组，每个活动工作单元只需要少量临时转折值数组和一个部分和。
Задача построена вокруг семейства функций
$$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.$$
Для каждой пары нечётных простых \(p<q\) из интервала \([m,n]\) величина \(F(p,q,p,2q-p)\) измеряет суммарное вертикальное “закручивание”, возникающее при сравнении кривых \(H(p,q,\cdot)\) и \(H(p,2q-p,\cdot)\). Итоговая сумма имеет вид
$$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p),$$
и реализации вычисляют именно целевой для Project Euler случай \(G(500,1000)\), печатая ответ с точностью до пяти знаков после запятой.
Главная идея решения состоит в том, что нужный функционал не вычисляется по плотной сетке значений. Каждая кривая сначала сжимается до последовательности своих настоящих значений в точках поворота, а затем весь winding path восстанавливается только по перекрытиям соседних монотонных интервалов по высоте.
Решение опирается на три конкретных для этой задачи факта: аналитическую форму \(H\), точное описание настоящих экстремумов и детерминированный проход по двум последовательностям turning-values.
Если \(a\) и \(b\) нечётны, кривая всегда стартует в высоте \(0\) и заканчивается в высоте \(1\):
$$H(a,b,0)=0,\qquad H(a,b,1)=1.$$
Это следует из равенств \(\cos 0=1\) и \(\cos(k\pi)=-1\) для нечётного \(k\). Значит, каждая кривая в задаче соединяет нижнюю и верхнюю границы диапазона высот.
Производная равна
$$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr).$$
Тождество \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\) даёт
$$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).$$
Следовательно, знак производной полностью определяется выражением
$$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right).$$
Нули этих двух тригонометрических множителей дают все кандидаты на критические точки внутри \((0,1)\):
$$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$$
и, если \(a\ne b\), дополнительно
$$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$$
Реализация собирает эти точки, добавляет концы \(0\) и \(1\), сортирует список и удаляет дубликаты. После этого известны все места, где \(H'\) может обращаться в ноль, но ещё не ясно, какие из них являются настоящими экстремумами.
Нуль производной важен только тогда, когда знак \(H'\) действительно меняется. Если оба тригонометрических множителя зануляются в одной точке, производная может коснуться нуля и не изменить направление движения. Такая точка не создаёт нового звена winding path.
Поэтому реализация проверяет знак \(H'\) на каждом открытом промежутке между соседними кандидатами и сохраняет только те точки, где знаки слева и справа различаются. Так аналитическое множество кандидатов сжимается до существенно меньшего множества настоящих точек поворота.
Пусть \(x_0<x_1<\dots<x_r\) - оставшиеся точки поворота функции \(H(a,b,\cdot)\). Хранятся только высоты
$$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$$
Каждая соседняя пара \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) задаёт один монотонный участок кривой. Для вычисления winding path точные координаты \(x\) внутри участка уже не нужны; достаточно знать только интервал достижимых высот
$$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].$$
Именно здесь происходит ключевое сжатие задачи. Как только известна последовательность turning-values, оставшаяся геометрия описывается исключительно перекрытиями этих интервалов высот.
Рассмотрим две кривые с последовательностями \(z^A\) и \(z^B\). Стартуем на первом сегменте каждой кривой при общей высоте \(z_{\mathrm{cur}}=0\), причём первый шаг направлен вверх. Если текущие активные интервалы по высоте равны \(I_i^A\) и \(I_j^B\), то следующая общая высота определяется только их перекрытием.
Для шага вверх траектория поднимается к верхней границе перекрытия:
$$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$$
Для шага вниз она опускается к нижней границе:
$$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$$
Вклад такого шага всегда равен
$$|z_{\mathrm{next}}-z_{\mathrm{cur}}|.$$
После достижения \(z_{\mathrm{next}}\) каждый сегмент, у которого был достигнут нужный конец, переходит к соседнему сегменту, а направление шага меняется на противоположное. Повторяя это правило до тех пор, пока обе кривые не исчерпают последний сегмент, мы получаем полную path variation \(F\).
Ключевой инвариант состоит в том, что текущая общая высота всегда лежит в перекрытии двух активных интервалов. Поэтому после фиксации двух последовательностей turning-values весь проход становится однозначным.
Для \((a,b)=(3,5)\) факторизация производной даёт внутренние кандидаты
$$x=\frac14,\ \frac12,\ \frac34.$$
Однако \(x=\frac12\) - лишь касательный нуль: знак производной там не меняется. Поэтому настоящими точками поворота остаются
$$0,\ \frac14,\ \frac34,\ 1,$$
а соответствующие значения равны
$$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1.$$
Этот пример наглядно показывает, почему недостаточно просто взять все нули \(H'\): структуру пути определяют только те нули, в которых действительно происходит смена направления.
Вторая кривая в этом контроле имеет последовательность значений поворота
$$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1.$$
Поэтому проход по перекрытиям начинается так:
$$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1.$$
Каждый раз, когда активный сегмент достигает верхней или нижней границы текущего перекрытия, направление разворачивается. Сумма всех абсолютных изменений высоты даёт
$$F(3,5,3,7)\approx 7.01772,$$
точно совпадая с численным контрольным значением из реализации.
Реализация вычисляет \(H\) в плавающей точке и возвращает к \([0,1]\) крошечные выходы за границы, возникающие из-за округления. Для каждой пары параметров строится набор кандидатов из двух множителей производной, он сортируется, очищается от почти одинаковых значений, после чего на промежутках между соседними кандидатами проверяется знак производной. Сохраняются только точки со сменой знака, а значения \(H\) в них образуют компактную последовательность turning-values.
Когда две последовательности turning-values уже построены, реализация хранит только индексы активных сегментов в обеих последовательностях, текущую высоту, текущее направление и накопленную сумму variation. На каждом шаге она вычисляет следующую общую высоту по формулам перекрытия, добавляет абсолютное изменение высоты, продвигает достигнутые концы сегментов и меняет направление. Если проход когда-либо выходит за допустимый диапазон сегментов, результат считается недопустимым; в нормальном случае процесс заканчивается на высоте \(1\).
Чтобы вычислить \(G(m,n)\), реализация сначала решетом строит список нечётных простых из \([m,n]\), а затем перебирает все пары \(p<q\). Каждая пара даёт вклад, равный variation пути между \(H(p,q,\cdot)\) и \(H(p,2q-p,\cdot)\). Программа на C++ делает это напрямую и умеет делить список пар между рабочими потоками. Версия на Java повторяет тот же численный алгоритм в своей среде выполнения. Python-вход представляет собой тонкую оболочку, которая собирает и запускает C++-решатель, так что все три варианта фактически опираются на одну и ту же численную схему. Контрольные значения \(F(7,17,9,19)\approx 26.79578\) и \(G(3,20)\approx 463.80866\) используются как проверка того, что извлечение поворотных точек и проход по перекрытиям согласованы.
Пусть \(P\) - число нечётных простых в интервале \([m,n]\). Тогда внешняя сумма содержит \(\binom{P}{2}\) пар. Для одной кривой \(H(a,b,\cdot)\) число кандидатов в критические точки равно \(O(a+b+|b-a|)=O(\max(a,b))\), а сам проход по перекрытиям линеен по суммарной длине двух последовательностей turning-values. Значит, стоимость одной пары линейна по масштабу её параметров.
Если верхняя граница интервала равна \(N\), безопасная верхняя оценка времени работы составляет \(O(P^2N)\). В фиксированном случае Project Euler параметры лежат между 500 и 1000, поэтому на практике доминирует квадратичный перебор пар простых, а не вычисление какого-то одного пути. Память расходуется умеренно: решето требует линейного массива по верхней границе, а каждому активному worker достаточно нескольких временных массивов turning-values и одной частичной суммы.
تدور المسألة حول عائلة الدوال
$$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.$$
لكل زوج من الأوليات الفردية \(p<q\) داخل المجال \([m,n]\)، تقيس الكمية \(F(p,q,p,2q-p)\) مقدار الالتفاف الرأسي المتراكم عند مقارنة المنحنى \(H(p,q,\cdot)\) بالمنحنى \(H(p,2q-p,\cdot)\). والمطلوب في النهاية هو
$$G(m,n)=\sum_{\substack{p<q\\ p,q\in\mathbb{P}_{\mathrm{odd}}\cap [m,n]}} F(p,q,p,2q-p),$$
وتحسب التطبيقات الحالة المطلوبة في Project Euler وهي \(G(500,1000)\)، ثم تطبع النتيجة بخمس خانات عشرية.
الفكرة الحاسمة هي أن هذه الكمية لا تُحسب عبر أخذ عينات كثيفة من المنحنيات. كل منحنى يُختزل أولاً إلى سلسلة قيمه عند نقاط الانعطاف الحقيقية، ثم يُعاد بناء functional المسار بالكامل من خلال تداخل مجالات الارتفاع الخاصة بالمقاطع الأحادية المتتالية.
يعتمد الحل على ثلاثة عناصر خاصة بهذه المسألة: البنية التحليلية للدالة \(H\)، والطريقة الصحيحة لاستخراج القيم القصوى الحقيقية، ثم مسار حتمي يتنقل بين سلسلتي قيم الانعطاف لمنحنيين.
إذا كان \(a\) و\(b\) عددين فرديين، فإن المنحنى يبدأ دائماً من الارتفاع \(0\) وينتهي عند الارتفاع \(1\):
$$H(a,b,0)=0,\qquad H(a,b,1)=1.$$
والسبب أن \(\cos 0=1\)، وأن \(\cos(k\pi)=-1\) عندما يكون \(k\) فردياً. ولذلك فكل منحنى في هذه المسألة يصل بين طرفي مجال الارتفاع \([0,1]\).
باشتقاق \(H\) نحصل على
$$H'(a,b,x)=\frac{\pi ab}{2(a+b)}\bigl(\sin(a\pi x)+\sin(b\pi x)\bigr).$$
وباستخدام الهوية \(\sin u+\sin v=2\sin\!\bigl(\frac{u+v}{2}\bigr)\cos\!\bigl(\frac{u-v}{2}\bigr)\) تصبح
$$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).$$
ولأن المعامل الأمامي موجب دائماً، فإن إشارة المشتقة يحكمها بالكامل العامل
$$\sin\!\left(\frac{(a+b)\pi x}{2}\right)\cos\!\left(\frac{|b-a|\pi x}{2}\right).$$
جذور العاملين المثلثيين تعطي جميع مرشحي النقاط الحرجة داخل \((0,1)\):
$$x=\frac{2k}{a+b}\quad\left(1\le k<\frac{a+b}{2}\right),$$
وإذا كان \(a\ne b\) نضيف أيضاً
$$x=\frac{2k+1}{|b-a|}\quad\left(0\le k<\frac{|b-a|}{2}\right).$$
تقوم الخوارزمية بجمع هذه النقاط، ثم تضيف النهايتين \(0\) و\(1\)، ثم ترتب القائمة وتحذف المكررات. عند هذه المرحلة نكون قد حددنا جميع المواضع التي قد تنعدم عندها المشتقة، لكننا لم نحدد بعد أيها يمثل انعطافاً حقيقياً.
ليس كل جذر للمشتقة مهماً. النقطة تكون مؤثرة فقط إذا تبدلت إشارة \(H'\) عند عبورها. فإذا انعدم العاملان المثلثيان معاً في الموضع نفسه، فقد تلامس المشتقة الصفر من دون أن تغيّر اتجاه الحركة، وعندئذ لا تتولد قطعة جديدة في مسار الالتفاف.
لهذا السبب تفحص الخوارزمية كل فترة مفتوحة بين مرشحين متجاورين، وتحدد إشارة \(H'\) داخلها، ثم تحتفظ فقط بالنقاط التي تختلف عندها الإشارة على اليسار عن الإشارة على اليمين. بهذه الطريقة ننتقل من مجموعة تحليلية كبيرة من المرشحين إلى مجموعة صغيرة من نقاط الانعطاف الحقيقية.
لنفترض أن نقاط الانعطاف المحتفَظ بها للدالة \(H(a,b,\cdot)\) هي \(x_0<x_1<\dots<x_r\). ما تحتفظ به الخوارزمية ليس الرسم الكامل للدالة، بل الارتفاعات فقط:
$$z_i^{(a,b)}=H(a,b,x_i)\qquad (0\le i\le r).$$
كل زوج متتالٍ \((z_i^{(a,b)},z_{i+1}^{(a,b)})\) يمثل مقطعاً أحادياً من المنحنى. وفي حساب الالتفاف لا تعود إحداثيات \(x\) الدقيقة داخل هذا المقطع مهمة؛ المهم فقط هو مجال الارتفاعات الذي يغطيه هذا المقطع:
$$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].$$
وهنا يحدث الاختزال الحاسم للمسألة: بمجرد معرفة سلسلة قيم الانعطاف يمكن وصف بقية الهندسة بالكامل من خلال تداخل هذه المجالات الرأسية.
لنأخذ منحنيين لهما سلسلتا قيم انعطاف \(z^A\) و\(z^B\). نبدأ من المقطع الأول في كل منهما عند الارتفاع المشترك \(z_{\mathrm{cur}}=0\)، ونفترض أن أول حركة تكون إلى الأعلى. إذا كان المجالان النشطان هما \(I_i^A\) و\(I_j^B\)، فإن الارتفاع المشترك التالي تحدده منطقة التداخل بينهما فقط.
في حالة الحركة الصاعدة، يصعد المسار إلى الحد الأعلى للتداخل:
$$z_{\mathrm{next}}=\min\!\bigl(\max I_i^A,\max I_j^B\bigr).$$
وفي حالة الحركة الهابطة، يهبط إلى الحد الأدنى للتداخل:
$$z_{\mathrm{next}}=\max\!\bigl(\min I_i^A,\min I_j^B\bigr).$$
ومساهمة هذه الخطوة هي دائماً
$$|z_{\mathrm{next}}-z_{\mathrm{cur}}|.$$
بعد الوصول إلى \(z_{\mathrm{next}}\)، فإن كل مقطع لامس طرفه المناسب يتقدم إلى المقطع المجاور، ثم تنعكس جهة الحركة. وبالاستمرار في هذه القاعدة الحتمية حتى يستنفد المنحنيان آخر مقطع فيهما نحصل على مقدار path variation الكلي \(F\).
الثابت الأساسي هنا هو أن الارتفاع الحالي يبقى دائماً داخل منطقة التداخل بين المجالين النشطين. ولهذا يكون المسار محدداً بشكل وحيد بمجرد تثبيت سلسلتي قيم الانعطاف.
عند \((a,b)=(3,5)\) يعطي تحليل المشتقة المرشحين الداخليين
$$x=\frac14,\ \frac12,\ \frac34.$$
لكن النقطة \(x=\frac12\) ليست إلا جذراً تماسياً؛ فإشارة المشتقة لا تتغير عندها. لذلك تصبح نقاط الانعطاف الحقيقية هي
$$0,\ \frac14,\ \frac34,\ 1,$$
وقيمها الموافقة
$$0,\ \frac12+\frac{\sqrt2}{4},\ \frac12-\frac{\sqrt2}{4},\ 1.$$
هذا المثال يوضح بدقة لماذا لا يكفي أخذ جميع جذور \(H'\) كما هي. البنية الفعلية للمسار يحددها فقط ما يغيّر اتجاه الحركة حقاً.
ينتج عن المنحنى الثاني في هذا الفحص تسلسل قيم الانعطاف
$$0,\ 0.6545084972,\ 0.6414213562,\ 0.9045084972,\ 0.0954915028,\ 0.3585786438,\ 0.3454915028,\ 1.$$
ومن ثم يبدأ مسار التداخل على الصورة
$$0\to 0.6545084972\to 0.6414213562\to 0.8535533906\to \cdots \to 1.$$
وفي كل مرة يبلغ فيها أحد المقطعين النشطين الحد الأعلى أو الأدنى للتداخل الحالي تنقلب جهة الحركة. وبجمع جميع التغيرات المطلقة في الارتفاع نحصل على
$$F(3,5,3,7)\approx 7.01772,$$
وهي نفس قيمة الفحص الرقمية المستخدمة في التطبيقات.
تقيّم الخوارزمية الدالة \(H\) باستخدام أعداد ذات فاصلة عائمة، ثم تعيد أي انحرافات عددية صغيرة جداً إلى داخل المجال \([0,1]\). ولكل زوج من المعاملات تبني مجموعة المرشحين من عاملي المشتقة، ثم ترتبها وتحذف القيم شبه المتطابقة، وتفحص إشارة المشتقة بين كل نقطتين متتاليتين، ولا تحتفظ إلا بالنقاط التي يحدث عندها انقلاب فعلي في الإشارة. وبعد ذلك تؤدي قيمة \(H\) في هذه النقاط إلى السلسلة المضغوطة لقيم الانعطاف.
عندما تصبح سلسلتا قيم الانعطاف جاهزتين، لا يحتاج التنفيذ إلا إلى فهارس المقطعين النشطين، والارتفاع الحالي، والاتجاه الحالي، والمجموع المتراكم للvariation. في كل خطوة يحسب الارتفاع المشترك التالي من صيغ التداخل، ويضيف مقدار التغير المطلق في الارتفاع، ثم يحرّك كل طرف من أطراف المقاطع التي تم الوصول إليها ويعكس الاتجاه. وإذا خرج المسار من نطاق المقاطع الصحيح تُعد النتيجة غير صالحة؛ وإلا فإنه ينتهي طبيعياً عند الارتفاع \(1\).
لحساب \(G(m,n)\)، يولد التنفيذ جميع الأوليات الفردية في \([m,n]\) باستخدام غربال، ثم يمر على كل زوج \(p<q\). كل زوج يضيف variation المسار بين \(H(p,q,\cdot)\) و\(H(p,2q-p,\cdot)\). ينفذ برنامج C++ هذه العملية مباشرة ويمكنه توزيع قائمة الأزواج على عدة خيوط عمل. ويعيد برنامج Java الخطوات العددية نفسها داخل بيئته الخاصة. أما مدخل Python فهو غلاف صغير يبني محلل C++ ويشغله، ولذلك يشترك معهما في السلوك العددي نفسه. وتُستخدم القيمتان \(F(7,17,9,19)\approx 26.79578\) و\(G(3,20)\approx 463.80866\) للتأكد من أن استخراج نقاط الانعطاف ومسار التداخل متوافقان تماماً.
إذا كان عدد الأوليات الفردية في \([m,n]\) هو \(P\)، فإن المجموع الخارجي يحتوي على \(\binom{P}{2}\) زوجاً. وبالنسبة إلى منحنى واحد \(H(a,b,\cdot)\)، فإن عدد المرشحين للنقاط الحرجة هو \(O(a+b+|b-a|)=O(\max(a,b))\)، كما أن مسار التداخل نفسه خطي في الطول المدمج لسلسلتي قيم الانعطاف. لذلك تكون كلفة الزوج الواحد خطية في حجم معامِلاته.
إذا كانت النهاية العليا للمجال هي \(N\)، فإن حداً أعلى آمناً للزمن هو \(O(P^2N)\). وفي حالة Project Euler الثابتة تقع المعاملات بين 500 و1000، ولذلك تهيمن عملياً عملية المسح التربيعي على أزواج الأوليات، لا حساب مسار منفرد بعينه. أما الذاكرة فتبقى متواضعة: الغربال يحتاج إلى مساحة خطية في حد المجال، وكل عامل نشط يحتاج فقط إلى بضع سلاسل مؤقتة لقيم الانعطاف ومجموع جزئي واحد.