# Integers which are the sum of both two and three consecutive squares

This is a math problem I came across the other day:

$365$ can be written as a sum of two and also three consecutive perfect squares: $$365=14^2+13^2=12^2+11^2+10^2$$ What is the next number with this property? Give the last 4 digits of the number. The perfect squares cannot be zero.

I would like to know what would be a good way (especially performance-wise) to check the first, let's say, $1000000$ natural numbers if they can be represented in the way above, using Mathematica?

• 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 Oct 1 '14 at 12:14
• @KConrad Why dont you put together an answer similar to this comment, just a little bit extended? It looks extremely interesting. – VividD Oct 1 '14 at 12:51
• @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. will be back soon Oct 22 '15 at 15:31
• Incidentally, sequences associated with this problem involve nice relations between second, third, and fifth powers. See math.stackexchange.com/questions/952216 – Tito Piezas III Jan 13 '16 at 2:12

You can just step through $i$ and $j$ while trying to simultaneously satisfy $$i^2+(i+1)^2=j^2+(j+1)^2+(j+2)^2$$

Just loop and if the inequality is too small on the left, increment $i$. If it's too small on the right, increment $j$. That looks like this:

Clear[f, g, i, j];
f[i_] = i^2 + (i + 1)^2;
g[j_] = j^2 + (j + 1)^2 + (j + 2)^2;

max = 10^6; i = 1; j = 1;
While[f[i] <= max && g[i] <= max,
If[f[i] == g[j], Print[{i, j, f[i]}]; i++;];
If[f[i] < g[j], i++];
If[f[i] > g[j], j++];
];

(*Output: {13, 10, 365}
{133, 108, 35645} *)


This executes almost instantaneously.

So $133^2+134^2 = 108^2 + 109^2 + 110^2 = 35645$. You can increase max to find more, like these:

{13, 10, 365}
{133, 108, 35645}
{1321, 1078, 3492725}
{13081, 10680, 342251285}
{129493, 105730, 33537133085}


That's up to $10^{12}$, which takes about 10 seconds.

Further discussion

Any useful algorithm here will focus on the $i$ and $j$, rather than the $n$, from $$i^2+(i+1)^2=j^2+(j+1)^2+(j+2)^2=n$$

If you are searching in a straightforward way, with all things equal, checking all possible $i$ and $j$ (keeping in mind that you iterate them together) takes about $\sqrt{n}$ time (whereas checking all possible $n$ takes, well, $n$ time).

You can try something using FindInstance, but even the following:

Timing[FindInstance[i^2 + (i + 1)^2 == j^2 + (j + 1)^2 + (j + 2)^2 &&
i > 0 && j > 0, {i, j}, Integers]]


will still take about ten times as long as the code above.

http://oeis.org/A007667

