20
$\begingroup$

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

enter image description here

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}
$\endgroup$
7
  • 4
    $\begingroup$ You are effectively solving for Pythagorean triples -- so you could do something like: n = 100; Reduce[(a/c)^2 + (b/c)^2 == 1 && 1 < a < n && 1 < b < n && 1 < c < n, {a, b, c} [Element] Integers]. This returns all the triples with fewer than 2 digits in the denominator. I suspect this might not scale well. $\endgroup$
    – bill s
    Commented Jan 6, 2023 at 19:08
  • 5
    $\begingroup$ There's a famous trick to solve this from the world of algebraic geometry where you take a line going through the point {-R,0} (where $R$ is the circle radius), and the line has a rational slope $m$ (which you can set to anything). You intersect the line with with the circle and you get coordinates of rational Pythagorean triples which you can then multiply up to get integer solutions. No need for optimization, brute force, or searching. $\endgroup$
    – flinty
    Commented Jan 6, 2023 at 19:56
  • 4
    $\begingroup$ ^ All the solutions end up looking like {x,y}={(1 - m^2)/(1 + m^2), (2 m)/(1 + m^2)} and since you picked a rational $m$, they are all rational pythag triples. - No need for sqrt or anything like that. $\endgroup$
    – flinty
    Commented Jan 6, 2023 at 20:05
  • 1
    $\begingroup$ @flinty Very cool. Why don't you post this as an answer? $\endgroup$ Commented Jan 6, 2023 at 21:11
  • 1
    $\begingroup$ Maybe obtain via TextRecognize from this 4000 year old spreadsheet. $\endgroup$ Commented Jan 7, 2023 at 12:08

4 Answers 4

9
$\begingroup$

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

enter image description here

$\endgroup$
1
  • 2
    $\begingroup$ Oh, that's pretty cool. I would not have expected that Solve is so fast! $\endgroup$ Commented Jan 7, 2023 at 10:40
19
$\begingroup$

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

enter image description here

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
$\endgroup$
2
  • 3
    $\begingroup$ In other words, the Weierstraß substitution $x=\cos \alpha=\frac{1-t^2}{1+t^2}$ and $y=\sin \alpha=\frac{2t}{1+t^2}$ (with $t=\tan\frac{\alpha}{2}$) gives rational points on the unit circle if the "half-slope" $t$ is rational. $\endgroup$
    – Roman
    Commented Jan 7, 2023 at 11:07
  • 2
    $\begingroup$ Yet another way of sampling, which is a rather natural sampling of the rational numbers, is the Calkin–Wilf tree: q[1] = 1; q[i_] := q[i] = 1/(2*Floor[q[i-1]] - q[i-1] + 1) and p[i_] := {(1-m^2)/(1+m^2), 2m/(1+m^2)} /. m -> q[i]. This sampling does not quite match the "number of digits in the denominator" prescription but it tends heavily in the right direction, in the sense that the size of the denominator tends to increase with i: try ListLogLogPlot[Table[{i, Max[Denominator[p[i]]]}, {i, 10^5}]] to see the trend. $\endgroup$
    – Roman
    Commented Jan 7, 2023 at 11:22
19
$\begingroup$

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 = Internal`Bag[Most[{0}]];
    Do[
     c = Compile`GetElement[cList, i];
     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]]}];
     , {i, start + 1, end}];
    Internal`BagPart[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.

$\endgroup$
9
  • $\begingroup$ +1 Henrik - thanks a lot! How do I make this work Compile`GetElement ? I am trying to run code but it fails. Could you please put explicitly how to load packages or anything else needed. $\endgroup$ Commented Jan 6, 2023 at 22:16
  • 1
    $\begingroup$ Oh sorry, I forgot to add also the function JobPointers. (It computes a (somewhat suboptimal) load balancing for cf.) Now it should work... $\endgroup$ Commented Jan 6, 2023 at 22:23
  • $\begingroup$ Henrik, the cod after Edit still does not run for me. Compile`GetElement looks blue and undefined and errors are: CompiledFunction::cfne: Numerical error encountered; proceeding with uncompiled evaluation. Do::iterb: Iterator {a,1,Quotient[c$17340,2]} does not have appropriate bounds. $\endgroup$ Commented Jan 7, 2023 at 4:15
  • 1
    $\begingroup$ Denominators (c) form a set corresponding to oeis.org/A008846 (Hypotenuses of primitive Pythagorean triangles) which is equal to oeis.org/A004613 (Numbers that are divisible only by primes congruent to 1 mod 4). No mod 12 needed! $\endgroup$
    – kirma
    Commented Jan 7, 2023 at 9:54
  • 2
    $\begingroup$ Ah, interesting. I think that is not related to mod 12, but rather to the fact that I used a too small search interval for a. I modified it to {a, 1, Floor[Sqrt[0.5 cc]]} now. Hopefully, it is correct. Anyways, the method is very clunky and way less efficient as flinty's solution. $\endgroup$ Commented Jan 7, 2023 at 10:36
6
$\begingroup$

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]

enter image description here

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]]
$\endgroup$
3
  • 3
    $\begingroup$ See here for an ultra-fast way of checking if an integer is a perfect square, much faster than your IntegerQ@*Sqrt. $\endgroup$
    – Roman
    Commented Jan 7, 2023 at 11:26
  • $\begingroup$ Thank's, I added this to my answer. $\endgroup$ Commented Jan 7, 2023 at 11:47
  • 2
    $\begingroup$ An additional (small) speedup: use Sow/Reap instead of Table/Nothing/Flatten: Reap[Do[If[CoprimeQ[n1, n2] && sQ[n1^2 + n2^2], Sow[{n1, n2}/Sqrt[n1^2 + n2^2]]], {n1, max}, {n2, n1 + 1, max}]][[2, 1]] $\endgroup$
    – Roman
    Commented Jan 7, 2023 at 12:19

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.