libmath-planepath-perl 117-1 / usr / share / perl5 / Math / PlanePath / Diagonals.pm

# Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.


# Leading diagonal 2,8,18 = 2*d^2
#
# cf A185787 lists numerous seqs for rows,columns,diagonals

package Math::PlanePath::Diagonals;
use 5.004;
use strict;
use Carp 'croak';
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 117;
use Math::PlanePath;
@ISA = ('Math::PlanePath');

use Math::PlanePath::Base::Generic
  'round_nearest';

# uncomment this to run the ### lines
#use Smart::Comments;

use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant n_frac_discontinuity => .5;

use constant parameter_info_array =>
  [ { name        => 'direction',
      share_key   => 'direction_downup',
      display     => 'Direction',
      type        => 'enum',
      default     => 'down',
      choices     => ['down','up'],
      choices_display => ['Down','Up'],
      description => 'Number points downwards or upwards along the diagonals.',
    },
    Math::PlanePath::Base::Generic::parameter_info_nstart1(),
    { name        => 'x_start',
      display     => 'X start',
      type        => 'integer',
      default     => 0,
      width       => 3,
      description => 'Starting X coordinate.',
    },
    { name        => 'y_start',
      display     => 'Y start',
      type        => 'integer',
      default     => 0,
      width       => 3,
      description => 'Starting Y coordinate.',
    },
  ];

sub x_minimum {
  my ($self) = @_;
  return $self->{'x_start'};
}
sub y_minimum {
  my ($self) = @_;
  return $self->{'y_start'};
}

sub dx_minimum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? undef  # down jumps back unlimited at bottom
          : -1);   # up at most -1 across
}
sub dx_maximum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? 1       # down at most +1 across
          : undef); # up jumps back across unlimited at top
}

sub dy_minimum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? -1      # down at most -1
          : undef); # up jumps down unlimited at top
}
sub dy_maximum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? undef  # down jumps up unlimited at bottom
          : 1);    # up at most +1
}

sub absdx_minimum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? 0   # N=1 dX=0,dY=1
          : 1); # otherwise always changes
}
sub absdy_minimum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? 1   # otherwise always changes
          : 0); # N=1 dX=1,dY=0
}

# within diagonal X+Y=k is dSum=0
# end of diagonal X=Xstart+k Y=Ystart
#             to  X=Xstart   Y=Ystart+k+1
# is (Xstart + Ystart+k+1) - (Xstart+k + Ystart) = 1 always, to next diagonal
#
use constant dsumxy_minimum => 0; # advancing diagonals
use constant dsumxy_maximum => 1;

sub ddiffxy_minimum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? undef  # "down" jumps back unlimited at bottom
          : -2);   # NW diagonal
}
sub ddiffxy_maximum {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? 2       # SE diagonal
          : undef); # "up" jumps down unlimited at top
}

sub dir_minimum_dxdy {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? (0,1)   # North, vertical at N=1
          : (1,0)); # East,  horiz at N=1
}
sub dir_maximum_dxdy {
  my ($self) = @_;
  return ($self->{'direction'} eq 'down'
          ? (1,-1)    # South-East at N=2
          : (2,-1));  # ESE        at N=3
}

# If Xstart>0 or Ystart>0 then the origin is not reached.
sub rsquared_minimum {
  my ($self) = @_;
  return ((  $self->{'x_start'} > 0 ? $self->{'x_start'}**2 : 0)
          + ($self->{'y_start'} > 0 ? $self->{'y_start'}**2 : 0));
}



#------------------------------------------------------------------------------

sub new {
  my $self = shift->SUPER::new(@_);
  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }

  my $direction = ($self->{'direction'} ||= 'down');
  if (! ($direction eq 'up' || $direction eq 'down')) {
    croak "Unrecognised direction option: ", $direction;
  }

  $self->{'x_start'} ||= 0;
  $self->{'y_start'} ||= 0;
  return $self;
}

# start each diagonal at 0.5 earlier than the integer point
#   d = [    0,   1,   2,   3,   4 ]
#   n = [ -0.5, 0.5, 2.5, 5.5, 9.5 ]
#             +1   +2   +3   +4
#                1    1    1
# N = (1/2 d^2 + 1/2 d - 1/2)
#   = (1/2*$d**2 + 1/2*$d - 1/2)
#   = ((1/2*$d + 1/2)*$d - 1/2)
# d = -1/2 + sqrt(2 * $n + 5/4)
#   = (sqrt(8*$n + 5) -1)/2

