4
$\begingroup$

I am a new to Mathematica. My goal is to find many (if not all) positive integer solutions to the equation:

$x^2+y^2=z^2+1$

using Mathematica. However the problem is that I can only find a particular solution not many solutions of the given equation using FindInstance.

Can someone help or give the codes in finding the solutions of the given equation?

Any help will be appreciated.

Improve this question
$\endgroup$
13
  • 1
    $\begingroup$ There are an infinite number of points that satisfy that equation. You can't list all of them. What are you trying to do? $\endgroup$
    – C. E. ♦
    Commented May 5, 2015 at 10:01
  • $\begingroup$ thanks for your comment @Pickett. I want to find many of them using Mathematica. $\endgroup$
    – Jr Antalan
    Commented May 5, 2015 at 10:02
  • $\begingroup$ You can pick values according to Reduce[x^2 + y^2 == z^2 + 1, {x, y, z}, Integers]. $\endgroup$
    – C. E. ♦
    Commented May 5, 2015 at 10:08
  • $\begingroup$ Thanks @Pickett will try it now. $\endgroup$
    – Jr Antalan
    Commented May 5, 2015 at 10:12
  • 1
    $\begingroup$ It seems obvious that triples of the form {1, n, n} and {n, 1, n} are solutions for any integer n. $\endgroup$
    – m_goldberg
    Commented May 5, 2015 at 12:06

5 Answers 5

Reset to default
5
$\begingroup$

Diophantus's approach

Diophantus (Book II, problem 9) gives parameterized solutions to x^2 + y^2 == z^2 + a^2, here parametrized by C[1], which may be a rational number (different than 1). We can use his method to find solutions to the OP's case, a == 1. Since Diophantus' method produces rational solutions, we have to clear denominators to get a solution in integers. This complicates searching for solutions with a == 1.

ClearAll[diophantusII9, opcase];
diophantusII9[x_, y_] := Module[{t, t0},
   t0 = Solve[(t + x)^2 + (C[1] t - y)^2 == x^2 + y^2, t];
   {x, y, t + x, C[1] t - y} /. Last@t0
   ];
