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


# 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 = 117;
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'],
    },
    Math::PlanePath::Base::Generic::parameter_info_nstart0(),
  ];

my %x_negative = (triangular => 1,
                  left       => 1,
                  right      => 0,
                  diagonal   => 0);
sub x_negative {
  my ($self) = @_;
  return $x_negative{$self->{'align'}};
}
sub x_negative_at_n {
  my ($self) = @_;
  return ($self->{'align'} eq 'triangular' || $self->{'align'} eq 'left'
          ? $self->n_start + 1
          : undef);
}
use constant class_y_negative => 0;
use constant default_n_start => 0;
use constant n_frac_discontinuity => .5;
use constant tree_num_children_list => (0,1,2);

sub x_maximum {
  my ($self) = @_;
  return ($self->{'align'} eq 'left'
          ? 0       # left all X<=0
          : undef); # others X unbounded
}
use constant sumxy_minimum => 0;  # triangular X>=-Y or all X>=0

# 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
}


#------------------------------------------------------------------------------
sub new {
  my $self = shift->SUPER::new(@_);
  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }
  my $align = ($self->{'align'} ||= 'triangular');
  if (! exists $x_negative{$align}) {
    croak "Unrecognised align option: ", $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)
    or return ($n,$n); # infinite

  ### $depthbits
  ### $ndepth

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

  # 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 $zero = $n * 0;
  my $x = digit_join_lowtohigh (\@xbits,    2, $zero);
  my $y = digit_join_lowtohigh ($depthbits, 2, $zero);

  ### final: "$x,$y"
  if ($self->{'align'} eq 'right') {
    return ($x, $y);
  } elsif ($self->{'align'} eq 'left') {
    return ($x-$y, $y);
  } elsif ($self->{'align'} eq 'diagonal') {
    return ($x, $y-$x);
  } else { # triangular
    return (-$y+2*$x, $y);
  }
}

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

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

  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() ...

  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;
  }
  my ($depthbits, $ndepth) = _n0_to_depthbits($n)
    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) = @_;

  $n = $n - $self->{'n_start'};   # N=0 basis
  if ($n < 0) {
    return;
  }
  my ($depthbits, $ndepth, $nwidth) = _n0_to_depthbits($n)
    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)
    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)
    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)  # 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 which give 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) = @_;

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

  my ($nwidth, $bitpos) = round_down_pow ($n, 3);
  ### $nwidth
  ### $bitpos

  my @depthbits;
  my $ndepth = 0;
  for (;;) {
    ### 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;
    }
    $bitpos--;
    last unless $bitpos >= 0;
    $nwidth /= 3;
  }

  # Nwidth = 2**count1bits(depth)
  ### @depthbits
  ### assert: $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__

=for stopwords eg Ryde Sierpinski Nlevel ie Ymin Ymax OEIS Online rowpoints Nleft Math-PlanePath Gould's Nend bitand Noffset Ndepth Nrem Dyck

=head1 NAME

Math::PlanePath::SierpinskiTriangle -- self-similar triangular path traversal

=head1 SYNOPSIS

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

=head1 DESCRIPTION

X<Sierpinski, Waclaw>This is an integer version of the Sierpinski triangle
with cells numbered horizontally across each row.

    65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80   15
      57      58      59      60      61      62      63      64     14
        49  50          51  52          53  54          55  56       13
          45              46              47              48         12
            37  38  39  40                  41  42  43  44           11
              33      34                      35      36             10
                29  30                          31  32                9
                  27                              28                  8
                    19  20  21  22  23  24  25  26                    7
                      15      16      17      18                      6
                        11  12          13  14                        5
                           9              10                          4
                             5   6   7   8                            3
                               3       4                              2
                                 1   2                                1
                                   0                             <- Y=0

         X= ... -9-8-7-6-5-4-3-2-1 0 1 2 3 4 5 6 7 8 9 ...

The base figure is the first two rows shape N=0 to N=2.  Notice the middle
"." position X=0,Y=1 is skipped

    1  .  2
       0

