NAME
Math::PlanePath::GcdRationals  rationals by triangular GCD
SYNOPSIS
use Math::PlanePath::GcdRationals;
my $path = Math::PlanePath::GcdRationals>new;
my ($x, $y) = $path>n_to_xy (123);
DESCRIPTION
This path enumerates X/Y rationals using a method by Lance Fortnow taking a greatest common divisor out of a triangular position.
The attraction of this approach is that it's both efficient to calculate and visits blocks of X/Y rationals using a modest range of N values, roughly a square N=2*max(num,den)^2 in the default rows style.
13  79 80 81 82 83 84 85 86 87 88 89 90
12  67 71 73 77 278
11  56 57 58 59 60 61 62 63 64 65 233 235
10  46 48 52 54 192 196
9  37 38 40 41 43 44 155 157 161
8  29 31 33 35 122 126 130
7  22 23 24 25 26 27 93 95 97 99 101 103
6  16 20 68 76 156
5  11 12 13 14 47 49 51 53 108 111 114
4  7 9 30 34 69 75 124
3  4 5 17 19 39 42 70 74 110
2  2 8 18 32 50 72 98
1  1 3 6 10 15 21 28 36 45 55 66 78 91
Y=0 

X=0 1 2 3 4 5 6 7 8 9 10 11 12 13
The mapping from N to rational is
N = i + j*(j1)/2 for upper triangle 1 <= i <= j
gcd = GCD(i,j)
rational = i/j + gcd1
which means X=numerator Y=denominator are
X = (i + j*(gcd1))/gcd = j + (ij)/gcd
Y = j/gcd
The i,j position is a numbering of points above the X=Y diagonal by rows in the style of Math::PlanePath::PyramidRows with step=1, but starting from i=1,j=1.
j=4  7 8 9 10
j=3  4 5 6
j=2  2 3
j=1  1
+
i=1 2 3 4
If GCD(i,j)=1 then X/Y is simply X=i,Y=j unchanged. This means fractions X/Y < 1 are numbered by rows with increasing numerator, but skipping positions where i,j have a common factor.
The skipped positions where i,j have a common factor become rationals X/Y>1, ie. below the X=Y diagonal. The integer part is GCD(i,j)1 so rational = gcd1 + i/j. For example
N=51 is at i=6,j=10 by rows
common factor gcd(6,10)=2
so rational R = 21 + 6/10 = 1+3/5 = 8/5
ie. X=8,Y=5
If j is prime then gcd(i,j)=1 and so X=i,Y=j. This means that in rows with prime Y are numbered by consecutive N across to the X=Y diagonal. For example in row Y=7 above N=22 to N=27.
Triangular Numbers
N=1,3,6,10,etc along the bottom Y=1 row is the triangular numbers N=k*(k1)/2. Such an N is at i=k,j=k and has gcd(i,j)=k which divides out to Y=1.
N=k*(k1)/2 i=k,j=k
Y = j/gcd
= 1 on the bottom row
X = (i + j*(gcd1)) / gcd
= (k + k*(k1)) / k
= k1 successive points on that bottom row
N=1,2,4,7,11,etc in the column at X=1 immediately follows each of those bottom row triangulars, ie. N+1.
N in X=1 column = Y*(Y1)/2 + 1
Primes
If N is prime then it's above the sloping line X=2*Y. If N is composite then it might be above or below, but the primes are always above. Here's the table with dots "..." marking the X=2*Y line.
primes and composites above

6  16 20 68
 .... X=2*Y
5  11 12 13 14 47 49 51 53 ....
 ....
4  7 9 30 34 .... 69
 ....
3  4 5 17 19 .... 39 42 70 only
 .... composite
2  2 8 .... 18 32 50 below
 ....
1  1 ..3. 6 10 15 21 28 36 45 55
 ....
Y=0  ....

