Math::PlanePath::GcdRationals(3pm)

Triangular Numbers

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.

    N=k*(k-1)/2  i=k,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 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 in X=1 column = Y*(Y-1)/2 + 1

Primes

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. The rest of the line X=2*k,Y=k is not visited since common factor k would mean 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 taking out 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.

Rows Reverse

Option "pairs_order => "rows_reverse"" reverses the order of points within the rows of i,j pairs,

    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" above, 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.

Diagonals

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.)

    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

X,Y to N Ω- Rows

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 (for "xy_to_n()" it is).

    X/Y = g-1 + (i/g)/(j/g)

The g-1 integer part is recovered by a division X divide Y,

    X = quot*Y + rem   division by Y rounded towards 0
      where 0 <= rem < Y
      unless Y=1 in which case use quot=X-1, rem=1
    g-1 = quot
    g = quot+1

The Y=1 special case can instead be left as the usual kind of division quot=X,rem=0, so 0<=rem<Y. This will give i=0 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 then it can be left rem=0.

Equating the denominators in the X/Y formula above gives j by

    Y = j/g          the definition above
    j = g*Y
      = (quot+1)*Y
    j = X+Y-rem      per the division X=quot*Y+rem

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 - quot*Y*g
      = X*g - (X-rem)*g     per the division X=quot*Y+rem
    i = rem*g
    i = rem*(quot+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 quot=2,rem=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 rem, by taking j=(quot+1)*Y instead of j=X+Y-rem, but usually a division operation gives the remainder at no extra cost, or a cost small enough that it's worth swapping a multiply for an add or two.

The gcd g can be recovered by rounding up in the division, instead of rounding down and then incrementing with g=quot+1.

   g = ceil(X/Y)
     = cquot for division X=cquot*Y - crem

But division in most programming languages is towards 0 or towards -infinity, not upwards towards +infinity.

X,Y to N Ω- Rows Reverse

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 rem=0 giving i=0 and j+1 since the reversal j-i+1 is then not correct. Either use rem=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 quot=5,rem=0 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).