This is replicated twice in the next row pair as N=3 to N=8.  Then the
resulting four-row shape is replicated twice again in the next four-row
group as N=9 to N=26, etc.

See the C<SierpinskiArrowheadCentres> path to traverse by a connected
sequence rather than rows jumping across gaps.

=head2 Row Ranges

The number of points in each row is always a power of 2.  The power is the
count of 1-bits in Y.  (This count is sometimes called Gould's sequence.)

    rowpoints(Y) = 2^count_1_bits(Y)

For example Y=13 is binary 1101 which has three 1-bits so in row Y=13 there
are 2^3=8 points.

Because the first point is N=0, the N at the left of each row is the
cumulative count of preceding points,

    Ndepth(Y) = rowpoints(0) + ... + rowpoints(Y-1)

Since the powers of 2 are always even except for 2^0=1 in row Y=0, this
Ndepth(Y) total is always odd.  The self-similar nature of the triangle
means the same is true of the sub-triangles, for example odd
N=31,35,41,47,etc on the left of the triangle at X=8,Y=8.  This means in
particular the primes (being odd) fall predominately on the left side of the
triangles and sub-triangles.

=head2 Replication Sizes

Counting the single point N=0 as level=0, then N=0,1,2 as level 1, each
replication level goes from

    Nstart = 0
    Nlevel = 3^level - 1     inclusive

For example level 2 is from N=0 to N=3^2-1=8.  Each level doubles in size,

               0  <= Y <= 2^level - 1
    - 2^level + 1 <= X <= 2^level - 1

=head2 Align Right

Optional C<align=E<gt>"right"> puts points to the right of the Y axis,
packed next to each other and so using an eighth of the plane.

=cut

# math-image --path=SierpinskiTriangle,align=right --all --output=numbers

=pod

    align => "right"

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

=head2 Align Left

Optional C<align=E<gt>"left"> puts points to the left of the Y axis,
ie. into negative X.  The rows are still numbered starting from the left, so
it's a shift across, not a negate of X.

=cut

# math-image --path=SierpinskiTriangle,align=left --all --output=numbers

=pod

    align => "left"

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

=head2 Align Diagonal

Optional C<align=E<gt>"diagonal"> puts rows on diagonals down from the Y
axis to the X axis.  This uses the whole of the first quadrant, with gaps
according to the pattern.

=cut

# math-image --expression='i<=80?i:0' --path=SierpinskiTriangle,align=diagonal --output=numbers

=pod

    align => "diagonal"

     15 | 65       ...
     14 | 57 66
     13 | 49    67
     12 | 45 50 58 68
     11 | 37          69
     10 | 33 38       59 70
      9 | 29    39    51    71
      8 | 27 30 34 40 46 52 60 72
      7 | 19                      73
      6 | 15 20                   61 74
      5 | 11    21                53    75
      4 |  9 12 16 22             47 54 62 76
      3 |  5          23          41          77       ...
      2 |  3  6       17 24       35 42       63 78
      1 |  1     7    13    25    31    43    55    79
    Y=0 |  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

This form visits all points X,Y where X and Y written in binary have no
1-bits in the same bit positions, ie. where S<X bitand Y> == 0.  For example
X=13,Y=3 is not visited because 13="1011" and 6="0110" both have bit "0010"
set.

This bit-and rule is an easy way to test for visited or not visited cells of
the pattern.  The visited cells can be calculated by this diagonal X,Y
bitand, but then plotted X,X+Y for the "right" align or X-Y,X+Y for
"triangular".

=head2 Cellular Automaton

The triangle arises in Stephen Wolfram's 1-D cellular automatons (per
L<Math::PlanePath::CellularRule> and L<Cellular::Automata::Wolfram>).

    align           rule
    -----           ----
    "triangular"    18,26,82,90,146,154,210,218
    "right"         60
    "left"          102

=over

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

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

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

=back

=cut

# rule 60 right hand octant
# rule 102 left hand octant
# math-image --path=CellularRule,rule=60 --all
# math-image --path=CellularRule,rule=102 --all

=pod

