Math-PlanePath-Toothpick

 view release on metacpan or  search on metacpan

lib/Math/PlanePath/OneOfEight.pm  view on Meta::CPAN

    # $exp == _log2_floor($depth+1) so at $depth==2*$pow-1 one less.
    $n -= 3*$depth + 2*$exp + 6;

  } elsif ($parts eq '3side') {
    @pending = ($depth+1, $depth, $depth-1);
    @mult = (1, 3, 2);
    # Duplicated diagonals, and no log2_extras on two outermost octants.
    # For plain depth each log2 at depth=2^k-2, so another log2 decrease
    # when depth=2^k-1.
    # For depth+1 block each log2 at depth=2^k-2, so another log2 decrease
    # when depth=2^k-2.
    # $exp == _log2_floor($depth+1) so at $depth==2*$pow-1 one less.
    $n -= 3*$depth + 2*$exp + ($depth == $pow-1 ? 3 : 4);

  } elsif ($parts eq 'side') {
    unshift @pending, $depth+1;
    @mult = (1, 1);
    # $exp == _log2_floor($depth+1)
    $n -= $depth + 1 + $exp;
  }

  while ($exp >= 0 && @pending) {
    ### at: "pow=$pow exp=$exp  n=$n"
    ### assert: $pow == 2 ** $exp
    ### pending: join(',',@pending)
    ### mult: join(',',@mult)

    my @new_pending;
    my @new_mult;
    my $oct_pow;
    foreach my $depth (@pending) {
      my $mult = shift @mult;
      ### assert: $depth >= 0

      if ($depth <= 1) {
        ### small depth: "depth=$depth mult=$mult * $oct_to_n[$depth]"
        $n += $mult * $depth;  # oct=0 at depth=0, oct=1 at depth=1
        next;
      }
      my $rem = $depth - $pow;
      if ($rem < 0) {
        push @new_pending, $depth;
        push @new_mult, $mult;
        next;
      }

      ### $depth
      ### $mult
      ### $rem
      ### assert: $rem >= 0 && $rem < $pow

      my $powmult = $mult;
      if ($rem <= 1) {
        if ($rem == 0) {
          ### rem=0, oct(pow) only ...
        } else { # $rem == 1
          ### rem=1, oct(pow)+1 ...
          $n += $mult;
        }
      } else {
        ### formula ...
        # oct(pow+rem) = oct(pow)
        #                + oct(rem+1)
        #                + 2*oct(rem)
        #                - floor(log2(rem+1))
        #                - rem - 3

        my $rem1 = $rem + 1;
        {
          my ($lpow,$lexp) = round_down_pow ($rem1, 2);
          $n -= ($lexp + $rem + 3)*$mult;
          ### sub also: ($lexp + $rem + 3). " *mult=$mult"
        }
        if ($rem1 == $pow) {
          ### rem+1 == pow, increase powmult ...
          $powmult *= 2;    # oct(pow)+oct(rem+1) is 2*oct(pow)
        } elsif (@new_pending && $new_pending[-1] == $rem1) {
          ### merge into previously pushed new_pending[] ...
          # print "rem+1=$rem1 ",join(',',@new_pending),"\n";
          $new_mult[-1] += $mult;
        } else {
          ### push: "depth=$rem1 mult=$mult"
          push @new_pending, $rem1;
          push @new_mult, $mult;
        }

        ### push: "depth=$rem mult=".2*$mult
        push @new_pending, $rem;
        push @new_mult, 2*$mult;
      }

      # oct(pow) = (2*pow*pow + 3*exp + 7)/9 + pow/2
      #          = ((4*pow+9)*pow + 6*exp + 14)/18
      #
      $oct_pow ||= ((4*$pow+9)*$pow + 6*$exp + 14)/18;
      $n += $oct_pow * $powmult;
      ### oct(pow): "pow=$pow is $oct_pow * powmult=$powmult"
    }
    @pending = @new_pending;
    @mult = @new_mult;

    $exp--;
    $pow /= 2;
  }

  ### return: $n
  return $n;
}


