Detailed mathematical approach

Problem Summary

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\).

Mathematical Approach

The key is to understand when two triangular numbers have a square product. That happens exactly when they belong to the same squarefree class.

Step 1: Reduce the Condition to a Pell-Type Equation

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.

Step 2: Describe One Family with Chebyshev Polynomials

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\).

Step 3: Turn a Family into Square Triangle Products

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\).

Step 4: Bound the Family Pairs and Sum over \(k\)

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\).

Step 5: Remove Non-Fundamental Bases

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.

Step 6: Worked Example

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.

Step 7: Final Summation Formula

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.

How the Code Works

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.

Complexity Analysis

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}\).

Footnotes and References

1. Project Euler Problem 833: https://projecteuler.net/problem=833
2. Triangular numbers: Wikipedia — Triangular number
3. Pell's equation: Wikipedia — Pell's equation
4. Chebyshev polynomials: Wikipedia — Chebyshev polynomials
5. Faulhaber's formula: Wikipedia — Faulhaber's formula
6. Lagrange interpolation: Wikipedia — Lagrange polynomial