In each row the rule 18 etc pattern turns a cell "on" in the next row if one
but not both its diagonal predecessors are "on".  This is a mod 2 sum giving
Pascal's triangle mod 2.

Some other cellular rules are variations on the triangle,

=over

=item *

Rule 22 is "triangular" but filling the gap between leaf points such as N=5
and N=6.

=item *

Rule 126 adds an extra point on the inward side of each visited.

=item *

Rule 182 fills in the big gaps leaving just a single-cell
empty border delimiting them.

=back

=head2 N Start

The default is to number points starting N=0 as shown above.  An optional
C<n_start> parameter can give a different start, with the same shape.  For
example starting at 1, which is the numbering of C<CellularRule> rule=60,

=cut

# math-image --path=SierpinskiTriangle,n_start=1 --expression='i<=27?i:0' --output=numbers

=pod

    n_start => 1

    20    21    22    23    24    25    26    27
       16          17          18          19
          12    13                14    15
             10                      11
                 6     7     8     9
                    4           5
                       2     3
                          1

=head1 FUNCTIONS

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

=over 4

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

=item C<$path = Math::PlanePath::SierpinskiTriangle-E<gt>new (align =E<gt> $str, n_start =E<gt> $n)>

Create and return a new path object.  C<align> is a string, one of the
following as described above.

    "triangular"    (the default)
    "right"
    "left"
    "diagonal"

=back

=head2 Descriptive Methods

=over

=item C<$n = $path-E<gt>n_start()>

Return the first N in the path.  This is 0 by default, or the given
C<n_start> parameter.

=back

=head2 Tree Methods

=over

=item C<@n_children = $path-E<gt>tree_n_children($n)>

Return the children of C<$n>, or an empty list if C<$n E<lt> n_start>
(ie. before the start of the path).

The children are the points diagonally up left and right on the next row
(Y+1).  There can be 0, 1 or 2 such points.  At even depth there's 2, on
depth=1mod4 there's 1.  On depth=3mod4 there's some 0s and some 1s.  See
L</N to Number of Children> below.

For example N=3 has two children N=5,N=6.  Then in turn N=5 has just one
child N=9 and N=6 has no children.  The way points are numbered across a row
means that when there's two children they're consecutive N values.

=item C<$n_parent = $path-E<gt>tree_n_parent($n)>

Return the parent node of C<$n>, or C<undef> if C<$n E<lt>= n_start> (the
top of the triangle).

=item C<$depth = $path-E<gt>tree_n_to_depth($n)>

Return the depth of node C<$n>, or C<undef> if there's no point C<$n>.  In
the "triangular", "right" and "left" alignments this is the same as the Y
coordinate from C<n_to_xy()>.  In the "diagonal" alignment it's X+Y.

=item C<$n = $path-E<gt>tree_depth_to_n($depth)>

=item C<$n = $path-E<gt>tree_depth_to_n_end($depth)>

Return the first or last N at tree level C<$depth>.  The start of the tree
is depth=0 at the origin X=0,Y=0.

This is the N at the left end of each row.  So in the default triangular
alignment it's the same as C<xy_to_n(-$depth,$depth)>.

=back

=head2 Level Methods

=over

=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>

Return C<(0, 3**$level - 1)>.

=back

=head1 FORMULAS

=head2 X,Y to N

For calculation it's convenient to turn the X,Y coordinates into the "right"
alignment style, so that Y is the depth and X is in the range
0E<lt>=XE<lt>=Y.

The starting position of each row of the triangle is given by turning 1-bits
of the depth into powers-of-3.

    Y = depth = 2^a + 2^b + 2^c + 2^d ...       a>b>c>d...

    Ndepth = first N at this depth
           =         3^a
             +   2 * 3^b
             + 2^2 * 3^c
             + 2^3 * 3^d
             + ...

For example depth=6=2^2+2^1 starts at Ndepth=3^2+2*3^1=15.  The powers-of-3
are the three parts of the triangle replication.  The power-of-2 doubling is
the doubling of the row Y when replicated.

Then the bits of X at the positions of the 1-bits of the depth become an N
offset into the row.

               a  b  c  d
    depth    = 10010010010     binary
    X        = m00n00p00q0
    Noffset  =        mnpq     binary

    N = Ndepth + Noffset

