File Coverage

blib/lib/Math/PlanePath/LCornerTree.pm

Criterion	Covered	Total	%
statement	113	487	23.2
branch	61	342	17.8
condition	16	136	11.7
subroutine	14	33	42.4
pod	21	21	100.0
total	225	1019	22.0

# 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.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath-Toothpick.  If not, see .


# parts=2 row total X = 3^count1bits(depth) = A048883
# parts=wedge row total X = (3^count1bits(depth) + 1) / 2


# A147562 U-W
# A160410
# A151920

#------------------------------------------------------------------------------
# cf A183126 triplet around each exposed end, starting from length=1 two ends
#    A183127 added  0,1,6,16,16,40,16,40,40,112,16
#    http://www.math.vt.edu/people/layman/sequences/A183126.htm
# added = 4 * 3^count1bits(depth) or thereabouts
# two diagonal halves
#
#------------------------------------------------------------------------------
# wedge odd/even

#     1  1
#  2  0  1
#  2  2

#  3  3        3  3
#  3  2  2  2  2  3
#  3  2  1  1  2
#  3  3  0  1  2
#        3  2  2  3
#        3  3  3  3

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

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


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

use vars '$VERSION', '@ISA';
$VERSION = 17;
use Math::PlanePath;
@ISA = ('Math::PlanePath');


use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits 119  # v.119 for round_up_pow()
  'round_up_pow',
  'round_down_pow',
  'bit_split_lowtohigh',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';

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


use constant n_start => 0;
use constant parameter_info_array =>
  [ { name            => 'parts',
      share_key       => 'parts_lcornertree',
      display         => 'Parts',
      type            => 'enum',
      default         => '4',
      choices         => ['4','3','2','1',
                          'octant','octant+1','octant_up','octant_up+1',
                          'wedge','wedge+1','diagonal','diagonal-1',
                         ],
      choices_display => ['4','3','2','1',
                          'Octant','Octant+1','Octant Up','Octant Up+1',
                          'Wedge','Wedge+1','Diagonal','Diagonal 1',
                         ],
      description     => 'Which parts of the plane to fill.',
    },
  ];

{
  my %x_negative = (4             => 1,
                    3             => 1,
                    2             => 1,
                    1             => 0,
                    octant        => 0,
                    'octant+1'    => 0,
                    octant_up     => 0,
                    'octant_up+1' => 0,
                    wedge         => 1,
                    'wedge+1'     => 1,
                    diagonal      => 1,
                    'diagonal-1'  => 1,
                   );
  sub x_negative {
    my ($self) = @_;
    return $x_negative{$self->{'parts'}};
  }
}
{
  my %y_negative = (4             => 1,
                    3             => 1,
                    2             => 0,
                    1             => 0,
                    octant        => 0,
                    'octant+1'    => 0,
                    octant_up     => 0,
                    'octant_up+1' => 0,
                    wedge         => 0,
                    'wedge+1'     => 0,
                    diagonal      => 1,
                    'diagonal-1'  => 1,
                   );
  sub y_negative {
    my ($self) = @_;
    return $y_negative{$self->{'parts'}};
  }
}

{
  my %x_negative_at_n = (4             => 1,
                         3             => 2,
                         2             => 1,
                         1             => undef,
                         octant        => undef,
                         'octant+1'    => undef,
                         octant_up     => undef,
                         'octant_up+1' => undef,
                         wedge         => 1,
                         'wedge+1'     => 1,
                         diagonal      => 2,
                         'diagonal-1'  => 11,
                        );
  sub x_negative_at_n {
    my ($self) = @_;
    return $x_negative_at_n{$self->{'parts'}};
  }
}
{
  my %y_negative_at_n = (4             => 2,
                         3             => 1,
                         2             => undef,
                         1             => undef,
                         octant        => undef,
                         'octant+1'    => undef,
                         octant_up     => undef,
                         'octant_up+1' => undef,
                         wedge         => undef,
                         'wedge+1'     => undef,
                         diagonal      => 0,
                         'diagonal-1'  => 4,
                        );
  sub y_negative_at_n {
    my ($self) = @_;
    return $y_negative_at_n{$self->{'parts'}};
  }
}