X=0 1 2 3 4 5 6 7 8 9 10
Values below X=2*Y such as 39 and 42 are always composite. Values above such as 19 and 30 are either prime or composite. Only X=2,Y=1 is exactly on the line, which is prime N=3 as it happens. The rest of the line X=2*k,Y=k is not visited since common factor k would mean X/Y is not a rational in least terms.
This pattern of primes and composites occurs because N is a multiple of gcd(i,j) when that gcd is odd, or a multiple of gcd/2 when that gcd is even.
N = i + j*(j1)/2
gcd = gcd(i,j)
N = gcd * (i/gcd + j/gcd * (j1)/2) when gcd odd
gcd/2 * (2i/gcd + j/gcd * (j1)) when gcd even
If gcd odd then either j/gcd or j1 is even, to take the "/2" divisor. If gcd even then only gcd/2 can come out as a factor since taking out the full gcd might leave both j/gcd and j1 odd and so the "/2" not an integer. That happens for example to N=70
N = 70
i = 4, j = 12 for 4 + 12*11/2 = 70 = N
gcd(i,j) = 4
but N is not a multiple of 4, only of 4/2=2
Of course knowing gcd or gcd/2 is a factor of N is only useful when that factor is 2 or more, so
odd gcd >= 2 means gcd >= 3
even gcd with gcd/2 >= 2 means gcd >= 4
so N composite when gcd(i,j) >= 3
If gcd<3 then the "factor" coming out is only 1 and says nothing about whether N is prime or composite. There are both prime and composite N with gcd<3, as can be seen among the values above the X=2*Y line in the table above.
Rows Reverse
Option pairs_order => "rows_reverse"
reverses the order of points within the rows of i,j pairs,
j=4  10 9 8 7
j=3  6 5 4
j=2  3 2
j=1  1
+
i=1 2 3 4
The X,Y numbering becomes
pairs_order => "rows_reverse"
11  66 65 64 63 62 61 60 59 58 57
10  55 53 49 47 209
9  45 44 42 41 39 38 170 168
8  36 34 32 30 135 131
7  28 27 26 25 24 23 104 102 100 98
6  21 17 77 69
5  15 14 13 12 54 52 50 48 118
4  10 8 35 31 76 70
3  6 5 20 18 43 40 75 71
2  3 9 19 33 51 73
1  1 2 4 7 11 16 22 29 37 46 56
Y=0 

X=0 1 2 3 4 5 6 7 8 9 10 11
The triangular numbers, per "Triangular Numbers" above, are now in the X=1 column, ie. at the left rather than at the Y=1 bottom row. That bottom row is now the next after each triangular, ie. T(X)+1.
Diagonals
Option pairs_order => "diagonals_down"
takes the i,j pairs by diagonals down from the Y axis. pairs_order => "diagonals_up"
likewise but upwards from the X=Y centre up to the Y axis. (These numberings are in the style of Math::PlanePath::DiagonalsOctant.)
diagonals_down diagonals_up
j=7  13 j=7  16
j=6  10 14 j=6  12 15
j=5  7 11 15 j=5  9 11 14
j=4  5 8 12 16 j=4  6 8 10 13
j=3  3 6 9 j=3  4 5 7
j=2  2 4 j=2  2 3
j=1  1 j=1  1
+ +
i=1 2 3 4 i=1 2 3 4
The resulting path becomes
pairs_order => "diagonals_down"
9  21 27 40 47 63 72
8  17 28 41 56 74
7  13 18 23 29 35 42 58 76
6  10 30 44
5  7 11 15 20 32 46 62 80
4  5 12 22 48 52
3  3 6 14 24 33 55
2  2 8 19 34 54
1  1 4 9 16 25 36 49 64 81
Y=0 

X=0 1 2 3 4 5 6 7 8 9
The Y=1 bottom row is the perfect squares which are at i=j in the DiagonalsOctant
and have gcd(i,j)=i thus becoming X=i,Y=1.
pairs_order => "diagonals_up"
9  25 29 39 45 58 65
8  20 28 38 50 80
7  16 19 23 27 32 37 63 78
6  12 26 48
5  9 11 14 17 35 46 59 74
4  6 10 24 44 54
3  4 5 15 22 34 51
2  2 8 18 33 52
1  1 3 7 13 21 31 43 57 73
Y=0 

X=0 1 2 3 4 5 6 7 8 9
N=1,2,4,6,9 etc in the X=1 column is the perfect squares k*k and the pronics k*(k+1) interleaved, also called the quartersquares. N=2,5,10,17,etc on Y=X+1 above the leading diagonal are the squares+1, and N=3,8,15,24,etc below on Y=X1 below the diagonal are the squares1.
The GCD division moves points downwards and shears them across horizontally. The effect on diagonal lines of i,j points is as follows
 1
 1 gcd=1 slope=1
 1
 1
 1
 1
 1
 1
 1
 . gcd=2 slope=0
 . 2
 . 2 3 gcd=3 slope=1
 . 2 3 gcd=4 slope=2
 . 2 3 4
 . 3 4 5 gcd=5 slope=3
 . 4 5
 . 4 5
 . 5
+
The line of "1"s is the diagonal with gcd=1 and thus X,Y=i,j unchanged.
The line of "2"s is when gcd=2 so X=(i+j)/2,Y=j/2. Since i+j=d is constant within the diagonal this makes X=d fixed, ie. vertical.
Then gcd=3 becomes X=(i+2j)/3 which slopes across by +1 for each i, or gcd=4 has X=(i+3j)/4 slope +2, etc.
Of course only some of the points in an i,j diagonal have a given gcd, but those which do are transformed this way. The effect is that for N up to a given diagonal row all the "*" points in the following are traversed, plus extras in wedge shaped arms out to the side.
 *
 * * up to a given diagonal points "*"
 * * * all visited, plus some wedges out
 * * * * to the right
 * * * * *
 * * * * * /
 * * * * * / 
 * * * * * 
 * * * * *
