I'm surprised no one bothered to pursue the method given in Simons and Alder's paper, as commented by KConrad. A *Mathematica* implementation of their idea is surprisingly short:

```
n = 10;
MapThread[With[{z = (#1 + 1)/2}, 3 #2 + 2 z - 2 + {z, 0}] &,
LinearRecurrence[{10, -1}, #, {1, n} + 1] & /@ {{1, 5}, {0, 2}}]
{{13, 10}, {133, 108}, {1321, 1078}, {13081, 10680}, {129493, 105730},
{1281853, 1046628}, {12689041, 10360558}, {125608561, 102558960},
{1243396573, 1015229050}, {12308357173, 10049731548}}
2 (# + 1) # + 1 & /@ %[[All, 1]]
{365, 35645, 3492725, 342251285, 33537133085, 3286296790925, 322023548377445,
31555021444198565, 3092070077983081805, 302991312620897818205}
```

As mentioned, the problem can be reduced to solving a certain Pell equation. I used `LinearRecurrence[]`

to generate the Pell solutions, after which they are transformed to the first terms of the sum of consecutive squares. To check, say, the second term of the sequence:

```
133^2 + 134^2 == 108^2 + 109^2 + 110^2 == 35645
True
```

I recommend reading that fascinating paper if you can.

Not being contented with what I've written above, I've decided to write a full *Mathematica* implementation of the Simons and Alder method. The function `consecutiveSquareSum[n, k]`

, given below, will list the first `k`

triples consisting of the first term of the sums of `n`

and `n + 1`

consecutive squares, and the sum itself:

```
SquarefreePart[n_Integer?Positive] :=
Times @@ Power @@@ MapAt[Mod[#, 2] &, FactorInteger[n], {All, 2}]
(* generate k solutions to the Diophantine equation x^2 - n y^2 == 1 *)
pellSolutions[n_Integer, k_Integer?Positive] := Module[{p, q, r, u, v},
r = Length[Last[ContinuedFraction[Sqrt[n]]]];
If[OddQ[r], r *= 2];
{p, q} =
Through[{Numerator, Denominator}[Last[Convergents[Sqrt[n], r]]]];
{u, v} = {2 p, n q^2 - p^2};
Transpose[LinearRecurrence[{u, v}, #, {2, k + 1}] & /@
{{1, p}, {0, q}}]]
consecutiveSquareSum[n_Integer?Positive, k_Integer?Positive] :=
Module[{a, a2b, b, cQ, d}, cQ = MatchQ[Mod[n, 4], 1 | 2];
a2b = If[cQ, n (n + 1)/2, n (n + 1)/4];
b = SquarefreePart[a2b]; a = Sqrt[a2b/b]; d = If[cQ, 2 b, b];
Block[{z = (#1 + 1)/2, y}, y = a b #2 + n (z - 1);
{y + z, y, ((n/3 + 1/2) n + 1/6) n + y (y + n) (n + 1)}] & @@@
pellSolutions[d, k]]
```

As an example, generate the first five sums of $100$ and $101$ consecutive squares:

```
consecutiveSquareSum[100, 5]
{{20201, 20100, 41008358350}, {8140601, 8100200, 6627019056398350},
{3272521201, 3256280300, 1070939533497408458350},
{1315545402001, 1309016600400, 173065970485621168960538350},
{528845979103001, 526221417100500, 27967806961346412612849840638350}}
```

Taking the first triple as an example, this says that $20201^2 + \cdots + 20300^2 = 20100^2 + \cdots + 20200^2 = 41008358350$.

7

This problem can be recast as the Pell equation $v^2-24w^2=1$, whose integral solutions $(v,w)$ come from coefficients of $(5+\sqrt{24})^k=v_k+w_k\sqrt{24}$ for integers $k$. This was done more generally for the problem of describing the integers that are a sum of $n$ and $n+1$ consecutive squares almost 50 years ago, by Alder and Simons. See "$n$ and $n + 1$ Consecutive Integers with Equal Sums of Squares," Amer. Math. Monthly 74 (1967), 28–30. Writing $n(n+1)=a^2b$ where $b$ is squarefree, the task amounts to solving $v^2-4bw^2=1$ in integers $v$ and $w$.

– KConrad – 2014-10-01T12:14:16.037@KConrad Why dont you put together an answer similar to this comment, just a little bit extended? It looks extremely interesting. – VividD – 2014-10-01T12:51:23.713

@KConrad, thank you for pointing that paper out; it was a fun read! Thanks to your suggestion, I was able to come up with an efficient routine. – J. M.'s ennui – 2015-10-22T15:31:31.203

Incidentally, sequences associated with this problem involve nice relations between

– Tito Piezas III – 2016-01-13T02:12:02.807second,third, andfifthpowers. See http://math.stackexchange.com/questions/952216/