*Crossposted: https://mathematica.stackexchange.com/q/278250/13 *
I’m not sure of the exact performance implications, but using CoprimeQ with just one pass might be better than creating all of the fractions and then selecting based on the length of the denominator. Using a helper function,
AllFractionsWithDenom[d_] := Select[Range@d, CoprimeQ[#, d] &]/d
and then timing this
Flatten[AllFractionsWithDenom /@ Range[1000, 9999]]
takes about 55 seconds on my machine.
On the other hand, timing this
Union[Divide @@@ Flatten[Table[{m, k}, {k, Range[10^4 - 1]}, {m, k - 1}], 1]]
(which doesn’t include the next step to filter down to just 4 digit denominators) takes 259 seconds on my machine.
Dear Vitaliy,
you always come up with interesting problems! My first guesses would be:
- making use of the build in function
FareySequence[]; - using
Select[]instead ofCases[]to avoid inefficient pattern matching;
So I end up with a simple one-liner:
ratUniFarey[n_] := Select[Sqrt[1 - #^2] & /@ FareySequence[10^n - 1], # \[Element] Rationals &]
Making “performance measurement” does not make much sense on my really old hardware. And I do not have an idea on using parallelization here.
Maybe that helps a bit! Best regards – Henrik
Here is a somewhat improved version (using Compile and the better performance of functions being Listable):
cTest = Compile[{{frac, _Integer, 1}},
With[{c = Sqrt[frac[[2]]^2 - frac[[1]]^2]},
If[IntegerPart[c] == c, frac, {-1, -1}]], RuntimeAttributes -> {Listable}];
ratUniFareyC[n_Integer] := Module[{fareyData},
fareyData = Function[frac, {Numerator[frac], Denominator[frac]}, Listable]@FareySequence[10^n - 1];
Rational @@@ Select[cTest@fareyData, # != {-1, -1} &]]
With this I get on my pretty old hardware (Linux makes it possible!):
ratUniFareyC[4] // Length // AbsoluteTiming
(* Out: {101.15, 3188} *)
StackOverflow, at the link above, has already provided a few very fast implementations; here is another one. It uses the CoprimeQ trick from @Eric Rimbey and the fact that if m/n is an x-coordinate of a rational point on a uni circle, then n^2-m^2 is an exact square. We can then hash all the exact squares once to speed up the calculation considerably.
ClearAll[RationalUni, DenominatorLength, PythagoreanTuples]
DenominatorLength[x_] :=
If[IntegerQ[x], 0, Ceiling[Log10[Denominator[x]]]]
PythagoreanTuples[n_, squares_] :=
{#/n, Sqrt[1 - (#/n)^2]} & /@
Select[Range[1, Ceiling[Sqrt[n^2/2]]],
CoprimeQ[n, #] && squares[n^2 - #^2] &]
RationalUni[n_] := Module[{squares},
squares[_] = False;
Scan[(squares[#^2] = True) &, Range[10^n]];
Union@Flatten@Select[
Flatten[
Parallelize[
PythagoreanTuples[#, squares] & /@
Range[10^(n - 1), 10^n - 1]], 1],
DenominatorLength[Last[#]] == n &]]
Let’s test it:
In[220]:= Sort@ratUni[3] == Sort@RationalUni[3]
Out[220]= True
In[221]:= ratUni[4] // Length // AbsoluteTiming
Out[221]= {970.987, 2870}
In[219]:= RationalUni[4] // Length // AbsoluteTiming
Out[219]= {7.80111, 2870}
In[213]:= RationalUni[5] // Length // AbsoluteTiming
Out[213]= {739.865, 28652}
We managed to speed up the calculation for n=4 by over 100 times, measuring on the same CPU, and calculate RationalUni[5] in just over 12 minutes, which is still faster than the original solution for n=4.
macOS Big Sur 11.7.2 /
Wolfram Desktop 13.2 /
3.6 GHz 8-Core Intel Core i9
@Henrik - wonderful idea - thank you very much!
@Victor - this is the best timing I’ve seen around – thank you! Do you want to post it this answer also on Stack Exchange?
I plan to post a bounty there and if no one beats your time I will give then the bounty to you 
…and @Victor – thanks a lot for computing RationalUni[5] – it was unattainable for me. I wonder what the next few points are for n=6,7,8,...etc. but we’d probably need a cluster for that. Time seems rapidly increasing.
Thanks for the kind words, but I don’t think my solution is that fast - e.g., the one using Solve is almost instantaneous for n=5.
E.g.
In[5]:= ParallelTable[{a, b}/c /. # & /@
Solve[a^2 + b^2 == c^2 && a <= b, PositiveIntegers], {c, 5,
100000, 4}] // Flatten[#, 1] & // SortBy[Denominator@*First] //
DeleteDuplicates // Length // AbsoluteTiming
Out[5]= {7.43017, 15919}
Which matches the total of the first five elements of your sequence below:
In[4]:= Total@{1, 15, 142, 1435, 14326}
Out[4]= 15919
I’m a bit late to the game, but…
For an octant there are roughly n solutions for denominators less or equal to 2 \pi n. See https://mathematica.stackexchange.com/a/278283/3056, and in particular A020882 - OEIS, which is the sequence of denominators of solutions to this problem. I don’t believe there is a closed form or a recurrence formulation for this.