view release on metacpan or  search on metacpan

lib/Math/PlanePath/GcdRationals.pm  view on Meta::CPAN

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

=head1 FUNCTIONS

See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.

=over

=item C<$path = Math::PlanePath::GcdRationals-E<gt>new ()>

=item C<$path = Math::PlanePath::GcdRationals-E<gt>new (pairs_order =E<gt> $str)>

Create and return a new path object.  The C<pairs_order> option can be

    "rows"               (default)
    "rows_reverse"
    "diagonals_down"
    "diagonals_up"

=back

=head1 FORMULAS

=head2 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 C<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 0E<lt>=remE<lt>Y.  This will give i=0 which is outside the
intended 1E<lt>=iE<lt>=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.

=head2 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
L</Rows Reverse> above).

=cut

# No, not quite
#
# =head2 Rectangle N Range -- Rows
#
# An over-estimate of the N range can be calculated just from the X,Y to N
# formula above.
#
# 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 of the
# rectangle, or the top of the next lower wedge, ie. smaller Y but bigger g.
# Conversely the minimum is either the bottom right of the rectangle, or the
# bottom of the next higher wedge, ie. smaller g but bigger Y.  (Is that
# right?)
#
# This is an over-estimate because it ignores which X,Y points are coprime and
# thus actually should have N values.
#
# =head2 Rectangle N Range -- Rows Reverse
#
# When row pairs are taken in reverse order increasing X is not increasing N,
# but rather the maximum N of a row is at the left end of the wedge.

=pod

=head1 OEIS

This enumeration of rationals is in Sloane's Online Encyclopedia of Integer
Sequences in the following forms

=over

L<http://oeis.org/A054531> (etc)

=back

    pairs_order="rows" (the default)
      A226314   X coordinate
      A054531   Y coordinate, being N/GCD(i,j)
      A000124   N in X=1 column, triangular+1
      A050873   ceil(X/Y), gcd by rows
      A050873-A023532  floor(X/Y)
                gcd by rows and subtract 1 unless i=j

    pairs_order="diagonals_down"
      A033638   N in X=1 column, quartersquares+1 and pronic+1
      A000290   N in Y=1 row, perfect squares

    pairs_order="diagonals_up"
      A002620   N in X=1 column, squares and pronics
      A002061   N in Y=1 row, central polygonals (extra initial 1)
      A002522   N at Y=X+1 above leading diagonal, squares+1

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::DiagonalRationals>,
L<Math::PlanePath::RationalsTree>,
L<Math::PlanePath::CoprimeColumns>,
L<Math::PlanePath::DiagonalsOctant>

=head1 HOME PAGE

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