19
$\begingroup$

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?

Improve this question
$\endgroup$
4
  • 7
    $\begingroup$ 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$. $\endgroup$
    – KCd
    Commented Oct 1, 2014 at 12:14
  • $\begingroup$ @KConrad Why dont you put together an answer similar to this comment, just a little bit extended? It looks extremely interesting. $\endgroup$
    – VividD
    Commented Oct 1, 2014 at 12:51
  • $\begingroup$ @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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Oct 22, 2015 at 15:31
  • $\begingroup$ Incidentally, sequences associated with this problem involve nice relations between second, third, and fifth powers. See math.stackexchange.com/questions/952216 $\endgroup$
    – Tito Piezas III
    Commented Jan 13, 2016 at 2:12

7 Answers 7

Reset to default
26
$\begingroup$

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

Improve this answer
$\endgroup$
16
$\begingroup$ 
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

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

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.

Improve this answer
$\endgroup$
0
8
$\begingroup$

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

Improve this answer
$\endgroup$
4
$\begingroup$

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

Improve this answer
$\endgroup$
3
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi ! Can you make this post self contained, e.g. add a few lines of code demonstrating the method outlined here. $\endgroup$
    – Sektor
    Commented Dec 24, 2014 at 8:22
  • $\begingroup$ @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. $\endgroup$
    – individ
    Commented Dec 24, 2014 at 8:38
  • $\begingroup$ 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. $\endgroup$
    – Sektor
    Commented Dec 24, 2014 at 8:40
3
$\begingroup$

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.