File Coverage

blib/lib/Math/PlanePath/HTree.pm

Criterion	Covered	Total	%
statement	35	214	16.3
branch	0	80	0.0
condition	0	33	0.0
subroutine	12	24	50.0
pod	12	12	100.0
total	59	363	16.2

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


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

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_join_lowtohigh';

use Math::PlanePath::LCornerTree;
*_divrem = \&Math::PlanePath::LCornerTree::_divrem;

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

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


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


use constant n_start => 1;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant tree_num_children_list => (0,1,2);

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

  if ($n < 1) { 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
  }

  my $zero = $n * 0;
  my @nbits = bit_split_lowtohigh($n);
  my $len = (2 + $zero) ** int($#nbits/2);   # half power
  ### $len

  # eg. N=9 "up" block, nbits even, $#nbits odd, dX=1 dY=0 East
  # High 1-bit of sub-tree position does rotate +90 to North for
  # initial step North N=8 to N=9.
  #
  # eg. N=17 "right" block, nbits odd, $#nbits even, dX=0 dY=-1 South
  # High 1-bit of sub-tree position does rotate +90 to East for
  # initial step East N=16 to N=17.
  #
  my $dx = ($#nbits % 2);
  my $dy = $dx - 1;
  ### initial direction: "$dx, $dy"

  my $x = my $y = $len;
  if ($dx) {
    $y *= 2;
  }
  $x -= 1;
  $y -= 1;
  ### initial xy: "$x, $y"

  ### assert: $nbits[-1] == 1
  pop @nbits;  # strip high 1-bit which is N=2^k spine bit
  ### strip high spine 1-bit to: @nbits
  # N=10001xxx
  #    ^^^^^^^ leaving these

  # Strip high 0-bits of sub-tree.
  # Could have @nbits all 0-bits if N=2^k spine point.
  while (@nbits && ! $nbits[-1]) {
    pop @nbits;
  }
  ### strip high zeros to: @nbits
  # N=10001xxx
  #       ^^^^ leaving these
  # or if an N=1000000 spine point then @nbits now empty and no move $x,$y.

  foreach my $bit (reverse @nbits) {   # high to low
    ### at: "$x,$y  bit=$bit len=$len dir=$dx,$dy"

    ($dx,$dy) = ($dy,-$dx);  # rotate -90 for $bit==0
    if ($bit) {
      $dx = -$dx;            # rotate 180 to give rotate +90 for $bit==1
      $dy = -$dy;
    }
    ### turn to: "dir=$dx,$dy"

    if ($dx) { $len /= 2; }  # halve when going horizontal
    $x += $dx * $len;
    $y += $dy * $len;
  }

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

#                                  |    [37]   |           34=10,0010
# 13            |                 4,13--5,13--6,13         35=10,0011
#               |                  |     |     |
# 12            |                        |                 36=10,0100
#               |          [33]          |                 37=10,0101
# 11       [35]1,7---------3,7----------5,7[34]
#               |           |            |
# 10  0,10      |           |     4,6    |
#      |        |    |      |      |     |     |
#  9  0,9-----1,9---2,9     |     4,5---5,5---6,5
#      |             |      |      |    [36]   |
#  8  0,8                   |     4,4
#                           |
#  7                   [32]3,7-------------------------------
#                           |
#  6  0,6            2,6    |      [30]        [29]          17=1,0001
#      |       [9]    |     |       |    [19]   |            18=1,0010
#  5  0,5------1,5---2,5    |     [23]---5,5---[22]          20=1,0100
#      |        |     |     |       |     |     |            24=1,1000
#  4  0,4       |    2,4    |      [31]   |    [28]
#     [15]      |           |             |
#  3        [8]1,3---------3,3-----------5,3[17]
#               |          [16]           |
#  2  0,2[3]    |    2,2[7]        [24]   |    [27]
#      |        |     |             |     |     |
#  1  0,1[2]---1,1---2,1[5]    [20]4,1---5,1---6,1[21]       21=1,0101
#      |        4     |             |    [18]   |
#  0  0,0[1]         2,0[6]        [25]        [26]
#
#      0        1     2      3      4     5      6

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

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

  my ($len,$exp);
  my $n;
  my ($xlen,$xexp) = round_down_pow($x+1, 2);
  my ($ylen,$yexp) = round_down_pow($y+1, 2);
  if ($yexp > $xexp) {
    ### Y bigger ...
    $len = $ylen/2;
    $exp = $yexp;
    $n = $len*$len*2;

    $y -= $ylen;
    ### to: "$x,$y  len=$len"
    if ($x == $len-1 && $y == -1) {
      ### spine ...
      return $n;
    }

  } else {
    ### X bigger, fake initial X bit ...
    $n = $xlen;
    $len = $xlen;
    $exp = $xexp;
    $n = $len*$len;

    if ($x == $len-1 && $y == $len-1) {
      ### spine ...
      return $n;
    }
  }

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

  ### $n
  my @nbits;  # high to low


  while ($exp-- >= 0) {
    ### X at: "$x,$y len=$len"
    ### assert: $len >= 1
    ### assert: $x >= 0
    ### assert: $y >= 0
    ### assert: $x <= 2*$len
    ### assert: $y <= 2*$len

    if ($x == $len-1 && $y == $len-1) {
      ### midpoint X ...

      # ### nbits HtoL: join('',@nbits)
      # ### shift off high nbits ...
      # shift @nbits;

      last;
    }
    if ($x >= $len) {
      ### move left, digit 0 ...
      $x -= $len;
      push @nbits, 0;
    } else {
      ### rotate 180, digit 0: "$x,$y"
      push @nbits, 1;
      $x = $len-2 - $x;
      $y = 2*$len-2 - $y;
      if ($x < 0 || $y < 0) {
        ### outside: "$x,$y"
        return undef;
      }
    }

    ### Y at: "$x,$y len=$len"
    ### assert: $x >= 0
    ### assert: $y >= 0
    ### assert: $x <= $len
    ### assert: $y <= 2*$len

    if ($y == $len-1 && $x == $len/2-1) {
      ### midpoint Y ...

      last;
    }
    if ($y >= $len) {
      ### move down only, digit 1 ...
      $y -= $len;
      push @nbits, 1;
    } else {
      ### rotate 180, digit 0 ...
      push @nbits, 0;
      $x = $len-2 - $x;
      $y = $len-2 - $y;
      if ($x < 0 || $y < 0) {
        ### outside: "$x,$y"
        return undef;
      }
    }

    $len /= 2;
  }

  if ($yexp > $xexp) {
  } else {
    ### nbits HtoL: join('',@nbits)
    ### shift off high nbits ...
    shift @nbits;
  }

  ### nbits HtoL: join('',@nbits)

  if ($yexp > $xexp) {
  } else {
  }


  @nbits = reverse @nbits;
  push @nbits, 1;

  ### nbits HtoL: join('',reverse @nbits)
  return $n + digit_join_lowtohigh(\@nbits,2);
}


#  7
#                |
#  6  |       |  |  |       |
#  5  *---*---*  |  *---*---*
#  4  |   |   |  |  |   |   |
#         |      |      |
#  3      *------*------*
#         |             |
#  2  |   |   |     |   |   |
#  1  *---*---*     *---*---*
#  0  |       |     |       |
#     0   1   2  3  4   5   6
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### HTree rect_to_n_range(): "$x1,$y1  $x2,$y2"

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

  $x2 = max($x1,$x2);
  $y2 = max($y1,$y2);
  if ($x2 < 0 || $y2 < 0) {
    ### all outside first quadrant ...
    return (1,0);
  }

  my ($pow) = round_down_pow(max(2*$x2, $y2) + 1,
                             2);
  return (1, $pow*$pow);

  # my ($xpow) = round_down_pow($x2+1, 2);
  # my ($ypow) = round_down_pow($y2+1, 2);
  # if ($xpow > $ypow) {
  #   return (1, 2*$xpow*$xpow);
  # } else {
  #   return (1, $ypow*$ypow);
  # }
}