{
  my %sumxy_minimum = (1             => 0,
                       octant        => 0,
                       'octant+1'    => 0,
                       octant_up     => 0,
                       'octant_up+1' => 0,
                       wedge         => -1,  # X>=-Y-1 so X+Y>=-1
                       'wedge+1'     => -2,  # X>=-Y-2 so X+Y>=-2
                       diagonal      => -1,
                       'diagonal-1'  => 0,   # X>=-Y
                      );
  sub sumxy_minimum {
    my ($self) = @_;
    return $sumxy_minimum{$self->{'parts'}};
  }
}
{
  my %diffxy_minimum = (octant     => 0,   # octant X>=Y so X-Y>=0
                        'octant+1' => -1,  # octant X>=Y-1 so X-Y>=-1
                       );
  sub diffxy_minimum {
    my ($self) = @_;
    return $diffxy_minimum{$self->{'parts'}};
  }
}
{
  my %diffxy_maximum = (octant_up     => 0,  # octant_up   X<=Y so X-Y<=0
                        'octant_up+1' => 1,  # octant_up+1 X<=Y+1 so X-Y<=1
                        wedge         => 0,  # wedge X<=Y so X-Y<=0
                       'wedge+1'      => 1,  # wedge+1 X>=Y+1 so X-Y>=1
                       );
  sub diffxy_maximum {
    my ($self) = @_;
    return $diffxy_maximum{$self->{'parts'}};
  }
}

# parts=1 Dir4 max 12,-11
#                 121,-110
#                 303,-213
#                1212,-1031
#               12121,-10310      -> 12,-10
#               30303,-21213      -> 3,-2
# parts=2  dX=big,dY=-1
# parts=3  dX=big,dY=-1
# parts=4  dx=0,dy=-1 at N=1
{
  my %dir_maximum_dxdy = (1             => [3,-2],  # supremum
                          2             => [0,0],   # supremum
                          3             => [0,0],   # supremum
                          4             => [0,-1],  # N=1  dX=0,dY=-1 South
                          octant        => [0,-2],  # N=4  dX=0,dY=-2 South
                          'octant+1'    => [1,-2],  # N=6  dX=1,dY=-2 SSE
                          octant_up     => [0,-1],  # N=8  dX=0,dY=-1 South
                          'octant_up+1' => [0,-1],  # N=11 dX=0,dY=-1 South
                          wedge         => [0,-1],  # N=13 dX=0,dY=-1 South
                          'wedge+1'     => [0,-1],  # N=6  dX=0,dY=-1 South
                          diagonal      => [2,-2],  # N=2  South-East
                          'diagonal-1'  => [3,-3],  # N=12 South-East
                         );
  sub dir_maximum_dxdy {
    my ($self) = @_;
    return @{$dir_maximum_dxdy{$self->{'parts'}}};
  }
}

{
  my %any_num_children_2
    = (octant        => 1,
       octant_up     => 1,
       wedge         => 1,
       diagonal      => 1,
      );
  sub tree_num_children_list {
    my ($self) = @_;
    return (0,
            ($any_num_children_2{$self->{'parts'}} ? (2) : ()),
            3);
  }
}

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

# how many toplevel root nodes in the tree of given $parts
my %parts_to_numroots = (4             => 4,
                         3             => 3,
                         2             => 2,
                         1             => 1,
                         octant        => 1,
                         'octant+1'    => 1,
                         octant_up     => 1,
                         'octant_up+1' => 1,
                         wedge         => 2,
                         'wedge+1'     => 2,
                         diagonal      => 3,
                         'diagonal-1'  => 1,
                        );
sub tree_num_roots {
  my ($self) = @_;
  return $parts_to_numroots{$self->{'parts'}};
}

sub new {
  my $self = shift->SUPER::new(@_);
  my $parts = ($self->{'parts'} ||= 4);
  if (! exists $parts_to_numroots{$parts}) {
    croak "Unrecognised parts: ",$parts;
  }
  return $self;
}

my @next_state = (0,12,0,4, 4,0,4,8, 8,4,8,12, 12,8,12,0);
my @digit_to_x = (0,1,1,0, 1,1,0,0, 1,0,0,1, 0,0,1,1);
my @digit_to_y = (0,0,1,1, 0,1,1,0, 1,1,0,0, 1,0,0,1);

