The problem asks for
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
with \(N=45000\) and \(L=10^{10}\). For one target angle \(\alpha\) measured in degrees, we consider integer right triangles coming from Pythagorean triples. Each admissible triangle determines a special angle \(\theta\), and \(f(\alpha,L)\) is the perimeter of the scaled triangle whose angle is closest to \(\alpha\), subject to the hypotenuse bound \(c\le L\) after scaling.
The implementation does not brute-force all triples. Instead it converts the geometry into a one-variable function of the ratio \(u=m/n\), inverts that function numerically, and then searches only for nearby rational values that can really occur in primitive Pythagorean triples.
The entire solution rests on two facts. First, every primitive integer right triangle has Euclid parameters \(m>n>0\). Second, the relevant angle depends only on the quotient \(u=m/n\), so the continuous optimization is one-dimensional before it is converted back to integer data.
Every primitive Pythagorean triple can be written as
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
with
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
Any non-primitive solution is just a scaled copy \((ka,kb,kc)\) of such a primitive triangle. Scaling does not change the associated angle, so once a primitive triple is chosen, the best allowed scale is simply
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$
The value returned by \(f(\alpha,L)\) is therefore
$$k(a+b+c).$$
So the mathematical task is: among primitive pairs \((m,n)\), find the one whose angle is closest to \(\alpha\), then multiply its perimeter by the largest admissible scale.
The problem-specific angle is encoded through
$$\tan\theta=\frac{3ab}{2c^2}.$$
Substituting Euclid's formulas and writing
$$u=\frac{m}{n}>1$$
gives
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$
This is the key collapse: the angle no longer depends on the absolute size of \(m\) and \(n\), only on their ratio.
For a target angle \(\alpha\), the code computes
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$
Using the tangent-difference identity, the angular error can be measured by
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$
Because \(\tan x\) is strictly increasing on the small angle range relevant here, minimizing the true angle difference is equivalent to minimizing \(\Delta\).
Differentiating yields
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$
For \(u>1\), there is a single critical point, obtained from \(u^4-6u^2+1=0\):
$$u_{\max}=1+\sqrt{2}.$$
At that point,
$$t(u_{\max})=\frac34,$$
so \(t(u)\) increases on \((1,1+\sqrt2)\) and decreases on \((1+\sqrt2,\infty)\). Therefore a target value \(t_\alpha<3/4\) has two inverse images:
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
Those are the two continuous ratios around which the search must be organized.
There is also a useful problem-specific observation: in the actual outer sum, \(\alpha=\sqrt[3]{i}\) with \(1\le i\le45000\), so
$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$
Hence every angle used in the final computation lies below the peak. The production run always has two branches, even though the implementations still handle the peak case for completeness.
A primitive triple corresponds to a reduced fraction
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
with size constraint
$$p^2+q^2\le L.$$
If \(p/q\) is close to one of the continuous roots \(u_\ast\), then \(p\approx u_\ast q\), so the bound above suggests the natural denominator scale
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
The implementations turn this into the Farey order
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$
They first find the Farey neighbors of order \(N_u\) that bracket \(u_\ast\). Then they walk outward in the Farey sequence until they encounter the nearest fractions on the lower and upper sides that satisfy all primitive-triple conditions. Those nearby admissible fractions are the concrete integer candidates that get evaluated.
This is why the solution is fast: it never scans the whole disk \(m^2+n^2\le L\). It solves the continuous inverse problem first and only then inspects a tiny local set of rational approximations around the relevant roots.
For the small checkpoint angle,
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735.$$
One excellent admissible fraction is
$$\frac{m}{n}=\frac92.$$
It generates
$$a=9^2-2^2=77,\qquad b=2\cdot 9\cdot 2=36,\qquad c=9^2+2^2=85.$$
Since \(c=85\le100\), the largest allowed scale is \(k=\lfloor 100/85\rfloor=1\). The associated tangent value is
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550,$$
so the error proxy is
$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875.$$
This candidate beats the nearby admissible alternatives, so
$$f(30,100)=77+36+85=198.$$
The implementations use this value as a validation checkpoint in the C++ version, and the same mathematics appears in all three languages.
For each \(\alpha\), the C++, Python, and Java implementations compute \(t_\alpha\). When \(t_\alpha<3/4\), they use bisection once on the increasing branch \((1,1+\sqrt2)\) and once on the decreasing branch \((1+\sqrt2,\infty)\), producing two continuous targets \(u_-\) and \(u_+\). If the target ever reached or exceeded the peak value \(3/4\), the search would collapse to the neighborhood of \(u=1+\sqrt2\).
For each continuous target \(u_\ast\), the implementation computes the corresponding Farey order \(N_u\), constructs the bracketing neighbors, and then walks outward to the closest admissible fractions on both sides. Every valid fraction is turned back into a primitive triple \((a,b,c)\), scaled by \(k=\lfloor L/c\rfloor\), and scored using the same quantity \(\Delta\).
If two candidates are numerically tied in angular error, the implementation prefers the one with larger scaled area. Since the scaled area is
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2},$$
it is enough to compare \(k^2ab\), which avoids any division.
After the single-angle routine is in place, the final answer is obtained by evaluating it at
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
and summing the resulting perimeters. The C++ and Java implementations split this outer loop across available processor threads and add partial sums at the end. The Python implementation performs the same computation serially.
For one angle, the inverse stage uses a fixed number of bisection iterations, and the local rational search advances through only nearby Farey neighbors. In the literal implementations, even the outward walk is guarded by a fixed iteration cap, so the coded per-angle cost is bounded by a constant. For the full problem this means the runtime scales essentially linearly with \(N\).
Memory usage is \(O(1)\) per worker for the single-angle search. The threaded implementations need \(O(T)\) additional storage for \(T\) partial sums, which is negligible compared with the arithmetic work.
Gesucht ist
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
mit \(N=45000\) und \(L=10^{10}\). Für einen Zielwinkel \(\alpha\) in Grad betrachtet man ganzzahlige rechtwinklige Dreiecke aus pythagoreischen Tripeln. Jedes zulässige Dreieck bestimmt einen speziellen Winkel \(\theta\), und \(f(\alpha,L)\) ist der Umfang des skalierten Dreiecks, dessen Winkel \(\alpha\) am nächsten kommt, wobei nach dem Skalieren die Hypotenuse höchstens \(L\) sein darf.
Die Implementierungen durchsuchen nicht alle Tripel direkt. Stattdessen wird das Problem auf eine eindimensionale Funktion des Quotienten \(u=m/n\) zurückgeführt, diese Funktion numerisch invertiert und anschließend nur in einer kleinen Umgebung nach tatsächlich zulässigen rationalen Werten gesucht.
Der Kern der Lösung besteht aus zwei Schritten. Zuerst werden primitive pythagoreische Tripel mit der Euclidischen Parametrisierung beschrieben. Danach wird aus dem Winkelproblem eine Optimierung über eine einzige reelle Variable \(u\), die später wieder auf ganzzahlige Kandidaten zurückgeführt wird.
Jedes primitive pythagoreische Tripel hat die Form
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
mit
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
Jede nicht-primitive Lösung ist nur eine skalierte Kopie \((ka,kb,kc)\) eines solchen Tripels. Die Skalierung verändert den zugehörigen Winkel nicht. Ist also ein primitives Tripel gewählt, dann ist der größte erlaubte Faktor
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$
Damit lautet der Rückgabewert von \(f(\alpha,L)\)
$$k(a+b+c).$$
Mathematisch muss man also unter allen primitiven Paaren \((m,n)\) dasjenige finden, dessen Winkel am nächsten an \(\alpha\) liegt, und anschließend seinen Umfang maximal skalieren.
Der problemspezifische Winkel wird durch
$$\tan\theta=\frac{3ab}{2c^2}$$
beschrieben. Setzt man Euclids Formeln ein und schreibt
$$u=\frac{m}{n}>1,$$
so erhält man
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$
Das ist die entscheidende Vereinfachung: Der Winkel hängt nicht von der absoluten Größe von \(m\) und \(n\) ab, sondern nur von ihrem Verhältnis.
Für einen Zielwinkel \(\alpha\) berechnet der Code
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$
Mit der Tangensformel für Winkeldifferenzen ergibt sich als Fehlermaß
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$
Da \(\tan x\) im hier relevanten kleinen Winkelbereich streng wächst, ist die Minimierung des echten Winkelabstands gleichbedeutend mit der Minimierung von \(\Delta\).
Die Ableitung lautet
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$
Für \(u>1\) gibt es genau einen kritischen Punkt, nämlich
$$u_{\max}=1+\sqrt{2}.$$
Dort gilt
$$t(u_{\max})=\frac34,$$
also steigt \(t(u)\) auf \((1,1+\sqrt2)\) und fällt auf \((1+\sqrt2,\infty)\). Deshalb besitzt jeder Zielwert \(t_\alpha<3/4\) zwei Urbilder:
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
Um diese beiden stetigen Lösungen herum organisiert sich die gesamte Suche.
Für das eigentliche Projekt-Euler-Summenproblem ist außerdem wichtig, dass \(\alpha=\sqrt[3]{i}\) mit \(1\le i\le45000\) gilt. Also
$$\alpha\le\sqrt[3]{45000}\approx 35{,}57^\circ\lt\arctan\!\left(\frac34\right)\approx 36{,}87^\circ.$$
Damit liegt jeder wirklich benötigte Winkel unter dem Maximum von \(t(u)\). Im produktiven Lauf treten also immer zwei Äste auf, auch wenn die Implementierungen den Spitzenfall vorsichtshalber ebenfalls behandeln.
Ein primitives Tripel entspricht einem vollständig gekürzten Bruch
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
mit der Größenbedingung
$$p^2+q^2\le L.$$
Liegt \(p/q\) nahe an einer stetigen Lösung \(u_\ast\), dann gilt näherungsweise \(p\approx u_\ast q\). Daraus folgt die natürliche Denominatorskala
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
Genau daraus wird in den Implementierungen die Farey-Ordnung
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor$$
gewonnen. Zuerst werden die Farey-Nachbarn dieser Ordnung bestimmt, die \(u_\ast\) einschließen. Danach wandert der Algorithmus in der Farey-Folge nach außen, bis auf beiden Seiten die nächstgelegenen zulässigen Brüche gefunden sind, die alle Bedingungen für primitive Tripel erfüllen.
Gerade dadurch wird das Verfahren schnell: Es durchsucht nicht die gesamte Kreisscheibe \(m^2+n^2\le L\), sondern lokalisiert zuerst das reelle Ziel und prüft dann nur wenige rationale Näherungen in dessen unmittelbarer Umgebung.
Für den kleinen Kontrollwert gilt
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0{,}57735.$$
Ein sehr guter zulässiger Bruch ist
$$\frac{m}{n}=\frac92.$$
Daraus entsteht
$$a=9^2-2^2=77,\qquad b=2\cdot9\cdot2=36,\qquad c=9^2+2^2=85.$$
Da \(c=85\le100\), ist der größte erlaubte Skalierungsfaktor \(k=\lfloor100/85\rfloor=1\). Der zugehörige Tangenswert ist
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0{,}57550,$$
also
$$\Delta=\frac{|0{,}57550-0{,}57735|}{1+0{,}57550\cdot0{,}57735}\approx 0{,}0013875.$$
Dieser Kandidat schlägt die nahen zulässigen Alternativen. Deshalb ist
$$f(30,100)=77+36+85=198.$$
Dieser Wert dient in der C++-Version als Prüfpunkt; dieselbe Mathematik steckt aber in allen drei Sprachfassungen.
Für jedes \(\alpha\) berechnen die C++-, Python- und Java-Implementierungen zunächst \(t_\alpha\). Falls \(t_\alpha<3/4\), wird \(t(u)=t_\alpha\) per Bisektion einmal auf dem steigenden Ast \((1,1+\sqrt2)\) und einmal auf dem fallenden Ast \((1+\sqrt2,\infty)\) gelöst. So entstehen zwei stetige Zielwerte \(u_-\) und \(u_+\). Würde ein Ziel den Spitzenwert \(3/4\) erreichen oder überschreiten, bliebe nur noch die Umgebung von \(u=1+\sqrt2\).
Für jeden Zielwert \(u_\ast\) bestimmt die Implementierung die passende Farey-Ordnung \(N_u\), konstruiert die einschließenden Nachbarn und läuft dann nach außen, bis die nächstgelegenen zulässigen Brüche auf beiden Seiten gefunden sind. Jeder gültige Bruch wird wieder in ein primitives Tripel \((a,b,c)\) übersetzt, mit \(k=\lfloor L/c\rfloor\) skaliert und über dasselbe Fehlermaß \(\Delta\) bewertet.
Sind zwei Kandidaten numerisch praktisch gleich gut, entscheidet die größere skalierte Fläche. Da diese Fläche gleich
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2}$$
ist, genügt der Vergleich von \(k^2ab\), ganz ohne zusätzliche Division.
Ist die Einzelfunktion implementiert, dann entsteht die Endantwort durch Auswertung bei
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
und anschließendes Aufsummieren aller Perimeter. C++ und Java verteilen diese äußere Schleife auf verfügbare Prozessorkerne und addieren danach die Teilsummen. Python führt dieselbe Rechnung seriell aus.
Für einen einzelnen Winkel besteht die inverse Phase aus einer festen Anzahl von Bisektionsschritten, und die rationale Suche bewegt sich nur über nahe Farey-Nachbarn. In den konkreten Implementierungen ist sogar die Auswärtsbewegung durch eine feste Iterationsgrenze abgesichert, sodass die kodierte Arbeit pro Winkel durch eine Konstante beschränkt ist. Für das Gesamtproblem skaliert die Laufzeit damit im Wesentlichen linear in \(N\).
Der Speicherbedarf beträgt \(O(1)\) pro Arbeiter für die Einzelsuche. Die parallelen Varianten benötigen zusätzlich \(O(T)\) Platz für \(T\) Teilsummen, was gegenüber der eigentlichen Rechenarbeit vernachlässigbar bleibt.
Hesaplanması istenen büyüklük
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
olup burada \(N=45000\) ve \(L=10^{10}\) alınır. Derece cinsinden verilen bir hedef açı \(\alpha\) için, Pisagor üçlülerinden gelen tam sayılı dik üçgenler incelenir. Her uygun üçgen belirli bir \(\theta\) açısı üretir ve \(f(\alpha,L)\), açısı \(\alpha\)'ya en yakın olan üçgenin, hipotenüsü ölçeklendikten sonra \(L\)'yi aşmayan en büyük kopyasının çevresidir.
Uygulamalar bütün üçlüleri kaba kuvvetle taramaz. Bunun yerine geometri, \(u=m/n\) oranının tek değişkenli bir fonksiyonuna indirgenir; sonra bu fonksiyon terslenir ve yalnızca gerçekten ilkel Pisagor üçlüsü üretebilen yakın rasyonel değerler test edilir.
Çözüm iki temel fikre dayanır. Birincisi, her ilkel tam sayılı dik üçgen Euclid parametreleriyle yazılabilir. İkincisi, problemdeki açı yalnızca \(m/n\) oranına bağlıdır; yani sürekli optimizasyon tek boyutludur ve daha sonra tekrar tam sayılı adaylara dönüştürülür.
Her ilkel Pisagor üçlüsü
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
şeklinde yazılır; burada
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
İlkel olmayan her çözüm, böyle bir üçlünün \((ka,kb,kc)\) biçimindeki ölçeklenmiş kopyasından ibarettir. Ölçekleme açıyı değiştirmez. Dolayısıyla ilkel üçlü belirlendiğinde izin verilen en büyük katsayı
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor$$
olur. Böylece \(f(\alpha,L)\) fonksiyonunun döndürdüğü değer
$$k(a+b+c)$$
yani ölçeklenmiş çevredir.
Matematiksel görev bu yüzden nettir: ilkel \((m,n)\) çiftleri arasında açısı hedefe en yakın olanı bul, sonra çevresini mümkün olan en büyük ölçekle çarp.
Problemde kullanılan özel açı şu bağıntıyla tanımlanır:
$$\tan\theta=\frac{3ab}{2c^2}.$$
Euclid formüllerini yerine koyup
$$u=\frac{m}{n}>1$$
yazarsak
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}$$
elde edilir. En kritik sadeleşme budur: açı artık \(m\) ve \(n\)'nin büyüklüklerine değil, yalnızca oranlarına bağlıdır.
Hedef açı \(\alpha\) için kod
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right)$$
değerini hesaplar. Tanjant farkı özdeşliği ile açısal hata
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|$$
şeklinde ölçülür. İlgili açı aralığında \(\tan x\) artan olduğundan, gerçek açı farkını küçültmek ile \(\Delta\)'yı küçültmek aynı şeydir.
Türev
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}$$
olur. \(u>1\) bölgesindeki tek kritik nokta
$$u_{\max}=1+\sqrt{2}$$
ve bu noktadaki en büyük değer
$$t(u_{\max})=\frac34$$
olduğundan, \(t(u)\) önce \((1,1+\sqrt2)\) aralığında artar, sonra \((1+\sqrt2,\infty)\) aralığında azalır. Bu nedenle \(t_\alpha<3/4\) olduğunda iki ters görüntü vardır:
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
Arama tam olarak bu iki sürekli hedef oranın etrafında örgütlenir.
Asıl Project Euler toplamında ayrıca şu gözlem çok değerlidir: \(\alpha=\sqrt[3]{i}\) ve \(1\le i\le45000\) olduğundan
$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$
Yani gerçek çalışmada kullanılan bütün açılar tepe değerin altındadır. Pratikte üretim koşusunda her zaman iki dal vardır; kodun tepe durumu için ayrı mantık içermesi ise sadece genel güvenlik içindir.
Bir ilkel üçlü, şu özellikleri sağlayan sade bir kesre karşılık gelir:
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
ve ayrıca
$$p^2+q^2\le L.$$
Eğer \(p/q\), sürekli çözümlerden biri olan \(u_\ast\)'a yakınsa yaklaşık olarak \(p\approx u_\ast q\) yazılabilir. Bu da payda için doğal ölçeği verir:
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
Uygulamalardaki Farey mertebesi tam olarak buradan gelir:
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$
Algoritma önce bu mertebedeki Farey komşularını kullanarak \(u_\ast\)'ı iki taraftan sarar. Daha sonra Farey dizisinde dışarı doğru yürüyerek, ilkel üçlü koşullarını gerçekten sağlayan en yakın alt ve üst kesirleri bulur. Sonunda yalnızca bu yerel geçerli adaylar açıkça değerlendirilir.
Hız buradan gelir: yöntem, \(m^2+n^2\le L\) diskindeki bütün noktalara bakmaz; önce sürekli hedefi bulur, sonra yalnızca o hedefin çevresindeki küçük bir rasyonel aday kümesini kontrol eder.
Küçük doğrulama örneğinde
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735$$
olur. Çok iyi bir geçerli aday
$$\frac{m}{n}=\frac92$$
kesridir. Bu oran
$$a=9^2-2^2=77,\qquad b=2\cdot9\cdot2=36,\qquad c=9^2+2^2=85$$
üçlüsünü üretir. \(c=85\le100\) olduğu için en büyük ölçek \(k=\lfloor100/85\rfloor=1\) olur. İlgili tanjant değeri
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550$$
ve hata ölçütü
$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875$$
çıkar. Bu aday yakındaki diğer geçerli seçeneklerden daha iyidir. Dolayısıyla
$$f(30,100)=77+36+85=198$$
elde edilir. C++ sürümü bu değeri doğrulama noktası olarak kullanır; aynı matematik Python ve Java sürümlerinde de vardır.
C++, Python ve Java uygulamaları her \(\alpha\) için önce \(t_\alpha\)'yı hesaplar. Eğer \(t_\alpha<3/4\) ise \(t(u)=t_\alpha\) denklemi, artan dal \((1,1+\sqrt2)\) üzerinde bir kez ve azalan dal \((1+\sqrt2,\infty)\) üzerinde bir kez ikili aramayla çözülür. Böylece iki sürekli hedef \(u_-\) ve \(u_+\) elde edilir. Hedef değer tepeye ulaşsa veya onu geçseydi arama \(u=1+\sqrt2\) çevresine sıkışacaktı.
Her sürekli hedef \(u_\ast\) için uygulama uygun Farey mertebesi \(N_u\)'yu hesaplar, hedefi saran komşuları kurar ve sonra dışarı doğru ilerleyerek iki taraftaki en yakın geçerli kesirleri bulur. Her geçerli kesir yeniden bir ilkel üçlü \((a,b,c)\) olarak kurulur, \(k=\lfloor L/c\rfloor\) ile ölçeklenir ve aynı \(\Delta\) ölçütüyle puanlanır.
İki aday açısal olarak neredeyse eşitse, daha büyük ölçeklenmiş alan seçilir. Çünkü ölçeklenmiş alan
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2}$$
olduğundan \(k^2ab\) karşılaştırması yeterlidir ve ek bölme gerektirmez.
Tek açı için rutin hazır olduktan sonra nihai cevap
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
değerlerinde hesaplanan çevrelerin toplanmasıyla oluşur. C++ ve Java sürümleri bu dış döngüyü mevcut işlemci iş parçacıklarına bölerek kısmi toplamlar üretir. Python sürümü aynı hesabı tek iş parçacığında yapar.
Tek bir açı için ters arama sabit sayıda ikili arama adımı kullanır ve rasyonel arama yalnızca yakın Farey komşuları arasında dolaşır. Üstelik somut uygulamalarda bu dışa yürüyüş de sabit bir iterasyon üst sınırı ile korunur. Bu yüzden kod düzeyinde tek açı başına maliyet sabitle sınırlıdır ve toplam çalışma süresi pratikte \(N\) ile yaklaşık doğrusal ölçeklenir.
Bellek kullanımı tek-açı araması için işçi başına \(O(1)\)'dir. Paralel sürümlerde yalnızca \(T\) adet kısmi toplam için ek \(O(T)\) alan gerekir; bu da aritmetik yükün yanında çok küçüktür.
La cantidad que se debe calcular es
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
con \(N=45000\) y \(L=10^{10}\). Para un ángulo objetivo \(\alpha\) medido en grados, se consideran triángulos rectángulos enteros procedentes de ternas pitagóricas. Cada triángulo admisible determina un ángulo especial \(\theta\), y \(f(\alpha,L)\) es el perímetro del triángulo escalado cuya \(\theta\) queda más cerca de \(\alpha\), siempre que la hipotenusa tras el escalado no supere \(L\).
La implementación no recorre todas las ternas. Primero convierte la geometría en una función de una sola variable, \(u=m/n\); luego invierte numéricamente esa función y, por último, busca sólo fracciones cercanas que de verdad puedan aparecer como parámetros de una terna pitagórica primitiva.
Todo el método se apoya en dos observaciones. La primera es la parametrización de Euclides para ternas pitagóricas primitivas. La segunda es que el ángulo relevante depende únicamente del cociente \(m/n\), así que la parte continua del problema es unidimensional antes de volver a los enteros.
Toda terna pitagórica primitiva puede escribirse como
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
con
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
Cualquier terna no primitiva es sólo una copia escalada \((ka,kb,kc)\) de una terna primitiva. El escalado no modifica el ángulo asociado, de modo que, una vez elegida la terna primitiva, el mayor factor permitido es
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$
Por tanto, el valor que devuelve \(f(\alpha,L)\) es
$$k(a+b+c).$$
La tarea matemática consiste entonces en encontrar, entre todos los pares primitivos \((m,n)\), el que produce el ángulo más cercano a \(\alpha\), y después tomar el perímetro de su mayor copia admisible.
El ángulo específico del problema viene dado por
$$\tan\theta=\frac{3ab}{2c^2}.$$
Al sustituir las fórmulas de Euclides y escribir
$$u=\frac{m}{n}>1,$$
se obtiene
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$
Esta reducción es decisiva: el ángulo ya no depende del tamaño absoluto de \(m\) y \(n\), sino sólo de su razón.
Para un ángulo objetivo \(\alpha\), el código calcula
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$
Usando la identidad del tangente de la diferencia, la discrepancia angular puede medirse con
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$
Como \(\tan x\) es creciente en el pequeño rango angular que aparece aquí, minimizar la diferencia angular real equivale a minimizar \(\Delta\).
La derivada es
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$
Para \(u>1\) existe un único punto crítico:
$$u_{\max}=1+\sqrt{2}.$$
En ese punto,
$$t(u_{\max})=\frac34,$$
de modo que \(t(u)\) crece en \((1,1+\sqrt2)\) y decrece en \((1+\sqrt2,\infty)\). Por eso, todo valor objetivo con \(t_\alpha<3/4\) tiene dos preimágenes:
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
Toda la búsqueda se organiza alrededor de esas dos soluciones continuas.
Hay además una observación concreta para este problema: en la suma final, \(\alpha=\sqrt[3]{i}\) con \(1\le i\le45000\), así que
$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$
Eso significa que todos los ángulos realmente usados quedan por debajo del máximo de \(t(u)\). En la ejecución principal siempre aparecen dos ramas, aunque las implementaciones también conservan el caso del pico por completitud.
Una terna primitiva corresponde a una fracción irreducible
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
sujeta además a
$$p^2+q^2\le L.$$
Si \(p/q\) está cerca de una raíz continua \(u_\ast\), entonces \(p\approx u_\ast q\), y la restricción circular sugiere la escala natural
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
De ahí sale el orden de Farey utilizado por el código:
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$
Primero se construyen los vecinos de Farey de orden \(N_u\) que encierran a \(u_\ast\). Después la búsqueda avanza hacia afuera en la sucesión de Farey hasta encontrar las fracciones válidas más cercanas por debajo y por encima, imponiendo todas las condiciones de primitividad.
Ésa es la razón de la eficiencia: el algoritmo no recorre todo el disco \(m^2+n^2\le L\), sino que localiza antes el objetivo continuo y luego inspecciona sólo una pequeña nube de aproximaciones racionales cercanas.
En el caso pequeño de validación,
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735.$$
Una fracción admisible muy buena es
$$\frac{m}{n}=\frac92.$$
Genera
$$a=9^2-2^2=77,\qquad b=2\cdot9\cdot2=36,\qquad c=9^2+2^2=85.$$
Como \(c=85\le100\), el mayor escalado posible es \(k=\lfloor100/85\rfloor=1\). El valor tangente asociado es
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550,$$
y por tanto
$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875.$$
Este candidato mejora a los competidores admisibles cercanos, así que
$$f(30,100)=77+36+85=198.$$
Ese valor se usa como punto de control en la versión C++, y la misma matemática aparece en las tres implementaciones.
Para cada \(\alpha\), las implementaciones en C++, Python y Java calculan primero \(t_\alpha\). Si \(t_\alpha<3/4\), resuelven \(t(u)=t_\alpha\) por bisección una vez en la rama creciente \((1,1+\sqrt2)\) y una vez en la rama decreciente \((1+\sqrt2,\infty)\). Así obtienen dos objetivos continuos, \(u_-\) y \(u_+\). Si el valor objetivo alcanzara o superara el máximo \(3/4\), la búsqueda se concentraría en el entorno de \(u=1+\sqrt2\).
Para cada objetivo continuo \(u_\ast\), la implementación calcula el orden de Farey \(N_u\), construye el encierro racional inicial y avanza hacia afuera hasta hallar las fracciones admisibles más próximas a ambos lados. Cada fracción válida se reconstruye como terna primitiva \((a,b,c)\), se escala con \(k=\lfloor L/c\rfloor\) y se puntúa con la misma cantidad \(\Delta\).
Si dos candidatos quedan prácticamente empatados en error angular, se elige el de mayor área escalada. Como esa área es
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2},$$
basta comparar \(k^2ab\), sin introducir divisiones adicionales.
Una vez resuelto el caso de un solo ángulo, la respuesta final se obtiene evaluando la rutina en
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
y sumando todos los perímetros resultantes. Las versiones de C++ y Java reparten este bucle exterior entre los hilos disponibles y luego agregan sumas parciales. La versión de Python realiza el mismo cálculo de forma secuencial.
Para un ángulo individual, la fase inversa usa un número fijo de iteraciones de bisección, y la búsqueda racional sólo recorre vecinos de Farey cercanos. En las implementaciones concretas, incluso ese recorrido hacia afuera está protegido por un tope fijo de iteraciones, así que el coste codificado por ángulo queda acotado por una constante. En consecuencia, el tiempo total crece esencialmente de manera lineal con \(N\).
El uso de memoria es \(O(1)\) por trabajador en la búsqueda de un ángulo. Las versiones paralelas sólo añaden \(O(T)\) espacio para \(T\) sumas parciales, una cantidad muy pequeña frente al trabajo aritmético.
题目要求计算
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
其中 \(N=45000\)，\(L=10^{10}\)。对每一个以“度”为单位的目标角 \(\alpha\)，我们考虑由勾股三元组得到的整数直角三角形。每个可行三角形都会对应一个特殊角 \(\theta\)，而 \(f(\alpha,L)\) 定义为：在允许缩放、且缩放后斜边不超过 \(L\) 的前提下，选择使 \(\theta\) 最接近 \(\alpha\) 的那个三角形，并返回其最大可行缩放副本的周长。
实现并没有暴力枚举所有三元组。它先把几何问题压缩成关于比值 \(u=m/n\) 的单变量函数，再对这个连续函数做反求，最后只在目标值附近搜索真正可能对应原始勾股三元组的有理数候选。
整个解法建立在两个核心事实上。第一，所有原始勾股三元组都可以用 Euclid 参数化表示。第二，题目中的角只依赖于 \(m/n\) 这个比值，因此连续优化阶段本质上是一维问题，之后再回到整数候选上。
每个原始勾股三元组都可写成
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
其中
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
任何非原始解都只是这样一个三元组的整数倍 \((ka,kb,kc)\)。缩放不会改变关联角，因此一旦原始三元组固定，允许的最大缩放倍数就是
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$
于是 \(f(\alpha,L)\) 返回的量就是
$$k(a+b+c).$$
换句话说，数学任务可以表述为：在所有原始参数 \((m,n)\) 中，找到使题目定义的角最接近 \(\alpha\) 的那一个，然后取它在限制 \(L\) 下的最大缩放周长。
题目中的特殊角通过下面的公式给出：
$$\tan\theta=\frac{3ab}{2c^2}.$$
把 Euclid 参数化代入，并记
$$u=\frac{m}{n}>1,$$
就得到
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$
这是整个问题最关键的降维步骤：角度不再依赖 \(m\) 和 \(n\) 的绝对大小，而只依赖于它们的比值。
对目标角 \(\alpha\)，代码先计算
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$
再利用正切差角公式，把角度误差改写成
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$
在本题涉及的角度范围内，\(\tan x\) 是严格递增的，因此最小化真实角差与最小化 \(\Delta\) 完全等价。
对 \(t(u)\) 求导可得
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$
在 \(u>1\) 的区域里，唯一的临界点是
$$u_{\max}=1+\sqrt{2}.$$
代回去可得峰值
$$t(u_{\max})=\frac34.$$
因此 \(t(u)\) 在 \((1,1+\sqrt2)\) 上单调递增，在 \((1+\sqrt2,\infty)\) 上单调递减。于是只要 \(t_\alpha<3/4\)，就会有两个连续解：
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
之后的有理搜索正是围绕这两个连续目标展开的。
这里还有一个很有价值的题目特定结论：在最终求和中，\(\alpha=\sqrt[3]{i}\)，且 \(1\le i\le45000\)，所以
$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$
这意味着实际大规模计算时，所有目标角都严格落在峰值以下。也就是说，真正的生产计算始终是“两分支”情形；代码保留峰值附近的处理，只是为了更完整和更稳健。
一个原始勾股三元组对应一个既约分数
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
并且还必须满足大小约束
$$p^2+q^2\le L.$$
如果 \(p/q\) 靠近连续目标 \(u_\ast\)，那么近似有 \(p\approx u_\ast q\)，于是圆盘约束自然给出分母尺度
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
实现中据此定义 Farey 阶数
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$
接下来程序先构造 Farey 序列中阶数为 \(N_u\) 且恰好把 \(u_\ast\) 夹住的一对邻分数，然后向左右两侧逐步扩展，直到找到满足所有原始条件的最近下侧分数和最近上侧分数。最后只需对这些局部候选做显式评估。
这正是算法快的原因：它并不扫描整个区域 \(m^2+n^2\le L\)，而是先解决连续逆问题，再只检查目标附近极少量的有理逼近。
在小型校验例子中，
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735.$$
一个非常好的可行候选是
$$\frac{m}{n}=\frac92.$$
它生成的三元组为
$$a=9^2-2^2=77,\qquad b=2\cdot9\cdot2=36,\qquad c=9^2+2^2=85.$$
由于 \(c=85\le100\)，最大缩放倍数是 \(k=\lfloor100/85\rfloor=1\)。对应的正切值为
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550,$$
因此误差代理量
$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875.$$
它优于附近其他可行候选，所以
$$f(30,100)=77+36+85=198.$$
这个值在 C++ 实现中被用作校验点，而三种语言的算法本质完全相同。
对于每个 \(\alpha\)，C++、Python 和 Java 实现都会先计算 \(t_\alpha\)。当 \(t_\alpha<3/4\) 时，它们分别在递增区间 \((1,1+\sqrt2)\) 与递减区间 \((1+\sqrt2,\infty)\) 上用二分法求解 \(t(u)=t_\alpha\)，从而得到两个连续目标 \(u_-\) 和 \(u_+\)。如果目标值达到或超过峰值 \(3/4\)，搜索就退化为围绕 \(u=1+\sqrt2\) 的局部寻找。
对于每个连续目标 \(u_\ast\)，实现先计算相应的 Farey 阶数 \(N_u\)，构造初始夹逼区间，再沿 Farey 邻分数向外行走，直到找到最近的可行下侧分数和上侧分数。每个有效分数都被重建成原始三元组 \((a,b,c)\)，再按 \(k=\lfloor L/c\rfloor\) 放缩，并用同一个 \(\Delta\) 指标评分。
如果两个候选在角度误差上几乎无法区分，程序会选择缩放后面积更大的那个。因为缩放后的面积是
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2},$$
所以只比较 \(k^2ab\) 就够了，不需要真的做除法。
当单个角度的子程序完成后，最终答案就是在
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
这些角度上逐个求值并把周长相加。C++ 和 Java 实现会把最外层循环分配到多个处理器线程上，再汇总部分和；Python 实现则按相同数学逻辑顺序执行。
对单个角而言，反解阶段只需要固定次数的二分迭代，而有理搜索只在附近的 Farey 邻分数之间移动。在实际实现里，这个向外扩展过程还带有固定的迭代上限，因此从代码层面看，单角成本被一个常数所界定。于是总运行时间在整体上近似与 \(N\) 成线性关系。
内存方面，单角搜索对每个工作线程只需 \(O(1)\) 额外空间。并行版本只再增加 \(O(T)\) 的部分和存储，其中 \(T\) 是线程数，这相对于主要运算量来说非常小。
Требуется вычислить
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
где \(N=45000\) и \(L=10^{10}\). Для каждого целевого угла \(\alpha\), измеряемого в градусах, рассматриваются целочисленные прямоугольные треугольники, получающиеся из пифагоровых троек. Каждый допустимый треугольник задаёт специальный угол \(\theta\), а \(f(\alpha,L)\) равно периметру наибольшей допустимой масштабированной копии того треугольника, чей угол \(\theta\) ближе всего к \(\alpha\).
Реализация не перебирает все тройки напрямую. Она сводит геометрию к функции одной переменной \(u=m/n\), численно обращает эту функцию и затем ищет только близкие рациональные значения, которые действительно могут соответствовать примитивным пифагоровым тройкам.
В основе решения лежат два ключевых наблюдения. Во-первых, любая примитивная пифагорова тройка имеет параметризацию Евклида. Во-вторых, нужный угол зависит только от отношения \(m/n\), поэтому непрерывная часть задачи одномерна, а затем уже переводится обратно в целочисленные кандидаты.
Каждая примитивная пифагорова тройка записывается в виде
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
где
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
Любая непримитивная тройка является просто масштабированной копией \((ka,kb,kc)\) некоторой примитивной. Масштабирование не меняет соответствующий угол, поэтому после выбора примитивной тройки максимальный допустимый коэффициент равен
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$
Следовательно, функция \(f(\alpha,L)\) возвращает величину
$$k(a+b+c).$$
Итак, математическая задача состоит в том, чтобы среди всех примитивных пар \((m,n)\) найти пару с углом, ближайшим к \(\alpha\), а затем взять периметр её наибольшей допустимой масштабированной копии.
Специальный угол задачи задаётся формулой
$$\tan\theta=\frac{3ab}{2c^2}.$$
Подставляя параметризацию Евклида и вводя
$$u=\frac{m}{n}>1,$$
получаем
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$
Это ключевое упрощение: угол больше не зависит от абсолютных размеров \(m\) и \(n\), а только от их отношения.
Для целевого угла \(\alpha\) код вычисляет
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$
По формуле тангенса разности углов угловую ошибку удобно измерять выражением
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$
Поскольку на рассматриваемом диапазоне \(\tan x\) строго возрастает, минимизация реальной угловой разницы эквивалентна минимизации \(\Delta\).
Производная равна
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$
На области \(u>1\) имеется единственная критическая точка
$$u_{\max}=1+\sqrt{2}.$$
В ней достигается максимум
$$t(u_{\max})=\frac34.$$
Значит, функция \(t(u)\) возрастает на \((1,1+\sqrt2)\) и убывает на \((1+\sqrt2,\infty)\). Поэтому для любого \(t_\alpha<3/4\) существуют две непрерывные прообразы:
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
Именно вокруг этих двух вещественных целей и строится дискретный поиск.
Для самой суммы Project Euler есть ещё одно важное наблюдение: \(\alpha=\sqrt[3]{i}\) при \(1\le i\le45000\), следовательно
$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$
То есть все реальные углы основной задачи лежат ниже вершины графика \(t(u)\). В рабочем запуске всегда присутствуют две ветви, хотя код на всякий случай умеет обрабатывать и граничный случай около максимума.
Примитивной тройке соответствует несократимая дробь
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
которая дополнительно должна удовлетворять ограничению
$$p^2+q^2\le L.$$
Если \(p/q\) близко к непрерывному корню \(u_\ast\), то \(p\approx u_\ast q\), и отсюда следует естественный масштаб для знаменателя:
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
Именно это приводит в реализации к выбору порядка Farey
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$
Сначала строятся соседи Farey этого порядка, между которыми лежит \(u_\ast\). Затем поиск идёт наружу по последовательности Farey, пока не будут найдены ближайшие допустимые дроби снизу и сверху, удовлетворяющие всем условиям примитивности.
В этом и состоит эффективность метода: он не просматривает весь диск \(m^2+n^2\le L\), а сначала локализует вещественную цель, после чего исследует только маленькое облако рациональных приближений вокруг неё.
Для малого проверочного случая
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735.$$
Очень хорошим допустимым кандидатом является дробь
$$\frac{m}{n}=\frac92.$$
Она даёт тройку
$$a=9^2-2^2=77,\qquad b=2\cdot9\cdot2=36,\qquad c=9^2+2^2=85.$$
Поскольку \(c=85\le100\), максимальный масштаб равен \(k=\lfloor100/85\rfloor=1\). Соответствующее значение тангенса равно
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550,$$
а значит
$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875.$$
Этот кандидат оказывается лучше близких допустимых альтернатив, поэтому
$$f(30,100)=77+36+85=198.$$
Это значение используется как контрольная точка в C++-версии, а та же математическая схема присутствует и в Python, и в Java.
Для каждого \(\alpha\) реализации на C++, Python и Java сначала вычисляют \(t_\alpha\). Если \(t_\alpha<3/4\), уравнение \(t(u)=t_\alpha\) решается методом бисекции дважды: на возрастающей ветви \((1,1+\sqrt2)\) и на убывающей ветви \((1+\sqrt2,\infty)\). Так получаются два непрерывных ориентира \(u_-\) и \(u_+\). Если бы целевое значение достигало максимума \(3/4\) или превосходило его, поиск сводился бы к окрестности \(u=1+\sqrt2\).
Для каждого непрерывного ориентира \(u_\ast\) вычисляется соответствующий порядок Farey \(N_u\), строится начальное рациональное зажатие, а затем выполняется продвижение наружу до ближайших допустимых дробей по обе стороны. Каждая валидная дробь восстанавливается как примитивная тройка \((a,b,c)\), масштабируется с коэффициентом \(k=\lfloor L/c\rfloor\) и оценивается тем же показателем \(\Delta\).
Если два кандидата практически совпадают по угловой ошибке, код выбирает тот, у которого больше масштабированная площадь. Поскольку она равна
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2},$$
достаточно сравнивать \(k^2ab\), не выполняя деления.
После реализации процедуры для одного угла итоговый ответ получается подстановкой
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
и суммированием всех полученных периметров. Реализации на C++ и Java делят внешний цикл между доступными потоками процессора и затем складывают частичные суммы. Версия на Python делает ту же работу последовательно.
Для одного угла обратная стадия использует фиксированное число шагов бисекции, а рациональный поиск проходит только по близким соседям Farey. Более того, в самих реализациях даже этот внешний проход защищён фиксированным верхним пределом числа итераций, так что закодированная стоимость на один угол ограничена константой. Поэтому общее время работы по существу растёт линейно по \(N\).
Память составляет \(O(1)\) на один рабочий поток в процедуре для одного угла. Параллельные варианты требуют лишь дополнительное \(O(T)\) место для \(T\) частичных сумм, что ничтожно по сравнению с объёмом вычислений.
المطلوب حساب الكمية
$$F(N,L)=\sum_{i=1}^{N} f\!\left(\sqrt[3]{i},L\right),$$
حيث \(N=45000\) و\(L=10^{10}\). لكل زاوية هدف \(\alpha\) مقاسة بالدرجات، ننظر إلى المثلثات القائمة ذات الأضلاع الصحيحة الناتجة من الثلاثيات الفيثاغورية. كل مثلث صالح يحدد زاوية خاصة \(\theta\)، وتكون \(f(\alpha,L)\) هي محيط أكبر نسخة مكبرة مسموح بها من المثلث الذي يجعل \(\theta\) أقرب ما يمكن إلى \(\alpha\)، مع شرط ألا تتجاوز الوتر بعد التكبير القيمة \(L\).
التنفيذ لا يفحص جميع الثلاثيات بالقوة الغاشمة. بدلاً من ذلك يحول المسألة الهندسية إلى دالة في متغير واحد هو النسبة \(u=m/n\)، ثم يعكس هذه الدالة عددياً، وبعد ذلك يبحث فقط في الكسور القريبة التي يمكن فعلاً أن تمثل ثلاثيات فيثاغورية أولية.
الحل مبني على حقيقتين أساسيتين. الأولى هي أن كل ثلاثية فيثاغورية أولية تكتب بتمثيل إقليدس. والثانية هي أن الزاوية المطلوبة تعتمد فقط على النسبة \(m/n\)، ولذلك يصبح الجزء المستمر من المسألة أحادي البعد قبل الرجوع إلى المرشحين الصحيحين.
كل ثلاثية فيثاغورية أولية يمكن كتابتها على الصورة
$$a=m^2-n^2,\qquad b=2mn,\qquad c=m^2+n^2,$$
مع الشروط
$$m>n>0,\qquad \gcd(m,n)=1,\qquad m+n\equiv 1 \pmod 2.$$
وأي ثلاثية غير أولية ليست إلا نسخة مكبرة \((ka,kb,kc)\) من ثلاثية أولية. التكبير لا يغير الزاوية المرتبطة، ولذلك بعد اختيار الثلاثية الأولية يكون أكبر عامل تكبير مسموح هو
$$k=\left\lfloor\frac{L}{c}\right\rfloor=\left\lfloor\frac{L}{m^2+n^2}\right\rfloor.$$
وعليه فإن القيمة التي تعيدها \(f(\alpha,L)\) هي
$$k(a+b+c).$$
إذن المهمة الرياضية هي اختيار الزوج الأولي \((m,n)\) الذي يعطي زاوية أقرب ما تكون إلى \(\alpha\)، ثم أخذ محيط أكبر نسخة مكبرة مسموح بها لذلك المثلث.
الزاوية الخاصة في هذه المسألة تُعطى بالعلاقة
$$\tan\theta=\frac{3ab}{2c^2}.$$
عند تعويض صيغة إقليدس وكتابة
$$u=\frac{m}{n}>1,$$
نحصل على
$$\tan\theta=t(u)=\frac{3u(u^2-1)}{(u^2+1)^2}.$$
هذا هو الاختزال الحاسم في الحل: الزاوية لم تعد تعتمد على الحجم المطلق لـ \(m\) و\(n\)، بل على نسبتهما فقط.
وبالنسبة إلى زاوية الهدف \(\alpha\)، يحسب الكود
$$t_\alpha=\tan\!\left(\frac{\alpha\pi}{180}\right).$$
وباستخدام هوية ظل فرق الزاويتين يمكن قياس الخطأ بالكمية
$$\Delta(u,\alpha)=\frac{|t(u)-t_\alpha|}{1+t(u)t_\alpha}=\tan|\theta-\alpha|.$$
ولأن \(\tan x\) دالة متزايدة على المجال الزاوي الصغير الذي يهمنا هنا، فإن تصغير فرق الزاوية الحقيقي يكافئ تماماً تصغير \(\Delta\).
المشتقة هي
$$t'(u)=-\frac{3(u^4-6u^2+1)}{(u^2+1)^3}.$$
وعلى المجال \(u>1\) توجد نقطة حرجة وحيدة هي
$$u_{\max}=1+\sqrt{2}.$$
وفيها تبلغ الدالة القيمة العظمى
$$t(u_{\max})=\frac34.$$
لذلك تكون \(t(u)\) متزايدة على \((1,1+\sqrt2)\) ومتناقصة على \((1+\sqrt2,\infty)\). ومن ثم إذا كان \(t_\alpha<3/4\) فهناك صورتان عكسيتان مستمرتان:
$$u_-\in(1,1+\sqrt2),\qquad u_+\in(1+\sqrt2,\infty).$$
كل البحث العددي في التنفيذ يتمحور حول هاتين القيمتين المستمرتين.
وهنا توجد ملاحظة خاصة بالمسألة نفسها: في المجموع النهائي لدينا \(\alpha=\sqrt[3]{i}\) حيث \(1\le i\le45000\)، ولذلك
$$\alpha\le\sqrt[3]{45000}\approx 35.57^\circ\lt\arctan\!\left(\frac34\right)\approx 36.87^\circ.$$
أي أن كل الزوايا المستخدمة فعلاً في الحساب الرئيسي تقع أسفل قمة الدالة. هذا يعني أن التشغيل الحقيقي يكون دائماً في حالة الفرعين، رغم أن التنفيذ يحتفظ أيضاً بحالة القمة من أجل الاكتمال والصلابة العددية.
الثلاثية الأولية تقابل كسراً مختزلاً من الصورة
$$u=\frac{p}{q},\qquad p>q>0,\qquad \gcd(p,q)=1,\qquad p+q\equiv 1 \pmod 2,$$
ويجب أيضاً أن يحقق الشرط
$$p^2+q^2\le L.$$
إذا كان \(p/q\) قريباً من الجذر المستمر \(u_\ast\)، فلدينا تقريباً \(p\approx u_\ast q\)، ومن ثم يوحي القيد الدائري بالمقياس الطبيعي للمقام:
$$q\lesssim \frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}.$$
ولهذا يختار التنفيذ رتبة Farey التالية:
$$N_u=\left\lfloor\frac{\sqrt{L}}{\sqrt{1+u_\ast^2}}\right\rfloor.$$
بعد ذلك يبني زوج الجوار في متتالية Farey من هذه الرتبة الذي يحصر \(u_\ast\)، ثم يتحرك إلى الخارج حتى يجد أقرب كسرين صالحين من الجهتين، مع فرض جميع شروط الأولية الخاصة بثلاثيات إقليدس.
ومن هنا تأتي السرعة: الخوارزمية لا تمسح كامل القرص \(m^2+n^2\le L\)، بل تحدد أولاً الهدف المستمر ثم تختبر فقط مجموعة محلية صغيرة جداً من التقريبات الكسرية القريبة منه.
في مثال التحقق الصغير نحصل على
$$t_\alpha=\tan 30^\circ=\frac{1}{\sqrt3}\approx 0.57735.$$
أحد أفضل المرشحين الصالحين هو
$$\frac{m}{n}=\frac92.$$
وهذا يعطي
$$a=9^2-2^2=77,\qquad b=2\cdot9\cdot2=36,\qquad c=9^2+2^2=85.$$
وبما أن \(c=85\le100\)، فإن أكبر عامل تكبير هو \(k=\lfloor100/85\rfloor=1\). وقيمة الظل المقابلة هي
$$t=\frac{3ab}{2c^2}=\frac{3\cdot77\cdot36}{2\cdot85^2}\approx 0.57550,$$
ومن ثم يكون مقياس الخطأ
$$\Delta=\frac{|0.57550-0.57735|}{1+0.57550\cdot0.57735}\approx 0.0013875.$$
هذا المرشح يتفوق على البدائل الصالحة القريبة، ولذلك نحصل على
$$f(30,100)=77+36+85=198.$$
وتستخدم نسخة C++ هذه القيمة كنقطة تحقق، بينما تعتمد اللغات الثلاث على الرياضيات نفسها.
لكل \(\alpha\)، تحسب تطبيقات C++ وPython وJava أولاً القيمة \(t_\alpha\). وإذا كان \(t_\alpha<3/4\)، فإنها تحل المعادلة \(t(u)=t_\alpha\) بالبحث الثنائي مرة على الفرع المتزايد \((1,1+\sqrt2)\) ومرة على الفرع المتناقص \((1+\sqrt2,\infty)\)، وبذلك تحصل على هدفين مستمرين \(u_-\) و\(u_+\). أما إذا بلغت القيمة الهدف القمة \(3/4\) أو تجاوزتها، فإن البحث يتركز قرب \(u=1+\sqrt2\).
لكل هدف مستمر \(u_\ast\)، يحسب التنفيذ رتبة Farey المناسبة \(N_u\)، ويبني الحصر الكسري الابتدائي، ثم يتحرك إلى الخارج حتى يصل إلى أقرب الكسور الصالحة من الجهتين. كل كسر صالح يعاد بناؤه كثلاثية أولية \((a,b,c)\)، ثم يكبر بالعامل \(k=\lfloor L/c\rfloor\)، ثم يقيم باستخدام المقدار نفسه \(\Delta\).
إذا تعادل مرشحان تقريباً في الخطأ الزاوي، يختار التنفيذ صاحب المساحة المكبرة الأكبر. وبما أن هذه المساحة تساوي
$$\frac{(ka)(kb)}{2}=\frac{k^2ab}{2},$$
فإن مقارنة \(k^2ab\) تكفي بالكامل ولا تحتاج إلى أي قسمة إضافية.
بعد اكتمال روتين الزاوية الواحدة، تُحسب الإجابة النهائية عند القيم
$$\alpha_i=\sqrt[3]{i}\qquad (1\le i\le N)$$
ثم تُجمع المحيطات الناتجة. نسختا C++ وJava توزعان هذه الحلقة الخارجية على الأنوية أو الخيوط المتاحة ثم تجمعان المجاميع الجزئية. أما نسخة Python فتؤدي الحساب نفسه بشكل تسلسلي.
لكل زاوية منفردة، تستخدم مرحلة العكس عدداً ثابتاً من خطوات البحث الثنائي، كما أن البحث الكسري يتحرك فقط بين جيران Farey القريبين. وفي التنفيذ الفعلي يوجد أيضاً حد ثابت لعدد خطوات التوسع إلى الخارج، لذلك تكون الكلفة البرمجية لكل زاوية محصورة بثابت. ونتيجة لذلك ينمو زمن التنفيذ الكلي تقريباً خطياً مع \(N\).
أما الذاكرة الإضافية فهي \(O(1)\) لكل عامل في بحث الزاوية الواحدة. النسخ المتوازية تحتاج فقط إلى \(O(T)\) لتخزين \(T\) من المجاميع الجزئية، وهو مقدار ضئيل جداً مقارنة بالعمل الحسابي نفسه.