File Coverage

blib/lib/Math/PlanePath/CellularRule190.pm

Criterion	Covered	Total	%
statement	101	110	91.8
branch	34	42	80.9
condition	2	3	66.6
subroutine	22	23	95.6
pod	5	5	100.0
total	164	183	89.6

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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://mathworld.wolfram.com/ElementaryCellularAutomaton.html
#
# Loeschian numbers strips on the right ...


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

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

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

use Math::PlanePath::CellularRule54;
*_rect_for_V = \&Math::PlanePath::CellularRule54::_rect_for_V;

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


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

use constant parameter_info_array =>
  [ { name        => 'mirror',
      display     => 'Mirror',
      type        => 'boolean',
      default     => 0,
      description => 'Mirror to "rule 246" instead.',
    },
    Math::PlanePath::Base::Generic::parameter_info_nstart1(),
  ];

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 => 2; # across gap
use constant dy_minimum => 0;
use constant dy_maximum => 1;
use constant absdx_minimum => 1;
use constant dsumxy_maximum => 2; # straight East dX=+2
use constant ddiffxy_maximum => 2; # straight East dX=+2
use constant dir_maximum_dxdy => (-1,0); # supremum, West except dY=+1


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

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

# 31 32 33   34 35 36    37 38 39    40
#    22 23 24   25 26 27    28 29 30
#       15 6 17    18 19 20    21
#          9 10 11    12 13 14
#             5  6  7     8
#                2  3  4
#                   1
#
# even  y = [ 0, 2, 4, 6 ]
#       N = [ 1, 5, 15, 31 ]
# Neven = (3/4 y^2 + 1/2 y + 1)
#       = (3y + 2)*y/4 + 1
#       = ((3y + 2)*y + 4) /4
#       = (3 (y/2)^2 + (y/2) + 1)
#       = (3*(y/2) + 1)*(y/2) + 1
#
# odd  y = [ 1, 3, 5,7 ]
#      N = [ 2,9,22,41 ]
# Nodd = (3/4 y^2 + 1/2 y + 3/4)
#      = ((3y+2)*y+ 3) / 4
#
# pair even d = [0,1,2,3]
#           N = [ 1, 5, 15, 31 ]
# Npair = (3 d^2 + d + 1)
# d = -1/6 + sqrt(1/3 * $n + -11/36)
#   = [ -1 + sqrt(1/3 * $n + -11/36)*6 ] / 6
#   = [ -1 + sqrt(1/3 * $n*36 + -11/36*36) ] / 6
#   = [ -1 + sqrt(12n-11) ] / 6
#
sub n_to_xy {
  my ($self, $n) = @_;
  ### CellularRule190 n_to_xy(): $n

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

  if ($n < 1) {
    return;
  }

  # d is the two-row number, ie. d=2*y, where n belongs
  # start of the two-row group is nbase = 3 d^2 + d + 1
  #
  my $d = int ((_sqrtint(12*$n-11) - 1) / 6);
  $n -= ((3*$d + 1)*$d + 1);   # remainder within two-row
  ### $d
  ### remainder: $n
  if ($n <= 3*$d) {
    # 3d+1 many points in the Y=0,2,4,6 etc even row
    $d *= 2;    # y=2*d
    return ($frac + $n + int(($n + ($self->{'mirror'} ? 2 : 0))/3) - $d,
            $d);
  } else {
    # 3*d many points in the Y=1,3,5,7 etc odd row, using 3 in 4 cells
    $n -= 3*$d+1;    # remainder 0 upwards into odd row
    $d = 2*$d+1;   # y=2*d+1
    return ($frac + $n + int($n/3) - $d,
            $d);
  }
}

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

  if ($y < 0 || $x > $y) {
    return undef;
  }

  $x += $y;  # move to have x=0 the start of the row
  if ($x < 0) {
    return undef;
  }

  ### x centred: $x
  if ($y % 2) {
    ### odd row, 3s from the start ...
    if (($x % 4) == 3) {
      return undef;
    }
    # 3y^2+2y-1 = (3y-1)*(y+1)
    return ($x
            - int($x/4)
            + ((3*$y+2)*$y-1)/4
            + $self->{'n_start'});
  } else {
    ### even row, 3s then 1 sep, or mirror 1 sep start ...
    my $mirror = $self->{'mirror'};
    if (($x % 4) == ($mirror ? 1 : 3)) {
      return undef;
    }
    return ($x
            - int(($x+($mirror ? 2 : 1))/4)
            + (3*$y+2)*$y/4
            + $self->{'n_start'});
  }
}

# exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### CellularRule190 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

  # # inherit bignum (before collapsing some y1 to x1 etc)
  # my $zero = ($x1 * 0 * $y1 * $x2 * $y2);

  my $mirror = $self->{'mirror'};
  my $unincremented_x1 = $x1;

  #    \+------+
  #     |      |     /
  #     |\     |    /
  #     | \    |   /
  #     |  \   |  /
  #  y1 +------+ /
  #    x1    \  /
  #           \/
  #
  if ($x1 < (my $neg_y1 = -$y1)) {
    ### bottom-left outside, move across to: "$neg_y1,$y1"
    $x1 = $neg_y1;

    # For the following blank checks a blank doesn't occur at the ends of a
    # row, so when on a blank it's always possible to increment or decrement
    # X to go to a non-blank -- as long as that adjacent space is within the
    # rectangle.
    #
  } elsif ((($mirror ? $y1-$x1 : $x1+$y1) % 4) == 3) {
    ### x1,y1 bottom left is on a blank: "x1+y1=".($x1+$y1)
    if ($x1 < $x2) {
      ### rect width >= 2, so increase x1 ...
      $x1 += 1;
    } else {
      ### rect is a single column width==1, increase y1 ...
      if (($y1 += 1) > $y2) {
        ### rect was a single blank square, contains no N ...
        return (1,0);
      }
    }
  }

  if ((($mirror ? $y2-$x2 : $x2+$y2) % 4) == 3) {
    ### x2,y2 top right is on a blank, decrement ...
    if ($x2 > $unincremented_x1) {
      ### rect width >= 2, so decrease x2 ...
      $x2 -= 1;
    } else {
      ### rect is a single column width==1, decrease y2 ...
      $y2 -= 1;

      # Can decrement without checking whether the rect is a single square.
      # If the rect was a single blank square then the x1+y1 bottom left
      # above detects and returns.  And the bottom left blank check
      # incremented y1 to leave a single square then that's a non-blank
      # because there's no vertical blank pairs (they go on the diagonal).
      ### assert $y2 >= $y1
    }
  }

  # At this point $x1,$y1 is a non-blank bottom left corner, and $x2,$y2
  # is a non-blank top right corner, being the N lo to hi range.

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


# old rect_to_n_range() of row endpoints
#
# # even right  y = [ 0, 2, 4, 6 ]
# #             N = [ 1,8,21,40 ]
# # Nright = (3/4 y^2 + 2 y + 1)
# #        = (3 y^2 + 8 y + 4) / 4
# #        = ((3y + 8)y + 4) / 4
# #
# # odd right  y = [ 1, 3, 5, 7 ]
# #            N = [ 4,14,30, 52 ]
# # Nright = (3/4 y^2 + 2 y + 5/4)
# #        = (3 y^2 + 8 y + 5) / 4
# #        = ((3y + 8)y + 5) / 4
# #
# # Nleft y even ((3y+2)*y + 4)/4
# # Nleft y odd  ((3y+2)*y + 3)/4
# # Nright even ((3(y+1)+2)*(y+1) + 3)/4 - 1
# #          = ((3y+3+2)*(y+1) + 3 - 4)/4
# #          = ((3y+5)*(y+1) - 1)/4
# #          = ((3y^2 + 8y + 5 - 1)/4
# #          = ((3y^2 + 8y + 4)/4
# #          = ((3y+8)y + 4)/4
# #          = ((3y+2)(y+2)/4
# #
# ### $y1
# ### $y2
# $y2 += $zero;
# $y1 += $zero;
# return (((3*$y1 + 2)*$y1 + 4 - ($y1%2)) / 4,    # even/odd Nleft
#         ((3*$y2 + 8)*$y2 + 4 + ($y2%2)) / 4);   # even/odd Nright

1;
__END__