# 5*2^n, 6*2^n successively
# being start of last two rows of each sub-tree
#   depth=13   160 10100000
#   depth=14   192 11000000
#
sub tree_depth_to_n {
  my ($self, $depth) = @_;
  ### HTree depth_to_n(): $depth
  $depth = int($depth);
  if ($depth < 3) {
    if ($depth < 0) { return undef; }
    return $depth + 1;
  }
  ($depth, my $rem) = _divrem ($depth-3, 2);
  return (5+$rem) * 2**$depth;
}

# spine 2^n
#
sub tree_depth_to_n_end {
  my ($self, $depth) = @_;
  ### HTree depth_to_n(): $depth
  if ($depth < 0) {
    return undef;
  }
  return 2**int($depth);
}

#  0  N=1
#      |
#  1  N=2--
#      |   \
#  2  N=3   N=4-------
#            |        \
#  3        N=5        N=8 ------------
#          /  \         |               \
#  4     N=6  N=7      N=9              N=16----------
#                     /    \             |            \
#  5              N=10      N=11        N=17           N=32------
#                 /  \      /  \        /   \           |        \
#  6           N=12 N=13  N=14 N=15   N=18 N=19        N=33       N=64---
#                                     /         \    /     \       |     \
#                                  N=20   [4of]      N=34 N=35     N=65   N=128
# 0  1  = 1
# 1  1  = 1
# 2  2  = 1 + 1
# 3  2  = 1 + 1
# 4  4  = 2 + 1 + 1
# 5  4  = 2 + 1 + 1
# 6  8  = 4 + 2 + 1 + 1
# 7  8
# 8  16
# 9  16
#
# 1 + sum i=0 to floor(depth/2)-1 of 2^i
#   = 2^floor(depth/2)

