Speeding code of plotting solutions to x^2+y^2= z^2 for integer x,y,z

Questions
,
1

“Can the code be improved?”. The first step takes a long time.

First find the values of z = Sqrt[x^2 + y^2] for x and y integers between 1 and 10,000

t = Flatten[Table[{{x, y}, Sqrt[x^2 + y^2]}, {x, 10^4}, {y, 10^4}], 
   1];

Next find the points with z an integer

In[2]:= Length[s = Select[t, IntegerQ[#[[2]]] &]]

Out[2]= 28948

ListPlot {x,y}. Many of the points are multiples of small value solutions.

fermat1.jpg
fermat1.jpg480×446 79.5 KB

Now select the points where x and y have no common factors

In[4]:= Length[ss = Select[s[[All, 1]], GCD[Sequence @@ #] == 1 &]]

Out[4]= 3576

Plot those points

fermat2jpg.jpg
fermat2jpg.jpg480×446 68.3 KB

2

Frank,

How about going from the other direction by generating only numbers that are squares of two other squared integers.

Original code on 2500 points:

In[1]:= AbsoluteTiming[
 t = Flatten[Table[{{x, y}, Sqrt[x^2 + y^2]}, {x, 2500}, {y, 2500}], 
    1];]

Out[1]= {24.1006, Null}

Generate all possible z’s that fit the pattern:

In[2]:= AbsoluteTiming[zs = Map[#^2 &, Range[1, 2500]]; 
 zss = Total /@ Permutations[zs, {2}];]

Out[2]= {5.23768, Null}

Then find the sums of squares that are in the original list of squares:

vals = Intersection[zss, zs];

Now you can sort out which ones you want. Does this help?

Regards,

Neil

3

I found this code on SE gives you instantaneous results (around half second to evaluate the entire notebook including plot and no compilation needed):

Open in Wolfram Cloud ↗
4

Thank you. I must admit I don’t understand how genPTunder works.

5

Related thread .

6

{2 m n, m^2 - n^2, m^2 + n^2}

is from this formula on wiki and you can easily verify that (2mn)^2 + (m^2-n^2)^2 = (m^2 + n^2)^2. The second part of the code is kind of clever way to generate

  • Same pythagorean triples up to common factor
  • Eleminate m^2 + n^2’s larger than the lim

Remove the last part of genPTunder gives:

    genPTunder2[lim_Integer?Positive] := 
     Module[{prim}, 
      prim = Join @@ 
        Table[If[
          CoprimeQ[m, n], {2 m n, m^2 - n^2, m^2 + n^2}, ## &[]], {m, 2, 
          Floor@Sqrt@lim}, {n, 1 + m~Mod~2, m, 2}]]

Notice that we need to add 2*{3,4,5}, 3*{3,4,5}, 4*{3,4,5} and remove {24,7,25} if we set lim to 20. And we also want to have each tuple sorted numertically:

In[9]:= data = genPTunder2[20]
Out[9]= {{4, 3, 5}, {12, 5, 13}, {8, 15, 17}, {24, 7, 25}}

Quotient function generate the common factors that can be used for multiplication of each triple:

In[11]:= {Range[20~Quotient~Max@#], Sort[#]} & /@ data
Out[11]= {{{1, 2, 3, 4}, {3, 4, 5}}, {{1}, {5, 12, 13}}, {{1}, {8, 15,
   17}}, {{}, {7, 24, 25}}}

20/5 = 4 therefore {3, 4, 5} can be multiplied up to four times so that the hypothenus is still no larger than 20. Similiarly \lfloor 20/13 \rfloor = 1, {5,12,13} is the only case; \lfloor 20/25 \rfloor = 0 so Range[20~Quotient~25] returns empty list.

KroneckerProduct works like map each multiplier to a given triple and the empty list is eaten by Union in the end:

In[13]:= KroneckerProduct[{1, 2, 3, 4}, {3, 4, 5}]
Out[13]= {{3, 4, 5}, {6, 8, 10}, {9, 12, 15}, {12, 16, 20}}
In[14]:= KroneckerProduct[{}, {7, 24, 25}]
Out[14]= {}

Because most of the funtion used here are either optimized for integer calculation, automatic multithreading or linear algebra operation, this SE implementation is very efficient for large number

7

Take a look at this, the function primpyths1, generates all primitive Pythagorean triangles with a given input value for x. The remainder of the code creates a table up to 10000, also creates a reversed copy and joins them and finally prints them out. On my pc it takes 1.88 seconds to complete.

The function can easily be altered to generate all Pythagorean triangles, just remove the GCD[x,#]==1 section and change the #<=x to #>x. and if you wish add , Sqrt[x^2 + m[[i]]^2] in the table section.

primpyths1[x_Integer] := (t = Select[Divisors[x^2], # <= x &];
  m = Sort[
    Select[(x^2 - t^2)/(2 t), 
     IntegerQ[#] && # <= x && GCD[x, #] == 1 &]];
  Table[{x, m[[i]]}, {i, 1, Length[m]}])


primp = Flatten[Table[primpyths1[i], {i, 3, 10000}], 1]; prmp1 = 
 RotateLeft /@ primp[[All, 1 ;; 2]]; allp = 
 Riffle[primp, prmp1]; ListPlot[allp, AspectRatio -> 1, 
 GridLines -> {{10000}, {10000}}]
8

Thank you everybody for your help.