1Math::PlanePath::GcdRatUisoenralCso(n3t)ributed Perl DocMuamtehn:t:aPtliaonnePath::GcdRationals(3)
2
3
4

NAME

6       Math::PlanePath::GcdRationals -- rationals by triangular GCD
7

SYNOPSIS

9        use Math::PlanePath::GcdRationals;
10        my $path = Math::PlanePath::GcdRationals->new;
11        my ($x, $y) = $path->n_to_xy (123);
12

DESCRIPTION

14       This path enumerates X/Y rationals using a method by Lance Fortnow
15       taking a greatest common divisor out of a triangular position.
16
17           <http://blog.computationalcomplexity.org/2004/03/counting-rationals-quickly.html>
18
19       The attraction of this approach is that it's both efficient to
20       calculate and visits blocks of X/Y rationals using a modest range of N
21       values, roughly a square N=2*max(num,den)^2 in the default rows style.
22
23           13  |      79  80  81  82  83  84  85  86  87  88  89  90
24           12  |      67              71      73              77     278
25           11  |      56  57  58  59  60  61  62  63  64  65     233 235
26           10  |      46      48              52      54     192     196
27            9  |      37  38      40  41      43  44     155 157     161
28            8  |      29      31      33      35     122     126     130
29            7  |      22  23  24  25  26  27      93  95  97  99 101 103
30            6  |      16              20      68              76     156
31            5  |      11  12  13  14      47  49  51  53     108 111 114
32            4  |       7       9      30      34      69      75     124
33            3  |       4   5      17  19      39  42      70  74     110
34            2  |       2       8      18      32      50      72      98
35            1  |       1   3   6  10  15  21  28  36  45  55  66  78  91
36           Y=0 |
37                --------------------------------------------------------
38                 X=0   1   2   3   4   5   6   7   8   9  10  11  12  13
39
40       The mapping from N to rational is
41
42           N = i + j*(j-1)/2     for upper triangle 1 <= i <= j
43           gcd = GCD(i,j)
44           rational = i/j + gcd-1
45
46       which means X=numerator Y=denominator are
47
48           X = (i + j*(gcd-1))/gcd  = j + (i-j)/gcd
49           Y = j/gcd
50
51       The i,j position is a numbering of points above the X=Y diagonal by
52       rows in the style of Math::PlanePath::PyramidRows with step=1, but
53       starting from i=1,j=1.
54
55           j=4  |  7  8  9 10
56           j=3  |  4  5  6
57           j=2  |  2  3
58           j=1  |  1
59                +-------------
60                 i=1  2  3  4
61
62       If GCD(i,j)=1 then X/Y is simply X=i,Y=j unchanged.  This means
63       fractions X/Y < 1 are numbered by rows with increasing numerator, but
64       skipping positions where i,j have a common factor.
65
66       The skipped positions where i,j have a common factor become rationals
67       X/Y>1, ie. below the X=Y diagonal.  The integer part is GCD(i,j)-1 so
68       rational = gcd-1 + i/j.  For example
69
70           N=51 is at i=6,j=10 by rows
71           common factor gcd(6,10)=2
72           so rational R = 2-1 + 6/10 = 1+3/5 = 8/5
73           ie. X=8,Y=5
74
75       If j is prime then gcd(i,j)=1 and so X=i,Y=j.  This means that in rows
76       with prime Y are numbered by consecutive N across to the X=Y diagonal.
77       For example in row Y=7 above N=22 to N=27.
78
79   Triangular Numbers
80       N=1,3,6,10,etc along the bottom Y=1 row is the triangular numbers
81       N=k*(k-1)/2.  Such an N is at i=k,j=k and has gcd(i,j)=k which divides
82       out to Y=1.
83
84           N=k*(k-1)/2  i=k,j=k
85
86           Y = j/gcd
87             = 1       on the bottom row
88
89           X = (i + j*(gcd-1)) / gcd
90             = (k + k*(k-1)) / k
91             = k-1     successive points on that bottom row
92
93       N=1,2,4,7,11,etc in the column at X=1 immediately follows each of those
94       bottom row triangulars, ie. N+1.
95
96           N in X=1 column = Y*(Y-1)/2 + 1
97
98   Primes
99       If N is prime then it's above the sloping line X=2*Y.  If N is
100       composite then it might be above or below, but the primes are always
101       above.  Here's the table with dots "..." marking the X=2*Y line.
102
103                  primes and composites above
104               |
105            6  |      16              20      68
106               |                                             .... X=2*Y
107            5  |      11  12  13  14      47  49  51  53 ....
108               |                                     ....
109            4  |       7       9      30      34 .... 69
110               |                             ....
111            3  |       4   5      17  19 .... 39  42      70   only
112               |                     ....                      composite
113            2  |       2       8 .... 18      32      50       below
114               |             ....
115            1  |       1 ..3.  6  10  15  21  28  36  45  55
116               |     ....
117           Y=0 | ....
118                ---------------------------------------------
119                 X=0   1   2   3   4   5   6   7   8   9  10
120
121       Values below X=2*Y such as 39 and 42 are always composite.  Values
122       above such as 19 and 30 are either prime or composite.  Only X=2,Y=1 is
123       exactly on the line, which is prime N=3 as it happens.  The rest of the
124       line X=2*k,Y=k is not visited since common factor k would mean X/Y is
125       not a rational in least terms.
126
127       This pattern of primes and composites occurs because N is a multiple of
128       gcd(i,j) when that gcd is odd, or a multiple of gcd/2 when that gcd is
129       even.
130
131           N = i + j*(j-1)/2
132           gcd = gcd(i,j)
133
134           N = gcd   * (i/gcd + j/gcd * (j-1)/2)  when gcd odd
135               gcd/2 * (2i/gcd + j/gcd * (j-1))   when gcd even
136
137       If gcd odd then either j/gcd or j-1 is even, to take the "/2" divisor.
138       If gcd even then only gcd/2 can come out as a factor since taking out
139       the full gcd might leave both j/gcd and j-1 odd and so the "/2" not an
140       integer.  That happens for example to N=70
141
142           N = 70
143           i = 4, j = 12     for 4 + 12*11/2 = 70 = N
144           gcd(i,j) = 4
145           but N is not a multiple of 4, only of 4/2=2
146
147       Of course knowing gcd or gcd/2 is a factor of N is only useful when
148       that factor is 2 or more, so
149
150           odd gcd >= 2                means gcd >= 3
151           even gcd with gcd/2 >= 2    means gcd >= 4
152
153           so N composite when gcd(i,j) >= 3
154
155       If gcd<3 then the "factor" coming out is only 1 and says nothing about
156       whether N is prime or composite.  There are both prime and composite N
157       with gcd<3, as can be seen among the values above the X=2*Y line in the
158       table above.
159
160   Rows Reverse
161       Option "pairs_order => "rows_reverse"" reverses the order of points
162       within the rows of i,j pairs,
163
164           j=4 | 10  9  8  7
165           j=3 |  6  5  4
166           j=2 |  3  2
167           j=1 |  1
168               +------------
169                i=1  2  3  4
170
171       The X,Y numbering becomes
172
173           pairs_order => "rows_reverse"
174
175           11  |      66  65  64  63  62  61  60  59  58  57
176           10  |      55      53              49      47     209
177            9  |      45  44      42  41      39  38     170 168
178            8  |      36      34      32      30     135     131
179            7  |      28  27  26  25  24  23     104 102 100  98
180            6  |      21              17      77              69
181            5  |      15  14  13  12      54  52  50  48     118
182            4  |      10       8      35      31      76      70
183            3  |       6   5      20  18      43  40      75  71
184            2  |       3       9      19      33      51      73
185            1  |       1   2   4   7  11  16  22  29  37  46  56
186           Y=0 |
187                ------------------------------------------------
188                 X=0   1   2   3   4   5   6   7   8   9  10  11
189
190       The triangular numbers, per "Triangular Numbers" above, are now in the
191       X=1 column, ie. at the left rather than at the Y=1 bottom row.  That
192       bottom row is now the next after each triangular, ie. T(X)+1.
193
194   Diagonals
195       Option "pairs_order => "diagonals_down"" takes the i,j pairs by
196       diagonals down from the Y axis.  "pairs_order => "diagonals_up""
197       likewise but upwards from the X=Y centre up to the Y axis.  (These
198       numberings are in the style of Math::PlanePath::DiagonalsOctant.)
199
200           diagonals_down            diagonals_up
201
202           j=7 | 13                   j=7 | 16
203           j=6 | 10 14                j=6 | 12 15
204           j=5 |  7 11 15             j=5 |  9 11 14
205           j=4 |  5  8 12 16          j=4 |  6  8 10 13
206           j=3 |  3  6  9             j=3 |  4  5  7
207           j=2 |  2  4                j=2 |  2  3
208           j=1 |  1                   j=1 |  1
209               +------------              +------------
210                i=1  2  3  4               i=1  2  3  4
211
212       The resulting path becomes
213
214           pairs_order => "diagonals_down"
215
216            9  |     21 27    40 47    63 72
217            8  |     17    28    41    56    74
218            7  |     13 18 23 29 35 42    58 76
219            6  |     10          30    44
220            5  |      7 11 15 20    32 46 62 80
221            4  |      5    12    22    48    52
222            3  |      3  6    14 24    33 55
223            2  |      2     8    19    34    54
224            1  |      1  4  9 16 25 36 49 64 81
225           Y=0 |
226                --------------------------------
227                 X=0  1  2  3  4  5  6  7  8  9
228
229       The Y=1 bottom row is the perfect squares which are at i=j in the
230       "DiagonalsOctant" and have gcd(i,j)=i thus becoming X=i,Y=1.
231
232           pairs_order => "diagonals_up"
233
234            9  |     25 29    39 45    58 65
235            8  |     20    28    38    50    80
236            7  |     16 19 23 27 32 37    63 78
237            6  |     12          26    48
238            5  |      9 11 14 17    35 46 59 74
239            4  |      6    10    24    44    54
240            3  |      4  5    15 22    34 51
241            2  |      2     8    18    33    52
242            1  |      1  3  7 13 21 31 43 57 73
243           Y=0 |
244                --------------------------------
245                 X=0  1  2  3  4  5  6  7  8  9
246
247       N=1,2,4,6,9 etc in the X=1 column is the perfect squares k*k and the
248       pronics k*(k+1) interleaved, also called the quarter-squares.
249       N=2,5,10,17,etc on Y=X+1 above the leading diagonal are the squares+1,
250       and N=3,8,15,24,etc below on Y=X-1 below the diagonal are the
251       squares-1.
252
253       The GCD division moves points downwards and shears them across
254       horizontally.  The effect on diagonal lines of i,j points is as follows
255
256             | 1
257             |   1     gcd=1 slope=-1
258             |     1
259             |       1
260             |         1
261             |           1
262             |             1
263             |               1
264             |                 1
265             |                 .    gcd=2 slope=0
266             |               .   2
267             |             .     2     3  gcd=3 slope=1
268             |           .       2   3           gcd=4 slope=2
269             |         .         2 3         4
270             |       .           3       4       5     gcd=5 slope=3
271             |     .                 4      5
272             |   .               4     5
273             | .                 5
274             +-------------------------------
275
276       The line of "1"s is the diagonal with gcd=1 and thus X,Y=i,j unchanged.
277
278       The line of "2"s is when gcd=2 so X=(i+j)/2,Y=j/2.  Since i+j=d is
279       constant within the diagonal this makes X=d fixed, ie. vertical.
280
281       Then gcd=3 becomes X=(i+2j)/3 which slopes across by +1 for each i, or
282       gcd=4 has X=(i+3j)/4 slope +2, etc.
283
284       Of course only some of the points in an i,j diagonal have a given gcd,
285       but those which do are transformed this way.  The effect is that for N
286       up to a given diagonal row all the "*" points in the following are
287       traversed, plus extras in wedge shaped arms out to the side.
288
289            | *
290            | * *                 up to a given diagonal points "*"
291            | * * *               all visited, plus some wedges out
292            | * * * *             to the right
293            | * * * * *
294            | * * * * *   /
295            | * * * * * /  --
296            | * * * * *  --
297            | * * * * *--
298            +--------------
299
300       In terms of the rationals X/Y the effect is that up to N=d^2 with
301       diagonal d=2j the fractions enumerated are
302
303           N=d^2
304           enumerates all num/den where num <= d and num+den <= 2*d
305

