Let \(T_n=\frac{n(n+1)}{2}\) be the \(n\)-th triangular number. The problem asks for the sum \(S(N)\) of all positive integers \(c\le N\) for which there exist two distinct triangular numbers \(T_a\) and \(T_b\) such that
$$T_aT_b=c^2.$$
For the actual input \(N=10^{35}\), direct enumeration of triangular pairs is impossible. The successful approach groups triangular numbers by squarefree part, turns each group into a Pell-type family, and then sums the corresponding square roots modulo \(136101521\).
The key is to understand when two triangular numbers have a square product. That happens exactly when they belong to the same squarefree class.
If \(T_aT_b\) is a square, then \(T_a\) and \(T_b\) have the same squarefree part. So we can write
$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$
where \(d\) is squarefree. Then automatically
$$c=d\,rs.$$
Thus, for each fixed squarefree \(d\), we need all triangular numbers of the form
$$T_n=d\,y^2.$$
Using \(T_n=\frac{n(n+1)}{2}\) and \(x=2n+1\), this becomes
$$x^2-8dy^2=1.$$
So every squarefree class of triangular numbers is controlled by a Pell equation.
Suppose a squarefree class has a fundamental representative \(T_k=d\). Then
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
The Pell family generated by this base produces all triangular numbers in the same squarefree class. The implementation expresses the square multipliers with Chebyshev polynomials of the second kind:
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
For this fixed base \(k\), the family of triangular numbers is
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
These are exactly the triangular numbers whose squarefree part is \(T_k\). The first term is \(A_0(k)=T_k\) because \(U_0=1\).
Take two distinct members of the same family, say \(A_i(k)\) and \(A_j(k)\) with \(i<j\). Their product is
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$
Therefore each admissible contribution to \(S(N)\) has the form
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$
So the entire problem becomes: enumerate all family pairs \((i,j)\), sum \(C_{i,j}(k)\) over all fundamental bases \(k\), and keep only the values not exceeding \(N\).
For fixed \(i\) and \(j\), the function \(C_{i,j}(k)\) is increasing in \(k\), because both \(T_k\) and \(U_m(2k+1)\) increase for \(k\ge 1\). The smallest base is \(k=1\), where \(2k+1=3\), so a necessary condition for the pair \((i,j)\) to contribute at all is
$$U_i(3)\,U_j(3)\le N.$$
This is why the implementations first build the short sequence
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
and keep only the pair indices whose minimum possible contribution is not already too large.
Once \((i,j)\) is fixed, there is a largest admissible base \(\kappa_{i,j}\) satisfying
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
Because \(U_m(2k+1)\) is a polynomial in \(k\) of degree \(m\), the whole expression \(C_{i,j}(k)\) is a polynomial in \(k\). Hence
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
can be reduced to a linear combination of power sums \(\sum_{k=1}^{\kappa_{i,j}} k^p\).
Not every base index \(k\) should start a new squarefree class. Some triangular numbers already occur inside an older family:
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
for a smaller base \(k_1\). Such \(k\) are non-fundamental. If we summed over all \(k\ge 1\) blindly, those classes would be counted more than once.
The duplicate base indices are generated by the same Pell mechanism. Starting from \(x_1=2k_1+1\), define
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$
Every \(k_r\) produced in this way belongs to the squarefree class already represented by \(k_1\), so its contribution must be subtracted from the all-\(k\) polynomial sum.
Take the smallest fundamental base \(k=1\). Then
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
The corresponding triangular family begins with
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
Now pair these family members:
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$
So \(6\), \(35\), and \(210\) are all square triangle products.
The same family also explains why some later bases are non-fundamental. The Pell recurrence from \(x_1=3\) gives
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
hence
$$k=0,1,8,49,\dots$$
and indeed
$$T_8=36,\qquad T_{49}=1225,$$
which were already present in the family of \(k=1\). Those later bases must therefore be excluded as new starting points.
If \(\mathcal{N}_{i,j}\) denotes the non-fundamental bases not exceeding the relevant bound for the pair \((i,j)\), then the target sum is
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
This is exactly the structure implemented by the programs.
The C++, Python, and Java implementations all follow the same pipeline. First they enumerate the admissible family pairs \((i,j)\) by growing the sequence \(U_m(3)\) until its products are too large to matter. For each surviving pair they find the largest admissible base by monotone binary search on \(C_{i,j}(k)\le N\).
Next they build the Chebyshev factors symbolically as polynomials in \(k\). Multiplying those polynomials by the triangular factor \(T_k=\frac{k(k+1)}{2}\) gives a polynomial representation of each contribution \(C_{i,j}(k)\). The full range sum for that pair is then evaluated through power sums, and those power sums are obtained modulo \(136101521\) by Lagrange-style interpolation of Faulhaber polynomials.
Finally the implementations generate every non-fundamental base only once, evaluate the corresponding pair contributions at those indices, and subtract them. The C++ implementation parallelizes this subtraction pass, while the Python and Java implementations perform the same arithmetic sequentially.
Let \(F\) be the number of family levels with \(U_F(3)\le N\). Because \(U_m(3)\) grows exponentially, \(F=O(\log N)\), so the number of family pairs is \(O(F^2)\). Let \(B\) be the largest admissible base among all pairs, let \(D\) be the maximum polynomial degree, and let \(L\) be the number of non-fundamental bases up to \(B\).
Enumerating pair indices costs \(O(F^2)\). Finding all pair bounds costs \(O(F^2\log B)\) evaluations, each using a short Chebyshev recurrence. Polynomial preprocessing and power-sum evaluation are low-degree operations with total cost polynomial in \(F\), and the correction phase costs \(O(F^2L)\) in the straightforward form used here. Memory is dominated by the stored pair metadata, polynomial coefficients, and the correction list, so it is \(O(F^2D+L)\).
The important practical fact is that the family depth grows only logarithmically, which keeps the entire computation manageable even at \(N=10^{35}\).
Sei \(T_n=\frac{n(n+1)}{2}\) die \(n\)-te Dreieckszahl. Gesucht ist die Summe \(S(N)\) aller positiven ganzen Zahlen \(c\le N\), für die es zwei verschiedene Dreieckszahlen \(T_a\) und \(T_b\) gibt mit
$$T_aT_b=c^2.$$
Für den eigentlichen Eingabewert \(N=10^{35}\) ist eine direkte Suche über Dreieckszahlpaare unmöglich. Deshalb gruppiert die Lösung die Dreieckszahlen nach ihrem quadratfreien Anteil, beschreibt jede Gruppe durch eine Pell-artige Familie und summiert dann die zugehörigen Quadratwurzeln modulo \(136101521\).
Entscheidend ist die Beobachtung, dass das Produkt zweier Dreieckszahlen genau dann ein Quadrat ist, wenn beide denselben quadratfreien Kern besitzen.
Ist \(T_aT_b\) ein Quadrat, so haben \(T_a\) und \(T_b\) denselben quadratfreien Anteil. Also können wir schreiben
$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$
wobei \(d\) quadratfrei ist. Dann gilt sofort
$$c=d\,rs.$$
Für festes \(d\) müssen wir also alle Dreieckszahlen der Form
$$T_n=d\,y^2$$
finden. Mit \(T_n=\frac{n(n+1)}{2}\) und \(x=2n+1\) wird daraus
$$x^2-8dy^2=1.$$
Jede quadratfreie Klasse von Dreieckszahlen wird also von einer Pell-Gleichung gesteuert.
Angenommen, eine quadratfreie Klasse besitzt den fundamentalen Repräsentanten \(T_k=d\). Dann ist
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
Die gesamte Pell-Familie zu dieser Basis liefert alle Dreieckszahlen derselben quadratfreien Klasse. Die Implementierung beschreibt die Quadratfaktoren mit Chebyshev-Polynomen zweiter Art:
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
Für dieses feste \(k\) lautet die Dreieckszahlfamilie
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
Genau das sind die Dreieckszahlen, deren quadratfreier Anteil \(T_k\) ist. Der erste Term ist \(A_0(k)=T_k\), weil \(U_0=1\).
Nimmt man zwei verschiedene Glieder derselben Familie, also \(A_i(k)\) und \(A_j(k)\) mit \(i<j\), so ist ihr Produkt
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$
Damit hat jeder Beitrag zu \(S(N)\) die Form
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$
Das Problem reduziert sich also darauf, alle Familienpaare \((i,j)\) zu betrachten, \(C_{i,j}(k)\) über alle fundamentalen Basen \(k\) zu summieren und nur Werte bis \(N\) zu behalten.
Für feste \(i\) und \(j\) wächst \(C_{i,j}(k)\) mit \(k\), weil sowohl \(T_k\) als auch \(U_m(2k+1)\) für \(k\ge 1\) wachsen. Die kleinste Basis ist \(k=1\), also \(2k+1=3\). Deshalb kann ein Paar \((i,j)\) nur dann überhaupt beitragen, wenn
$$U_i(3)\,U_j(3)\le N.$$
Genau deswegen bauen die Implementierungen zuerst die kurze Folge
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
auf und behalten nur die Indexpaare, deren minimal möglicher Beitrag nicht schon zu groß ist.
Ist \((i,j)\) fixiert, so gibt es eine größte zulässige Basis \(\kappa_{i,j}\) mit
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
Weil \(U_m(2k+1)\) ein Polynom in \(k\) vom Grad \(m\) ist, ist auch \(C_{i,j}(k)\) ein Polynom in \(k\). Daher kann
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
auf Potenzsummen \(\sum_{k=1}^{\kappa_{i,j}} k^p\) zurückgeführt werden.
Nicht jeder Index \(k\) darf eine neue quadratfreie Klasse beginnen. Manche Dreieckszahlen tauchen bereits in einer älteren Familie auf:
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
für ein kleineres \(k_1\). Solche \(k\) sind nicht-fundamental. Würde man blind über alle \(k\ge 1\) summieren, entstünden Mehrfachzählungen.
Diese doppelten Basisindizes entstehen durch dieselbe Pell-Mechanik. Mit \(x_1=2k_1+1\) definiert man
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$
Jeder so erzeugte Wert \(k_r\) gehört zur bereits von \(k_1\) repräsentierten quadratfreien Klasse und muss daher vom vollständigen Polynom-Summenwert subtrahiert werden.
Nehmen wir die kleinste fundamentale Basis \(k=1\). Dann gilt
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
Die zugehörige Dreieckszahlfamilie beginnt also mit
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
Bildet man daraus Paarprodukte, erhält man
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$
Damit sind \(6\), \(35\) und \(210\) allesamt Quadrat-Dreiecks-Produkte.
Dieselbe Familie erklärt auch, warum manche spätere Basen nicht fundamental sind. Die Pell-Rekurrenz zu \(x_1=3\) liefert
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
und damit
$$k=0,1,8,49,\dots$$
sowie tatsächlich
$$T_8=36,\qquad T_{49}=1225.$$
Diese Zahlen standen bereits in der Familie von \(k=1\), also dürfen \(k=8\) und \(k=49\) nicht noch einmal als neue Startpunkte gezählt werden.
Bezeichnet \(\mathcal{N}_{i,j}\) die nicht-fundamentalen Basen bis zur jeweiligen Schranke des Paares \((i,j)\), dann ist die Zielsumme
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
Genau diese Struktur setzen die Programme um.
Die C++-, Python- und Java-Implementierungen folgen derselben Abfolge. Zuerst werden alle relevanten Familienpaare \((i,j)\) erzeugt, indem die Folge \(U_m(3)\) so lange fortgesetzt wird, bis ihre Produkte sicher zu groß sind. Für jedes verbleibende Paar wird die größte zulässige Basis mit einer monotonen Binärsuche über die Bedingung \(C_{i,j}(k)\le N\) bestimmt.
Danach werden die Chebyshev-Faktoren symbolisch als Polynome in \(k\) aufgebaut. Durch Multiplikation mit dem Dreieckszahlfaktor \(T_k=\frac{k(k+1)}{2}\) entsteht eine Polynomdarstellung jedes Beitrags \(C_{i,j}(k)\). Die Summation über den vollständigen Bereich erfolgt dann über Potenzsummen; diese werden modulo \(136101521\) mit einer Lagrange-Interpolation der Faulhaber-Polynome ausgewertet.
Zum Schluss erzeugen die Implementierungen die Liste aller nicht-fundamentalen Basen nur einmal, werten die Paarbeiträge an diesen Stellen aus und subtrahieren sie. Die C++-Version parallelisiert diese Korrekturphase, während Python und Java dieselbe Rechnung sequentiell durchführen.
Sei \(F\) die Zahl der Familienstufen mit \(U_F(3)\le N\). Wegen des exponentiellen Wachstums von \(U_m(3)\) gilt \(F=O(\log N)\), also gibt es \(O(F^2)\) relevante Familienpaare. Sei \(B\) die größte zulässige Basis über alle Paare, \(D\) der maximale Polynomgrad und \(L\) die Anzahl der nicht-fundamentalen Basen bis \(B\).
Die Erzeugung der Indexpaare kostet \(O(F^2)\). Das Finden aller Paargrenzen benötigt \(O(F^2\log B)\) Auswertungen, jeweils mit einer kurzen Chebyshev-Rekurrenz. Polynom-Vorbereitung und Potenzsummen sind nur von kleinem Grad abhängig und haben daher eine niedrige Gesamtkomplexität in \(F\). Die Korrekturphase kostet in der hier verwendeten direkten Form \(O(F^2L)\). Der Speicherbedarf wird von Paarmetadaten, Polynomkoeffizienten und der Korrekturliste dominiert und ist daher \(O(F^2D+L)\).
Praktisch entscheidend ist, dass die Familientiefe nur logarithmisch wächst. Genau das macht die Rechnung auch bei \(N=10^{35}\) handhabbar.
\(T_n=\frac{n(n+1)}{2}\) \(n\). üçgensel sayı olsun. Problem,
$$T_aT_b=c^2$$
eşitliğini sağlayan iki farklı üçgensel sayı \(T_a\) ve \(T_b\) için ortaya çıkan tüm pozitif \(c\le N\) değerlerinin toplamı \(S(N)\)'i istiyor. Gerçek sınır \(N=10^{35}\) olduğundan, üçgensel sayı çiftlerini doğrudan taramak mümkün değildir. Bu yüzden çözüm, üçgensel sayıları karekökten arındırılmış ortak sınıflara ayırır, her sınıfı Pell tipinde bir aile olarak üretir ve elde edilen karekökleri \(136101521\) modunda toplar.
Ana gözlem şudur: İki üçgensel sayının çarpımı ancak ve ancak karekökten arındırılmış kısımları aynıysa tam kare olur.
Eğer \(T_aT_b\) bir kare ise, \(T_a\) ve \(T_b\) aynı karekökten arındırılmış çekirdeğe sahiptir. Dolayısıyla
$$T_a=d\,r^2,\qquad T_b=d\,s^2$$
şeklinde yazabiliriz; burada \(d\) karekökten arındırılmıştır. O zaman
$$c=d\,rs$$
olur. Demek ki sabit bir \(d\) için çözmemiz gereken şey
$$T_n=d\,y^2$$
şeklindeki üçgensel sayılardır. \(T_n=\frac{n(n+1)}{2}\) ve \(x=2n+1\) yazınca denklem
$$x^2-8dy^2=1$$
haline gelir. Yani her karekökten arındırılmış sınıf bir Pell denklemiyle kontrol edilir.
Bir sınıfın temel temsilcisi \(T_k=d\) olsun. O zaman
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
Bu tabandan türeyen Pell ailesi, aynı karekökten arındırılmış sınıftaki tüm üçgensel sayıları verir. Uygulama kare çarpanlarını ikinci tür Chebyshev polinomlarıyla ifade eder:
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
Böylece sabit bir \(k\) için aile
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
şeklindedir. Bunlar tam olarak karekökten arındırılmış kısmı \(T_k\) olan üçgensel sayılardır. İlk terim \(U_0=1\) olduğu için \(A_0(k)=T_k\) olur.
Aynı aileden iki farklı terim seçelim: \(A_i(k)\) ve \(A_j(k)\), burada \(i<j\). Çarpımları
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2$$
olur. Dolayısıyla \(S(N)\)'e katkı yapan her değer
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j$$
biçimindedir. Problem artık şuna dönüşür: tüm aile çiftlerini \((i,j)\) dolaş, tüm temel tabanlar \(k\) üzerinde \(C_{i,j}(k)\) değerlerini topla ve yalnızca \(N\)'yi aşmayanları tut.
Sabit \(i\) ve \(j\) için \(C_{i,j}(k)\), \(k\) arttıkça büyür; çünkü hem \(T_k\) hem de \(U_m(2k+1)\) \(k\ge 1\) için artandır. En küçük taban \(k=1\) olduğundan \(2k+1=3\) olur. Bu nedenle bir \((i,j)\) çiftinin herhangi bir katkı yapabilmesi için gerekli koşul
$$U_i(3)\,U_j(3)\le N$$
olmalıdır. Bu yüzden uygulamalar önce
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
dizisini üretir ve yalnızca en küçük olası katkısı yeterince küçük olan indeks çiftlerini saklar.
\((i,j)\) sabitlenince,
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
vardır. \(U_m(2k+1)\), \(k\) cinsinden \(m\). dereceden bir polinom olduğundan \(C_{i,j}(k)\) de \(k\) cinsinden bir polinomdur. Bu nedenle
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
kuvvet toplamlarının \(\sum_{k=1}^{\kappa_{i,j}} k^p\) doğrusal birleşimine indirgenebilir.
Her \(k\) yeni bir karekökten arındırılmış sınıf başlatmaz. Bazı üçgensel sayılar daha eski bir ailenin içinde zaten vardır:
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
ve burada \(k_1<k\) olur. Böyle \(k\) değerleri temel olmayan tabanlardır. Eğer tüm \(k\ge 1\) için körlemesine toplarsak aynı sınıfları tekrar saymış oluruz.
Bu tekrar eden tabanlar yine aynı Pell mekanizmasından gelir. \(x_1=2k_1+1\) ile başlayıp
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}$$
tanımını yaparız. Böyle oluşan her \(k_r\), zaten \(k_1\) tarafından temsil edilen sınıfa aittir; dolayısıyla tam aralık polinom toplamından çıkarılmalıdır.
En küçük temel taban \(k=1\) olsun. O zaman
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
Buna karşılık gelen aile şöyle başlar:
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
Şimdi bu aile üyelerini ikili çarparsak
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2$$
elde edilir. Yani \(6\), \(35\) ve \(210\) gerçekten kare üçgensel çarpım kökleridir.
Aynı aile, daha sonraki bazı tabanların neden temel olmadığını da gösterir. \(x_1=3\) için Pell yinelemesi
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
verir; buradan
$$k=0,1,8,49,\dots$$
çıkar ve gerçekten
$$T_8=36,\qquad T_{49}=1225$$
olduğundan, \(k=8\) ile \(k=49\) zaten \(k=1\) ailesinin içindedir. Bunlar yeni başlangıç olarak tekrar sayılmamalıdır.
\(\mathcal{N}_{i,j}\), \((i,j)\) çifti için ilgili sınırın altındaki temel olmayan tabanlar kümesi olsun. O zaman hedef toplam
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
Programların uyguladığı yapı tam olarak budur.
C++, Python ve Java uygulamaları aynı işlem hattını izler. Önce \(U_m(3)\) dizisi, ürünleri artık işe yaramayacak kadar büyüyene kadar genişletilir ve buradan anlamlı aile çiftleri \((i,j)\) elde edilir. Her kalan çift için, \(C_{i,j}(k)\le N\) koşulunu sağlayan en büyük taban monoton bir ikili aramayla bulunur.
Sonra Chebyshev çarpanları \(k\) cinsinden sembolik polinomlar olarak kurulur. Bunlar üçgensel sayı çarpanı \(T_k=\frac{k(k+1)}{2}\) ile çarpılarak her katkı için polinom katsayıları elde edilir. Tam aralık toplamı, kuvvet toplamları üzerinden hesaplanır; kuvvet toplamları da Faulhaber polinomlarının Lagrange türü enterpolasyonu ile \(136101521\) modunda değerlendirilir.
Son aşamada tüm temel olmayan tabanlar tek bir liste halinde üretilir, bu indekslerdeki katkılar hesaplanır ve tam aralık toplamından çıkarılır. C++ sürümü bu düzeltme aşamasını paralelleştirir; Python ve Java ise aynı hesabı sıralı şekilde yapar.
\(F\), \(U_F(3)\le N\) koşulunu sağlayan aile derinliği olsun. \(U_m(3)\) üstel büyüdüğü için \(F=O(\log N)\) olur; dolayısıyla aile çifti sayısı \(O(F^2)\)'dir. Tüm çiftler içindeki en büyük uygun taban \(B\), en yüksek polinom derecesi \(D\) ve \(B\)'ye kadar olan temel olmayan taban sayısı \(L\) olsun.
İndeks çiftlerini üretmek \(O(F^2)\) zaman ister. Tüm çift sınırlarını bulmak, kısa bir Chebyshev yinelemesi kullanan \(O(F^2\log B)\) kadar değerlendirme gerektirir. Polinom hazırlığı ve kuvvet toplamı aşaması, dereceler küçük olduğu için \(F\) cinsinden düşük maliyetlidir. Düzeltme fazı burada kullanılan doğrudan biçimde \(O(F^2L)\) zaman alır. Bellek kullanımı ise çift bilgileri, polinom katsayıları ve düzeltme listesinden dolayı \(O(F^2D+L)\) mertebesindedir.
Pratikte belirleyici nokta, aile derinliğinin yalnızca logaritmik büyümesidir. Bu sayede \(N=10^{35}\) için bile hesap yapılabilir.
Sea \(T_n=\frac{n(n+1)}{2}\) el \(n\)-ésimo número triangular. El problema pide la suma \(S(N)\) de todos los enteros positivos \(c\le N\) para los que existen dos números triangulares distintos \(T_a\) y \(T_b\) tales que
$$T_aT_b=c^2.$$
Para el valor real \(N=10^{35}\), es imposible recorrer directamente pares de números triangulares. Por eso la solución agrupa los triangulares por su parte libre de cuadrados, describe cada grupo como una familia de tipo Pell y luego suma las raíces cuadradas correspondientes módulo \(136101521\).
La observación central es que el producto de dos números triangulares es un cuadrado exactamente cuando ambos comparten la misma parte libre de cuadrados.
Si \(T_aT_b\) es un cuadrado, entonces \(T_a\) y \(T_b\) tienen la misma parte libre de cuadrados. Podemos escribir
$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$
donde \(d\) es libre de cuadrados. Entonces
$$c=d\,rs.$$
Así, para cada \(d\) fijo, debemos encontrar todos los números triangulares de la forma
$$T_n=d\,y^2.$$
Usando \(T_n=\frac{n(n+1)}{2}\) y \(x=2n+1\), esto se transforma en
$$x^2-8dy^2=1.$$
Por tanto, cada clase libre de cuadrados queda gobernada por una ecuación de Pell.
Supongamos que una clase tiene como representante fundamental \(T_k=d\). Entonces
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
La familia de Pell generada por esa base produce todos los triangulares de la misma clase libre de cuadrados. La implementación expresa los factores cuadrados mediante polinomios de Chebyshev de segunda especie:
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
Para este \(k\) fijo, la familia triangular es
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
Estos son exactamente los números triangulares cuya parte libre de cuadrados es \(T_k\). El primer término es \(A_0(k)=T_k\) porque \(U_0=1\).
Tomemos dos miembros distintos de la misma familia, \(A_i(k)\) y \(A_j(k)\), con \(i<j\). Su producto vale
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$
Por consiguiente, cada contribución a \(S(N)\) tiene la forma
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$
Todo el problema pasa a ser: enumerar los pares de familia \((i,j)\), sumar \(C_{i,j}(k)\) sobre todas las bases fundamentales \(k\) y conservar solo los valores que no superan \(N\).
Para \(i\) y \(j\) fijos, la función \(C_{i,j}(k)\) crece con \(k\), porque tanto \(T_k\) como \(U_m(2k+1)\) crecen para \(k\ge 1\). La base mínima es \(k=1\), donde \(2k+1=3\). Así, una condición necesaria para que el par \((i,j)\) aporte algo es
$$U_i(3)\,U_j(3)\le N.$$
Por eso las implementaciones primero construyen la breve sucesión
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
y solo conservan los pares de índices cuyo aporte mínimo posible todavía cabe dentro del límite.
Una vez fijado \((i,j)\), existe una base máxima \(\kappa_{i,j}\) tal que
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
Como \(U_m(2k+1)\) es un polinomio en \(k\) de grado \(m\), toda la expresión \(C_{i,j}(k)\) también es un polinomio en \(k\). Entonces
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
se reduce a una combinación lineal de sumas de potencias \(\sum_{k=1}^{\kappa_{i,j}} k^p\).
No todo índice \(k\) debe iniciar una nueva clase libre de cuadrados. Algunos triangulares ya aparecen dentro de una familia anterior:
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
para cierto \(k_1\) menor. Esos \(k\) son bases no fundamentales. Si sumáramos ingenuamente sobre todos los \(k\ge 1\), contaríamos varias veces la misma clase.
Los índices duplicados se generan con la misma mecánica de Pell. Partiendo de \(x_1=2k_1+1\), definimos
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$
Cada \(k_r\) obtenido así pertenece a la clase ya representada por \(k_1\), por lo que su contribución debe restarse del total polinómico completo.
Tomemos la base fundamental más pequeña, \(k=1\). Entonces
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
La familia triangular correspondiente empieza con
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
Al emparejar esos términos, obtenemos
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$
Así, \(6\), \(35\) y \(210\) son productos triangulares cuadrados válidos.
La misma familia muestra por qué algunas bases posteriores no son fundamentales. La recurrencia de Pell para \(x_1=3\) produce
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
y por tanto
$$k=0,1,8,49,\dots$$
con
$$T_8=36,\qquad T_{49}=1225.$$
Esos valores ya estaban en la familia de \(k=1\), de modo que \(k=8\) y \(k=49\) no pueden volver a contarse como nuevas bases.
Si \(\mathcal{N}_{i,j}\) denota el conjunto de bases no fundamentales que no superan el límite relevante del par \((i,j)\), entonces la suma buscada es
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
Esa es exactamente la estructura que implementan los programas.
Las implementaciones en C++, Python y Java siguen la misma secuencia. Primero enumeran los pares de familia \((i,j)\) relevantes prolongando la sucesión \(U_m(3)\) hasta que sus productos ya no pueden contribuir. Para cada par superviviente, hallan la base máxima admisible mediante una búsqueda binaria monótona sobre la condición \(C_{i,j}(k)\le N\).
Después construyen simbólicamente los factores de Chebyshev como polinomios en \(k\). Al multiplicarlos por el factor triangular \(T_k=\frac{k(k+1)}{2}\), obtienen una representación polinómica de cada contribución \(C_{i,j}(k)\). La suma completa del rango se evalúa entonces mediante sumas de potencias, y esas sumas de potencias se calculan módulo \(136101521\) usando interpolación de tipo Lagrange sobre los polinomios de Faulhaber.
Por último, las implementaciones generan una sola vez todas las bases no fundamentales, evalúan allí las contribuciones de cada par y las restan. La versión en C++ paraleliza esta fase de corrección, mientras que Python y Java realizan la misma aritmética de forma secuencial.
Sea \(F\) el número de niveles de familia con \(U_F(3)\le N\). Como \(U_m(3)\) crece exponencialmente, se tiene \(F=O(\log N)\), y por tanto el número de pares de familia es \(O(F^2)\). Sea \(B\) la mayor base admisible entre todos los pares, \(D\) el grado polinómico máximo y \(L\) el número de bases no fundamentales hasta \(B\).
Enumerar los pares cuesta \(O(F^2)\). Encontrar todos los límites de par requiere \(O(F^2\log B)\) evaluaciones, cada una con una recurrencia corta de Chebyshev. El preprocesamiento polinómico y la evaluación de sumas de potencias dependen de grados pequeños, así que su coste total es moderado en función de \(F\). La fase de corrección cuesta \(O(F^2L)\) en la forma directa usada aquí. La memoria está dominada por los metadatos de pares, los coeficientes polinómicos y la lista de corrección, de modo que es \(O(F^2D+L)\).
En la práctica, lo decisivo es que la profundidad de familia solo crece logarítmicamente. Eso mantiene el cálculo viable incluso cuando \(N=10^{35}\).
设 \(T_n=\frac{n(n+1)}{2}\) 为第 \(n\) 个三角数。题目要求求出所有满足下面条件的正整数 \(c\le N\) 的总和 \(S(N)\)：存在两个不同的三角数 \(T_a\) 与 \(T_b\)，使得
$$T_aT_b=c^2.$$
真正的上界是 \(N=10^{35}\)，因此绝不可能直接枚举三角数对。可行的方法是先按“平方因子去掉之后的核心部分”对三角数分组，再把每一组写成 Pell 型递推族，最后把对应的平方根贡献在模 \(136101521\) 下累加。
整个问题的核心在于：两个三角数的乘积什么时候会是完全平方数？答案是，当且仅当它们具有相同的平方自由部分。
如果 \(T_aT_b\) 是一个平方，那么 \(T_a\) 和 \(T_b\) 的平方自由部分相同。于是可以写成
$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$
其中 \(d\) 是平方自由整数。这样一来
$$c=d\,rs$$
就自动成立。所以对每一个固定的平方自由 \(d\)，我们只需要找出所有满足
$$T_n=d\,y^2$$
的三角数。利用 \(T_n=\frac{n(n+1)}{2}\)，再令
$$x=2n+1,$$
便得到
$$x^2-8dy^2=1.$$
这说明每一个平方自由类最终都对应到一个 Pell 方程。
设某个平方自由类的最小基本代表是 \(T_k=d\)。则有
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
从这个基本解出发，Pell 递推会生成同一平方自由类中的全部三角数。实现里把平方因子写成第二类 Chebyshev 多项式：
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
于是对这个固定的 \(k\)，整条三角数家族可以写成
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
它们恰好就是平方自由部分等于 \(T_k\) 的所有三角数。因为 \(U_0=1\)，所以首项就是 \(A_0(k)=T_k\)。
现在从同一家族中选两个不同成员 \(A_i(k)\) 与 \(A_j(k)\)，其中 \(i<j\)。它们的乘积满足
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$
因此，每一个对 \(S(N)\) 有贡献的数都具有形式
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$
也就是说，原问题可以重写为：枚举所有家族层数对 \((i,j)\)，对所有基本底数 \(k\) 累加 \(C_{i,j}(k)\)，同时只保留不超过 \(N\) 的贡献。
对固定的 \(i,j\) 而言，\(C_{i,j}(k)\) 随 \(k\) 单调增大，因为 \(T_k\) 和 \(U_m(2k+1)\) 在 \(k\ge 1\) 时都增大。最小的底数是 \(k=1\)，此时 \(2k+1=3\)。所以一个层数组合 \((i,j)\) 想要产生任何有效贡献，必要条件是
$$U_i(3)\,U_j(3)\le N.$$
这正是实现首先构造序列
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
并仅保留最小可能贡献仍不超过上界的配对的原因。
一旦固定 \((i,j)\)，就存在最大的可行底数 \(\kappa_{i,j}\)，满足
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
由于 \(U_m(2k+1)\) 是关于 \(k\) 的 \(m\) 次多项式，\(C_{i,j}(k)\) 也就是关于 \(k\) 的多项式。因此
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
可以化成若干幂和 \(\sum_{k=1}^{\kappa_{i,j}} k^p\) 的线性组合。
并不是每个 \(k\) 都应该被当作新的平方自由类起点。有些三角数其实早已出现在更早的家族中：
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
其中 \(k_1<k\)。这样的 \(k\) 就是非基本底数。如果把所有 \(k\ge 1\) 都直接相加，同一个平方自由类就会被重复统计。
这些重复底数仍然由同样的 Pell 机制生成。设 \(x_1=2k_1+1\)，再定义
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$
这样产生出的每个 \(k_r\)，都属于已经由 \(k_1\) 代表过的平方自由类，因此必须从“全范围多项式求和”的结果中扣掉。
取最小的基本底数 \(k=1\)。这时
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
对应的三角数家族前几项就是
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
于是有
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$
这说明 \(6\)、\(35\)、\(210\) 都是合法贡献。
同一个家族还能解释为什么后面的某些底数并不基本。由 \(x_1=3\) 出发的 Pell 递推给出
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
从而得到
$$k=0,1,8,49,\dots$$
并且确实有
$$T_8=36,\qquad T_{49}=1225.$$
这些值本来就已经出现在 \(k=1\) 的家族中，所以 \(k=8\)、\(k=49\) 不能再次被当作新起点。
记 \(\mathcal{N}_{i,j}\) 为配对 \((i,j)\) 对应界限内的所有非基本底数集合，则目标和为
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
程序实现的正是这个结构。
C++、Python 和 Java 三份实现走的是同一条路线。首先，它们通过不断扩展 \(U_m(3)\) 来枚举所有可能有贡献的层数组合 \((i,j)\)。因为一旦最小可能值都超过上界，该组合就不可能再贡献任何结果。对每一个保留下来的配对，再用单调二分搜索找出满足 \(C_{i,j}(k)\le N\) 的最大底数。
接着，实现把 Chebyshev 因子显式构造成关于 \(k\) 的多项式，再乘上三角数因子 \(T_k=\frac{k(k+1)}{2}\)，从而得到每个 \(C_{i,j}(k)\) 的多项式系数。完整区间上的求和于是转化为若干幂和，而这些幂和是在模 \(136101521\) 下通过 Faulhaber 型多项式的 Lagrange 插值来计算的。
最后，程序只生成一次所有非基本底数，并在这些底数上重新计算配对贡献，然后从完整区间总和中减去。C++ 版本把这一步的减法修正做了并行化；Python 和 Java 版本则按相同的数学过程顺序执行。
设 \(F\) 为满足 \(U_F(3)\le N\) 的家族深度。由于 \(U_m(3)\) 呈指数增长，因此 \(F=O(\log N)\)，家族配对数量就是 \(O(F^2)\)。再设所有配对中最大的可行底数为 \(B\)，最高多项式次数为 \(D\)，直到 \(B\) 为止的非基本底数数量为 \(L\)。
配对枚举代价为 \(O(F^2)\)。所有配对上界的搜索需要 \(O(F^2\log B)\) 次评估，每次评估只运行一段很短的 Chebyshev 递推。多项式预处理与幂和计算都只涉及很小的次数，因此总成本对 \(F\) 来说很低。修正阶段在当前直接实现下是 \(O(F^2L)\)。空间主要消耗在配对信息、多项式系数和修正列表上，因此为 \(O(F^2D+L)\)。
真正关键的事实是：家族层数只按对数增长。正因为这一点，\(N=10^{35}\) 这样的规模仍然可以处理。
Пусть \(T_n=\frac{n(n+1)}{2}\) обозначает \(n\)-е треугольное число. Требуется найти сумму \(S(N)\) всех положительных чисел \(c\le N\), для которых существуют два различных треугольных числа \(T_a\) и \(T_b\), удовлетворяющие равенству
$$T_aT_b=c^2.$$
При реальном значении \(N=10^{35}\) прямой перебор пар треугольных чисел невозможен. Поэтому решение группирует треугольные числа по общему квадратсвободному ядру, описывает каждую группу как Pell-подобное семейство и затем суммирует соответствующие квадратные корни по модулю \(136101521\).
Главная идея состоит в том, что произведение двух треугольных чисел является квадратом тогда и только тогда, когда их квадратсвободные части совпадают.
Если \(T_aT_b\) является квадратом, то \(T_a\) и \(T_b\) имеют одну и ту же квадратсвободную часть. Поэтому можно написать
$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$
где \(d\) квадратсвободно. Тогда автоматически
$$c=d\,rs.$$
Значит, для фиксированного \(d\) нужно описать все треугольные числа вида
$$T_n=d\,y^2.$$
Подставляя \(T_n=\frac{n(n+1)}{2}\) и вводя
$$x=2n+1,$$
получаем уравнение
$$x^2-8dy^2=1.$$
Итак, каждая квадратсвободная категория треугольных чисел задается уравнением Пелля.
Предположим, что у некоторой квадратсвободной категории есть фундаментальный представитель \(T_k=d\). Тогда
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
Pell-семейство, порождённое этой базой, дает все треугольные числа с той же квадратсвободной частью. В реализации квадратные множители выражаются через полиномы Чебышёва второго рода:
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
Для фиксированного \(k\) семейство имеет вид
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
Это в точности все треугольные числа, квадратсвободная часть которых равна \(T_k\). Первый член равен \(A_0(k)=T_k\), поскольку \(U_0=1\).
Возьмем два различных члена одного семейства: \(A_i(k)\) и \(A_j(k)\) при \(i<j\). Их произведение равно
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$
Следовательно, любой вклад в сумму \(S(N)\) имеет форму
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$
Задача тем самым превращается в перечисление всех пар уровней \((i,j)\), суммирование \(C_{i,j}(k)\) по всем фундаментальным базам \(k\) и отсечение вкладов, превышающих \(N\).
При фиксированных \(i\) и \(j\) функция \(C_{i,j}(k)\) возрастает по \(k\), потому что возрастают и \(T_k\), и \(U_m(2k+1)\) при \(k\ge 1\). Минимальная база равна \(k=1\), то есть \(2k+1=3\). Поэтому необходимое условие существования вклада от пары \((i,j)\) таково:
$$U_i(3)\,U_j(3)\le N.$$
Именно поэтому реализации сначала строят короткую последовательность
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
и сохраняют только те пары индексов, у которых минимально возможный вклад ещё не слишком велик.
После фиксации \((i,j)\) существует максимальная допустимая база \(\kappa_{i,j}\), для которой
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
Поскольку \(U_m(2k+1)\) является полиномом по \(k\) степени \(m\), выражение \(C_{i,j}(k)\) тоже является полиномом по \(k\). Поэтому сумма
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
сводится к линейной комбинации степенных сумм \(\sum_{k=1}^{\kappa_{i,j}} k^p\).
Не каждый индекс \(k\) должен начинать новую квадратсвободную категорию. Некоторые треугольные числа уже встречаются внутри более раннего семейства:
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
для меньшего \(k_1\). Такие \(k\) называются нефундаментальными. Если суммировать по всем \(k\ge 1\) без поправки, одни и те же категории будут учтены несколько раз.
Повторяющиеся базы порождаются той же Pell-рекурсией. Если \(x_1=2k_1+1\), то определяем
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$
Каждый полученный таким образом индекс \(k_r\) уже принадлежит классу, представленному базой \(k_1\), поэтому его вклад нужно вычесть из полной полиномиальной суммы.
Возьмем наименьшую фундаментальную базу \(k=1\). Тогда
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
Соответствующее семейство треугольных чисел начинается так:
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
Теперь попарно перемножим эти элементы:
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$
Значит, \(6\), \(35\) и \(210\) действительно входят в искомую сумму.
Это же семейство показывает, почему некоторые более поздние базы не являются фундаментальными. Pell-рекурсия для \(x_1=3\) дает
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
а значит
$$k=0,1,8,49,\dots$$
и действительно
$$T_8=36,\qquad T_{49}=1225.$$
Эти значения уже были в семействе базы \(k=1\), поэтому \(k=8\) и \(k=49\) нельзя снова считать новыми начальными точками.
Если \(\mathcal{N}_{i,j}\) обозначает множество нефундаментальных баз, не превосходящих нужной границы для пары \((i,j)\), то искомая сумма равна
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
Именно эту схему и реализуют программы.
Реализации на C++, Python и Java следуют одному и тому же конвейеру. Сначала они перечисляют все существенные пары уровней \((i,j)\), постепенно наращивая последовательность \(U_m(3)\) до тех пор, пока её произведения не становятся заведомо слишком большими. Для каждой оставшейся пары максимальная допустимая база находится монотонным бинарным поиском по условию \(C_{i,j}(k)\le N\).
Затем Chebyshev-множители строятся символически как полиномы по \(k\). После умножения на треугольный множитель \(T_k=\frac{k(k+1)}{2}\) получается полиномиальное представление каждого вклада \(C_{i,j}(k)\). Сумма по полному диапазону вычисляется через степенные суммы, а сами степенные суммы находятся по модулю \(136101521\) с помощью лагранжевой интерполяции полиномов Фаульхабера.
Наконец, реализации один раз генерируют список всех нефундаментальных баз, пересчитывают значения пар на этих индексах и вычитают их из полной суммы. C++-версия распараллеливает эту корректирующую фазу, а Python и Java выполняют те же вычисления последовательно.
Пусть \(F\) обозначает глубину семейства, то есть максимальный индекс, для которого \(U_F(3)\le N\). Так как \(U_m(3)\) растет экспоненциально, имеем \(F=O(\log N)\), а число пар уровней равно \(O(F^2)\). Пусть \(B\) — наибольшая допустимая база среди всех пар, \(D\) — максимальная степень полинома, а \(L\) — число нефундаментальных баз до \(B\).
Перечисление пар стоит \(O(F^2)\). Поиск всех границ для пар требует \(O(F^2\log B)\) вычислений, каждое из которых использует короткую рекурсию Чебышёва. Полиномиальная подготовка и вычисление степенных сумм зависят только от небольших степеней, поэтому их стоимость по \(F\) невелика. Корректирующая фаза в прямой реализации имеет стоимость \(O(F^2L)\). Память определяется в основном данными о парах, полиномиальными коэффициентами и списком коррекции, то есть имеет порядок \(O(F^2D+L)\).
Практически решающим фактом является логарифмический рост глубины семейства. Именно это делает вычисление выполнимым даже при \(N=10^{35}\).
ليكن \(T_n=\frac{n(n+1)}{2}\) هو العدد المثلثي رقم \(n\). المطلوب هو حساب مجموع جميع الأعداد الصحيحة الموجبة \(c\le N\) التي توجد لها قيمتان مثلثيتان مختلفتان \(T_a\) و\(T_b\) بحيث
$$T_aT_b=c^2.$$
عندما يكون الحد الحقيقي \(N=10^{35}\)، يصبح الفحص المباشر لكل أزواج الأعداد المثلثية مستحيلًا عمليًا. لذلك تجمع الفكرة الرياضية الأعداد المثلثية بحسب الجزء الخالي من المربعات، ثم تصف كل مجموعة على أنها عائلة من نوع Pell، وبعد ذلك تجمع الجذور التربيعية الناتجة بترديد \(136101521\).
الفكرة الحاسمة هي أن حاصل ضرب عددين مثلثيين يكون مربعًا كاملًا إذا وفقط إذا كان لهما الجزء نفسه الخالي من المربعات.
إذا كان \(T_aT_b\) مربعًا كاملًا، فهذا يعني أن \(T_a\) و\(T_b\) يملكان الجزء نفسه الخالي من المربعات. لذلك يمكن كتابة
$$T_a=d\,r^2,\qquad T_b=d\,s^2,$$
حيث إن \(d\) عدد خالٍ من المربعات. وعندها نحصل مباشرة على
$$c=d\,rs.$$
إذن، لكل قيمة ثابتة من \(d\)، نحتاج إلى جميع الأعداد المثلثية التي تأخذ الصورة
$$T_n=d\,y^2.$$
وباستعمال \(T_n=\frac{n(n+1)}{2}\) ووضع
$$x=2n+1,$$
تصبح المعادلة
$$x^2-8dy^2=1.$$
وهكذا نرى أن كل فئة خالية من المربعات تخضع في النهاية لمعادلة Pell.
افترض أن لهذه الفئة ممثلًا أساسيًا هو \(T_k=d\). عندئذٍ
$$d=T_k=\frac{k(k+1)}{2},\qquad x_1=2k+1.$$
العائلة المولدة من هذه القاعدة عبر معادلة Pell تعطي جميع الأعداد المثلثية التي تمتلك الجزء الخالي من المربعات نفسه. وتكتب الشيفرات العوامل التربيعية باستعمال متعددات حدود Chebyshev من النوع الثاني:
$$U_0(x)=1,\qquad U_1(x)=2x,\qquad U_{m+1}(x)=2x\,U_m(x)-U_{m-1}(x).$$
وعند تثبيت \(k\)، تصبح العائلة
$$A_m(k)=T_k\,U_m(2k+1)^2,\qquad m=0,1,2,\dots$$
وهذه هي بالضبط الأعداد المثلثية التي يساوي جزؤها الخالي من المربعات \(T_k\). وأول حد هو \(A_0(k)=T_k\) لأن \(U_0=1\).
خذ حدين مختلفين من العائلة نفسها، وليكنا \(A_i(k)\) و\(A_j(k)\) مع \(i<j\). عندها يكون حاصل ضربهما
$$A_i(k)A_j(k)=\left(T_k\,U_i(2k+1)\,U_j(2k+1)\right)^2.$$
وبالتالي فإن كل مساهمة تدخل في \(S(N)\) تأخذ الصورة
$$C_{i,j}(k)=T_k\,U_i(2k+1)\,U_j(2k+1),\qquad i<j.$$
إذن تتحول المسألة كلها إلى تعداد أزواج مستويات العائلة \((i,j)\)، ثم جمع \(C_{i,j}(k)\) على جميع القواعد الأساسية \(k\)، مع الاحتفاظ فقط بما لا يتجاوز \(N\).
لأجل \(i\) و\(j\) ثابتين، فإن \(C_{i,j}(k)\) دالة متزايدة في \(k\)، لأن كلًّا من \(T_k\) و\(U_m(2k+1)\) يزداد عندما \(k\ge 1\). أصغر قاعدة ممكنة هي \(k=1\)، وعندها \(2k+1=3\). لذلك فشرطًا لازمًا لكي يسهم الزوج \((i,j)\) بشيء هو
$$U_i(3)\,U_j(3)\le N.$$
ولهذا تبدأ التطبيقات ببناء السلسلة القصيرة
$$U_0(3),U_1(3),U_2(3),\dots=1,6,35,204,\dots$$
ثم تحتفظ فقط بأزواج الفهارس التي لا يكون أصغر إسهام ممكن لها أكبر من الحد المطلوب.
وبعد تثبيت \((i,j)\)، توجد أكبر قاعدة مسموح بها \(\kappa_{i,j}\) تحقق
$$C_{i,j}(k)\le N,\qquad 1\le k\le \kappa_{i,j}.$$
ولأن \(U_m(2k+1)\) كثير حدود في \(k\) من الدرجة \(m\)، فإن \(C_{i,j}(k)\) هو أيضًا كثير حدود في \(k\). لذلك يمكن اختزال
$$\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)$$
إلى تركيب خطي من مجاميع القوى \(\sum_{k=1}^{\kappa_{i,j}} k^p\).
ليس كل فهرس \(k\) صالحًا كبداية لفئة جديدة خالية من المربعات. فبعض الأعداد المثلثية يظهر مسبقًا داخل عائلة أقدم:
$$T_k=T_{k_1}\,U_r(2k_1+1)^2$$
حيث \(k_1\) أصغر. هذه القيم من \(k\) تسمى قواعد غير أساسية. ولو جمعنا على جميع \(k\ge 1\) من دون تصحيح، فسنعيد عدّ الفئة نفسها أكثر من مرة.
وتولد هذه القواعد المكررة من آلية Pell نفسها. إذا كان \(x_1=2k_1+1\)، فنعرّف
$$x_0=1,\qquad x_{r+1}=2x_1x_r-x_{r-1},\qquad k_r=\frac{x_r-1}{2}.$$
كل قيمة \(k_r\) تنتج بهذه الطريقة تنتمي إلى الفئة التي سبق أن مثلتها القاعدة \(k_1\)، ولذلك يجب طرح مساهمتها من مجموع كثيرات الحدود الكامل.
لنأخذ أصغر قاعدة أساسية وهي \(k=1\). عندها
$$T_1=1,\qquad 2k+1=3,\qquad U_0(3)=1,\qquad U_1(3)=6,\qquad U_2(3)=35.$$
فتبدأ العائلة المقابلة كما يلي:
$$A_0(1)=1,\qquad A_1(1)=1\cdot 6^2=36,\qquad A_2(1)=1\cdot 35^2=1225.$$
وبأخذ الأزواج داخل هذه العائلة نحصل على
$$1\cdot 36=6^2,\qquad 1\cdot 1225=35^2,\qquad 36\cdot 1225=210^2.$$
إذًا فإن \(6\) و\(35\) و\(210\) كلها مساهمات صحيحة.
والعائلة نفسها تشرح أيضًا سبب كون بعض القواعد اللاحقة غير أساسية. فالتكرار Pell عند \(x_1=3\) يعطي
$$x_0=1,\qquad x_1=3,\qquad x_2=17,\qquad x_3=99,\dots$$
ومن ثم
$$k=0,1,8,49,\dots$$
وبالفعل
$$T_8=36,\qquad T_{49}=1225.$$
وهاتان القيمتان ظهرتا أصلًا داخل عائلة \(k=1\)، لذا لا يجوز عدّ \(k=8\) و\(k=49\) كبدايات جديدة.
إذا رمزنا بـ \(\mathcal{N}_{i,j}\) إلى مجموعة القواعد غير الأساسية الواقعة تحت الحد المناسب للزوج \((i,j)\)، فإن المجموع المطلوب يساوي
$$\boxed{S(N)=\sum_{i<j}\left(\sum_{k=1}^{\kappa_{i,j}} C_{i,j}(k)-\sum_{k\in \mathcal{N}_{i,j}} C_{i,j}(k)\right)\pmod{136101521}.}$$
وهذا هو البناء نفسه الذي تطبقه البرامج.
تتبع تطبيقات C++ وPython وJava المسار الحسابي نفسه. فهي تبدأ بتعداد أزواج مستويات العائلة \((i,j)\) الممكنة من خلال توسيع السلسلة \(U_m(3)\) إلى أن تصبح جداءات حدودها كبيرة أكثر من اللازم. ولكل زوج باقٍ، يُحسب أكبر أساس مسموح به عبر بحث ثنائي أحادي على الشرط \(C_{i,j}(k)\le N\).
بعد ذلك تُبنى عوامل Chebyshev رمزيًا على هيئة كثيرات حدود في \(k\). وبضربها في العامل المثلثي \(T_k=\frac{k(k+1)}{2}\)، نحصل على تمثيل كثير حدود لكل مساهمة \(C_{i,j}(k)\). ثم يجري جمع المجال الكامل عبر مجاميع القوى، وتحسب مجاميع القوى هذه بترديد \(136101521\) باستخدام استيفاء من نوع Lagrange على كثيرات حدود Faulhaber.
وفي النهاية، تُولَّد جميع القواعد غير الأساسية مرة واحدة فقط، ثم تُقيَّم مساهمات الأزواج عند هذه القواعد وتُطرح من المجموع الكامل. نسخة C++ تجعل مرحلة الطرح هذه متوازية، بينما تنفذ نسختا Python وJava العملية نفسها بصورة متسلسلة.
ليكن \(F\) هو عمق العائلة بحيث \(U_F(3)\le N\). وبما أن \(U_m(3)\) ينمو نموًا أسيًا، فإن \(F=O(\log N)\)، وبالتالي يكون عدد أزواج المستويات \(O(F^2)\). ولتكن \(B\) أكبر قاعدة مسموح بها بين جميع الأزواج، و\(D\) أعلى درجة لكثيرات الحدود، و\(L\) عدد القواعد غير الأساسية حتى \(B\).
تعداد الأزواج يكلف \(O(F^2)\). أما إيجاد جميع الحدود القصوى للأزواج فيحتاج إلى \(O(F^2\log B)\) عملية تقييم، وكل تقييم يعتمد على تكرار قصير من نوع Chebyshev. وتجهيز كثيرات الحدود ومجاميع القوى يعتمد فقط على درجات صغيرة، لذلك تبقى كلفته الكلية منخفضة بالنسبة إلى \(F\). مرحلة التصحيح تكلف \(O(F^2L)\) في الصورة المباشرة المستخدمة هنا. وتسيطر على الذاكرة بيانات الأزواج، ومعاملات كثيرات الحدود، وقائمة التصحيح، ولذلك يكون الاستهلاك \(O(F^2D+L)\).
والحقيقة العملية الأهم هي أن عمق العائلة ينمو لوغاريتميًا فقط. وهذا بالذات ما يجعل الحساب ممكنًا حتى عند \(N=10^{35}\).