File Coverage

blib/lib/Math/PlanePath/UlamWarburtonQuarter.pm

Criterion	Covered	Total	%
statement	176	237	74.2
branch	67	112	59.8
condition	16	26	61.5
subroutine	21	28	75.0
pod	15	15	100.0
total	295	418	70.5

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Kevin Ryde

# This file is part of Math-PlanePath.
#
# Math-PlanePath is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3, or (at your option) any later
# version.
#
# Math-PlanePath is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
# for more details.
#
# You should have received a copy of the GNU General Public License along
# with Math-PlanePath.  If not, see .


package Math::PlanePath::UlamWarburtonQuarter;
use 5.004;
use strict;
use Carp 'croak';
use List::Util 'sum';

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

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

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


use constant parameter_info_array =>
  [
   { name            => 'parts',
     share_key       => 'parts_ulamwarburton_quarter',
     display         => 'Parts',
     type            => 'enum',
     default         => '1',
     choices         => ['1','octant','octant_up' ],
     choices_display => ['1','Octant','Octant Up' ],
     description     => 'Which parts of the plane to fill.',
   },
   Math::PlanePath::Base::Generic::parameter_info_nstart1(),
  ];

use constant class_x_negative => 0;
use constant class_y_negative => 0;

sub diffxy_minimum {
  my ($self) = @_;
  return ($self->{'parts'} eq 'octant' ? 0 : undef);
}
sub diffxy_maximum {
  my ($self) = @_;
  return ($self->{'parts'} eq 'octant_up' ? 0 : undef);
}

# Minimum dir=0 at N=13 dX=2,dY=0.
# Maximum dir seems dX=13,dY=-9 at N=149 going top-left part to new bottom
# right diagonal.
my %dir_maximum_dxdy = (1         => [13,-9],
                        octant    => [1,-1],  # South-East
                        octant_up => [0,-1],  # South
                       );
sub dir_maximum_dxdy {
  my ($self) = @_;
  return @{$dir_maximum_dxdy{$self->{'parts'}}};
}

sub tree_num_children_list {
  my ($self) = @_;
  return ($self->{'parts'} =~ /octant/
          ? (0, 1, 2, 3)
          : (0, 1,    3));
}

#------------------------------------------------------------------------------
sub new {
  my $self = shift->SUPER::new(@_);
  if (! defined $self->{'n_start'}) {
    $self->{'n_start'} = $self->default_n_start;
  }
  my $parts = ($self->{'parts'} ||= '1');
  if (! exists $dir_maximum_dxdy{$parts}) {
    croak "Unrecognised parts option: ", $parts;
  }
  return $self;
}

# 7   7   7   7
#   6       6
# 7   5   5   7
#       4
# 3   3   5   7
#   2       6
# 1   3   7   7
#
# 1+1+3=5
# 5+1+3*5=21
# 1+3 = 4
# 1+3+3+9 = 16
#
#       0
# 1  0 +1
# 2  1 +1       <- 1
# 3  2 +3
# 4  5 +1       <- 1 + 4 = 5
# 5  6 +3
# 6  9 +3
# 7  12 +9
# 8  21         <- 1 + 4 + 16 = 21