my @diagonal1_n_to_xy = ([0,0],
                         [1,0],
                         [1,1],
                         [0,1]);

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

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$n); }
  {
    my $int = int($n);
    ### $int
    ### $n
    if ($n != $int) {
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      my $frac = $n - $int;  # inherit possible BigFloat
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
    $n = $int;       # BigFloat int() gives BigInt, use that
  }

  if (is_infinite($n)) { return ($n,$n); }
  my $zero = ($n * 0); # inherit bignum 0
  my $parts = $self->{'parts'};

  if ($parts eq 'diagonal-1') {
    if ($n <= 3) {
      return @{$diagonal1_n_to_xy[$n]};
    }
  }

  my ($depthbits, $ndepth, $nwidth) = _n0_to_depthbits($n, $parts);

  ### $n
  ### $ndepth
  ### $nwidth
  ### $parts
  ### $depthbits
  ### assert: $nwidth == $self->tree_depth_to_n(1+digit_join_lowtohigh($depthbits,2,$n*0)) - $self->tree_depth_to_n(digit_join_lowtohigh($depthbits,2,$n*0))

  $n -= $ndepth;
  ### N remainder offset into row: $n
  ### assert: $n >= 0
  ### assert: $n < $nwidth

  my $quad;
  my $x = 0;
  my $y = 0;
  if ($parts eq 'wedge') {
    ### assert: $nwidth % 2 == 0
    my $noct = $nwidth/2;
    if ($n < $noct) {
      $n += $noct - 1;
      $quad = 0;
    } else {
      $n -= $noct;
      $quad = 1;
    }
  } elsif ($parts eq 'wedge+1') {
    ### assert: $nwidth % 2 == 0
    my $noct = $nwidth/2;
    if ($n < $noct) {
      ### first half, add to N: $noct - 3
      $n += $noct - 3;
      $quad = 0;
    } else {
      ### second half ...
      $n -= $noct;
      $quad = 1;
    }

  } elsif ($parts eq 'diagonal') {
    ### assert: ($nwidth+1)%4 == 0
    my $noct = ($nwidth+1)/4;
    ### $noct
    if ($n < $noct) {
      ### first oct is quad=3 ...
      $n += $noct - 1;
      $quad = 3;
    } elsif ($n >= $nwidth - $noct) {
      ### last oct is quad=1 ...
      $n -= ($nwidth - $noct);
      $quad = 1;
    } else {
      $n -= $noct;
      $quad = 0;
    }

  } elsif ($parts eq 'diagonal-1') {
    ### assert: ($nwidth+3) % 4 == 0
    my $noct = ($nwidth+3)/4;
    ### $noct
    if ($n < $noct) {
      ### first oct is quad=3 ...
      $n += $noct - 3;
      $quad = 3;
      $x = -1;
      $y = 1;
    } elsif ($n >= $nwidth - $noct) {
      ### last oct is quad=1 ...
      $n -= ($nwidth - $noct);
      $quad = 1;
      $x = 1;
      $y = -1;
    } else {
      $n -= $noct;
      $quad = 0;
      $x = 1;
      $y = 1;
    }

  } elsif ((my $numroots = $parts_to_numroots{$parts}) > 1) {
    # like a mixed-radix high digit radix $numroots then rest radix 3
    $nwidth /= $numroots;
    ($quad, $n) = _divrem($n,$nwidth);
    ### $quad
    ### assert: $quad >= 0
    ### assert: $quad < $numroots
    if ($parts eq '3') {
      if ($quad == 1) { $quad = 3; }   # quad=1 -> 3
      if ($quad == 2) { $quad = 1; }   # quad=2 -> 1
    }
  } else {
    $quad = 0;
    if ($parts eq 'octant_up') {
      $n += $nwidth - 1;
    } elsif ($parts eq 'octant_up+1') {
      $n += $nwidth - 3;
    }
  }
  ### $quad
  ### $n

  my @nternary = digit_split_lowtohigh($n, 3);
  ### @nternary

  # Ternary digits for triple parts of Noffset mapped out to base4 digits in
  # the style of LCornerReplicate.
  # Where there's a 0-bit in the depth is a 0-digit for Nbase4.
  # Where there's a 1-bit in the depth takes a ternary+1 for Nbase4.
  # Small Noffset has less trits than the depth 1s, hence "nternary || 0".
  #
  my @nbase4 = map {$_ && (1 + (shift @nternary || 0))} @$depthbits;
  ### @nbase4

  my $state = 0;
  my (@xbits, @ybits);
  foreach my $i (reverse 0 .. $#nbase4) {    # digits high to low
    $state += $nbase4[$i];
    $xbits[$i] = $digit_to_x[$state];
    $ybits[$i] = $digit_to_y[$state];
    $state = $next_state[$state];
  }

  ### xbits join: digit_join_lowtohigh (\@xbits, 2, $zero)
  ### ybits join: digit_join_lowtohigh (\@ybits, 2, $zero)
  $x += digit_join_lowtohigh (\@xbits, 2, $zero);
  $y += digit_join_lowtohigh (\@ybits, 2, $zero);

  if ($quad & 1) {
    ($x,$y) = (-1-$y,$x); # rotate +90
  }
  if ($quad & 2) {
    $x = -1-$x; # rotate +180
    $y = -1-$y;
  }
  ### final: "$x,$y"
  return $x,$y;
}

