Extending the triangular square numbers sequence​
​by Bill Gosper
An elaboration of a brief email exchange with noted mathematician Professor Mike Hirschhorn, concerning the sequence 0, 1, 36, 1225, 41616, 1413721, 48024900, 1631432881, 55420693056, 1882672131025, 63955431761796, ..., the numbers which are simultaneously triangular and square.
Dear Bill,
If the p_k/q_k are the convergents to sqrt{2}, then (p_kq_k)^2 is a triangular number, and these are the only numbers that are both squares and triangles. They satisfy the four-term recurrence
x_k=35x_{k-1}-35x_{k-2}+x_{k-3}.​

Hi Mike, alternatively the three term inhomogeneous:
x(n)=2+34x(n-1)-x(n-2)
Best regards, Mike
As starkly as possible,
In[]:=
2
n

2
-1
+
n

2
+1
2
2
==
2
2n

2
+1
-
2n

2
-1
4
2
In[]:=
Simplify[ReleaseHold[%]]
Out[]=
True
But my point was, you needn’t know all this. You can just take a 4◊4 determinant of any six consecutive values, and extend the sequence quodlibet. And it’s easy to get started with the conveniently bunched values on both sides of n = 0, choosing the first few squares that happen to be triangular. But not by choosing the first few triangulars that happen to be square! Even simpler: Fred Lunnon points out you only need four consecutive values and 3◊3 determinants to extend the sequence of signed square roots xHnL = ..., - 6, - 1, 0, 1, 6, 35, 204, 1189, 6930, ... .
Comparably superlative regards,
Bill
First, let' s test Mike' s claim about the convergents to
2
.
In[]:=
Convergents[
2
,9]
Out[]=
1,
3
2
,
7
5
,
17
12
,
41
29
,
99
70
,
239
169
,
577
408
,
1393
985

These are the approximations to
2
which are best for their size:
In[]:=
%^2
Out[]=
1,
9
4
,
49
25
,
289
144
,
1681
841
,
9801
4900
,
57121
28561
,
332929
166464
,
1940449
970225

In[]:=
2-%
Out[]=
1,-
1
4
,
1
25
,-
1
144
,
1
841
,-
1
4900
,
1
28561
,-
1
166464
,
1
970225