# 1+3 = 4  power 2
# 1+3+3+9 = 16    power 3
# 1+3+3+9+3+9+9+27 = 64    power 4
#
# (1+4+16+...+4^(l-1)) = (4^l-1)/3
#    l=1 total=(4-1)/3 = 1
#    l=2 total=(16-1)/3 = 5
#    l=3 total=(64-1)/3=63/3 = 21
#
# n = 1 + (4^l-1)/3
# n-1 = (4^l-1)/3
# 3n-3 = (4^l-1)
# 3n-2 = 4^l
#
# 3^0+3^1+3^1+3^2 = 1+3+3+9=16
# x+3x+3x+9x = 16x = 256
#
#               22
# 20  19  18  17
#   12      11
# 21   9   8  16
#        6
#  5   4   7  15
#    2      10
#  1   3  13  14
#

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

  if ($n < $self->{'n_start'}) { 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
  }

  $n = $n - $self->{'n_start'} + 1;  # N=1 basis
  if ($n == 1) { return (0,0); }

  my ($depthsum, $nrem, $rowwidth) = _n1_to_depthsum_rem_width($self,$n)
    or return ($n,$n); # N==nan or N==+inf

  ### assert: $nrem >= 0
  ### assert: $nrem < $width
  if ($self->{'parts'} eq 'octant_up') {
    $nrem += ($rowwidth-1)/2;
    ### assert: $nrem < $width
  }

  my @ndigits = digit_split_lowtohigh($nrem,3);
  my $dhigh = shift(@$depthsum) - 1;  # highest term
  my $x = 0;
  my $y = 0;
  foreach my $depthsum (reverse @$depthsum) { # depth terms low to high
    my $ndigit = shift @ndigits;              # N digits low to high
    ### $depthsum
    ### $ndigit

    $x += $depthsum;
    $y += $depthsum;
    ### depthsum to xy: "$x,$y"

    if ($ndigit) {
      if ($ndigit == 2) {
        ($x,$y) = (-$y,$x);   # rotate +90
      }
    } else {
      # digit==0 (or undef when run out of @ndigits)
      ($x,$y) = ($y,-$x);   # rotate -90
    }
    ### rotate to: "$x,$y"
  }

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

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

  $x = round_nearest ($x);
  $y = round_nearest ($y);
  my $parts = $self->{'parts'};
  if ($y < 0
      || $x < ($parts eq 'octant' ? $y : 0)
      || ($parts eq 'octant_up' && $x > $y)) {
    return undef;
  }
  if ($x == 0 && $y == 0) {
    return $self->{'n_start'};
  }
  $x += 1;  # pushed away by 1 ...
  $y += 1;

  my ($len, $exp) = round_down_pow ($x + $y, 2);
  if (is_infinite($exp)) { return $exp; }

  my $depth
    = my $n
      = ($x * 0 * $y);  # inherit bignum 0
  my $rowwidth = $depth + 1;

  while ($exp-- >= 0) {
    ### at: "$x,$y  n=$n len=$len"

    # first quadrant square
    ### assert: $x >= 0
    ### assert: $y >= 0
    # ### assert: $x < 2*$len
    # ### assert: $y < 2*$len

    if ($x >= $len || $y >= $len) {
      # one of three quarters away from origin
      #     +---+---+
      #     | 2 | 1 |
      #     +---+---+
      #     |   | 0 |
      #     +---+---+

      $x -= $len;
      $y -= $len;
      ### shift to: "$x,$y"

      if ($x) {
        unless ($y) {
          return undef;  # x==0, y!=0, nothing
        }
      } else {
        if ($y) {
          return undef;  # x!=0, y-=0, nothing
        }
      }

      $depth += $len;
      if ($x || $y) {
        $rowwidth *= 3;
        $n *= 3;
        if ($y < 0) {
          ### bottom right, digit 0 ...
          ($x,$y) = (-$y,$x);  # rotate +90
        } elsif ($x >= 0) {
          ### top right, digit 1 ...
          $n += 1;
        } else {
          ### top left, digit 2 ...
          ($x,$y) = ($y,-$x);  # rotate -90
          $n += 2;
        }
      }
    }

    $len /= 2;
  }

  ### $n
  ### $depth

  if ($self->{'parts'} eq 'octant_up') {
    $n -= ($rowwidth-1)/2;
  }

  return $n + $self->tree_depth_to_n($depth-1);
}

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

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

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

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

  if ($x1 < 0) { $x1 *= 0; }
  if ($y1 < 0) { $y1 *= 0; }

  # level numbers
  my $dlo = ($x1 > $y1 ? $x1 : $y1)+1;
  my $dhi = ($x2 > $y2 ? $x2 : $y2);
  ### $dlo
  ### $dhi

  # round down to level=2^k numbers
  if ($dlo) {
    ($dlo) = round_down_pow ($dlo,2);
  }
  ($dhi) = round_down_pow ($dhi,2);

  ### rounded to pow2: "$dlo  ".(2*$dhi)

  return ($self->tree_depth_to_n($dlo-1),
          $self->tree_depth_to_n(2*$dhi-1));
}

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

# ENHANCE-ME: step by the bits, not by X,Y
sub tree_n_children {
  my ($self, $n) = @_;
  if ($n < $self->{'n_start'}) {
    return;
  }
  my ($x,$y) = $self->n_to_xy($n);
  my @ret;
  my $dx = 1;
  my $dy = 1;
  foreach (1 .. 4) {
    if (defined (my $n_child = $self->xy_to_n($x+$dx,$y+$dy))) {
      if ($n_child > $n) {
        push @ret, $n_child;
      }
    }
    ($dx,$dy) = (-$dy,$dx); # rotate +90
  }
  return sort {$a<=>$b} @ret;
}
sub tree_n_parent {
  my ($self, $n) = @_;
  if ($n <= $self->{'n_start'}) {
    return undef;
  }
  my ($x,$y) = $self->n_to_xy($n);
  my $dx = 1;
  my $dy = 1;
  foreach (1 .. 4) {
    if (defined (my $n_parent = $self->xy_to_n($x+$dx,$y+$dy))) {
      if ($n_parent < $n) {
        return $n_parent;
      }
    }
    ($dx,$dy) = (-$dy,$dx); # rotate +90
  }
  return undef;
}