FUNCTIONS

307       See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
308       classes.
309
310       "$path = Math::PlanePath::GcdRationals->new ()"
311       "$path = Math::PlanePath::GcdRationals->new (pairs_order => $str)"
312           Create and return a new path object.  The "pairs_order" option can
313           be
314
315               "rows"               (default)
316               "rows_reverse"
317               "diagonals_down"
318               "diagonals_up"
319

FORMULAS

321   X,Y to N -- Rows
322       The defining formula above for X,Y can be inverted to give i,j and N.
323       This calculation doesn't notice if X,Y have a common factor, so a
324       coprime(X,Y) test must be made separately if necessary (for "xy_to_n()"
325       it is).
326
327           X/Y = g-1 + (i/g)/(j/g)
328
329       The g-1 integer part is recovered by a division X divide Y,
330
331           X = quot*Y + rem   division by Y rounded towards 0
332             where 0 <= rem < Y
333             unless Y=1 in which case use quot=X-1, rem=1
334           g-1 = quot
335           g = quot+1
336
337       The Y=1 special case can instead be left as the usual kind of division
338       quot=X,rem=0, so 0<=rem<Y.  This will give i=0 which is outside the
339       intended 1<=i<=j range, but j is 1 bigger and the combination still
340       gives the correct N.  It's as if the i=g,j=g point at the end of a row
341       is moved to i=0,j=g+1 just before the start of the next row.  If only N
342       is of interest not the i,j then it can be left rem=0.
343
344       Equating the denominators in the X/Y formula above gives j by
345
346           Y = j/g          the definition above
347
348           j = g*Y
349             = (quot+1)*Y
350           j = X+Y-rem      per the division X=quot*Y+rem
351
352       And equating the numerators gives i by
353
354           X = (g-1)*Y + i/g     the definition above
355
356           i = X*g - (g-1)*Y*g
357             = X*g - quot*Y*g
358             = X*g - (X-rem)*g     per the division X=quot*Y+rem
359           i = rem*g
360           i = rem*(quot+1)
361
362       Then N from i,j by the definition above
363
364           N = i + j*(j-1)/2
365
366       For example X=11,Y=4 divides X/Y as 11=4*2+3 for quot=2,rem=3 so
367       i=3*(2+1)=9 j=11+4-3=12 and so N=9+12*11/2=75 (as shown in the first
368       table above).
369
370       It's possible to use only the quotient p, not the remainder rem, by
371       taking j=(quot+1)*Y instead of j=X+Y-rem, but usually a division
372       operation gives the remainder at no extra cost, or a cost small enough
373       that it's worth swapping a multiply for an add or two.
374
375       The gcd g can be recovered by rounding up in the division, instead of
376       rounding down and then incrementing with g=quot+1.
377
378          g = ceil(X/Y)
379            = cquot for division X=cquot*Y - crem
380
381       But division in most programming languages is towards 0 or towards
382       -infinity, not upwards towards +infinity.
383
384   X,Y to N -- Rows Reverse
385       For pairs_order="rows_reverse", the horizontal i is reversed to j-i+1.
386       This can be worked into the triangular part of the N formula as
387
388           Nrrev = (j-i+1) + j*(j-1)/2        for 1<=i<=j
389                 = j*(j+1)/2 - i + 1
390
391       The Y=1 case described above cannot be left to go through with rem=0
392       giving i=0 and j+1 since the reversal j-i+1 is then not correct.
393       Either use rem=1 as described, or if not then compensate at the end,
394
395           if r=0 then j -= 2            adjust
396           Nrrev = j*(j+1)/2 - i + 1     same Nrrev as above
397
398       For example X=5,Y=1 is quot=5,rem=0 gives i=0*(5+1)=0 j=5+1-0=6.
399       Without adjustment it would be Nrrev=6*7/2-0+1=22 which is wrong.  But
400       adjusting j-=2 so that j=6-2=4 gives the desired Nrrev=4*5/2-0+1=11
401       (per the table in "Rows Reverse" above).
402