sub n_to_xy {
  my ($self, $n) = @_;
  ### Diagonals n_to_xy(): "$n   ".(ref $n || '')

  # adjust to N=0 at origin X=0,Y=0
  $n = $n - $self->{'n_start'};

  my $d;
  {
    my $r = 8*$n + 5;
    if ($r < 1) {
      ### which is N < -0.5 ...
      return;
    }
    $d = int((sqrt(int($r)) - 1) / 2);
    ### assert: $d >= 0
  }

  # subtract for offset into diagonal, range -0.5 <= $n < $d+0.5
  $n -= $d*($d+1)/2;

  my $y = -$n + $d;  # $n first so BigFloat not BigInt from $d
  # and X=$n

  if ($self->{'direction'} eq 'up') {
    ($n,$y) = ($y,$n);
  }
  return ($n + $self->{'x_start'},
          $y + $self->{'y_start'});
}

# round y on an 0.5 downwards so that x=-0.5,y=0.5 gives n=1 which is the
# inverse of n_to_xy() ... or is that inconsistent with other classes doing
# floor() always?
#
# d(d+1)/2+1
#   = (d^2 + d + 2) / 2
#
sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### xy_to_n(): $x, $y
  $x = $x - $self->{'x_start'};   # "-" operator to provoke warning if x==undef
  $y = $y - $self->{'y_start'};
  if ($self->{'direction'} eq 'up') {
    ($x,$y) = ($y,$x);
  }
  $x = round_nearest ($x);
  $y = round_nearest (- $y);
  ### rounded
  ### $x
  ### $y
  if ($x < 0 || $y > 0) {
    return undef;  # outside
  }

  my $d = $x - $y;
  ### $d
  return $d*($d+1)/2 + $x + $self->{'n_start'};
}

# bottom-left to top-right, used by DiagonalsAlternating too
# exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;

  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  if ($x2 < $self->{'x_start'} || $y2 < $self->{'y_start'}) {
    return (1, 0); # rect all negative, no N
  }

  $x1 = max ($x1, $self->{'x_start'});
  $y1 = max ($y1, $self->{'y_start'});

  # exact range bottom left to top right
  return ($self->xy_to_n ($x1,$y1),
          $self->xy_to_n ($x2,$y2));
}

1;
__END__

=for stopwords PlanePath Ryde Math-PlanePath OEIS triangulars sqrt

=head1 NAME

Math::PlanePath::Diagonals -- points in diagonal stripes

=head1 SYNOPSIS

 use Math::PlanePath::Diagonals;
 my $path = Math::PlanePath::Diagonals->new;
 my ($x, $y) = $path->n_to_xy (123);

=head1 DESCRIPTION

This path follows successive diagonals going from the Y axis down to the X
axis.

      6  |  22
      5  |  16  23
      4  |  11  17  24
      3  |   7  12  18  ...
      2  |   4   8  13  19
      1  |   2   5   9  14  20
    Y=0  |   1   3   6  10  15  21
         +-------------------------
           X=0   1   2   3   4   5

X<Triangular numbers>N=1,3,6,10,etc on the X axis is the triangular numbers.
N=1,2,4,7,11,etc on the Y axis is the triangular plus 1, the next point
visited after the X axis.

=head2 Direction

Option C<direction =E<gt> 'up'> reverses the order within each diagonal to
count upward from the X axis.

=cut

# math-image --path=Diagonals,direction=up  --all --output=numbers

=pod

    direction => "up"

      5  |  21
      4  |  15  20
      3  |  10  14  19 ...
      2  |   6   9  13  18  24
      1  |   3   5   8  12  17  23
    Y=0  |   1   2   4   7  11  16  22
         +-----------------------------
           X=0   1   2   3   4   5   6

This is merely a transpose changing X,Y to Y,X, but it's the same as in
C<DiagonalsOctant> and can be handy to control the direction when combining
C<Diagonals> with some other path or calculation.

=head2 N Start

The default is to number points starting N=1 as shown above.  An optional
C<n_start> can give a different start, in the same diagonals sequence.  For
example to start at 0,

=cut

# math-image --path=Diagonals,n_start=0 --all --output=numbers --size=35x5
# math-image --path=Diagonals,n_start=0,direction=up --all --output=numbers --size=35x5