# my @next_state = (0, 1, 3, 2,
# my @yx_to_digit = (0, 1, 3, 2,
#                    0, 1, 3, 2,     # rot +90
#                   );

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

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

  my $parts = $self->{'parts'};
  my $quad = 0;

  if ($parts eq '3') {
    if ($x < 0) {
      if ($y < 0) {
        return undef;
      }
      ($x,$y) = ($y,-1-$x); # rotate -90 and offset
      $quad = 2;
    } else {
      if ($y < 0) {
        ($x,$y) = (-1-$y,$x); # rotate +90 and offset
        $quad = 1;
      }
    }

  } elsif ($parts eq 'octant') {
    if ($y < 0 || $y > $x) { return undef; }
  } elsif ($parts eq 'octant+1') {
    if ($x < 0 || $y < 0 || $y > $x+1) { return undef; }

  } elsif ($parts eq 'octant_up') {
    if ($x < 0 || $x > $y) { return undef; }
  } elsif ($parts eq 'octant_up+1') {
    if ($y < 0 || $x < 0 || $x > $y+1) { return undef; }

  } elsif ($parts eq 'wedge') {
    if ($x < -1-$y || $x > $y) { return undef; }
    if ($x < 0) {
      ($x,$y) = ($y,-1-$x); # rotate -90 and offset
      $quad = 1;
    }
  } elsif ($parts eq 'wedge+1') {
    if ($y < 0 || $x < -2-$y || $x > $y+1) { return undef; }
    if ($x < 0) {
      ($x,$y) = ($y,-1-$x); # rotate -90 and offset
      $quad = 1;
    }

  } elsif ($parts eq 'diagonal') {
    if ($x < -1-$y) { return undef; }  # must have X+Y>=-1
    if ($y < 0) {
      ### lower triangular Y neg ...
      ($x,$y) = (-1-$y,$x); # rotate +90 and offset
    } elsif ($x < 0) {
      ### left triangular X neg ...
      ($x,$y) = ($y,-1-$x); # rotate -90 and offset
      $quad = 2;
    } else {
      ### first quad ...
      $quad = 1;
    }

  } elsif ($parts eq 'diagonal-1') {
    if ($x == 0 && $y == 0) {
      return 0;
    }
    if ($x < -$y) { return undef; }  # must have X+Y>=0
    if ($y <= 0) {
      ### lower triangular Y<=0 ...
      ($x,$y) = (-$y,$x-1); # rotate +90 and offset
    } elsif ($x <= 0) {
      ### left triangular X<=0 ...
      ($x,$y) = ($y-1,-$x); # rotate -90 and offset
      $quad = 2;
    } else {
      ### first quad ...
      $x -= 1;
      $y -= 1;
      $quad = 1;
    }

  } else {
    # parts=1,2,4

    if ($y < 0) {
      if ($parts ne '4') {
        return undef;
      }
      $x = -1-$x; # rotate +180
      $y = -1-$y;
      $quad = 2;
    }

    if ($x < 0) {
      if ($parts eq '1') {
        return undef;
      }
      ($x,$y) = ($y,-1-$x); # rotate +90 and offset
      $quad++;
    }
  }
  ### $quad
  ### xy rotated into first quad: "$x,$y"

  if (is_infinite($x)) {
    return $x;
  }
  if (is_infinite($y)) {
    return $y;
  }

  my $zero = ($x * 0 * $y); # inherit bignum 0
  # my @xbits = bit_split_lowtohigh($x);
  # my @ybits = bit_split_lowtohigh($y);
  # my $exp = max($#xbits, $#ybits);
  # my $len = 2**$exp;

  my ($len,$exp) = round_down_pow(max($x,$y), 2);
  my @depthbits;
  my @ndigits;  # high to low

  foreach my $i (reverse 0 .. $exp) {
    ### at: "x=$x,y=$y  ndigits=".join(',',@ndigits)." len=$len"

    ### assert: $x >= 0
    ### assert: $y >= 0
    ### assert: $x < 2 * $len
    ### assert: $y < 2 * $len
    ### assert: $len == int($len)

    if ($depthbits[$i] = ($x >= $len || $y >= $len ? 1 : 0)) {
      # one of the three parts away from the origin

      if ($y < $len) {
        ### lower right, digit 0 ...
        ($x,$y) = ($len-1-$y,$x-$len);  # rotate +90 and offset
        push @ndigits, 0;
      } elsif ($x >= $len) {
        ### diagonal, digit 1 ...
        ### right, digit 1 ...
        $x -= $len;
        $y -= $len;
        push @ndigits, 1;
      } else {
        ### top left, digit 2 ...
        ($x,$y) = ($y-$len,$len-1-$x);  # rotate -90 and offset
        push @ndigits, 2;
      }
    }

    $len /= 2;
  }

  @ndigits = reverse @ndigits;
  my $n = digit_join_lowtohigh(\@ndigits,3,$zero);
  ### $n
  ### $quad

  if ($quad) {
    ### npower: 3**scalar(@ndigits)
    ### quad npower: $quad * 3**scalar(@ndigits)
    $n += $quad * 3**scalar(@ndigits);
  }
  if ($parts eq 'octant_up' || $parts eq 'octant_up+1'
      || $parts eq 'wedge' || $parts eq 'wedge+1'
      || $parts eq 'diagonal' || $parts eq 'diagonal-1') {
    $n -= (3**scalar(@ndigits) - 1) / 2;
  }

  my $depth = digit_join_lowtohigh(\@depthbits,2,$zero);
  if ($parts eq 'diagonal-1') {
    if ($depth) { $n += 1; }
    $depth += 1;
  }
  if ($parts eq 'octant_up+1') {
    if ($depth) { $n += 1; }
  } elsif ($parts eq 'wedge+1') {
    if ($depth) { $n += 1; }
  }

  ### @depthbits
  ### $depth
  $n += $self->tree_depth_to_n($depth);

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


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

  $x1 = round_nearest ($x1);
  $x2 = round_nearest ($x2);
  $y1 = round_nearest ($y1);
  $y2 = round_nearest ($y2);

  ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
  ($y1,$y2) = ($y2,$y1) if $y1 > $y2;

  my $parts = $self->{'parts'};
  my $xymax;
  if ($parts eq 'octant') {
    if ($y2 < 0 || $x2 < $y1) { return (1,0); }
    $xymax = $x2;
  } elsif ($parts eq 'octant+1') {
    if ($x2 < 0 || $y2 < 0 || $y1 > $x2+1) { return (1,0); }
    $xymax = $x2+1;

  } elsif ($parts eq 'octant_up') {
    if ($x2 < 0 || $y2 < $x1) { return (1,0); }
    $xymax = $y2;
  } elsif ($parts eq 'octant_up+1') {
    if ($x2 < 0 || $y2 < 0 || $x1 > $y2+1) { return (1,0); }
    $xymax = $y2+1;

  } elsif ($parts eq 'wedge') {
    if ($x2 < -1-$y2 || $x1 > $y2) { return (1,0); }
    $xymax = $y2;
  } elsif ($parts eq 'wedge+1') {
    if ($x2 < -2-$y2 || $x1 > $y2+1) { return (1,0); }
    $xymax = $y2+1;

  } else {
    $xymax = max($x2,$y2);
    if ($parts eq 'diagonal') {
      if ($x2 < -1-$y2) { return (1,0); }

    } elsif ($parts eq 'diagonal-1') {
      if ($x2 < -$y2) { return (1,0); }

    } elsif ($parts eq '1') {
      if ($x2 < 0 || $y2 < 0) { return (1,0); }

    } else {
      # parts=2,3,4
      $xymax = max($xymax, -1-$x1);

      if ($parts eq '2') {
        if ($y2 < 0) { return (1,0); }
      } else {
        # parts=3,4
        $xymax = max($xymax, -1-$y1);

        if ($parts eq '3') {
          if ($x2 < 0 && $y2 < 0)  { return (1,0); }
        }
      }
    }
  }
  ### $xymax

  my ($depth) = round_down_pow($xymax,2);
  $depth *= 2;
  if ($parts eq 'diagonal-1') {
    $depth += 1;
  }
  ### depth: 2*$xymax
  return (0,
          $self->tree_depth_to_n($depth));
}

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

