# On finding all the positive integral solutions of $x^2+y^2=z^2+1$

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.

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

## 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 *)  • May I ask where can I find a English translation of "Diophantus (Book II, problem 9)"? – Kattern May 8 '15 at 2:49 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}}  • exactly what I am looking for, thanks @shrx – Jr Antalan May 5 '15 at 12:56 • 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]. – 2012rcampion May 5 '15 at 15:07 • 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. – george2079 May 5 '15 at 19:21 ## 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} *)  • I've no idea, how I came to this, but as you're obviously right, I'll take down my answer. – LLlAMnYP May 5 '15 at 15:04 • 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$) – Michael E2 May 5 '15 at 15:16 • 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. – Michael Seifert May 5 '15 at 15:16 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]]]}]  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];  • Thanks @george2079. Doing your codes in my mathematica using findinstance does not give an output, why is that? – Jr Antalan May 5 '15 at 22:36 • version issue? Im using 10.1 – george2079 May 6 '15 at 3:49 • That maybe is the problem since Im using version 7. – Jr Antalan May 6 '15 at 8:57 • I like your way of simplifying the job for FindInstance; the graphic is a nice insight too, thank you. – djp May 7 '15 at 14:13 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.