# _depth_to_octant_added() returns the number of cells added at a given
# $depth level in parts=octant.  This is the same as
#     $added = tree_depth_to_n(depth+1) - tree_depth_to_n(depth)
#
# @$depth_aref is a list of depth values.
# @$mult_aref is the multiple of oct(depth) desired for each @depth_aref.
#
# On input @$depth_aref must have $depth_aref->[0] as the highest value.
#
# Within the code the depth list is mostly high to low and growing by one
# extra depth value at each $exp level.  But sometimes it grows a bit more

lib/Math/PlanePath/OneOfEight.pm  view on Meta::CPAN

  ### assert: scalar(@$depth_aref) == scalar(@$mult_aref)

  # $depth_aref->[0] must be the biggest depth, to make the $pow finding easy
  ### assert: scalar(@$depth_aref) >= 1
  ### assert: max(@$depth_aref) == $depth_aref->[0]

  my ($pow,$exp) = round_down_pow ($depth_aref->[0], 2);
  if (is_infinite($exp)) {
    return $exp;
  }
  ### $pow
  ### $exp

  my $added = $zero;

  # running $pow down to 2 (inclusive)
  while ($exp >= 0 && @$depth_aref) {
    ### at: "pow=$pow exp=$exp"
    ### assert: $pow == 2 ** $exp

    ### depth: join(',',@$depth_aref)
    ### mult: join(',',@$mult_aref)
    my @new_depth;
    my @new_mult;
    foreach my $depth (@$depth_aref) {
      my $mult = shift @$mult_aref;
      ### assert: $depth >= 0

      if ($depth <= $#_depth_to_octant_added) {
        ### small depth: "depth=$depth mult=$mult * $_depth_to_octant_added[$depth]"
        $added += $mult * $_depth_to_octant_added[$depth];
        next;
      }
      if ($depth < $pow) {
        push @new_depth, $depth;
        push @new_mult, $mult;
        next;
      }

      my $rem = $depth - $pow;

      ### $depth
      ### $mult
      ### $rem
      ### assert: $rem >= 0 && $rem < $pow

      if ($rem <= 1) {
        if ($rem == 0) {
          ### rem=0, grow 1 ...
          $added += $mult;
        } else {
          ### rem=1, grow 3 ...
          $added += 3 * $mult;
        }
      } else {
        my $rem1 = $rem + 1;
        if ($rem1 == $pow) {
          ### rem+1=pow, no lower part, 3/2 of pow ...
          $added += ($pow/2) * (3*$mult);
        } else {
          ### formula ...
          # oadd(pow+rem) = oadd(rem+1) + 2*oadd(rem)
          #                 + (is_pow2($rem+2) ? -2 : -1)

          # upper/lower diagonal overlap, and no log2_extras in lower
          $added -= (_is_pow2($rem+2) ? 2*$mult : $mult);

          if (@new_depth && $new_depth[-1] == $rem1) {
            ### merge into previously pushed new_depth ...
            # print "rem=$rem ",join(',',@new_depth),"\n";
            $new_mult[-1] += $mult;
          } else {
            ### push: "rem+1  depth=$rem1 mult=$mult"
            push @new_depth, $rem1;
            push @new_mult, $mult;
          }

          ### push: "rem    depth=$rem mult=".2*$mult
          push @new_depth, $rem;
          push @new_mult, 2*$mult;
        }
      }
    }
    $depth_aref = \@new_depth;
    $mult_aref = \@new_mult;

    $exp--;
    $pow /= 2;
  }

  ### return: $added
  return $added;
}


#------------------------------------------------------------------------------
# tree_n_to_subheight()

#use Smart::Comments;