For example in depth=6 binary "110" then at X=4="100" take the bits of X
where depth has 1-bits, which is X="10_" so Noffset="10" binary and
N=15+2=17, as per the "right" table above at X=4,Y=6.

If X has any 1-bits which are a 0-bits in the depth depth then that X,Y is
not visited.  For example if depth=6="110" then X=3="11" is not visited
because the low bit X="__1" has depth="__0" at that position.

=head2 N to Depth

The row containing N can be found by working down the Ndepth formula shown
above.  The "a" term is the highest 3^a E<lt>= N, thus giving a bit 2^a for
the depth.  Then for the remaining Nrem = N - 3^a find the highest "b" where
2*3^b E<lt>= Nrem.  And so on until reaching an Nrem which is too small to
subtract any more terms.

It's convenient to go by bits high to low of the prospective depth, deciding
at each bit whether Nrem is big enough to give the depth a 1-bit there, or
whether it must be a 0-bit.

    a = floor(log3(N))     round down to power-of-3
    pow = 3^a
    Nrem = N - pow

    depth = high 1-bit at bit position "a" (counting from 0)

    factor = 2
    loop bitpos a-1 down to 0
      pow /= 3
      if pow*factor <= Nrem
      then depth 0-bit, factor *= 2
      else depth 1-bit

    factor is 2^count1bits(depth)
    Noffset = Nrem     offset into row
    0 <= Noffset < factor

=head2 N to X,Y

N is turned into depth and Noffset as per above.  X in "right" alignment
style is formed by spreading the bits of Noffset out according to the 1-bits
of the depth.

    depth   = 100110  binary
    Noffset =    abc  binary
    Xright  = a00bc0

For example in depth=5 this spreads an Noffset=0to3 to make X=000, 001, 100,
101 in binary and in "right" alignment style.

From an X,Y in "right" alignment the other alignments are formed

    alignment   from "right" X,Y
    ---------   ----------------
    triangular     2*X-Y, Y       so -Y <= X < Y
    right          X,     Y       unchanged
    left           X-Y,   Y       so -Y <= X <= 0
    diagonal       X,   Y-X       downwards sloping

=head2 N to Number of Children

The number of children follows a pattern based on the depth.

    depth      number of children
    -----      ------------------

     12    2       2       2       2   
     11     1 0 0 1         1 0 0 1
     10      2   2           2   2
      9       1 1             1 1
      8        2               2
      7         1 0 0 0 0 0 0 1   
      6          2   2   2   2 
      5           1 1     1 1  
      4            2       2   
      3             1 0 0 1   
      2              2   2
      1               1 1
      0                2   

If depth is even then all points have 2 children.  For example row depth=6
has 4 points and all have 2 children each.

At odd depth the number of children is either 1 or 0 according to the
Noffset position in the row masked down by the trailing 1-bits of the depth.

    depth  = ...011111 in binary, its trailing 1s

    Noffset = ...00000   \ num children = 1
            = ...11111   /
            = ...other   num children = 0

For example depth=11 is binary "1011" which has low 1-bits "11".  If those
two low bits of Noffset are "00" or "11" then 1 child.  Any other bit
pattern in Noffset ("01" or "10" in this case) is 0 children.  Hence the
pattern 1,0,0,1,1,0,0,1 reading across the depth=11 row.

In general when the depth doubles the triangle is replicated twice and the
number of children is carried with the replications, except the middle two
points are 0 children.  For example the triangle of depth=0to3 is repeated
twice to make depth=4to7, but the depth=7 row is not children 10011001 of a
plain doubling from the depth=3 row, but instead 10000001 which is the
middle two points becoming 0.

=head2 N to Number of Siblings

The number of siblings of a given node is determined by its depth,

    depth      number of siblings
    -----      ------------------

      4            0       0   
      3             1 1 1 1   
      2              0   0
      1               1 1
      0                0   

    depth     number of siblings
    -----     ------------------
     odd             1
     even            0