# level = depth+1 = 2^a + 2^b + 2^c + 2^d ...       a>b>c>d...
# Ndepth = 1 + (-1
#               +       4^a
#               +   3 * 4^b
#               + 3^2 * 4^c
#               + 3^3 * 4^d + ...) / 3
sub tree_depth_to_n {
  my ($self, $depth) = @_;
  ### tree_depth_to_n(): $depth
  if (is_infinite($depth)) {
    return $depth;
  }
  unless ($depth >= 0) {
    return undef;
  }
  my $n = $depth*0;        # inherit bignum 0
  my $pow3 = 1 + $n;       # inherit bignum 1
  foreach my $bit (reverse bit_split_lowtohigh($depth+1)) {  # high to low
    $n *= 4;
    if ($bit) {
      $n += $pow3;
      $pow3 *= 3;
    }
  }
  if ($self->{'parts'} =~ /octant/) {
    $n = ($n + (3*$depth-1))/6;
  } else {
    $n = ($n-1)/3;
  }
  return $n + $self->{'n_start'};
}

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

  $n = int($n - $self->{'n_start'} + 1);  # N=1 basis
  if ($n < 1) {
    return undef;
  }
  (my $depthsum, $n) = _n1_to_depthsum_rem_width($self,$n)
    or return $n;  # N==nan or N==+infinity
  return sum(-1, @$depthsum);
}

# Return ($aref, $remaining_n).
# sum(@$aref) = depth starting depth=1
#
# depth+1 = 2^k
# Ndepth(depth) = (4^k+2)/3
#   3N-2 = 4^k
# NdepthOct(depth) = ((4^k+2)/3 + 2^k)/2
#   6N-2 = 4^k + 3*2^k
#
sub _n1_to_depthsum_rem_width {
  my ($self, $n) = @_;
  ### _n1_to_depthsum_rem_width(): $n

  my $octant = ($self->{'parts'} =~ /octant/);
  my ($power, $exp) = round_down_pow (($octant ? 6 : 3)*$n - 2, 4);
  if (is_infinite($exp)) {
    return;
  }

  ### $power
  ### $exp
  ### pow base: ($power - 1)/3 + 1

  {
    my $sub = ($power + 2)/3;  # (power-1)/3 + 1
    if ($octant) {
      $sub = ($sub + 2**$exp) / 2;
      ### prospective sub: $sub
      ### assert: $sub == ($power + 3 * 2 ** $exp + 2)/6

      if ($sub > $n) {
        $exp -= 1;
        $power /= 4;
        $sub = ($power + 3*2**$exp + 2)/6;
      }
    }
    ### assert: $sub <= $n
    $n -= $sub;
  }
  ### n less pow base: $n

  my @depthsum = (2**$exp);

  # find the cumulative levelpoints total <= $n, being the start of the
  # level containing $n
  #
  my $factor = 1;
  while (--$exp >= 0) {
    $power /= 4;
    my $sub = $power * $factor;
    if ($octant) {
      $sub = ($sub + 2**$exp)/2;
    }
    ### $sub
    my $rem = $n - $sub;

    ### $n
    ### $power
    ### $factor
    ### consider subtract: $sub
    ### $rem

    if ($rem >= 0) {
      $n = $rem;
      push @depthsum, 2**$exp;
      $factor *= 3;
    }
  }

  ### _n1_to_depthsum_rem_width() result ...
  ### @depthsum
  ### remaining n: $n
  ### assert: $n >= 0
  ### assert: $n < $factor

  return (\@depthsum, $n, $factor);
}


# at 0,2 turn and new height limit
# at 1 keep existing depth limit
# N=30 rem=1 = 0,1 depth=11=8+2+1=1011 width=9
#
sub tree_n_to_subheight {
  my ($self, $n) = @_;
  ### tree_n_to_subheight(): $n

  $n = int($n - $self->{'n_start'} + 1);  # N=1 basis
  if ($n < 1) {
    return undef;
  }
  my ($depthsum, $nrem, $rowwidth) = _n1_to_depthsum_rem_width($self,$n)
    or return $n;  # N==nan or N==+infinity
  ### $depthsum
  ### $nrem

  if ($self->{'parts'} eq 'octant_up') {
    $nrem += ($rowwidth-1)/2;
  }

  my $sub = pop @$depthsum;
  while (@$depthsum && _divrem_mutate($nrem,3) == 1) {
    $sub += pop @$depthsum;
  }
  if (@$depthsum) {
    return $depthsum->[-1] - 1 - $sub;
  } else {
    return undef; # $nrem all 1-digits
  }
}

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

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

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