He claims that the numerators times the denominators of these approximations to 2 are exactly the triangular squares:
In[]:=
Numerator[%%]*Denominator[%%]
Out[]=
{1,36,1225,41616,1413721,48024900,1631432881,55420693056,1882672131025}
Familiar? See http://oeis.org/A001110 for lots of formulas and references. It will turn out that the square roots of these,
In[]:=
%
Out[]=
{1,6,35,204,1189,6930,40391,235416,1372105}
i.e., those numbers which become triangular when squared, obey a simpler recurrence formula. (See http://oeis.org/A001109 .)
Note my “starkly as possible” answer. It’s pretty stiff algebra to confirm that the two sides are algebraically equivalent. The equation is between the binomial of something and the square of something else. It is worth your while to see why these somethings are always integers. Mathematica directly (and rather impressively) finds the general formula for which triangular number equals which square by using Pell’s equation:
In[]:=
Reduce[j^2==k*(k+1)/2,{k,j},Integers]//TraditionalForm
Out[]//TraditionalForm=
c
1
∈∧
c
1
≥0∧k
1
2
-1+
1
2

c
1
(3-2
2
)
+
c
1
(3+2
2
)
∧j
c
1
(3-2
2
)
-
c
1
(3+2
2
)
4
2
∨
c
1
∈∧
c
1
≥0∧k
1
2
-1+
1
2
-
c
1
(3-2
2
)
-
c
1
(3+2
2
)
∧j
c
1
(3-2
2
)
-
c
1
(3+2
2
)
4
2
∨
c
1
∈∧
c
1
≥0∧k
1
2
-1+
1
2

c
1
(3-2
2
)
+
c
1
(3+2
2
)
∧j-
c
1
(3-2
2
)
-
c
1
(3+2
2
)
4
2
∨
c
1
∈∧
c
1
≥0∧k
1
2
-1+
1
2
-
c
1
(3-2
2
)
-
c
1
(3+2
2
)
∧j-
c
1
(3-2
2
)
-
c
1
(3+2
2
)
4
2
c
1
is an arbitrary integer constant
...,-1,0,1,2,...
which we shall rename
n
.
In[]:=
FullSimplify[%]/.a_==b_Solve[a==b,{j,k}]/.C[1]->n
Out[]=
n∈&&j-
n
(3-2
2
)
-
n
(3+2
2
)
4
2
||j
n
(3-2
2
)
-
n
(3+2
2
)
4
2
&&k
1
4
-2-
n
(3-2
2
)
-
n
(3+2
2
)
||k
1
4
-2+
n
(3-2
2
)
+
n
(3+2
2
)
&&n≥0
If FullSimplify were really smart, it would discard the &&n³0 and one of the two j solutions and one of the two k solutions, which are merely swapped if n->-n. The purpose of this note is not Pell’s equation, but rather to give a more elementary solution, once we propose that the triangular squares obey a linear recurrence. The recurrence I gave needs only two consecutive values, say 0 and 1, to produce all the rest:
In[]:=
Nest[Append[#,2+34*#〚-1〛-#〚-2〛]&,{0,1},9]
Out[]=
{0,1,36,1225,41616,1413721,48024900,1631432881,55420693056,1882672131025,63955431761796}
But, as I mentioned, their square roots
In[]:=
recurrenceSquared=
%
Out[]=
{0,1,6,35,204,1189,6930,40391,235416,1372105,7997214}
obey a simpler recurrence. If we knew what it was, we’d only need 0 and 1 to get the rest. Or better yet, solve the recurrence, getting them all at once. We can do all this knowing only the fragment
...,0,1,6,...
because we can extend it to negative:
...,-6,-1,0,1,6,...,
and, for any second order recurrence (of this type), four consecutive values determine the next, because the determinant formed from five consecutive terms vanishes:
In[]:=
Solve[0==Det[{{-6,-1,0},{-1,0,1},{0,1,x}}]]
Out[]=
{{x6}}
In[]:=
Solve[0==Det[{{-1,0,1},{0,1,6},{1,6,x}}]]
Out[]=
{{x35}}
Given a fairly short burst of values, Mathematica (correctly) guesses the general formula.
In[]:=
FindSequenceFunction[{1,6,35,204,1189,6930,40391},n]
Out[]=
-
(4+3
2
)
n
(3-2
2
)
-
n
(3+2
2
)

8(3+2
2
)
In[]:=
FullSimplify[%,n∈Integers]
Out[]=
-
n
(3-2
2
)
+
n
(3+2
2
)
4
2
Instead of the “stiff” algebra of directly verifying the equality of the triangular and square expressions, it is much easier to verify that they individually satisfy the aforementioned second order recurrence:
In[]:=
Simplifyx[n]==2+34x[n-1]-x[n-2]/.x[n_]->
2
n

2
+1
+
n

2
-1
2
2
//ReleaseHold
Out[]=
True
In[]:=
Simplifyx[n]==2+34x[n-1]-x[n-2]/.x[n_]->
2
2n

2
+1
-
2n

2
-1
4
2
//ReleaseHold
Out[]=
True
Then we need only check that both expressions start out 0 and 1 for n = 0 and 1:
In[]:=
Table
2
n

2
+1
+
n

2
-1
2
2
->
2
2n

2
+1
-
2n

2
-1
4
2
,{n,0,2}
Out[]=
00,1
1
32
2
-
2
(-1+
2
)
+
2
(1+
2
)

,
1
8
2

2
(-1+
2
)
+
2
(1+
2
)

-1+
1
4
2

2
(-1+
2
)
+
2
(1+
2
)


1
32
2
-
4
(-1+
2
)
+
4
(1+
2
)


In[]:=
Simplify[%]
Out[]=
{00,11,3636}
The recurrence then constrains them to agree forever. But how does FindSequence function guess recurrence formulæ? (It has to be guessing, for nothing constrains us from throwing in crazy values later in the data.) Suppose we group the triplets of consecutive terms of 0, 1, 6, 35, ... into 2◊2 matrices:
In[]:=
Table[MatrixForm[HankelMatrix[Take[recurrenceSquared,{k,k+1}],Take[recurrenceSquared,{k+1,k+2}]]],{k,3}]
Out[]=

0
1
1
6
,
1
6
6
35
,
6
35
35
204

The recurrence can be mechanized by a magic matrix M which advances each sequence matrix to the next. E.g., for the first step, at least,
In[]:=
M.%〚1〛==%〚2〛
Out[]=
M.
0
1
1
6

1
6
6
35

What M does this? Right-multiply both sides by the inverse of

0
1
1
6

:
In[]:=
MapAt[Inverse,%〚1,2〛,1]
Out[]//MatrixForm=

-6
1
1
0

In[]:=
%%/.MatrixForm[m_]:>MatrixForm[m.%]
Out[]=
M.
1
0
0
1

0
1
-1
6

I.e., the advancer matrix claims to be
In[]:=
%/.Dot(#&)
Out[]=
M
0
1
-1
6

Does it repeatedly advance

0
1
1
6

?
In[]:=
MatrixForm/@NestList
0
1
-1
6
.#&,
0
1
,6
Out[]=

0
1
,
1
6
,
6
35
,
35
204
,
204
1189
,
1189
6930
,
6930
40391

Thus if f(n=0,1,2,...) = 0, 1, 6, ..., advancement looks like
f[1+n]
f[2+n]
f[2+n]
f[3+n]

0
1
-1
6
.
f[n]
f[1+n]
f[1+n]
f[2+n]
Equating lower left elements gives the recurrence
Solving such (linear, constant coefficient) recurrences is routine:
(http://oeis.org/A055997)