view release on metacpan or  search on metacpan

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

=item C<$n = $path-E<gt>xy_to_n ($x,$y)>

Return the point number for coordinates C<$x,$y>.  C<$x> and C<$y> are
each rounded to the nearest integer, which has the effect of treating each
point C<$n> as a square of side 1, so the quadrant x>=-0.5, y>=-0.5 is
entirely covered.

=item C<($n_lo, $n_hi) = $path-E<gt>rect_to_n_range ($x1,$y1, $x2,$y2)>

The returned range is exact, meaning C<$n_lo> and C<$n_hi> are the smallest
and biggest in the rectangle.

=back

=head1 FORMULAS

=head2 X,Y to N

The sum d=X+Y numbers each diagonal from d=0 upwards, corresponding to the Y
coordinate where the diagonal starts (or X if direction=up).

    d=2
        \
    d=1  \
        \ \
    d=0  \ \
        \ \ \

N is then given by

    d = X+Y
    N = d*(d+1)/2 + X + Nstart

The d*(d+1)/2 shows how the triangular numbers fall on the Y axis when X=0
and Nstart=0.  For the default Nstart=1 it's 1 more than the triangulars, as
noted above.

=cut

# N = (X+Y)*(X+Y+1)/2 + X + Nstart
#   = [ (X+Y)*(X+Y+1) + 2X ]/2 + Nstart
#   = [ X^2 + XY + X + XY + Y^2 + Y + 2X ]/2 + Nstart
#   = [ X^2 + 3X + 2XY + Y + Y^2 ]/2 + Nstart
# GP-Test  my(d='x+'y); d*(d+1)/2 + 'x == (('x+'y)^2 + 3*'x + 'y)/2

=pod

d can be expanded out to the following quite symmetric form.  This almost
suggests something parabolic but is still the straight line diagonals.

        X^2 + 3X + 2XY + Y + Y^2
    N = ------------------------ + Nstart
                   2

        (X+Y)^2 + 3X + Y
      = ---------------- + Nstart       (using one square)
                2

=head2 N to X,Y

The above formula N=d*(d+1)/2 can be solved for d as

    d = floor( (sqrt(8*N+1) - 1)/2 )
    # with n_start=0

For example N=12 is d=floor((sqrt(8*12+1)-1)/2)=4 as that N falls in the
fifth diagonal.  Then the offset from the Y axis NY=d*(d-1)/2 is the X
position,

    X = N - d*(d+1)/2
    Y = X - d

In the code, fractional N is handled by imagining each diagonal beginning
0.5 back from the Y axis.  That's handled by adding 0.5 into the sqrt, which
is +4 onto the 8*N.

    d = floor( (sqrt(8*N+5) - 1)/2 )
    # N>=-0.5

The X and Y formulas are unchanged, since N=d*(d-1)/2 is still the Y axis.
But each diagonal d begins up to 0.5 before that and therefore X extends
back to -0.5.

=cut

# GP-DEFINE  \\ diagonal length
# GP-DEFINE  d(n) = floor( (sqrt(8*n+1) - 1)/2 );
# GP-Test  vector(16,n,n--; d(n)) == \
# GP-Test    [0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5]
# GP-DEFINE  X(n) = my(d=d(n)); n - d*(d+1)/2;
# GP-DEFINE  Y(n) = my(d=d(n)); d - X(n);
# GP-Test  vector(16,n,n--; X(n)) == \
# GP-Test    [0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0]
# GP-Test  vector(16,n,n--; Y(n)) == \
# GP-Test    [0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0, 5]
# GP-Test  my(l=List([])); \
# GP-Test  for(d=0,3, for(x=0,d, my(y=d-x); x+y==d || error(); \
# GP-Test     listput(l,[x,y]))); \
# GP-Test  vector(#l,n,n--; [X(n),Y(n)]) == Vec(l)

=pod

=head2 Rectangle to N Range

Within each row increasing X is increasing N, and in each column increasing
Y is increasing N.  So in a rectangle the lower left corner is minimum N and
the upper right is maximum N.

    |            \     \ N max
    |       \ ----------+
    |        |     \    |\
    |        |\     \   |
    |       \| \     \  |
    |        +----------
    |  N min  \  \     \
    +-------------------------

=head1 OEIS

Entries in Sloane's Online Encyclopedia of Integer Sequences related to this
path include

=over

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

=back

    direction=down (the default)
      A002262    X coordinate, runs 0 to k
      A025581  	 Y coordinate, runs k to 0
      A003056  	 X+Y coordinate sum, k repeated k+1 times
      A114327  	 Y-X coordinate diff
      A101080    HammingDist(X,Y)

      A127949    dY, change in Y coordinate

      A000124    N on Y axis, triangular numbers + 1
      A001844    N on X=Y diagonal