File Coverage

blib/lib/Math/PlanePath/CellularRule54.pm

Criterion	Covered	Total	%
statement	102	119	85.7
branch	25	36	69.4
condition	6	9	66.6
subroutine	21	23	91.3
pod	5	5	100.0
total	159	192	82.8

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 .

# 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 = 129;
use Math::PlanePath;
@ISA = ('Math::PlanePath');
*_divrem = \&Math::PlanePath::_divrem;
*_sqrtint = \&Math::PlanePath::_sqrtint;

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((_sqrtint(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__