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

# Copyright 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/>.

# math-image --path=CellularRule54 --all --scale=10
# math-image --path=CellularRule54 --all --output=numbers --size=132x50

package Math::PlanePath::CellularRule54;
use 5.004;
use strict;

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

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

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


use constant parameter_info_array =>
  [ Math::PlanePath::Base::Generic::parameter_info_nstart1(),
  ];

use constant class_y_negative => 0;
use constant n_frac_discontinuity => .5;
sub x_negative_at_n {
  my ($self) = @_;
  return $self->n_start + 1;
}
use constant sumxy_minimum => 0;  # triangular X>=-Y so X+Y>=0
use constant diffxy_maximum => 0; # triangular X<=Y so X-Y<=0
use constant dx_maximum => 4;
use constant dy_minimum => 0;
use constant dy_maximum => 1;
use constant absdx_minimum => 1;
use constant dsumxy_maximum => 4;  # straight East dX=+4
use constant ddiffxy_maximum => 4; # straight East dX=+4
use constant dir_maximum_dxdy => (-1,0); # supremum, West and dY=+1 up


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

sub new {
  my $self = shift->SUPER::new (@_);
  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }
  return $self;
}
#            left   add
# even  y=0    0     1
#         2    1     2
#         4    3     3
#         6    6     4
# left = y/2*(y/2+1)/2
#      = y*(y+2)/8   of 4-cell figures
# inverse y = -1 + sqrt(2 * $n + -1)
#
#            left   add
# odd   y=1    0     3
#         3    3     6
#         5    9     9
#         7   18    12
# left = 3*(y-1)/2*((y-1)/2+1)/2
#      = 3*(y-1)*(y+1)/8     of 4-cell figures
#
# nbase y even = y*(y+2)/8 + 3*((y+1)-1)*((y+1)+1)/8
#              = [ y*(y+2) + 3*y*(y+2) ] / 8
#              = y*(y+2)/2
# y=0  nbase=0
# y=2  nbase=4
# y=4  nbase=12
# y=6  nbase=24
#
# nbase y odd = 3*(y-1)*(y+1)/8  + (y+1)*(y+3)/8
#             = (y+1) * (3y-3 + y+3)/8
#             = (y+1)*4y/8
#             = y*(y+1)/2
# y=1  nbase=1
# y=3  nbase=6
# y=5  nbase=15
# y=7  nbase=28
# inverse y = -1/2 + sqrt(2 * $n + -7/4)
#           = sqrt(2n-7/4) - 1/2
#           = (2*sqrt(2n-7/4) - 1)/2
#           = (sqrt(4n-7)-1)/2
#
# dual
# d = [ 0, 1,  2,  3 ]
# N = [ 1, 5, 13, 25 ]
# N = (2 d^2 + 2 d + 1)
#   = ((2*$d + 2)*$d + 1)
# d = -1/2 + sqrt(1/2 * $n + -1/4)
#   = sqrt(1/2 * $n + -1/4) - 1/2
#   = [ 2*sqrt(1/2 * $n + -1/4) - 1 ] / 2
#   = [ sqrt(4/2 * $n + -4/4) - 1 ] / 2
#   = [ sqrt(2*$n - 1) - 1 ] / 2
#

