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

18

10

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?

VividD

Posted 2014-09-30T07:15:06.977

Reputation: 3 500

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 second, third, and fifth powers. See http://math.stackexchange.com/questions/952216/

– Tito Piezas III – 2016-01-13T02:12:02.807

Answers

25

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.

See also

http://oeis.org/A007667

Kellen Myers

Posted 2014-09-30T07:15:06.977

Reputation: 2 558

16

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  

enter image description here

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 

enter image description here


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 

enter image description here

kglr

Posted 2014-09-30T07:15:06.977

Reputation: 302 076

Where is the next set? – Kellen Myers – 2014-09-30T08:53:43.657

@kguler very nice and instructive – ubpdqn – 2014-09-30T12:20:42.213

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 – 2014-09-30T12:31:20.317

@GuntramBlohm, good observation. – kglr – 2014-09-30T15:00:16.543

@GuntramBlohm, ... more data added to check such patterns:) – kglr – 2014-09-30T15:01:28.730

2Limit[genOEISa007667[n + 1]/genOEISa007667[n], n -> Infinity] gives 49 + 20 Sqrt[6] or 97.98979485566356 – Bob Hanlon – 2014-09-30T15:52:40.850

Thank you @Bob. It is amazing that GunthramBlohm got something that close based on the first 10 numbers. – kglr – 2014-09-30T15:57:30.083

hmmm I think 0.005 sec would be good enough :) – VividD – 2014-09-30T18:31:43.863

@BobHanlon: Specifically, it is a square of the fundamental unit $U_6$. Hence $49+20\sqrt{6} = (5+2\sqrt{6})^2$. – Tito Piezas III – 2016-01-13T02:17:56.153

13

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.

ybeltukov

Posted 2014-09-30T07:15:06.977

Reputation: 41 907

7

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

J. M.'s ennui

Posted 2014-09-30T07:15:06.977

Reputation: 115 520

4

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:

enter image description here

ubpdqn

Posted 2014-09-30T07:15:06.977

Reputation: 53 491

3

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.

individ

Posted 2014-09-30T07:15:06.977

Reputation: 141

Hi ! Can you make this post self contained, e.g. add a few lines of code demonstrating the method outlined here. – Sektor – 2014-12-24T08:22:24.753

@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 – 2014-12-24T08:38:30.607

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 – 2014-12-24T08:40:43.713

3

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}}

Niki Estner

Posted 2014-09-30T07:15:06.977

Reputation: 34 978