File Coverage

blib/lib/Math/PlanePath/SierpinskiTriangle.pm

Criterion	Covered	Total	%
statement	163	252	64.6
branch	61	130	46.9
condition	16	38	42.1
subroutine	25	41	60.9
pod	23	23	100.0
total	288	484	59.5

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


# Maybe:
#
# rule 22 includes the midpoint between adjacent leaf points.
# math-image --path=CellularRule,rule=22 --all --text
#
# rule 126 extra cell to the inward side of each
# math-image --path=CellularRule,rule=60 --all --text
#
# cf rule 150 double ups, something base 2 instead
# math-image --path=CellularRule,rule=150 --all
#
# cf rule 182 filled gaps
# math-image --path=CellularRule,rule=182 --all

# math-image --path=SierpinskiTriangle --all --scale=5
# math-image --path=SierpinskiTriangle --all --output=numbers
# math-image --path=SierpinskiTriangle --all --text --size=80

# Number of cells in a row:
#    numerator of (2^k)/k!
#
# cf A067771  vertices of sierpinski graph, joins up replications
#             so 1 less each giving 3*(3^k-1)/2
#




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

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

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_down_pow',
  'bit_split_lowtohigh',
  'digit_join_lowtohigh';

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

use constant parameter_info_array =>
  [ { name      => 'align',
      share_key => 'align_trld',
      display   => 'Align',
      type      => 'enum',
      default   => 'triangular',
      choices   => ['triangular', 'right', 'left','diagonal'],
      choices_display => ['Triangular', 'Right', 'Left','Diagonal'],
    },
    # { name      => 'parts',
    #   share_key => 'parts_alr',
    #   display   => 'Parts',
    #   type      => 'enum',
    #   default   => 'all',
    #   choices   => ['all', 'left', 'right'],
    #   choices_display => ['All', 'Left', 'Right'],
    # },
    Math::PlanePath::Base::Generic::parameter_info_nstart0(),
  ];

use constant default_n_start => 0;
use constant class_y_negative => 0;
use constant n_frac_discontinuity => .5;
use constant tree_num_children_list => (0,1,2);

sub x_negative {
  my ($self) = @_;
  return ($self->{'align'} eq 'left'
          || ($self->{'align'} eq 'triangular'
              && $self->{'parts'} ne 'right'));
}
sub x_negative_at_n {
  my ($self) = @_;
  return ($self->{'align'} eq 'left'
          || ($self->{'align'} eq 'triangular'
              && $self->{'parts'} ne 'right')
          ? $self->n_start + 1
          : undef);
}

# Note: this method shared by SierpinskiArrowhead
sub x_maximum {
  my ($self) = @_;
  return ($self->{'align'} eq 'left'
          || ($self->{'align'} eq 'triangular'
              && ($self->{'parts'}||'all') eq 'left')
          ? 0       # left all X<=0
          : undef); # others X to +infinity
}
use constant sumxy_minimum => 0;  # triangular X>=-Y or all X>=0

sub diffxy_minimum {
  my ($self) = @_;
  return ($self->{'align'} eq 'right'
          && $self->{'parts'} eq 'right'
          ? 0       # X>=Y so X-Y>=0
          : undef);
}

# Note: this method shared by SierpinskiArrowhead, SierpinskiArrowheadCentres
sub diffxy_maximum {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal'
          ? undef
          : 0);    # triangular X<=Y so X-Y<=0
}

sub dy_minimum {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal' ? undef : 0);
}
sub dy_maximum {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal' ? undef : 1);
}
{
  my %absdx_minimum = (triangular => 1,
                       left       => 1,
                       right      => 0,  # at N=0
                       diagonal   => 0); # at N=0
  sub absdx_minimum {
    my ($self) = @_;
    return $absdx_minimum{$self->{'align'}};
  }
}
{
  my %absdy_minimum = (triangular => 0,  # rows
                       left       => 0,  # rows
                       right      => 0,  # rows
                       diagonal   => 1); # diagonal always moves
  sub absdy_minimum {
    my ($self) = @_;
    return $absdy_minimum{$self->{'align'}};
  }
}

sub dsumxy_minimum {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal'
          ? 0         # X+Y constant along diagonals
          : undef);
}
sub dsumxy_maximum {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal'
          ? 1         # X+Y increase by 1 to next diagonal
          : undef);
}