opcase[seedx_, seedy_] :=
  Module[{p, sol},
   sol = diophantusII9[seedx, seedy];
   Select[
    LCM @@ Denominator[#] # & /@ Abs[sol /. Flatten[
        Solve[# == 1, C[1], Rationals] & /@ Numerator@Together@sol,
        1]],
    MemberQ[#, 1] &]
   ];

Example:

diophantusII9[5, 7]
opcase[5, 7]
(*
  {5, 7, 5 + (2 (-5 + 7 C[1]))/(1 + C[1]^2), -7 + (2 C[1] (-5 + 7 C[1]))/(1 + C[1]^2)}

  {{25, 35, 43, 1}}
*)

Now, it is not the case that every opcase[x, y] produces solutions. A necessary condition is for there to be a solution with a == 1 is that

(x + Sqrt[x^2 + y + y^2])/y

is an integer. Searching up to x or y equal to 50 yields 16 solutions:

With[{$max = 50},
     Flatten[
      opcase @@@ (Position[
          Table[IntegerQ[(x + Sqrt[x^2 + y + y^2])/y], {x, 2, $max}, {y, 
            2, $max}], True] + 1),
  1]
 ]
Subtract @@@ Total[Partition[%^2, {Length[%], 2}], {4}]   (* verify *)
Length@First[%]
(*
  {{10, 15, 18, 1}, {25, 35, 43, 1}, {40, 55, 68, 1}, {55, 75, 93, 1},
   {221, 119, 251, 1}, {70, 95, 118, 1}, {85, 115, 143, 1}, {100, 135, 168, 1},
   {115, 155, 193, 1}, {130, 175, 218, 1}, {476, 255, 540, 1}, {145, 195, 243, 1},
   {1184, 407, 1252, 1}, {160, 215, 268, 1}, {175, 235, 293, 1}, {731, 391, 829, 1}}

  {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}

  16
*)

Improvement to Diophantus' approach

Note that the approach above misses solutions, such as {4, 7, 8} or {14, 17, 22}. Indeed, I do not think it can produce {14, 17, 22}. Applying a little thought to Diophantus's method yields the following improvement to finding nontrivial solutions. (A trivial solution is one in which x or y is 1; there are no nontrivial solutions with x and y less than 4.)

ClearAll[nontrivial];
nontrivial[a_] /; a >= 4 := Sort /@ Join[
  With[{d = Divisors[a + 1]},
   Select[{a, (# (a - 1) - (a + 1)/#)/2} & /@ d[[2 ;; -2]], VectorQ[#, IntegerQ] &]
   ],
  With[{d = Divisors[a - 1]},
   Select[{a, (# (a + 1) - (a - 1)/#)/2} & /@ d[[2 ;; -2]], VectorQ[#, IntegerQ] &]
   ]
  ]

Example:

Join @@ Table[nontrivial[n], {n, 4, 10}] // DeleteDuplicates
(*  {{5, 5}, {4, 7}, {7, 11}, {8, 9}, {9, 19}, {10, 15}}  *)

Here we find several thousand nontrivial solutions:

(sols = Join @@ Table[nontrivial[n], {n, 4, 1000}] // DeleteDuplicates) //
   Length // AbsoluteTiming
(*  {0.160005, 6602}  *)

AllTrue[sols, IntegerQ@Sqrt[#.# - 1] &]
(*  True  *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ May I ask where can I find a English translation of "Diophantus (Book II, problem 9)"? $\endgroup$
    – Kattern
    Commented May 8, 2015 at 2:49
5
$\begingroup$

FindInstance supports an optional argument that defines how many instances it should find. For example, this code will find 10 distinct solutions to your equation:

{x, y, z} /. 
 FindInstance[
  x^2 + y^2 == z^2 + 1 && x > 0 && y > 0 && z > 0, {x, y, z}, 
  Integers, 10]

{{2, 1, 2}, {600, 1, 600}, {1, 69, 69}, {114, 1, 114}, {92, 1, 92},
{558, 1, 558}, {1, 932, 932}, {1, 1, 1}, {996, 1, 996}, {1, 132, 132}}
Improve this answer
$\endgroup$
3
  • $\begingroup$ exactly what I am looking for, thanks @shrx $\endgroup$
    – Jr Antalan
    Commented May 5, 2015 at 12:56
  • 2
    $\begingroup$ To find the nontrivial solutions: FindInstance[ x^2 + y^2 == z^2 + 1 && ! (x == 1 && y == z) && ! (x == z && y == 1) && 0 < x && 0 < y && 0 < z, {x, y, z}, Integers]. $\endgroup$
    – 2012rcampion
    Commented May 5, 2015 at 15:07
  • $\begingroup$ even simpler, FindInstance[ x^2 + y^2 == z^2 + 1 && 1 < x && 1 < y && 1 < z, {x, y, z}, Integers] but note this and @2012rcampion form can only apparently find the one {4,7,8} solution. $\endgroup$
    – george2079
    Commented May 5, 2015 at 19:21
2
$\begingroup$

Solution that finds pairs where both x,y > 1

The other solutions miss values where both x>1 and y>1. For example, x=11, y=7, z=13.

First, iterate through the pairs of integers. I order the pairs {1,1}, {2,1}, {2,2}, {3,1}, {3,2}, ...:

pair[1] = {1, 1}
pair[n_] /; MatchQ[pair[n - 1], {x_, x_}] := pair[n] = {1, 0} pair[n - 1] + {1, 1}
pair[n_] := pair[n] = pair[n - 1] + {0, 1}

Table[pair[i], {i, 10}]
(* {{1, 1}, {2, 1}, {2, 2}, {3, 1}, {3, 2}, {3, 3}, {4, 1}, {4, 2}, {4, 3}, {4, 4}} *)

Now test a pair for whether it satisfies x^2 + y^2 == z^2 + 1:

testXY[{x_, y_}] := IntegerQ[Sqrt[x^2 + y^2 - 1]]

testXY[{11,7}]
(* True *)

Now generate a bunch of results:

Last@Reap[
  Do[ 
     If[testXY[pair[i]], Sow[pair[i]]]
  , {i, 1, 100}]
]

(* {{{1, 1}, {2, 1}, {3, 1}, {4, 1}, {5, 1}, {5, 5}, {6, 1}, {7, 1}, 
    {7, 4}, {8, 1}, {9, 1}, {9, 8}, {10, 1}, {11, 1}, {11, 7}, {12, 1}, 
    {13, 1}, {13, 11}, {14, 1}}}
*)

Obviously you'd have to run the algorithm forever to find all values.

FindInstance doesn't work

If I run FindInstance but require both x>1 and y>1 then it gives an error:

{x, y, z} /. 
 FindInstance[
   x^2 + y^2 == z^2 + 1 && x > 1 && y > 1 && z > 1, {x, y, z}, 
   Integers, 5]

(* The methods available to FindInstance are insufficient to find the requested instances or prove they do not exist. >> *)

I'd love to know why. Seems like a bug.

Diophantus' version

Michael E2 points out that Diophantus was interested in x^2 + y^2 == z^2 + a^2. Here's some code:

triple[1] = {1, 1, 1};
triple[n_] /; MatchQ[triple[n - 1], {x_, x_, x_}] := 
      triple[n] = {1, 0, 0} triple[n - 1] + {1, 1, 1}
triple[n_] /; MatchQ[triple[n - 1], {x_, y_, y_}] := 
      triple[n] = {1, 1, 0} triple[n - 1] + {0, 1, 1}
triple[n_] := triple[n] = triple[n - 1] + {0, 0, 1}

testXYZ[{x_, y_, z_}] := IntegerQ[Sqrt[x^2 + y^2 - z^2]]

Run (ignoring trivial cases where y==z):

Last@Reap[
   Do[
       If[
          testXYZ[triple[i] && Length@Union@triple[i] == 3], 
          Sow[triple[i]]],
     {i, 1, 1000}
   ]]

(* e.g. {17, 9, 3} *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ I've no idea, how I came to this, but as you're obviously right, I'll take down my answer. $\endgroup$
    – LLlAMnYP
    Commented May 5, 2015 at 15:04
  • 1
    $\begingroup$ Nice, Diophantus' solution to his problem Book II, 9, is $10^2+15^2 =18^2 + 1$. (His problem is to solve $x^2+y^2 = z^2+a^2$ and his method does not always yield $a=1$) $\endgroup$
    – Michael E2
    Commented May 5, 2015 at 15:16
  • $\begingroup$ Note that @2012rcampion has a use of FindInstance above that does find non-trivial solutions. Not sure why theirs doesn't work and yours doesn't, though. $\endgroup$
    – Michael Seifert
    Commented May 5, 2015 at 15:16
2
$\begingroup$

the simple approach works ok if you just need the first few thousand..

 allxyz = Append[##, Sqrt[Total[#^2] - 1]] & /@ 
            Select[ Tuples[ Range[2, 1000], {2}] , 
            IntegerQ[Sqrt[Total[#^2] - 1]] & ];


 Show[{ListPlot[allxy[[All, 1 ;; 2]]]}]

enter image description here

Better...FindInstance is much happier if we give it a z value and find x,y: This generates 45,000 solutions in about 30 seconds.

 allxyz = Flatten[Table[  
            {x, y, z} /. # & /@ 
             FindInstance[x^2 + y^2 == z^2 + 1 && 1 < x && 1 < y,
                   {x, y}, Integers, 100], {z, 10000}], 1];
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks @george2079. Doing your codes in my mathematica using findinstance does not give an output, why is that? $\endgroup$
    – Jr Antalan
    Commented May 5, 2015 at 22:36
  • $\begingroup$ version issue? Im using 10.1 $\endgroup$
    – george2079
    Commented May 6, 2015 at 3:49
  • $\begingroup$ That maybe is the problem since Im using version 7. $\endgroup$
    – Jr Antalan
    Commented May 6, 2015 at 8:57
  • $\begingroup$ I like your way of simplifying the job for FindInstance; the graphic is a nice insight too, thank you. $\endgroup$
    – djp
    Commented May 7, 2015 at 14:13
2
$\begingroup$

As this is a PE question, my answer is deliberately vague and incomplete. Consider rewriting the equation as $z^2 - y^2 = x^2 - 1$. Expand the left-hand side as $(z-y)(z+y)=p*q=x^2-1$, and identify $p=z-y$ and $q=z+y$. Thus, $z=(q+p)/2$ and $y=(q-p)/2$. The restriction of non-decreasing $x\le y\le z$ implies a lower bound on $q$ from $y=(q-p)/2\ge x$.

Use the following Mr.Wizard function DivisorPairs to get $\{p,q\}$ from $n=x^2-1$.

DivisorPairs[n_] := 
    Thread[{#, Reverse[#]}][[ ;; Ceiling[Length[#]/2]]] &[Divisors[n]]

Keep only those divisor pairs at or above the lower limit on $q$, then form triples $\{x,y,z\}=\{2x,q-p,q+p\}/2$. There are different criteria for selecting divisor pairs which give integer solutions, depending on the parity of $x$. Still, the final codes for even and odd $x$ are only two or three lines each.

Functions like Reduce, FindInstance, and Solve are amazing, but can be slow. Code based on divisor pairs will return all solutions as requested in this question's title, without a guess at the number of solutions such as required by FindInstance.

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.