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. It has the attraction of being both efficient to calculate (a GCD) and completing X/Y blocks with a much smaller N range than the tree based rationals.
http://blog.computationalcomplexity.org/2004/03/counting-rationals-quickly.html
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 13The mapping from N to rational is
N = i + j*(j-1)/2 upper triangle 1 <= i <= j
gcd = GCD(i,j)
rational = i/j + gcd-1which means
X = (i + j*(gcd-1)) / gcd
Y = j/gcdThe i,j position is a numbering of points above the X=Y diagonal by rows, in the style of PyramidRows step=1.
j=4 7 8 9 10
j=3 4 5 6
j=2 2 3
j=1 1
i=1 2 3 4If 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. GCD(i,j)-1 is the integer part as R = i/j+(gcd-1). For example N=51 is at i=6,j=10 by rows and that i,j has common factor gcd(6,10)=2 so becomes rational R = 6/10+(2-1) = 3/5+1 = 8/5, ie. X=8,Y=5.
Triangular Numbers
The bottom row Y=1 is the triangular numbers N=1,3,6,10,etc, k*(k-1)/2. Such an N is at i=k,j=k and thus gcd(i,j)=k which divides out to Y=1.
Y = j/gcd
= 1 on the bottom row
X = (i + j*(gcd-1)) / gcd
= (k + k*(k-1)) / k
= k-1 successive points on that bottom rowN=1,2,4,7,11,etc in the vertical at X=1 immediately following those triangulars on the bottom row, ie.
N on X=1 column = Y*(Y-1)/2 + 1Primes
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 "..." for 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 always
| .... 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 10Values 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. Other X=2*k,Y=k are not an X/Y rational in least terms because it has common factor k.
This pattern of primes and composites occurs because N is a multiple of gcd(i,j) when gcd odd, or a multiple of gcd/2 when gcd even.
N = i + j*(j-1)/2
gcd = gcd(i,j)
N = gcd * (i/gcd + j/gcd * (j-1)/2) when gcd odd
gcd/2 * (2i/gcd + j/gcd * (j-1)) when gcd evenIf gcd odd then either j/gcd or j-1 is even, taking the "/2". If gcd even then only gcd/2 can come out as a factor since the full gcd might leave both j/gcd and j-1 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=2Of course knowing gcd or gcd/2 is a factor is only useful when that factor is 2 or more, so only
odd gcd with gcd >= 2 means gcd >= 3
even gcd with gcd/2 >= 2 means gcd >= 4
so N composite when gcd(i,j) >= 3If 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 for 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 4The point 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 11The triangular numbers per "Triangular Numbers" above are now in the X=1 column, ie. at the left rather than the bottom. The Y=1 bottom row is 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. This is in the style of the DiagonalsOctant path.
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 4The 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
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 9For "diagonals_down" 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.
The gcd shears moves points downwards and shears them across horizontally.
| 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. a vertical.
Then gcd=3 becomes X=(i+2j)/3 which slopes across by +1 for each i, or gcd=4 X=(i+3j)/4 slope +2, etc.
Of course only some of the points in a 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 diagonal d=2j the fractions enumerated are
N=d^2
enumerates to num <= d and num+den <= 2*dFUNCTIONS
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_orderoption can be"rows" (default) "rows_reverse" "diagonals_down" "diagonals_up"
FORMULAS
X,Y to N
The defining formula above for X,Y can be reversed
X/Y = i/j + g-1
g-1 = floor(X/Y)
Y = j/g
X = ((g-1)*j + i)/gso
j = Y*g
i = X*g - (g-1)*Y*gSo N = i + j*(j-1)/2 = X*g - (g-1)*Y*g + Y*g*(Y*g-1)/2 = X*g + ((Y-2)*g + 1)*Y*g/2 = (((Y-2)*g + 1)*Y + 2X)*g/2
The /2 division is exact. If Y and g are both odd and so don't take that divisor then the term (Y-2)*g+1 is odd*odd+1 so even.
Y*g in the formulas is the first multiple of Y which is strictly greater than X. It can be formed from the g-1=floor(X/Y) division,
X = Y*quot + rem division
g = quot+1
Y*g = Y*(q+1) = X+Y-rem
= X+Y-remIf a division gives quotient and remainder for the same price then X+Y-rem instead of Y*g might reduce a multiply to instead an add or subtract.
OEIS
This enumeration of rationals is in Sloane's Online Encyclopedia of Integer Sequences in the following forms
http://oeis.org/A054531 (etc)
A054531 - Y coordinate, ie. denominatorsSEE ALSO
Math::PlanePath, Math::PlanePath::DiagonalRationals, Math::PlanePath::RationalsTree, Math::PlanePath::CoprimeColumns, Math::PlanePath::DiagonalsOctant
HOME PAGE
http://user42.tuxfamily.org/math-planepath/index.html
LICENSE
Copyright 2011, 2012 Kevin Ryde
This file is part of Math-PlanePath.
Math-PlanePath 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.
Math-PlanePath 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 Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.