sub dir_minimum_dxdy {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal'
          ? (0,1)   # North
          : (1,0)); # East
}
sub dir_maximum_dxdy {
  my ($self) = @_;
  return ($self->{'align'} eq 'diagonal'
          ? (1,-1)   # South-Eest
          : (-1,0)); # supremum, West and 1 up
}


#------------------------------------------------------------------------------
my %align_known = (triangular => 1,
                   left       => 1,
                   right      => 1,
                   diagonal   => 1);

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

  my $align = ($self->{'align'} ||= 'triangular');
  if (! $align_known{$align}) {
    croak "Unrecognised align option: ", $align;
  }
  ### $align

  return $self;
}

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

  # written as $n - n_start() rather than "-=" so as to provoke an
  # uninitialized value warning if $n==undef
  $n = $n - $self->{'n_start'};   # N=0 basis

  # this frac behaviour slightly unspecified yet
  my $frac;
  {
    my $int = int($n);
    $frac = $n - $int;
    if (2*$frac >= 1) {        # $frac>=0.5 and BigInt friendly
      $frac -= 1;
      $int += 1;
    } elsif (2*$frac < -1) {   # $frac<0.5 and BigInt friendly
      $frac += 1;
      $int -= 1;
    }
    $n = $int;
  }
  ### $n
  ### $frac

  if ($n < 0) {
    return;
  }
  if ($n == 0) {
    return ($n,$n);
  }

  my ($depthbits, $ndepth) = _n0_to_depthbits($n, $self->{'parts'})
    or return ($n,$n); # infinite

  ### $depthbits
  ### $ndepth

  my $zero = $n * 0;
  $n -= $ndepth;  # offset into row
  my @nbits = bit_split_lowtohigh($n);

  # Where there's a 0-bit in the depth remains a 0-bit.
  # Where there's a 1-bit in the depth takes a bit from Noffset.
  # Small Noffset has less bits than the depth 1s, hence "|| 0".
  #
  my @xbits = map {$_ && (shift @nbits || 0)} @$depthbits;
  ### @xbits

  my $x = digit_join_lowtohigh (\@xbits,    2, $zero);
  my $y = digit_join_lowtohigh ($depthbits, 2, $zero);

  ### n_to_xy as right: "$x,$y"

  # $x,$y is in the style of align=right, transform to others
  if ($self->{'align'} eq 'left') {
    $x -= $y;
  } elsif ($self->{'align'} eq 'diagonal') {
    $y -= $x;
  } elsif ($self->{'align'} eq 'triangular') {
    $x = 2*$x - $y;
  }

  ### n_to_xy final: "$x,$y"
  return ($x, $y);
}

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

  $y = round_nearest ($y);
  $x = round_nearest($x);

  # transform $x,$y to the style of align=right
  if ($self->{'align'} eq 'diagonal') {
    $y += $x;
  } elsif ($self->{'align'} eq 'left') {
    $x += $y;
  } elsif ($self->{'align'} eq 'triangular') {
    $x += $y;
    if (_divrem_mutate ($x, 2)) {
      # if odd point
      return undef;
    }
  }
  ### adjusted xy: "$x,$y"

  return _right_xy_to_n ($self, $x, $y);
}

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

  unless ($x >= 0 && $x <= $y && $y >= 0) {
    ### outside horizontal row range ...
    return undef;
  }
  if (is_infinite($y)) {
    return $y;
  }

  my $zero = ($y * 0);
  my $n = $zero;          # inherit bignum 0
  my $npower = $zero+1;   # inherit bignum 1

  my @xbits = bit_split_lowtohigh($x);
  my @depthbits = bit_split_lowtohigh($y);

  my @nbits;  # N offset into row
  foreach my $i (0 .. $#depthbits) {      # x,y bits low to high
    if ($depthbits[$i]) {
      $n = 2*$n + $npower;
      push @nbits, $xbits[$i] || 0;   # low to high
    } else {
      if ($xbits[$i]) {
        return undef;
      }
    }
    $npower *= 3;
  }

  ### n at left end of y row: $n
  ### n offset for x: @nbits
  ### total: $n + digit_join_lowtohigh(\@nbits,2,$zero) + $self->{'n_start'}

  return $n + digit_join_lowtohigh(\@nbits,2,$zero) + $self->{'n_start'};
}

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

  $y1 = round_nearest ($y1);
  $y2 = round_nearest ($y2);
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1) }

  $x1 = round_nearest ($x1);
  $x2 = round_nearest ($x2);
  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1) }

  # $y1 to $y2 is the depth range for "triangular", "right" and "left".
  # For "diagonal" must use X+Y to reckon by anti-diagonals.
  #
  if ($self->{'align'} eq 'diagonal') {
    $y2 += $x2;
    $y1 += $x1;
  }

  if ($y2 < 0) {
    return (1, 0);
  }
  if ($y1 < 0) {
    $y1 *= 0;  # preserve any bignum $y1
  }
  return ($self->tree_depth_to_n($y1),
          $self->tree_depth_to_n_end($y2));
}