# quad(d) = sum i=0tod 3^count1bits(i)
# quad(d) = 2*oct(d) + d
# oct(d) = (quad(d) + d) / 2
# oct(d) = sum i=0tod (3^count1bits(d) + 1)/2
# quad add   1,3,3, 9, 3, 9, 9,27, 3, 9, 9,27,9,27,27,81   A048883
# oct add    1,2,2, 5, 2, 5, 5,14, 2, 5, 5,14,5,14,14,41,  A162784
# oct total  1,3,5,10,12,17,22,36,38,43,48,
#
sub tree_depth_to_n {
  my ($self, $depth) = @_;
  ### tree_depth_to_n(): "depth=$depth"

  if ($depth < 0) {          return undef; }
  if (is_infinite($depth)) { return $depth; }

  my $parts = $self->{'parts'};
  if ($parts eq 'diagonal-1') {
    if ($depth <= 1) {
      return $depth;
    }
    $depth -= 1;
  }

  # pow3 = 3^count1bits(depth)
  my $n = ($depth*0);   # inherit bignum 0
  my $pow3 = $n + 1;    # inherit bignum 1

  foreach my $bit (reverse bit_split_lowtohigh($depth)) {  # high to low
    $n *= 4;
    if ($bit) {
      $n += $pow3;
      $pow3 *= 3;
    }
  }

  if ($parts eq 'octant' || $parts eq 'octant_up') {
    $n = ($n + $depth) / 2;
  } elsif ($parts eq 'octant+1' || $parts eq 'octant_up+1') {
    $n = ($n + 3*$depth) / 2;
    if ($depth) { $n -= 1; }
  } elsif ($parts eq 'wedge') {
    $n += $depth;
  } elsif ($parts eq 'wedge+1') {
    $n += 3*$depth;
    if ($depth) { $n -= 2; }
  } elsif ($parts eq 'diagonal') {
    $n = 2*$n + $depth;
  } elsif ($parts eq 'diagonal-1') {
    $n = 2*$n + 3*$depth - 1;
  } else {
    $n *= $parts;
  }

  return $n;
}

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

  if ($n < 0) { return undef; }
  $n = int($n);
  if (is_infinite($n)) {
    return $n;
  }

  if ($self->{'parts'} eq 'diagonal-1') {
    if ($n == 0) { return 0; }
  }

  my ($depthbits) = _n0_to_depthbits ($n, $self->{'parts'});
  my $depth = digit_join_lowtohigh ($depthbits, 2, $n*0);
  if ($self->{'parts'} eq 'diagonal-1') {
    $depth += 1;
  }
  return $depth;
}

