NAME
Math::PlanePath::GcdRationals -- rationals by prime factorization
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, ie. 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 in 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. Only X=2,Y=1 is exactly on the line (which is prime N=3 as it happens) because X=2*k,Y=k is not an X/Y rational in least terms (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.
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 10This occurs because N is a multiple of gcd(i,j) when the 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 formula only shows N is composite when odd gcd>=3 or even gcd/2>=2, so gcd>=3, since otherwise the factor coming out is only 1. When gcd<3 there are in fact both prime and composite N, the values above the X=2*Y line in the sample above.
FUNCTIONS
See "FUNCTIONS" in Math::PlanePath for the 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 = g * ( X + Y*((Y-2)*g + 1)/2 )
with g = floor(X/Y) + 1Either Y or (Y-2)*g+1 is even to take the /2. If Y is odd then g must be odd which makes (Y-2)*g odd and (Y-2)*g+1 even.
Y*g is the next multiple of Y which is strictly greater than X. It can be formed from the floor(X/Y) division
X = Y*q + r division
g = q+1
Y*g = Y*(q+1) = X+Y-r
= X+Y-rUsing X+Y-r instead of Y*g in the N formula might swap a multiply for an add or subtract if you get the remainder for free with the X/Y division.
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. higher g. And conversely the minimum is either the bottom right, or the start of the next higher wedge, ie. lower g. (That's right is it?)
SEE 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 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/>.