# To get N within a triangle row, based on the X range ...
#
# use Math::PlanePath::CellularRule54;
# *_rect_for_V = \&Math::PlanePath::CellularRule54::_rect_for_V;
#
# if ($self->{'align'} eq 'diagonal') {
#   if ($x2 < 0 || $y2 < 0) {
#     return (1,0);
#   }
#   if ($x1 < 0) { $x1 *= 0; }
#   if ($y1 < 0) { $y1 *= 0; }
#
#   return ($self->xy_to_n(0, $x1+$y1),
#           $self->xy_to_n($x2+$y2, 0));
# }
#
# ($x1,$y1, $x2,$y2) = _rect_for_V ($x1,$y1, $x2,$y2)
#   or return (1,0); # rect outside pyramid
#
# return ($self->xy_to_n($self->{'align'} eq 'right' ? 0 : -$y1,
#                        $y1),
#         $self->xy_to_n($self->{'align'} eq 'left' ? 0 : $y2,
#                        $y2));


#------------------------------------------------------------------------------
use constant tree_num_roots => 1;

sub tree_n_num_children {
  my ($self, $n) = @_;

  $n = $n - $self->{'n_start'};   # N=0 basis
  if ($n < 0) {
    return undef;
  }
  if ($n == 0 && $self->{'parts'} ne 'all') {
    # parts=left or parts=right have only 1 child under the root n=0
    return 1;
  }
  my ($depthbits, $ndepth) = _n0_to_depthbits($n, $self->{'parts'})
    or return 1;  # infinite

  unless (shift @$depthbits) {  # low bit
    # Depth even (incl zero), two children under every point.
    return 2;
  }

  # Depth odd, either 0 or 1 child.
  # If depth==1mod4 then 1-child.
  # If depth==3mod4 so two or more trailing 1-bits then some 0-child and
  # some 1-child.
  #
  $n -= $ndepth;  # Noffset into row
  my $repbit = _divrem_mutate($n,2); # low bit of $n
  while (shift @$depthbits) {               # bits of depth low to high
    if (_divrem_mutate($n,2) != $repbit) {  # bits of $n offset low to high
      return 0;
    }
  }
  return 1;
}

sub tree_n_children {
  my ($self, $n) = @_;
  ### tree_n_num_children(): $n

  $n = $n - $self->{'n_start'};   # N=0 basis
  if ($n < 0) {
    return;
  }
  if ($n == 0 && $self->{'parts'} ne 'all') {
    # parts=left or parts=right have only 1 child under the root n=0
    return ($n+1 + $self->{'n_start'});
  }
  my ($depthbits, $ndepth, $nwidth) = _n0_to_depthbits($n, $self->{'parts'})
    or return $n;  # infinite

  $n -= $ndepth;  # Noffset into row

  if (shift @$depthbits) {
    # Depth odd, single child under some or all points.
    # When depth==1mod4 it's all points, when depth has more than one
    # trailing 1-bit then it's only some points.
    while (shift @$depthbits) {  # depth==3mod4 or more low 1s
      my $repbit = _divrem_mutate($n,2);
      if (($n % 2) != $repbit) {
        return;
      }
    }
    return $n + $ndepth+$nwidth + $self->{'n_start'};

  } else {
    # Depth even (or zero), two children under every point.
    $n = 2*$n + $ndepth+$nwidth + $self->{'n_start'};
    return ($n,$n+1);
  }
}
sub tree_n_parent {
  my ($self, $n) = @_;

  my ($x,$y) = $self->n_to_xy($n)
    or return undef;

  if ($self->{'align'} eq 'diagonal') {
    my $n_parent = $self->xy_to_n($x-1, $y);
    if (defined $n_parent) {
      return $n_parent;
    } else {
      return $self->xy_to_n($x,$y-1);
    }
  }

  $y -= 1;
  my $n_parent = $self->xy_to_n($x-($self->{'align'} ne 'left'), $y);
  if (defined $n_parent) {
    return $n_parent;
  }
  return $self->xy_to_n($x+($self->{'align'} ne 'right'),$y);
}