Clear[fa, ga];
fa = Total[{#, # + 1}^2] &;
ga = Total[{#, # + 1, # + 2}^2] &;


Update: A closed form function

soln = Assuming[{C[1] ∈ Integers && C[1] >= 0 && x > 0 && y > 0},
Simplify@ Reduce[Total[{x, x + 1}^2] == Total[{y, y + 1, y + 2}^2], {x, y}, Integers]] /.
C[1] -> n;

fa@{ToRules[soln]}[[All, All, -1]][[1, 1]] // FullSimplify


ClearAll[genOEISa007667];
genOEISa007667[n_Integer] := 1/8 (10 + 3 (5 - 2 Sqrt[6])^(2 n) (5 + 2 Sqrt[6]) +
(15 - 6 Sqrt[6]) (5 + 2 Sqrt[6])^(2 n))


Finding all in the range $1$ to $10^{197}$ in less than 0.005 seconds:

genOEISa007667/@ Range[100]; // AbsoluteTiming
(* {0.005990,Null} *)


The first 50 numbers:

genOEISa007667/@ Range[50] // Simplify // Column


Previous version:

ClearSystemCache[]
(res = Reduce[{fa[x] == ga[y] && 1 < x < 100000000000000000 &&
1 < y < 100000000000000000}, {x, y}, Integers]); // AbsoluteTiming
(*{0.102073, Null} *)

initials = (res /. Or | And -> List)[[All, All, -1]];
{Row[#, ","], fa[First@#]} & /@ initials // TableForm


• Where is the next set? – Kellen Myers Sep 30 '14 at 8:53
• @kguler very nice and instructive – ubpdqn Sep 30 '14 at 12:20
• Is there any reason, mathematically easy to explain, why the factor between i_n and i_(n+1) seems to converge to 9.8989..? – Guntram Blohm Sep 30 '14 at 12:31
• @GuntramBlohm, good observation. – kglr Sep 30 '14 at 15:00
• Limit[genOEISa007667[n + 1]/genOEISa007667[n], n -> Infinity] gives 49 + 20 Sqrt[6] or 97.98979485566356 – Bob Hanlon Sep 30 '14 at 15:52

Update: the lhs of $i^2+(i+1)^2=j^2+(j+1)^2+(j+2)^2$ is always odd so $j$ should be even.

The fastest solution I found:

max = 10^12;

{(Sqrt[3 + 6 (# + 1)^2] - 1)/2, #, #^2 + (# + 1)^2 + (# + 2)^2} & @@@
(2 Position[#, 0]) &@UnitStep[Abs[# - Round[#]] - 1.*^-10] &[
Sqrt[0.75 + 1.5 #^2] - 0.5] &@ Range[3., Sqrt[max/3.], 2] // AbsoluteTiming

(* {0.023773,
{{13, 10, 365},
{133, 108, 35645},
{1321, 1078, 3492725},
{13081, 10680, 342251285},
{129493, 105730, 33537133085}}} *)


Only 0.02 sec for max = 10^12! It fast because it checks all number simultaneously and uses PackedArrays. I use floating point numbers, but it is enough to obtain the correct result.

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$.

This is not particularly fast but works:

f = (-1 + Sqrt[3 + 6 #^2])/2 &;
q = IntegerQ /@ (f /@ Range[1000000]);
answer = Pick[Range[1000000], q]


yields the triples : {1, 11, 109, 1079, 10681, 105731}

fanswer = f /@ answer


{1, 13, 133, 1321, 13081, 129493}

Dropping first case which is 0^2+1^2+2^2=1^2+2^2:

trip = Rest[{(# - 1), #, # + 1} & /@ answer];
ftrip = Rest[{#, (# + 1)} & /@ fanswer];
s = StringTemplate["\!$$\*SuperscriptBox[\($$, $$2$$]\)"];
func = Function[x,
StringJoin @@ Riffle[TemplateApply[s, #] & /@ x, "+"]];
Grid[MapThread[{#1, "=", #2, "=", ToExpression[#2]} &, {func /@ trip,
func /@ ftrip}]]


yields:

Solutions can be greatly simplified simply computing the next number.

For the equation:

$$n^2+(n+1)^2=k^2+(k+1)^2+(k+2)^2$$

Using the first number. $(p_1 ; s_1) - (1 ; 0 )$

Let's use these numbers. Which are the sequence. The following is found using the previous value according to the formula.

$$p_2=5p_1+12s_1$$

$$s_2=2p_1+5s_1$$

Using the numbers $p,s$ - you can find solutions on formulas.

$$n=p^2\pm6ps+12s^2$$

$$k=6s(2s\pm{p})$$

$$***$$

$$n=-(2p^2\pm6ps+6s^2)$$

$$k=-2p(p\pm3s)$$

So the formula looked compact, you can do the replacement.

$$x=p^2\pm6ps+6s^2$$

$$y=p^2\pm4ps+6s^2$$

Then the solutions are.

$$n=4x^2\pm6xy+3y^2$$

$$k=2x(2x\pm3y)$$

$$***$$

$$n=-(2x^2\pm6xy+6y^2)$$

$$k=-6y(y\pm{x})$$

Or such replacement.

$$z=2ps$$

$$q=p^2+6s^2$$

Then the solutions are.

$$n=-(2q^2\pm6qz+6z^2)$$

$$k=-2q(q\pm3z)$$

$$***$$

$$n=q^2\pm6qz+12z^2$$

$$k=6z(2z\pm{q})$$

Interestingly, all this variety of formulas give the same solution. So that one can restrict the upper formula. The rest of the formula was drawn to show what interesting patterns there.

• Hi ! Can you make this post self contained, e.g. add a few lines of code demonstrating the method outlined here. – Sektor Dec 24 '14 at 8:22
• @Sektor I tried it as easy as possible to write. So you can even use a calculator. I wrote not that interesting calculation, and to show how different calculations can lead to the same result. – individ Dec 24 '14 at 8:38
• I can see that you are trying to write in a clear fashion. I just wanted you to demonstrate that this method actually solves the problem at hand. – Sektor Dec 24 '14 at 8:40

You could always use brute force:

max = 10^12;

(
sq2 = Range[Floor[Sqrt[max/2]]]^2;
sq2 = sq2[[2 ;;]] + sq2[[ ;; -2]];

sq3 = Range[Floor[Sqrt[max/3]]]^2;
sq3 = sq3[[3 ;;]] + sq3[[2 ;; -2]] + sq3[[;; -3]];
Intersection[sq2, sq3]
) // Timing


{0.046800, {365, 35645, 3492725, 342251285, 33537133085}}

Not as elegant as the others, but quite fast: Finding all numbers < 10^18 takes about 22 seconds. (It needs a lot of RAM, though!)

{21.793340, {365, 35645, 3492725, 342251285, 33537133085, 3286296790925, 322023548377445, 31555021444198565}}