sub tree_depth_to_width {
  my ($self, $depth) = @_;
  ### HTree tree_n_to_subheight(): $depth

  if ($depth < 0) {
    return undef;
  }
  $depth = int($depth/2);
  return 2**$depth;
}


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

  if ($n < 1) { return undef; }
  $n = int($n);
  if (is_infinite($n)) { return $n; }
  my ($pow,$depth) = round_down_pow($n,2);
  if ($n -= $pow) {
    ($pow, my $exp) = round_down_pow($n,2);
    $depth += $exp+1;
  }
  return $depth;
}


# (n-pow)*2 + pow
#   = 2*n-2*pow+pow
#   = 2*n-pow
# (n-pow)*2 < pow
#   2*n-2*pow < pow
#   2*n-pow < 2*pow
# 1011 -> 011 -> 110 -> 1110
sub tree_n_children {
  my ($self, $n) = @_;
  ### HTree tree_n_children(): $n

  if ($n < 1) {
    return;
  }

  my ($pow) = round_down_pow($n,2);
  if ($pow == 1) {
    return $n+1;
  }
  if ($n == $pow) {
    return ($n+1, 2*$n);
  }

  $n *= 2;
  $n -= $pow;
  if ($n < 2*$pow) {
    return ($n, $n+1);
  } else {
    return;
  }
}

# 1
# 2 3
# 4 5,   6,7
#   101  11x
# 8 9,   10,11,    12,13,14,15
#   1001 1010,1011  1100,1101
#
# (n-pow)/2 + pow
#   = (n-pow+2*pow)/2
#   = (n+pow)/2
#
sub tree_n_parent {
  my ($self, $n) = @_;
  ### HTree tree_n_parent(): $n

  if ($n < 2) {
    return undef;
  }
  my ($pow) = round_down_pow($n,2);
  if ($n == $pow) {
    return $n/2;
  }
  return ($n + $pow - ($n%2)) / 2;
}

# length of the run of 0s immediately below the high 1-bit
#   10001111
#    ^^^------ 3 0-bits
#
sub tree_n_to_subheight {
  my ($self, $n) = @_;
  ### HTree tree_n_to_subheight(): $n

  if ($n < 2) {
    return undef;
  }
  my ($pow,$exp) = round_down_pow($n,2);
  if ($n == $pow) {
    return undef;  # spine, infinite
  }
  $n -= $pow;
  ($pow, my $exp2) = round_down_pow($n,2);
  return $exp - $exp2 - 1;
}


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

sub level_to_n_range {
  my ($self, $level) = @_;
  return (1, 2**(2*$level+1) - 1);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 1) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  _divrem_mutate ($n, 2);
  my ($pow, $exp) = round_up_pow ($n+1, 4);
  return $exp;
}

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