sub tree_n_to_depth {
  my ($self, $n) = @_;
  ### SierpinskiTriangle n_to_depth(): $n
  $n = $n - $self->{'n_start'};
  if ($n < 0) {
    return undef;
  }
  my ($depthbits) = _n0_to_depthbits($n, $self->{'parts'})
    or return $n;  # infinite
  return digit_join_lowtohigh ($depthbits, 2, $n*0);
}
sub tree_depth_to_n {
  my ($self, $depth) = @_;
  return ($depth >= 0 ? _right_xy_to_n($self,0,$depth) : undef);
}

# sub _NOTWORKING__tree_depth_to_n_range {
#   my ($self, $depth) = @_;
#   if (is_infinite($depth)) {
#     return $depth;
#   }
#   if ($depth < 0) {
#     return undef;
#   }
#
#   my $zero = my $n = ($depth * 0);    # inherit bignum 0
#   my $width = my $npower = $zero+1;   # inherit bignum 1
#
#   foreach my $dbit (bit_split_lowtohigh($depth)) {
#     if ($dbit) {
#       $n = 2*$n + $npower;
#       $width *= 2;
#     }
#     $npower *= 3;
#   }
#   $n += $self->{'n_start'};
#
#   return ($n, $n+$width-1);
# }


#------------------------------------------------------------------------------
# In align=diagonal style, height is following a straight line X increment
# until hit bit in common with Y, meaning the end of Y low 0s.  Or follow
# straight line Y until hit bit in common with X, meaning end of X low 0s.
#
# If X,Y both even then X or Y lines are the same.
# If X odd then follow X to limit of Y low 0s.
# If Y odd then follow Y to limit of X low 0s.
#
#  | 65       ...
#  | 57 66
#  | 49    67
#  | 45 50 58 68
#  | 37          69
#  | 33 38       59 70
#  | 29    39    51    71
#  | 27 30 34 40 46 52 60 72
#  | 19                      73
#  |  |                       |
#  | 15-20                   61-74
#  |  |                       |
#  | 11    21                53    75
#  |  |     |                 |     |
#  |  9-12-16-22             47-54-62-76
#  |  |                       |
#  |  5          23          41          77
#  |  |           |           |           |
#  |  3--6       17-24       35-42       63-78
#  |  |           |           |           |
#  |  1     7    13    25    31    43    55    79
#  |  |     |     |     |     |     |     |     |
#  |  0--2--4--8-10-14-18-26-28-32-36-44-48-56-64-80
#  +-------------------------------------------------
#   X=0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15
#
# depthbits   1 0 0 0 1   Y of "right"
# nbits             n n
# xbits       n 0 0 0 n
# ybits      1-n     1-n  of Y-X for "diagonal"
#
# Y odd when ylow==1,nlow==0
#       follow its X low 0s by nbit==0 and invert of ybits==1
# X odd when ylow==1,nlow==1
#       follow its Y low 0s by nbit==1 and invert of xbits=nbits==1
#
# At a given depth<=2^k can go at most to its 2^k-1 limit, which means
# height = 2^k-1 - depth which is depth with bits flipped.
# Then bits of Noffset may put it in the middle of somewhere which limits
# the height to a sub-part 2^j < 2^k.
#
sub tree_n_to_subheight {
  my ($self, $n) = @_;
  ### SierpinskiTriangle tree_n_to_subheight(): $n

  $n = $n - $self->{'n_start'};
  if ($n < 0) {
    return undef;
  }
  my ($depthbits, $ndepth) = _n0_to_depthbits($n, $self->{'parts'})
    or return $n;  # infinite
  $n -= $ndepth;     # offset into row
  my @nbits = bit_split_lowtohigh($n);

  my $target = $nbits[0] || 0;
  foreach my $i (0 .. $#$depthbits) {
    unless ($depthbits->[$i] ^= 1) {  # flip 0<->1, at original==1 take nbit
      if ((shift @nbits || 0) != $target) {
        $#$depthbits = $i-1;
        return digit_join_lowtohigh($depthbits, 2, $n*0);
      }
    }
  }
  return undef; # first or last of row, infinite
}


