Optimize search for rational numbers on unit circle?

Question

Crossposted: https://community.wolfram.com/groups/-/m/t/2763509

I'd like to classify rational numbers on a unit circle by the following property: number of digits in the denominator. It relates to an article I read:

https://www.quantamagazine.org/secret-link-uncovered-between-pure-math-and-physics-20171201

I built a function for the search, but it runs terribly slow - any ideas on speeding this up are highly appreciated:

ClearAll[ratUni];
ratUni[n_Integer]:=
(filter out on rational numbers)
Cases[
(find counterpart Pythagorean coordinate on unit square)
 Map[Sqrt[1-#^2]&,
(also remove fractions that simplify to less-than-n digits in denominator)
  Parallelize@Select[
(*build all fractions with denominators and numerators of exactly n digits and remove duplicates*)
Union[Divide@@@Flatten[Table[{m,k},{k,Range[10^n-1]},{m,k-1}],1]],


Length[IntegerDigits[Denominator[#]]]>n-1&]],
_Rational]

Simplest example of rational numbers on unit circle with denominator of 1 digit:

In[]:= ratUni[1]
Out[]= {4/5, 3/5}

Works fine for rational numbers on unit circle with denominator of 2 digits:

In[]:= ratUni[2]
% // Length
Out[]= {84/85, 60/61, 40/41, 63/65, 24/25, 35/37, 12/13, 77/85,
80/89, 15/17, 56/65, 45/53, 55/73, 72/97, 21/29, 20/29, 65/97, 48/73, 28/53, 33/65, 8/17, 39/89, 36/85, 5/13, 12/37, 7/25, 16/65, 9/41, 11/61, 13/85}
Out[]= 30

Still OK for 3 digits but already takes 7 seconds to compute:

In[]:= ratUni[3] // Length // AbsoluteTiming
Out[]= {6.47679, 284}

For 4 digits it is a disastrous 1237 seconds and for 5 digits it hangs. Any help on optimizing it is highly appreciated.

In[]:= ratUni[4] // Length // AbsoluteTiming
Out[]= {1237.74, 2870}

Answer 1

This is a compiled program to find all Pythagorean triples $(a,b,c)$ with $a^2 + b^2 = c^2$ and $a \leq b$ and where c lies between cMin and cMax.

cf = Compile[{{cMin, _Integer}, {cMax, _Integer}},
   Module[{cc, aa, b, bag},
    bag = Internal`Bag[Most[{0}]];
    Do[
     cc = c c;
     Do[
      aa = a a;
      b = Round[Sqrt[N[cc - aa]]];
      If[aa + b b == cc ,
       Internal`StuffBag[bag, a];
       Internal`StuffBag[bag, b];
       Internal`StuffBag[bag, c];
       ];
      , {a, 1, Floor[Sqrt[0.5 cc]]}]
     , {c, cMin, cMax}];
    Internal`BagPart[bag, All]
    ],
   CompilationTarget -> "C"
   ];

So finding the triples with $c \leq 100000$ can be done with this:

{a, b, c} = Transpose[Partition[cf[1, 100000], 3]]; // AbsoluteTiming // First

71.5202

Generating fractional numbers from this and discarding duplicates:

pts = DeleteDuplicates[Transpose[{a/c, b/c}]]; // 
  AbsoluteTiming // First

0.942499

So pts contains all points where the denominator has up to 5 digits.

Note that cf[1,cMax] has quadratical runtime in cMax. So this is not a very efficient method...

Edit

I have found out experimentally, that it might suffice to check only those $c$ that have remainder either $1$ or $5$ after division by $12$. (The margin is to narrow to give proof... ;)

Here a refined, parallelized version of the code above, that allows us to reduce the set of candidates before running the brute-force method:

JobPointers[jobCount_Integer?Positive, threadCount_Integer?Positive] :=
   Ceiling[Subdivide[0, jobCount, Min[threadCount, jobCount]]];
cf = Compile[{{cList, _Integer, 1}, {start, _Integer}, {end, _Integer}},
   Module[{c, cc, aa, b, bag},
    bag = InternalBag[Most[{0}]]; Do[ c = CompileGetElement[cList, i];
     cc = c c;
     Do[
      aa = a a;
      b = Round[Sqrt[N[cc - aa]]];
      If[aa + b b == cc ,
       InternalStuffBag[bag, a]; InternalStuffBag[bag, b];
       InternalStuffBag[bag, c]; ]; , {a, 1, Floor[Sqrt[0.5 cc]]}]; , {i, start + 1, end}]; InternalBagPart[bag, All]
    ],
   CompilationTarget -> "C",
   Parallelization -> True,
   RuntimeAttributes -> {Listable},
   RuntimeOptions -> "Quality"
   ];

It finds the same points in just 2.5 seconds:

cMax = 100000;
threadCount = 8;
candidates = Join[Range[1, cMax, 12], Range[5, cMax, 12]];
perm = PermutationList[RandomPermutation[Length[candidates]]];
candidates = candidates[[perm]];
jobptr = JobPointers[Length[candidates], threadCount];
data = Join @@ cf[candidates, Most[jobptr], Rest[jobptr]];
{a, b, c} = DeleteDuplicates[Transpose[Partition[data, 3]]];
pts = Union[Transpose[{a/c, b/c}]];

If you wonder what the random permutation perm is good for: I use it to improve the load balancing for the parallel computations.

Answer 2

This method is called the Stereographic approach and can be used to generate arbitary rational triples without needing brute force searching. I was fortunate to have learned this trick off Minhyong Kim himself when I met him over a decade ago.

Take the unit circle $x^2 + y^2 = 1$ and a line with a rational slope $m$ which goes through the point $(-1, 0)$. If you intersect this line with the circle, the point of intersection will give you a rational Pythagorean triple which you can multiply up to get integer Pythagorean triples too.

Here's a picture in Mathematica:

Manipulate[
 Plot[m x + m, {x, -1, 3}, PlotRange -> {{-2, 2}, {-2, 2}}, 
  Epilog -> Circle[], AspectRatio -> 1],
 {m, -1, 1}
 ]

The equation of the line is $y = m x + c$ and because we chose $(-1, 0)$ you know that $0 = -m + c$, so we can write the line equation as $y = m x + m$.

If we substitute that into the circle equation we get: $$ x^2 + (mx+m)^2 =1 $$ This has two solutions for $x$, and we already know $-1$ is a solution:

Solve[x^2 + (m x + m )^2 == 1, x]
(* {{x -> -1}, {x -> (1 - m^2)/(1 + m^2)}} *)

Substitute the $x$ solution back in to get $y$

Simplify[(m x + m) /. x -> (1 - m^2)/(1 + m^2)]
(* (2 m)/(1 + m^2) *)

Therefore our line intersects the circle at two points $(-1,0)$ and: $$ \left(\frac{1-m^2}{m^2+1},\frac{2 m}{m^2+1}\right) $$ Since we chose $m$ to be rational, the above point is in $\mathbb{Q}^2$. We can now pick some arbitrary $m$ and try it out, generating as many rational Pythagorean triples as we like:

p = {(1 - m^2)/(1 + m^2), (2 m)/(1 + m^2)} /. m -> 54/93
(* {637/1285,1116/1285} *)
(* verify it works *)
p[[1]]^2 + p[[2]]^2

... and if you multiply by the common denominator here you can get the integer Pythagorean triple too:

637^2 + 1116^2 == 1285^2
(* True *)

With that out of the way, we can now study the length of the denominators. I'm not sure if any care needs to be taken to ensure uniform sampling of the slope $m$. Here I just vary from -1 to 1 and I use IntegerLength instead of Length[IntegerDigits[...]] to avoid allocating the list:

With[{m = Rationalize[#, 0.0001]}, 
    IntegerLength[
      Denominator[{(1 - m^2)/(1 + m^2), (2 m)/(1 + m^2)}]][[1]]] & /@ 
  RandomReal[{-1, 1}, 100000] // Histogram

Another sampling strategy might be to use an angle between -Pi/2 and Pi/2, and Rationalize[Sin[#],0.0001] to get the slope, but I haven't noticed a difference. Yet another strategy would be to 'bounce' the intersecting line around like a ball on a circular pool table, using the last rational intersection point and a random rational slope to compute the next.

If we vary the Rationalize tolerance we have a log relationship between the tolerance denominator length and the denominator length of the triples:

Table[{tol, 
   Mean[With[{m = Rationalize[#, 10^-tol]}, 
       IntegerLength[
         Denominator[{(1 - m^2)/(1 + m^2), (2 m)/(1 + m^2)}]][[
        1]]] & /@ RandomReal[{-1, 1}, 100000]]}, {tol, 1, 
   7}] // ListLinePlot

Answer 3

This is not the fastest solution in the world for finding rational points on the unit circle (takes 4.3 seconds on my laptop for denominators up to 100000), but at least it's short:

ParallelTable[
    {a, b}/c /. # & /@
     Solve[a^2 + b^2 == c^2 && a <= b, PositiveIntegers],
    {c, 5, 100000, 4}] //
   Flatten[#, 1] & // SortBy[Denominator@*First] // DeleteDuplicates
(* {{3/5, 4/5}, {5/13, 12/13}, {8/17, 15/17}, {7/25, 24/25}, ...} *)

This relies on the Diophantine solver used by Mathematica, and the fact that per OEIS A008846 we don't need to check more than every fourth value of c, which is the denominator.

EDIT: Actually the runtime can be cut to about 1.9 seconds on my laptop by using the definition of A004613 for the sequence and thus testing even smaller set of values of c:

ParallelTable[
    {a, b}/c /. # & /@
     Solve[a^2 + b^2 == c^2 && a <= b, PositiveIntegers],
    {c,
     Select[Range[5, 100000, 4],
      AllTrue[FactorInteger[#][[All, 1]], Mod[#, 4] == 1 &] &]}] //
   Flatten[#, 1] & // SortBy[Denominator@*First] // DeleteDuplicates

When looking at the distribution of these unique solutions on basis of their magnitude of denominator, the distribution seems uniform on large scales. It would actually appear that on average there's one such point per $2\pi$ denominator values (see A020882), as can be seen on this plot for ten million first denominators:

Show[
 ListLogLogPlot[MapIndexed[{First[#2], Denominator[First[#1]]} &, %],
  PlotStyle -> PointSize[Large],
  Ticks -> Table[10^n (2 Pi)^d, {d, 0, 1}, {n, 0, 7}],
  GridLines -> Table[10^n (2 Pi)^d, {d, 0, 1}, {n, 0, 7}],
  ImageSize -> 1000],
 LogLogPlot[2 Pi n, {n, 1, Length[%]},
  PlotStyle -> 
   Directive[Thickness[1/300], ColorData[97, "ColorList"][[2]]]]]

Answer 4

As mentioned above, we may search for Phytogorean triples: n1^2+n2^2==n3^2. To simplify, we may restrict the search to the first quadrant, as all the other solutions are mirror images, mirrored at the axis and the origin. We may further restrict the search by noting, that the quarter circles is symmetric to the line x==y, that means if {n1,n2} is a solution, so is {n2,n1}. Therefore we may request: n2>n1 (as n1==n2 is no solution). Further, we may note that if n1 and n2 have a common factor, say a, then we get the same solution {n1/a,n2/a}/Sqrt[(n1/a)^2+(n2/a)^2]. Therefore, we may request that n1 and n2 are coprime. With this we may, e.g. write:

max = 100;
Flatten[Table[
  If[CoprimeQ[n1, n2] && IntegerQ[tmp = Sqrt[n1^2 + n2^2]], {n1/tmp, 
    n2/tmp}, Nothing], {n1, max}, {n2, n1 + 1, max}], 1]

The ultimate step would be to compile the above, but I leave this to soembody else.

Addendum

According to Romans note, an additional speed up can by achieved by using:

sQ[n_] := FractionalPart@Sqrt[n + 0``1] == 0;

with this:

Reap[Do[If[CoprimeQ[n1, n2] && sQ[n1^2 + n2^2], 
    tmp = Sqrt[n1^2 + n2^2]; Sow[{n1/tmp, n2/tmp}]], {n1, max}, {n2, n1 + 1, max}]][[2, 1]]