Math::PlanePath::GcdRationals -- rationals by triangular GCD
use Math::PlanePath::GcdRationals; my $path = Math::PlanePath::GcdRationals->new; my ($x, $y) = $path->n_to_xy (123);
This path enumerates X/Y rationals using a method by Lance Fortnow taking a greatest common divisor out of a triangular position.
http://blog.computationalcomplexity.org/2004/03/counting-rationals-quickly.html
The attraction of this approach is that it's both efficient to calculate and it traverses blocks of X/Y rationals using a modest range of N values, roughly a square N=2*max(num,den)^2 in the default i,j by rows.
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*(j-1)/2 for upper triangle 1 <= i <= j gcd = GCD(i,j) rational = i/j + gcd-1
which means
X = (i + j*(gcd-1)) / 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).
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 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 = gcd-1 + i/j. For example
N=51 is at i=6,j=10 by rows common factor gcd(6,10)=2 so rational R = 2-1 + 6/10 = 1+3/5 = 8/5 ie. X=8,Y=5
N=1,3,6,10,etc along the bottom Y=1 row is the triangular numbers N=k*(k-1)/2. Such an N is at i=k,j=k and has 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 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 on X=1 column = Y*(Y-1)/2 + 1
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. Others on the line X=2*k,Y=k are not visited since common factor k means 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*(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 even
If gcd odd then either j/gcd or j-1 is even, to take the "/2" divisor. 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=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.
Option pairs_order => "rows_reverse" reverses the order of points within the rows of i,j pairs,
pairs_order => "rows_reverse"
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" 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.
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.)
pairs_order => "diagonals_down"
pairs_order => "diagonals_up"
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 quarter-squares. 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=X-1 below the diagonal are the squares-1.
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
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
pairs_order
"rows" (default) "rows_reverse" "diagonals_down" "diagonals_up"
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 (which for xy_to_n() it is).
xy_to_n()
X/Y = g-1 + (i/g)/(j/g)
The g-1 integer part is recovered by a division X divide Y,
X = Y*p + r division by Y to p=quotient r=remainder where 0<=r<Y unless Y=1 in which case use r=1, p=X-1 g-1 = p g = p+1
The Y=1 special case can instead be left as the usual kind of division r=0,p=X, so 0<=r<Y. This ends up i=0 below 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 as such then it can be left r=0.
Equating the denominators in the X/Y formula above gives j by
Y = j/g the definition above j = g*Y = (p+1)*Y j = X+Y-r per the division X=p*Y+r
And equating the numerators gives i by
X = (g-1)*Y + i/g the definition above i = X*g - (g-1)*Y*g = X*g - p*Y*g = X*g - (X-r)*g per the division X=p*Y+r i = r*g i = r*(p+1)
Then N from i,j by the definition above
N = i + j*(j-1)/2
For example X=11,Y=4 divides X/Y as 11=4*2+3 for p=2,r=3 so i=3*(2+1)=9 j=11+4-3=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 r, by taking j=(p+1)*Y instead of j=X+Y-r, but usually a division operation gives the remainder at no extra cost, or a cost small enough that it's worth it to swap a multiply for an add or two.
For pairs_order="rows_reverse", the horizontal i is reversed to j-i+1. This can be worked into the triangular part of the N formula as
Nrrev = (j-i+1) + j*(j-1)/2 for 1<=i<=j = j*(j+1)/2 - i + 1
The Y=1 case described above cannot be left to go through with r=0 giving i=0 and j+1 since the reversal j-i+1 is then not correct. Either use r=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 r=0,p=5 gives i=0*(5+1)=0 j=5+1-0=6. Without adjustment it would be Nrrev=6*7/2-0+1=22 which is wrong. But adjusting j-=2 so that j=6-2=4 gives the desired Nrrev=4*5/2-0+1=11 (per the table in "Rows Reverse" above).
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) A054531 Y coordinate, being N/GCD(i,j) A000124 N in X=1 column, triangular+1 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
Math::PlanePath, Math::PlanePath::DiagonalRationals, Math::PlanePath::RationalsTree, Math::PlanePath::CoprimeColumns, Math::PlanePath::DiagonalsOctant
http://user42.tuxfamily.org/math-planepath/index.html
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/>.
To install Math::PlanePath, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Math::PlanePath
CPAN shell
perl -MCPAN -e shell install Math::PlanePath
For more information on module installation, please visit the detailed CPAN module installation guide.