sub n_to_xy {
  my ($self, $n) = @_;
  ### CellularRule54 n_to_xy(): $n

  $n = $n - $self->{'n_start'}; # to N=0 basis
  my $frac;
  {
    my $int = int($n);
    $frac = $n - $int;
    $n = $int;       # BigFloat int() gives BigInt, use that
    if (2*$frac >= 1) {  # $frac>=0.5 and BigInt friendly
      $frac -= 1;
      $n += 1;
    }
    # -0.5 <= $frac < 0.5
    ### assert: $frac >= -0.5
    ### assert: $frac < 0.5
  }

  if ($n < 0) {
    return;
  }

  # d is the two-row group number, d=2*y, where n belongs
  # start of the two-row group is nbase = 2 d^2 + 2 d starting from N=0 
  #
  my $d = int((sqrt(2*$n+1) - 1) / 2);
  $n -= (2*$d + 2)*$d;   # remainder within two-row
  ### $d
  ### remainder: $n
  if ($n <= $d) {
    # d+1 many points in the Y=0,2,4,6 etc even row, spaced 4*n apart
    $d *= 2;    # y=2*d
    return ($frac + 4*$n - $d,
            $d);
  } else {
    # 3*d many points in the Y=1,3,5,7 etc odd row, using 3 in 4 cells
    $n -= $d+1;    # remainder 0 upwards into odd row
    $d = 2*$d+1;   # y=2*d+1
    my ($q) = _divrem($n,3);
    return ($frac + $n + $q - $d,
            $d);
  }
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  $x = round_nearest ($x);
  $y = round_nearest ($y);
  ### CellularRule54 xy_to_n(): "$x,$y"

  if ($y < 0
      || $x < -$y
      || $x > $y) {
    return undef;
  }
  $x += $y;
  ### x centred: $x
  if ($y % 2) {
    ### odd row, 3 in 4 ...
    if (($x % 4) == 3) {
      return undef;
    }
    return $x - int($x/4) + $y*($y+1)/2 + $self->{'n_start'};
  } else {
    ## even row, sparse ...
    if ($x % 4) {
      return undef;
    }
    return $x/4 + $y*($y+2)/2 + $self->{'n_start'};
  }
}

# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### CellularRule54 rect_to_n_range(): "$x1,$y1, $x2,$y2"

  ($x1,$y1, $x2,$y2) = _rect_for_V ($x1,$y1, $x2,$y2)
    or return (1,0); # rect outside pyramid

  my $zero = ($x1 * 0 * $y1 * $x2 * $y2);  # inherit bignum

  # nbase y even y*(y+2)/2
  # nbase y odd  y*(y+1)/2
  # y even end (y+1)*(y+2)/2
  # y odd end  (y+1)*(y+3)/2

  $y2 += 1;
  return (# even/odd left end
          $zero + $y1*($y1 + 2-($y1%2))/2 + $self->{'n_start'},

          # even/odd right end
          $zero + $y2*($y2 + 2-($y2%2))/2 + $self->{'n_start'} - 1);
}

# Return ($x1,$y1, $x2,$y2) which is the rectangle part chopped to the top
# row entirely within the pyramid V and the bottom row partly within.
#
sub _rect_for_V {
  my ($x1,$y1, $x2,$y2) = @_;
  ### _rect_for_V(): "$x1,$y1, $x2,$y2"

  $y1 = round_nearest ($y1);
  $y2 = round_nearest ($y2);
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); } # swap to y1<=y2

  unless ($y2 >= 0) {
    ### rect all negative, no N ...
    return;
  }
  unless ($y1 >= 0) {
    # increase y1 to zero, including negative infinity discarded
    $y1 = 0;
  }

  $x1 = round_nearest ($x1);
  $x2 = round_nearest ($x2);
  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); } # swap to x1<=x2
  my $neg_y2 = -$y2;

  #     \        /
  #   y2 \      / +-----
  #       \    /  |
  #        \  /
  #         \/    x1
  #
  #        \        /
  #   ----+ \      /  y2
  #       |  \    /
  #           \  /
  #       x2   \/
  #
  if ($x1 > $y2            # off to the right
      || $x2 < $neg_y2) {  # off to the left
    ### rect all off to the left or right, no N
    return;
  }

  #     \        /  x2
  #      \   +------+ y2
  #       \  | /    |
  #        \ +------+
  #         \/
  #
  if ($x2 > $y2) {
    ### top-right beyond pyramid, reduce ...
    $x2 = $y2;
  }

  #
  #    x1  \        /
  # y2 +--------+  /  y2
  #    |     \  | /
  #    +--------+/
  #            \/
  #
  if ($x1 < $neg_y2) {
    ### top-left beyond pyramid, increase ...
    $x1 = $neg_y2;
  }

  #     \       | /
  #      \      |/
  #       \    /|       |
  #    y1  \  / +-------+
  #         \/  x1
  #
  #        \|       /
  #         \      /
  #         |\    /
  #  -------+ \  /   y1
  #        x2  \/
  #
  # in both of the following y1=x2 or y1=-x2 leaves y1<=y2 because have
  # already established some part of the rectangle is in the V shape
  #
  if ($x1 > $y1) {
    ### x1 off to the right, so y1 row is outside, increase y1 ...
    $y1 = $x1;

  } elsif ((my $neg_x2 = -$x2) > $y1) {
    ### x2 off to the left, so y1 row is outside, increase y1 ...
    $y1 = $neg_x2;
  }

  # values ordered
  ### assert: $x1 <= $x2
  ### assert: $y1 <= $y2

  # top row x1..x2 entirely within pyramid
  ### assert: $x1 >= -$y2
  ### assert: $x2 <= $y2

  # bottom row x1..x2 some part within pyramid
  ### assert: $x1 <= $y1
  ### assert: $x2 >= -$y1

  return ($x1,$y1, $x2,$y2);
}