#------------------------------------------------------------------------------
#   \                             /
#    4   0   0   0   0   0   0   4
#     \ /     \ /     \ /     \ /
#      1       1       1       1
#       \     /         \     /
#        2   2           2   2
#         \ /             \ /
#          3               3
#           \             /
#            4   0   0   4
#             \ /     \ /
#              1       1
#               \     /
#                2   2
#                 \ /
#                  3

# sub _EXPERIMENTAL__tree_n_to_leafdist {
#   my ($self, $n) = @_;
#   ### _EXPERIMENTAL__tree_n_to_leafdist() ...
#   my $d = $self->tree_n_to_depth($n);
#   if (defined $d) {
#     $d = 3 - ($d % 4);
#     if ($d == 0 && $self->tree_n_num_children($n) != 0) {
#       $d = 4;
#     }
#   }
#   return $d;
# }
sub _EXPERIMENTAL__tree_n_to_leafdist {
  my ($self, $n) = @_;
  ### _EXPERIMENTAL__tree_n_to_leafdist(): $n

  $n = $n - $self->{'n_start'};   # N=0 basis
  if ($n < 0) {
    return undef;
  }
  my ($depthbits, $ndepth) = _n0_to_depthbits($n, $self->{'parts'})  # low to high
    or return 1;  # infinite
  ### $depthbits

  # depth bits leafdist
  #   0    0,0    3
  #   1    0,1    2
  #   2    1,0    1
  #   3    1,1    0 or 4
  #
  my $ret = 3 - ((shift @$depthbits)||0);
  if (shift @$depthbits) { $ret -= 2; }
  ### $ret
  if ($ret) {
    return $ret;
  }

  $n -= $ndepth;
  ### Noffset into row: $n

  # Low bits of Nrem unchanging while trailing 1-bits in @depthbits,
  # to distinguish between leaf or non-leaf.  Same as tree_n_children().
  #
  my $repbit = _divrem_mutate($n,2); # low bit of $n
  ### $repbit
  do {
    ### next bit: $n%2
    if (_divrem_mutate($n,2) != $repbit) {  # bits of $n offset low to high
      return 0;  # is a leaf
    }
  } while (shift @$depthbits);
  return 4; # is a non-leaf
}

#------------------------------------------------------------------------------
# Return ($depthbits, $ndepth, $nwidth).
# $depthbits is an arrayref of bits low to high which are the tree depth of $n.
# $ndepth is first N of the row.
# $nwidth is the number of points in the row.
#
sub _n0_to_depthbits {
  my ($n, $parts) = @_;
  ### _n0_to_depthbits(): "$n  $parts"

  if (is_infinite($n)) {
    return;
  }
  if ($n == 0) {
    return ([], 0, 1);
  }

  my ($nwidth, $bitpos) = round_down_pow ($parts eq 'all' ? $n : 2*$n-1,
                                          3);
  ### $nwidth
  ### $bitpos

  my @depthbits;
  $depthbits[$bitpos] = 1;
  my $ndepth = ($parts eq 'all' ? $nwidth : ($nwidth + 1)/2);
  $nwidth *= 2;

  while (--$bitpos >= 0) {
    $nwidth /= 3;
    ### at: "n=$n nwidth=$nwidth bitpos=$bitpos depthbits=".join(',',map{$_||0}@depthbits)

    if ($n >= $ndepth + $nwidth) {
      $depthbits[$bitpos] = 1;
      $ndepth += $nwidth;
      $nwidth *= 2;
    } else {
      $depthbits[$bitpos] = 0;
    }
  }

  # Nwidth = 2**count1bits(depth), when parts=all
  ### @depthbits
  ### assert: $parts ne 'all' || $nwidth == (1 << scalar(grep{$_}@depthbits))

  return (\@depthbits, $ndepth, $nwidth);
}

#------------------------------------------------------------------------------
# levels

use Math::PlanePath::SierpinskiArrowheadCentres;
*level_to_n_range = \&Math::PlanePath::SierpinskiArrowheadCentres::level_to_n_range;
*n_to_level       = \&Math::PlanePath::SierpinskiArrowheadCentres::n_to_level;

#-----------------------------------------------------------------------------
1;
__END__