# nwidth = 4^k next 4^(k-1) or to 3*4^(k-1)
#
# octant nwidth = (4^k + 2^k)/2
#               = 2^k*(2^k+1)/2
# next (4^(k-1) + 2^(k-1))/2
#      = 2^(k-1)*(2^(k-1) + 1)/2
#
#       0 1 2   4      6     8          12            16
# oct   0,1,4,7,13,16,22,28,43,46,52,58,73,79,94,109,151,154,160,166,
# oct+1 0,1,3,5,10,12,17,22,36,38,43,48,62,67,81, 95,136,138,143,148,
#
# wedge+1 0,2,8,14,26,32,44,56,86,92,104,116,146,158,188,218,302,308,320,332,
#
sub _n0_to_depthbits {
  my ($n, $parts) = @_;
  ### _n0_to_depthbits(): $n
  ### $parts

  my $numroots = $parts_to_numroots{$parts};
  if ($n < $numroots) {
    ### root point, depth=0 ...
    return ([], 0, $numroots); # $n is in row depth=0
  }

  my $ndepth = 0;
  my ($nmore, $nhalf, $bitpos);
  if ($parts eq 'octant' || $parts eq 'octant_up') {
    ($nmore, $bitpos) = round_down_pow (2*$n, 4);
    $nhalf = 2**$bitpos;
  } elsif ($parts eq 'octant+1' || $parts eq 'octant_up+1') {
    ($nmore, $bitpos) = round_down_pow (2*$n, 4);
    $nhalf = 2**$bitpos;
    $ndepth = -1;

  } elsif ($parts eq 'wedge') {
    ($nmore, $bitpos) = round_down_pow ($n, 4);
    $nhalf = 2**$bitpos;
  } elsif ($parts eq 'wedge+1') {
    ($nmore, $bitpos) = round_down_pow ($n, 4);
    $nhalf = 2**$bitpos;
    $ndepth = -2;

  } elsif ($parts eq 'diagonal') {
    ($nmore, $bitpos) = round_down_pow ($n/2, 4);
    $nmore *= 2;
    $nhalf = 2**$bitpos;

  } elsif ($parts eq 'diagonal-1') {
    if ($n <= 3) {
      ### n in row depth=1 points ...
      return ([], 1, 3);
    }
    ($nmore, $bitpos) = round_down_pow ($n, 4);
    $nmore *= 2;
    $nhalf = 2**$bitpos;
    $ndepth = -1;

  } else {
    ($nmore, $bitpos) = round_down_pow ($n/$numroots, 4);
    $nmore *= $parts;
    $nhalf = 0;
  }
  ### $nmore
  ### $nhalf
  ### $bitpos

  my @depthbits;
  for (;;) {
    ### at: "n=$n ndepth=$ndepth nmore=$nmore nhalf=$nhalf bitpos=$bitpos depthbits=".join(',',map{$_||0}@depthbits)

    my $ncmp;
    if ($parts eq 'wedge' || $parts eq 'diagonal') {
      $ncmp = $ndepth + $nmore + $nhalf;
    } elsif ($parts eq 'wedge+1') {
      $ncmp = $ndepth + $nmore + 3*$nhalf;
    } elsif ($parts eq 'diagonal-1') {
      $ncmp = $ndepth + $nmore + 3*$nhalf;
    } elsif ($nhalf) { # octant, octant_up
      $ncmp = $ndepth + ($nmore + $nhalf)/2;
      if ($parts eq 'octant+1' || $parts eq 'octant_up+1') {
        $ncmp += $nhalf;
      }
    } else {
      $ncmp = $ndepth + $nmore;
    }
    ### $ncmp

    if ($n >= $ncmp) {
      ### use this depthbit ...
      $depthbits[$bitpos] = 1;
      $ndepth = $ncmp;
      $nmore *= 3;
    } else {
      $depthbits[$bitpos] = 0;
    }
    $bitpos--;
    last unless $bitpos >= 0;

    $nmore /= 4;
    $nhalf /= 2;
  }

  # Nwidth = 3**count1bits(depth)
  ### final ...
  ### $nmore
  ### $nhalf
  ### @depthbits
  # except for parts=diagonal
  # ### assert: $nmore == $numroots * 3 ** (scalar(grep{$_}@depthbits))

  if ($parts eq 'wedge' || $parts eq 'diagonal') {
    $nmore += 1;
  } elsif ($parts eq 'wedge+1') {
    $nmore += 3;
  } elsif ($parts eq 'diagonal-1') {
    $nmore += 3;
  } elsif ($parts eq 'octant+1' || $parts eq 'octant_up+1') {
    $nmore = ($nmore + 3) / 2;
  } elsif ($nhalf) {
    ### assert: $nmore % 2 == 1
    $nmore = ($nmore + 1) / 2;
  }

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

# ENHANCE-ME: step by the bits, not by X,Y
# ENHANCE-ME: tree_n_to_depth() by probe?
my @surround8_dx = (1, 0, -1, 0, 1, -1, 1, -1);
my @surround8_dy = (0, 1, 0, -1, 1, 1, -1, -1);
sub tree_n_children {
  my ($self, $n) = @_;
  ### LCornerTree tree_n_children(): $n

  if ($n < 0) {
    ### before n_start ...
    return;
  }
  my ($x,$y) = $self->n_to_xy($n);
  my @n_children;
  foreach my $i (0 .. 7) {
    if (defined (my $n_surround = $self->xy_to_n($x + $surround8_dx[$i],
                                                 $y + $surround8_dy[$i]))) {
      ### $n_surround
      if ($n_surround > $n) {
        my $n_parent = $self->tree_n_parent($n_surround);
        ### $n_parent
        if (defined $n_parent && $n_parent == $n) {
          push @n_children, $n_surround;
        }
      }
    }
  }
  ### @n_children
  # ### assert: scalar(@n_children) == 0 || scalar(@n_children) == 3
  return sort {$a<=>$b} @n_children;
}

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

  my $want_depth = $self->tree_n_to_depth($n);
  if (! defined $want_depth || ($want_depth -= 1) < 0) {
    return undef;
  }
  my ($x,$y) = $self->n_to_xy($n);
  ### $want_depth

  foreach my $i (0 .. 7) {
    if (defined (my $n_surround = $self->xy_to_n($x + $surround8_dx[$i],
                                                 $y + $surround8_dy[$i]))) {
      my $depth_surround = $self->tree_n_to_depth($n_surround);
      ### $n_surround
      ### $depth_surround
      if ($depth_surround == $want_depth) {
        return $n_surround;
      }
    }
  }
  ### no parent ...
  return undef;
}

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

  if ($n < 0)          { return undef; }
  if (is_infinite($n)) { return $n; }

  my ($depthbits, $ndepth, $nwidth) = _n0_to_depthbits($n, $self->{'parts'});
  $n -= $ndepth;      # remaining offset into row

  my $parts = $self->{'parts'};
  if ($parts eq 'octant+1') {
    if ($ndepth > 0 && $n == $nwidth-1) {
      return 0;  # last point in each row doesn't grow, except N=0
    }

  } elsif ($parts eq 'octant_up') {
    $n += $nwidth - 1;   # add to be second half of parts=1 row

  } elsif ($parts eq 'octant_up+1') {
    if ($n == 0) {
      return ($ndepth == 0
              ? undef  # N=0 is infinite
              : 0);    # first of any other row doesn't grow at all
    }
    $n += $nwidth - 3;    # add to be second half of parts=1 row

  } elsif ($parts eq 'wedge') {
    # swap row halves into style of parts=1
    my $noct = $nwidth/2;
    if ($n < $noct) {
      $n += $noct-1;
    } else {
      $n -= $noct;
    }

  } elsif ($parts eq 'wedge+1') {
    if ($ndepth > 0 && ($n == 0 || $n == $nwidth-1)) {
      return 0;  # first and last points don't grow, except N=0,N=1
    }
    # swap row halves into style of parts=1
    my $noct = $nwidth/2;
    if ($n < $noct) {
      $n += $noct-3;
    } else {
      $n -= $noct;
    }

  } elsif ($parts eq 'diagonal') {
    # swap row ends into style of parts=1
    my $noct = ($nwidth+1)/4;
    if ($n < $noct) {
      $n += $noct-1;
    } elsif ($n >= $nwidth - $noct) {
      $n -= ($nwidth - $noct);
    } else {
      $n -= $noct;
    }

  } elsif ($parts eq 'diagonal-1') {
    # swap row ends into style of parts=1
    if (@$depthbits < 2) {
      # N=0to3 depth=0,depth=1 on diagonals so infinite subtree
      return undef;
    }
    if ($n == 0 || $n == $nwidth-1) {
      # first and last of row are extras which never grow
      return 0;
    }
    my $noct = ($nwidth+3)/4;
    if ($n < $noct) {
      $n += $noct - 3;
    } elsif ($n >= $nwidth - $noct) {
      $n -= ($nwidth - $noct);
    } else {
      $n -= $noct;
    }

  } elsif ((my $numroots = $parts_to_numroots{$parts}) > 1) {
    # parts=2,3,4 reduce to parts=1 style
    ### assert: $nwidth % $numroots == 0
    $nwidth /= $numroots;
    $n %= $nwidth;    # Nrem in level, as per n_to_xy()
  }

  ### $depthbits
  ### $n

  foreach my $i (0 .. $#$depthbits) {
    ### $i
    ### N ternary digit: $n%3
    unless ($depthbits->[$i] ^= 1) {  # invert, taken Nrem digit at bit=1
      if (_divrem_mutate($n,3) != 1) {  # stop at lowest non-"1" ternary digit
        $#$depthbits = $i;  # truncate
        return digit_join_lowtohigh($depthbits, 2, $n*0);
      }
    }
  }
  return undef;  # Nrem all 1-digits, so on central infinite spine
}

