General formula to obtain triangular-square numbers
Since all triangular numbers are of the form $\frac{n(n+1)}{2}$, we must have the condition cited $$ m^2=\frac{n(n+1)}{2}\tag{1} $$ Equation $(1)$ is equivalent to $$ 2 = \frac{(2n+1)^2-1}{(2m)^2}\tag{2} $$ According to standard continued fraction theory, a rational approximation to $\sqrt{2}$ as good as $(2)$ must also be an approximant for the continued fraction of $\sqrt{2}$. So we must find an overestimate for $\sqrt{2}$ with an even denominator (continued fraction approximants alternate between over- and under-estimates). The continued fraction for $\sqrt{2}$ is $\{1;2,2,2,\dots\}$, so the sequences of numerators and denominators for the approximants satisfy $$ \begin{array}{}a_k=2a_{k-1}+a_{k-2}&\text{ and }&b_k=2b_{k-1}+b_{k-2}\end{array}\tag{3} $$ where $a_0=a_1=b_1=1$ and $b_0=0$. Computing the first few approximants, we get $$ \frac{1}{1},\frac{3}{2},\frac{7}{5},\frac{17}{12},\frac{41}{29},\frac{99}{70},\dots\tag{4} $$ It follows from $(3)$ that every other denominator is even, and that they correspond to the overestimates. Therefore, every other approximant yields a square triangular number.
Every other term of a sequence that satisfies $(3)$, also satisfies $$ \begin{array}{}a_{2k}=6a_{2k-2}-a_{2k-4}\text{ and }b_{2k}=6b_{2k-2}-b_{2k-4}\end{array}\tag{5} $$ where $a_0=1$, $a_2=3$, $b_0=0$, and $b_2=2$. Let $2n_k+1=a_{2k}$ and $2m_k=b_{2k}$. Applying $(5)$ yields $$ \begin{array}{}n_k=6n_{k-1}-n_{k-2}+2\text{ and }m_k=6m_{k-1}-m_{k-2}\end{array}\tag{6} $$ where $n_0=m_0=0$ and $n_1=m_1=1$.
Using standard recurrence methods, we can solve $(6)$ for $m_k$ to get $$ \begin{align} m_k^2 &=\left(\frac{(3+2\sqrt{2})^k-(3-2\sqrt{2})^k}{4\sqrt{2}}\right)^2\\ &=\frac{(17+12\sqrt{2})^k+(17-12\sqrt{2})^k-2}{32}\tag{7} \end{align} $$ Thus, the sequence $TS_k=m_k^2$.
See this page: Square Triangular Numbers. In particular:
In 1778 Leonhard Euler determined the explicit formula: $$ \left( \frac{(3 + 2\sqrt{2})^k - (3 - 2\sqrt{2})^k}{4\sqrt{2}} \right)^2.$$
The wikipedia page above is very useful. It explains the connection to Pells Equation.
I will post just a sketch of different solution, which I quite like. (With the hope that OP can do the rest of work by himself and/or someone will be able to provide some links to materials where similar approach is used and the details are not omitted.)
I like this approach because a) it has many pictures; b) it can be explained to someone who knows nothing about Pell's equation. (In fact, similar approach can be used for some special cases of Pell's equation, you can have a look here. These notes are written in Slovak language, but I guess that the equations and pictures there are sufficient to understand it.)
I have learned about this from diploma thesis of my classmate some years ago. (I do not know the exact title - I only have a ps file without title pages.) In fact, the notation I am using here and some of the pictures are taken from her thesis. I apologize for the quality of the pictures, but I did not know a good way how to convert them from ps format.
Let us start with the first picture.
$\triangle_m = \square_n$ $\Rightarrow$ $\triangle_{2n-m-1}=2\triangle_{m-n}$
(I guess this is self-explaining - we basically just gave the intersection away. One only has to be careful with one thing - we must not forget to subtract one. We are not working with areas of triangles and squares, but with the number of stones that are placed in them.)
Now we have to solve a new problem - when is one triangular number twice as large as another one. We can try our luck with a picture again.
We see that $2\triangle_p=\triangle_q$ $\Rightarrow$ $\square_{q-p}=\triangle_{2p-q}$.
Combining the above, with some algebraic manipulation, we get that from each solution $(m,n)$ we get a new solution of the form $(3m-4n+1,3n-2m-1)$.
Some details to fill in (I am omitting the proofs):
- If we start with a solution in positive integers different from (1,1), then the new solution is still positive and smaller.
- The only case in which the new solution is not in positive integers is $m=n=1$.
If we prove the above facts then we can argue that from each solution we can descend in finitely many steps to the solution (1,1) and, the other way round, each solution can be generated from (1,1) using the recurrence: $$ \begin{gather*} P_{n+1}=3P_n+4Q_n+1\\ Q_{n+1}=2P_n+3Q_n+1 \end{gather*} $$
There are many ways how to solve this, I will sketch one of them. First, let us substitute $P_n=p_n+\frac12$, $Q_n=q_n$, which leads us to $$ \begin{gather*} p_{n+1}=3p_n+4q_n\\ q_{n+1}=2p_n+3q_n \end{gather*} $$ Now we can notice that $p_{n+1}\pm\sqrt{2}q_{n+1}=(3\pm2\sqrt2)p_n + (4\pm3\sqrt2)q_n=(3\pm2\sqrt2)(p_n\pm\sqrt2q_n)$. Hence both $p_n\pm\sqrt2q_n$ are geometric progressions and with some work we can express both $p_n$ and $q_n$ as a linear combination of $(3\pm2\sqrt2)^n$.
(Side note: Anyone who knows about eigenvectors, eigenvalues and diagonalization can see where the above numbers come from - so there is no need of guessing. The numbers $3\pm2\sqrt2$ are precisely the eigenvalues of the matrix $\begin{pmatrix}3&4\\2&4\end{pmatrix}$.)
P.S. If you have seen this "graphical approach" somewhere, I will be glad to learn about such references. I know only about these papers, where similar approach is used to show that $x^2=2y^2$ and $x^2=3y^2$ have no solutions.
- S. J. Miller, D. Montague: Rational irrationality proofs (Wayback Machine)
- S. J. Miller, D. Montague: Irrationality From The Book