=pod

    n_start => 0,                    n_start=>0
    direction=>"down"                direction=>"up"

      4  |  10                       |  14
      3  |   6 11                    |   9 13
      2  |   3  7 12                 |   5  8 12
      1  |   1  4  8 13              |   2  4  7 11
    Y=0  |   0  2  5  9 14           |   0  1  3  6 10
         +-----------------          +-----------------
           X=0  1  2  3  4             X=0  1  2  3  4

X<Triangular numbers>N=0,1,3,6,10,etc on the Y axis of "down" or the X axis
of "up" is the triangular numbers Y*(Y+1)/2.

=head2 X,Y Start

Options C<x_start =E<gt> $x> and C<y_start =E<gt> $y> give a starting
position for the diagonals.  For example to start at X=1,Y=1

      7  |   22               x_start => 1,
      6  |   16 23            y_start => 1
      5  |   11 17 24
      4  |    7 12 18 ...
      3  |    4  8 13 19
      2  |    2  5  9 14 20
      1  |    1  3  6 10 15 21
    Y=0  |
         +------------------
         X=0  1  2  3  4  5

The effect is merely to add a fixed offset to all X,Y values taken and
returned, but it can be handy to have the path do that to step through
non-negatives or similar.

=head1 FUNCTIONS

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

=over 4

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

=item C<$path = Math::PlanePath::Diagonals-E<gt>new (direction =E<gt> $str, n_start =E<gt> $n, x_start =E<gt> $x, y_start =E<gt> $y)>

Create and return a new path object.  The C<direction> option (a string) can
be

    direction => "down"       the default
    direction => "up"         number upwards from the X axis

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

Return the X,Y coordinates of point number C<$n> on the path.

For C<$n E<lt> 0.5> the return is an empty list, it being considered the
path begins at 1.

=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

=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

=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 = d - X

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 therefor X extends back
to -0.5.

=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 the minimum N
and the upper right is the 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

      A185787    total N in row to X=Y diagonal
      A185788    total N in row to X=Y-1
      A100182    total N in column to Y=X diagonal
      A101165    total N in column to Y=X-1
      A185506    total N in rectangle 0,0 to X,Y

    direction=down, n_start=0
      A023531    dSum = dX+dY, being 1 at N=triangular+1 (and 0)
      A000096    N on X axis, X*(X+3)/2
      A000217    N on Y axis, the triangular numbers
      A129184    turn 1=left,0=right
      A103451    turn 1=left or right,0=straight, but extra initial 1
      A103452    turn 1=left,0=straight,-1=right, but extra initial 1
    direction=up, n_start=0
      A129184    turn 0=left,1=right

    direction=up, n_start=-1
      A023531    turn 1=left,0=right
    direction=down, n_start=-1
      A023531    turn 0=left,1=right

    in direction=up the X,Y coordinate forms are the same but swap X,Y

    either direction, n_start=1
      A038722    permutation N at transpose Y,X
                   which is direction=down <-> direction=up

    n_start=1, x_start=1, y_start=1, either direction
      A003991    X*Y coordinate product
      A003989    GCD(X,Y) greatest common divisor starting (1,1)
      A003983    min(X,Y)
      A051125    max(X,Y)
    n_start=1, x_start=1, y_start=1, direction=down
      A057046    X for N=2^k
      A057047    Y for N=2^k

    n_start=0 (either direction)
      A049581    abs(X-Y) coordinate diff
      A004197    min(X,Y)
      A003984    max(X,Y)
      A004247    X*Y coordinate product
      A048147    X^2+Y^2
      A109004    GCD(X,Y) greatest common divisor starting (0,0)
      A004198    X bit-and Y
      A003986    X bit-or Y
      A003987    X bit-xor Y
      A156319    turn 0=straight,1=left,2=right

      A061579    permutation N at transpose Y,X
                   which is direction=down <-> direction=up

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::DiagonalsAlternating>,
L<Math::PlanePath::DiagonalsOctant>,
L<Math::PlanePath::Corner>,
L<Math::PlanePath::Rows>,
L<Math::PlanePath::Columns>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2010, 2011, 2012, 2013, 2014 Kevin Ryde

This file is part of Math-PlanePath.

Math-PlanePath is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation; either version 3, or (at your option) any later
version.

Math-PlanePath is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
more details.

You should have received a copy of the GNU General Public License along with
Math-PlanePath.  If not, see <http://www.gnu.org/licenses/>.

=cut