In an even row the points are all spread apart so there are no siblings.
The points in such a row are cousins or second cousins, etc, but none share
a parent.

In an odd row each parent node (an even row) has 2 children and so each of
those points has 1 sibling.

The effect is to conflate the NumChildren=1 and NumChildren=0 cases in the
picture above, those two becoming a single sibling.

    num children of N      num siblings of N
    -----------------      -----------------
          0 or 1                   1
            2                      0

=head2 Rectangle to N Range

An easy range can be had just from the Y range by noting the diagonals X=Y
and X=-Y are always visited, so just take the Ndepth of Ymin and Nend of
Ymax,

    # align="triangular"
    Nmin = N at X=-Ymin,Y=Ymin
    Nmax = N at X=Ymax,Y=Ymax

Or in "right" style the left end is at X=0 instead,

    # align="right"
    Nmin = N at X=0,Ymin
    Nmax = N at Ymax,Ymax

For less work but a bigger over-estimate, invert the Nlevel formulas given
in L</Row Ranges> above to round up to the end of a depth=2^k replication,

    level = floor(log2(Ymax)) + 1
    Nmax = 3^level - 1

For example Y=11, level=floor(log2(11))+1=4, so Nmax=3^4-1=80, which is the
end of the Y=15 row, ie. rounded up to the top of the replication block Y=8
to Y=15.

=head1 OEIS

The Sierpinski triangle is in Sloane's Online Encyclopedia of Integer
Sequences in various forms,

=over

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

=back

    A001316   number of cells in each row (Gould's sequence)
    A001317   rows encoded as numbers with bits 0,1
    A006046   total cells to depth, being tree_depth_to_n(), 
    A074330   Nend, right hand end of each row (starting Y=1)

A001316 is the "rowpoints" described above.  A006046 is the cumulative total
of that sequence which is the "Ndepth", and A074330 is 1 less for "Nend".

    align="triangular" (the default)
      A047999   0,1 cells by rows
      A106344   0,1 cells by upwards sloping dX=3,dY=1
      A130047   0,1 cells of half X<=0 by rows

A047999 etc is every second point in the default triangular lattice, or all
points in align="right" or "left".

    align="triangular" (the default)
      A002487   count points along dX=3,dY=1 slopes
                  is the Stern diatomic sequence
      A106345   count points along dX=5,dY=1 slopes

dX=3,dY=1 sloping lines are equivalent to opposite-diagonals dX=-1,dY=1 in
align="right".

    align="right"
      A075438   0,1 cells by rows including 0 blanks at left of pyramid

    align="right", n_start=0
      A006046   N on Y axis, being Ndepth
      A074330   N on Diagonal starting from Y=1, being Nend
    align="left", n_start=0
      A006046   N on NW diagonal, being Ndepth
      A074330   N on Y axis starting from Y=1, being Nend

    A080263   Dyck encoding of the tree structure
    A080264     same in binary
    A080265     position in list of all balanced binary

    A080268   Dyck encoding breadth-first
    A080269     same in binary
    A080270     position in list of all balanced binary

    A080318   Dyck encoding breadth-first of branch-reduced
                (duplicate each bit)
    A080319     same in binary
    A080320     position in list of all balanced binary

For the Dyck encoding see for example L<Math::NumSeq::BalancedBinary/Binary
Trees>.  The position in all balanced binary which is A080265 etc
corresponds to C<value_to_i()> in that C<NumSeq>.

A branch-reduced tree has any single-child node collapsed out, so that all
remaining nodes are either a leaf node or have 2 (or more) children.  The
effect of this on the Sierpinski triangle in breadth-first encoding is to
duplicate each bit, so A080269 with each bit repeated gives the
branch-reduced A080319.

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::SierpinskiArrowhead>,
L<Math::PlanePath::SierpinskiArrowheadCentres>,
L<Math::PlanePath::CellularRule>,
L<Math::PlanePath::ToothpickUpist>

L<Math::NumSeq::SternDiatomic>,
L<Math::NumSeq::BalancedBinary>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2011, 2012, 2013, 2014 Kevin Ryde

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