1;
__END__

=for stopwords straight-ish Ryde Math-PlanePath ie hexagonals 18-gonal Xmax-Xmin Nleft Nright OEIS

=head1 NAME

Math::PlanePath::CellularRule54 -- cellular automaton points

=head1 SYNOPSIS

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

=head1 DESCRIPTION

X<Wolfram, Stephen>This is the pattern of Stephen Wolfram's "rule 54"
cellular automaton

=over

L<http://mathworld.wolfram.com/Rule54.html>

=back

arranged as rows,

    29 30 31  . 32 33 34  . 35 36 37  . 38 39 40     7
       25  .  .  . 26  .  .  . 27  .  .  . 28        6
          16 17 18  . 19 20 21  . 22 23 24           5
             13  .  .  . 14  .  .  . 15              4
                 7  8  9  . 10 11 12                 3
                    5  .  .  .  6                    2
                       2  3  4                       1
                          1                      <- Y=0

    -7 -6 -5 -4 -3 -2 -1 X=0 1  2  3  4  5  6  7

The initial figure N=1,2,3,4 repeats in two-row groups with 1 cell gap
between figures.  Each two-row group has one extra figure, for a step of 4
more points than the previous two-row.

X<Hexagonal numbers>The rightmost N on the even rows Y=0,2,4,6 etc is the
hexagonal numbers N=1,6,15,28, etc k*(2k-1).  The hexagonal numbers of the
"second kind" 1, 3, 10, 21, 36, etc j*(2j+1) are a steep sloping line
upwards in the middle too.  Those two taken together are the
X<Triangular numbers>triangular numbers 1,3,6,10,15 etc, k*(k+1)/2.

The 18-gonal numbers 18,51,100,etc are the vertical line at X=-3 on every
fourth row Y=5,9,13,etc.

=head2 Row Ranges

The left end of each row is

    Nleft = Y*(Y+2)/2 + 1     if Y even
            Y*(Y+1)/2 + 1     if Y odd

The right end is

    Nright = (Y+1)*(Y+2)/2    if Y even
             (Y+1)*(Y+3)/2    if Y odd

           = Nleft(Y+1) - 1   ie. 1 before next Nleft

The row width Xmax-Xmin is 2*Y but with the gaps the number of visited
points in a row is less than that, being either about 1/4 or 3/4 of the
width on even or odd rows.

    rowpoints = Y/2 + 1        if Y even
                3*(Y+1)/2      if Y odd

For any Y of course the Nleft to Nright difference is the number of points
in the row too

    rowpoints = Nright - Nleft + 1

=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 pattern.  For example to
start at 0,

=cut

# math-image --path=CellularRule54,n_start=0 --all --output=numbers --size=75x6

=pod

    n_start => 0

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

    -5 -4 -3 -2 -1 X=0 1  2  3  4  5

=head1 FUNCTIONS

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

=over 4

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

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

Create and return a new path object.

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

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

=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 cell
as a square of side 1.  If C<$x,$y> is outside the pyramid or on a skipped
cell the return is C<undef>.

=back

=head1 OEIS

This pattern is in Sloane's Online Encyclopedia of Integer Sequences in a
couple of forms,

=over

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

=back

    A118108    whole-row used cells as bits of a bignum
    A118109    1/0 used and unused cells across rows

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::CellularRule>,
L<Math::PlanePath::CellularRule57>,
L<Math::PlanePath::CellularRule190>,
L<Math::PlanePath::PyramidRows>

L<Cellular::Automata::Wolfram>

L<http://mathworld.wolfram.com/Rule54.html>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 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