+
In terms of the rationals X/Y the effect is that up to N=d^2 with diagonal d=2j the fractions enumerated are
N=d^2
enumerates all num/den where num <= d and num+den <= 2*d
FUNCTIONS
See "FUNCTIONS" in Math::PlanePath for behaviour common to all path classes.
$path = Math::PlanePath::GcdRationals>new ()
$path = Math::PlanePath::GcdRationals>new (pairs_order => $str)

Create and return a new path object. The
pairs_order
option can be"rows" (default) "rows_reverse" "diagonals_down" "diagonals_up"
FORMULAS
X,Y to N  Rows
The defining formula above for X,Y can be inverted to give i,j and N. This calculation doesn't notice if X,Y have a common factor, so a coprime(X,Y) test must be made separately if necessary (for xy_to_n()
it is).
X/Y = g1 + (i/g)/(j/g)
The g1 integer part is recovered by a division X divide Y,
X = quot*Y + rem division by Y rounded towards 0
where 0 <= rem < Y
unless Y=1 in which case use quot=X1, rem=1
g1 = quot
g = quot+1
The Y=1 special case can instead be left as the usual kind of division quot=X,rem=0, so 0<=rem<Y. This will give i=0 which is outside the intended 1<=i<=j range, but j is 1 bigger and the combination still gives the correct N. It's as if the i=g,j=g point at the end of a row is moved to i=0,j=g+1 just before the start of the next row. If only N is of interest not the i,j then it can be left rem=0.
Equating the denominators in the X/Y formula above gives j by
Y = j/g the definition above
j = g*Y
= (quot+1)*Y
j = X+Yrem per the division X=quot*Y+rem
And equating the numerators gives i by
X = (g1)*Y + i/g the definition above
i = X*g  (g1)*Y*g
= X*g  quot*Y*g
= X*g  (Xrem)*g per the division X=quot*Y+rem
i = rem*g
i = rem*(quot+1)
Then N from i,j by the definition above
N = i + j*(j1)/2
For example X=11,Y=4 divides X/Y as 11=4*2+3 for quot=2,rem=3 so i=3*(2+1)=9 j=11+43=12 and so N=9+12*11/2=75 (as shown in the first table above).
It's possible to use only the quotient p, not the remainder rem, by taking j=(quot+1)*Y instead of j=X+Yrem, but usually a division operation gives the remainder at no extra cost, or a cost small enough that it's worth swapping a multiply for an add or two.
The gcd g can be recovered by rounding up in the division, instead of rounding down and then incrementing with g=quot+1.
g = ceil(X/Y)
= cquot for division X=cquot*Y  crem
But division in most programming languages is towards 0 or towards infinity, not upwards towards +infinity.
X,Y to N  Rows Reverse
For pairs_order="rows_reverse", the horizontal i is reversed to ji+1. This can be worked into the triangular part of the N formula as
Nrrev = (ji+1) + j*(j1)/2 for 1<=i<=j
= j*(j+1)/2  i + 1
The Y=1 case described above cannot be left to go through with rem=0 giving i=0 and j+1 since the reversal ji+1 is then not correct. Either use rem=1 as described, or if not then compensate at the end,
if r=0 then j = 2 adjust
Nrrev = j*(j+1)/2  i + 1 same Nrrev as above
For example X=5,Y=1 is quot=5,rem=0 gives i=0*(5+1)=0 j=5+10=6. Without adjustment it would be Nrrev=6*7/20+1=22 which is wrong. But adjusting j=2 so that j=62=4 gives the desired Nrrev=4*5/20+1=11 (per the table in "Rows Reverse" above).
OEIS
This enumeration of rationals is in Sloane's Online Encyclopedia of Integer Sequences in the following forms
http://oeis.org/A054531 (etc)
pairs_order="rows" (the default)
A226314 X coordinate
A054531 Y coordinate, being N/GCD(i,j)
A000124 N in X=1 column, triangular+1
A050873 ceil(X/Y), gcd by rows
A050873A023532 floor(X/Y)
gcd by rows and subtract 1 unless i=j
pairs_order="diagonals_down"
A033638 N in X=1 column, quartersquares+1 and pronic+1
A000290 N in Y=1 row, perfect squares
pairs_order="diagonals_up"
A002620 N in X=1 column, squares and pronics
A002061 N in Y=1 row, central polygonals (extra initial 1)
A002522 N at Y=X+1 above leading diagonal, squares+1
SEE ALSO
Math::PlanePath, Math::PlanePath::DiagonalRationals, Math::PlanePath::RationalsTree, Math::PlanePath::CoprimeColumns, Math::PlanePath::DiagonalsOctant
HOME PAGE
http://user42.tuxfamily.org/mathplanepath/index.html
LICENSE
Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Kevin Ryde
This file is part of MathPlanePath.
MathPlanePath is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version.
MathPlanePath is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with MathPlanePath. If not, see <http://www.gnu.org/licenses/>.