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

@Victor - this is the best timing I've seen around -- thank you! Do you want to post it this answer also on Stack Exchange?

https://mathematica.stackexchange.com/q/278250/13

I plan to post a bounty there and if no one beats your time I will give then the bounty to you :-)

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

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 of Cases[] 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