1Math::PlanePath::GcdRatUisoenralCso(n3t)ributed Perl DocMuamtehn:t:aPtliaonnePath::GcdRationals(3)
2
6 Math::PlanePath::GcdRationals -- rationals by triangular GCD
7
9 use Math::PlanePath::GcdRationals;
10 my $path = Math::PlanePath::GcdRationals->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
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 <http://blog.computationalcomplexity.org/2004/03/counting-rationals-quickly.html>
18 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 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 The mapping from N to rational is
41 N = i + j*(j-1)/2 for upper triangle 1 <= i <= j
43 gcd = GCD(i,j)
44 rational = i/j + gcd-1
45 which means X=numerator Y=denominator are
47 X = (i + j*(gcd-1))/gcd = j + (i-j)/gcd
49 Y = j/gcd
50 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 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 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 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 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 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 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 N=k*(k-1)/2 i=k,j=k
85 Y = j/gcd
87 = 1 on the bottom row
88 X = (i + j*(gcd-1)) / gcd
90 = (k + k*(k-1)) / k
91 = k-1 successive points on that bottom row
92 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 N in X=1 column = Y*(Y-1)/2 + 1
97 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 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 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 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 N = i + j*(j-1)/2
132 gcd = gcd(i,j)
133 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 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 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 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 odd gcd >= 2 means gcd >= 3
151 even gcd with gcd/2 >= 2 means gcd >= 4
152 so N composite when gcd(i,j) >= 3
154 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 Rows Reverse
161 Option "pairs_order => "rows_reverse"" reverses the order of points
162 within the rows of i,j pairs,
163 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 The X,Y numbering becomes
172 pairs_order => "rows_reverse"
174 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 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 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 diagonals_down diagonals_up
201 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 The resulting path becomes
213 pairs_order => "diagonals_down"
215 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 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 pairs_order => "diagonals_up"
233 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 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 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 | 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 The line of "1"s is the diagonal with gcd=1 and thus X,Y=i,j unchanged.
277 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 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 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 | *
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 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 N=d^2
304 enumerates all num/den where num <= d and num+den <= 2*d
305
307 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
308 classes.
309 "$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 "rows" (default)
316 "rows_reverse"
317 "diagonals_down"
318 "diagonals_up"
319
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 X/Y = g-1 + (i/g)/(j/g)
328 The g-1 integer part is recovered by a division X divide Y,
330 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 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 Equating the denominators in the X/Y formula above gives j by
345 Y = j/g the definition above
347 j = g*Y
349 = (quot+1)*Y
350 j = X+Y-rem per the division X=quot*Y+rem
351 And equating the numerators gives i by
353 X = (g-1)*Y + i/g the definition above
355 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 Then N from i,j by the definition above
363 N = i + j*(j-1)/2
365 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 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 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 g = ceil(X/Y)
379 = cquot for division X=cquot*Y - crem
380 But division in most programming languages is towards 0 or towards
382 -infinity, not upwards towards +infinity.
383 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 Nrrev = (j-i+1) + j*(j-1)/2 for 1<=i<=j
389 = j*(j+1)/2 - i + 1
390 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 if r=0 then j -= 2 adjust
396 Nrrev = j*(j+1)/2 - i + 1 same Nrrev as above
397 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
404 This enumeration of rationals is in Sloane's Online Encyclopedia of
405 Integer Sequences in the following forms
406 <http://oeis.org/A054531> (etc)
408 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 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 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
427 Math::PlanePath, Math::PlanePath::DiagonalRationals,
428 Math::PlanePath::RationalsTree, Math::PlanePath::CoprimeColumns,
429 Math::PlanePath::DiagonalsOctant
430
432 <http://user42.tuxfamily.org/math-planepath/index.html>
433
435 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020
436 Kevin Ryde
437 This file is part of Math-PlanePath.
439 Math-PlanePath is free software; you can redistribute it and/or modify
441 it under the terms of the GNU General Public License as published by
442 the Free Software Foundation; either version 3, or (at your option) any
443 later version.
444 Math-PlanePath is distributed in the hope that it will be useful, but
446 WITHOUT ANY WARRANTY; without even the implied warranty of
447 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
448 General Public License for more details.
449 You should have received a copy of the GNU General Public License along
451 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
452perl v5.38.0 2023-07-20 Math::PlanePath::GcdRationals(3)