# return ($quotient, $remainder)
sub _divrem {
  my ($n, $d) = @_;
  if (ref $n && $n->isa('Math::BigInt')) {
    my ($quot,$rem) = $n->copy->bdiv($d);
    if (! ref $d || $d < 1_000_000) {
      $rem = $rem->numify;  # plain remainder if fits
    }
    return ($quot, $rem);
  }
  my $rem = $n % $d;
  return (int(($n-$rem)/$d), # exact division stays in UV
          $rem);
}

# return $remainder, modify $n
# the scalar $_[0] is modified, but if it's a BigInt then a new BigInt is made
# and stored there, the bigint value is not changed
sub _divrem_mutate {
  my $d = $_[1];
  my $rem;
  if (ref $_[0] && $_[0]->isa('Math::BigInt')) {
    ($_[0], $rem) = $_[0]->copy->bdiv($d);  # quot,rem in array context
    if (! ref $d || $d < 1_000_000) {
      return $rem->numify;  # plain remainder if fits
    }
  } else {
    $rem = $_[0] % $d;
    $_[0] = int(($_[0]-$rem)/$d); # exact division stays in UV
  }
  return $rem;
}

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

sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, $self->tree_depth_to_n_end(2**$level-1));
}
sub n_to_level {
  my ($self, $n) = @_;
  my $depth = $self->tree_n_to_depth($n);
  if (! defined $depth) { return undef; }
  my ($pow, $exp) = round_up_pow ($depth+1, 2);
  return $exp;
}

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