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 rationals X/Y using a method by Lance Fortnow taking a GCD out of a triangular position.
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 X/Y is
N = i + j*(j-1)/2 upper triangle 1 <= i <= j
gcd = GCD(i,j)
rational = i/j + gcd-1
X = (i + j*(gcd-1)) / gcd
Y = j/gcdThe i,j position is a numbering of points above the X=Y diagonal by rows,
j=4 7 8 9 10
j=3 4 5 6
j=2 2 3
j=1 1
i=1 2 3 4When gcd=1, X/Y is simply X=i,Y=j. So the fractions X/Y < 1 above the X=Y diagonal are numbered by rows, so increasing numerator, skipping positions where X,Y have a common factor.
The skipped positions where i,j have a common factor become rationals R>1 below the X=Y diagonal. gcd(i,j)-1 is made the integer part as R = i/j+(gcd-1). For example N=51 is at i=6,j=10 by rows, but they have common factor gcd=2 so 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 has i=k,j=k and thus gcd(i,j)=k,
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 rowPrimes
All prime N values are above the sloping line X=2*Y. There's composites both above and below, but the primes are all above. Here's the table with "..." 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 composites
| .... below
2 | 2 8 .... 18 32 50
| ....
1 | 1 ..3. 6 10 15 21 28 36 45 55
| ....
Y=0 | ....
---------------------------------------------
X=0 1 2 3 4 5 6 7 8 9 10Only X=2,Y=1 is exactly on the line, which is prime N=3 as it happens. Other X=2*k,Y=k is not an X/Y rational in least terms because it has common factor k. Values below X=2*Y such as 39 and 42 are all composites, values above like 19 and 30 are either prime or composite.
This 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 evenFor gcd odd either j/gcd or j-1 is even to take the "/2". When gcd is even 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 with N=70
N = 70
i = 4, j = 12 for 4 + 12*11/2 = 70
gcd = 4
but N is not a multiple of 4, only of 4/2=2Of course the calculation only shows N is composite for odd gcd>=3 or even gcd/2>=2, together meaning gcd>=3. If gcd<3 then the "factor" coming out is only 1 and so says nothing about whether N is prime or composite. There are in fact both prime and composite N for gcd<3, as per the values above the X=2*Y line in the sample above.
FUNCTIONS
See "FUNCTIONS" in Math::PlanePath for behaviour common to all path classes.
FORMULAS
X,Y to N
The defining formula above can be reversed
X/Y = i/j + g-1
g-1 = floor(X/Y)
Y = j/g
X = ((g-1)*j + i)/g
j = Y*g
i = X*g - (g-1)*Y*g
N = i + j*(j-1)/2So
N = 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/2The /2 division is exact. If Y and g are both odd then (Y-2)*g+1 is odd*odd+1 so even.
Y*g in the formulas is the next multiple of Y strictly greater than X. It can be formed from the g-1=floor(X/Y) division,
X = Y*q + r division
g = q+1
Y*g = Y*(q+1) = X+Y-r
= X+Y-rIf you get the remainder for free in the X/Y division then using X+Y-r instead of Y*g in the N formula might reduce a multiply to an add or subtract.
Rectangle N Range
An over-estimate of the N range can be calculated just from the X,Y to N formula above ignoring which X,Y points are coprime and thus actually should have N values.
Within a row N increases with increasing X, so for a rectangle the minimum is in the left column and the maximum in the right column.
Within a column N values increase until reaching the end of a "g" wedge, then drop down a bit. So the maximum is either the top-right corner or the top of the next lower wedge, ie. smaller Y but bigger g. And conversely the minimum is either the bottom right, or the bottom of the next higher wedge, ie. smaller g but bigger Y. (That's right is it?)
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
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/>.