OEIS

404       This enumeration of rationals is in Sloane's Online Encyclopedia of
405       Integer Sequences in the following forms
406
407           <http://oeis.org/A054531> (etc)
408
409           pairs_order="rows" (the default)
410             A226314   X coordinate
411             A054531   Y coordinate, being N/GCD(i,j)
412             A000124   N in X=1 column, triangular+1
413             A050873   ceil(X/Y), gcd by rows
414             A050873-A023532  floor(X/Y)
415                       gcd by rows and subtract 1 unless i=j
416
417           pairs_order="diagonals_down"
418             A033638   N in X=1 column, quartersquares+1 and pronic+1
419             A000290   N in Y=1 row, perfect squares
420
421           pairs_order="diagonals_up"
422             A002620   N in X=1 column, squares and pronics
423             A002061   N in Y=1 row, central polygonals (extra initial 1)
424             A002522   N at Y=X+1 above leading diagonal, squares+1
425

SEE ALSO

427       Math::PlanePath, Math::PlanePath::DiagonalRationals,
428       Math::PlanePath::RationalsTree, Math::PlanePath::CoprimeColumns,
429       Math::PlanePath::DiagonalsOctant
430

HOME PAGE

432       <http://user42.tuxfamily.org/math-planepath/index.html>
433

LICENSE

435       Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Kevin Ryde
436
437       This file is part of Math-PlanePath.
438
439       Math-PlanePath is free software; you can redistribute it and/or modify
440       it under the terms of the GNU General Public License as published by
441       the Free Software Foundation; either version 3, or (at your option) any
442       later version.
443
444       Math-PlanePath is distributed in the hope that it will be useful, but
445       WITHOUT ANY WARRANTY; without even the implied warranty of
446       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
447       General Public License for more details.
448
449       You should have received a copy of the GNU General Public License along
450       with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.
451
452
453
454perl v5.30.0                      2019-08-17  Math::PlanePath::GcdRationals(3)
Impressum