{
  my %tree_n_to_subheight
    = do {
      my $depth0 = [ ]; # depth=0
      (wedge   => [ $depth0,
                    [ undef, 0 ], # depth=1
                  ],
       '3mid'  => [ $depth0,
                    [ undef, 0, undef, 0 ], # depth=1
                  ],
       '3side' => [ $depth0,
                    [ undef, 0, undef ],           # depth=1
                    [ 0, undef, undef, 0 ], # depth=2 N=4to8
                  ],
      )
    };

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

lib/Math/PlanePath/OneOfEight.pm  view on Meta::CPAN


parts=4,1,3mid have 5 children growing out of the 1-child of the X=2^k,Y=2^k
corner.  In an parts=octant, octant_up, and wedge there's only 3 children
around that point since that pattern doesn't go above the X=Y diagonal.

=back

=head2 Level Methods

=over

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

Return C<(0, tree_depth_to_n_end(2**$level - 1)>, or for parts=3side
C<tree_depth_to_n_end(2**$level)>.

parts=3side

=back

=head1 FORMULAS

=head2 Depth to N

The first point is N=0 so C<tree_depth_to_n($depth)> is the total number of
points up to and not including C<$depth>.  For the full pattern this
total(depth) follows a recurrence

    total(0)         = 0
    total(pow)       = (16*pow^2 + 24*exp - 7) / 9
    total(pow + rem) = total(pow) + 2*total(rem) + total(rem+1)
                         - 8*floor(log2(rem+1)) + 1
    where depth = pow + rem
      with pow=2^k the biggest power-of-2 <= depth
      and rem the remainder

For parts=octant the equivalent total points is

    oct(0)         = 0
    oct(pow)       = (4*pow^2 + 9*pow + 6*exp + 14) / 18
    oct(pow + rem) = oct(pow) + 2*oct(rem) + oct(rem+1)
                       - floor(log2(rem+1)) - rem - 3

The way this recurrence works can be seen from the self-similar pattern
described in L</One Octant> above.

    oct(pow)                # "base"
    + oct(rem)              # "extend"
    + oct(rem)              # "upper"
    + oct(rem+1)            # "lower"
    - floor(log2(rem+1))    # no pow2 points in lower
    - rem                   # unduplicate diagonal upper/lower
    - 3                     # unduplicate centre points

oct(rem)+oct(rem+1) of upper and lower would count their common diagonal
twice, hence "-rem" being the length of that diagonal.  The "centre" point
at X=pow,Y=pow is repeated by each of extend, upper, lower so "-2" to count
just once, and the X=pow+1,Y=pow point is repeated by extend and upper, so
"-1" to count it just once.

The 2*total(rem)+total(rem+1) in the formula is the same recurrence as the
toothpick pattern and the approach there can calculate it as a set of
pending depths and pow subtractions.  See
L<Math::PlanePath::ToothpickTree/Depth to N>.

The other patterns can be expressed as combinations of octants,

    parts=4 total   = 8*oct(n) - 4*n - 7
    parts=1 total   = 2*oct(n) - n
    3mid V2 total   = 2*oct(n+1) + 4*oct(n)
                        - 3n - 2*floor(log(n+1)) - 6
    3side V1 total  =   oct(n+1) + 3*oct(n) + 2*oct(n-1)
                        - 3n - floor(log(n+1)) - floor(log(n)) - 4

The depth offsets n,n+1,etc in these formulas become initial pending depth
for the toothpick style depth to N algorithm (and with respective initial
multipliers).

From the V1,V2 formulas it can be seen that V2(n)+V2(n+1) gives the same
combination of 1,3,2 times oct n-1,n,n+1 which is in V1, and that therefore
as noted in the Ndepth part of L</Three Side> above

    V1(n) = (V2(n) + V2(n-1) + 1) / 2

=head1 OEIS

This cellular automaton is in Sloane's Online Encyclopedia of Integer
Sequences as

=over

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

=back

    parts=4 (the default)
      A151725   total cells "V", tree_depth_to_n()
      A151726   added cells "v"

    parts=1
      A151735   total cells, tree_depth_to_n()
      A151737   added cells

    parts=3mid
      A170880   total cells, tree_depth_to_n()
      A151728   added cells "v2"
      A151727   added cells "v2" * 4
      A151729   (added cells - 1) / 2

    parts=3side
      A170879   total cells, tree_depth_to_n()
      A151747   added cells "v1"

=head1 SEE ALSO

L<Math::PlanePath>,
L<Math::PlanePath::ToothpickTree>,
L<Math::PlanePath::UlamWarburton>

=head1 HOME PAGE

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

=head1 LICENSE

Copyright 2012, 2013, 2014, 2015 Kevin Ryde

This file is part of Math-PlanePath-Toothpick.

Math-PlanePath-Toothpick 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-Toothpick 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.



( run in 0.532 second using v1.01-cache-2.11-cpan-39bf76dae61 )