File Coverage

blib/lib/Math/PlanePath/UlamWarburton.pm

Criterion	Covered	Total	%
statement	205	278	73.7
branch	83	140	59.2
condition	13	20	65.0
subroutine	24	34	70.5
pod	20	20	100.0
total	345	492	70.1

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 .


#------------------------------------------------------------------------------
# cf
# Ulam/Warburton with cells turning off too
# A079315 cells OFF -> ON
# A079317 cells ON at stage n
# A079316 cells ON at stage n, in first quadrant
# A151921 net gain ON cells


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

package Math::PlanePath::UlamWarburton;
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 = \&Math::PlanePath::_divrem;
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

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

use Math::PlanePath::UlamWarburtonQuarter;

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


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

# octant_up goes up the Y axis spine, dX=0
# all others always have dX!=0
sub absdx_minimum {
  my ($self) = @_;
  return ($self->{'parts'} eq 'octant_up' ? 0 : 1);
}

# used also to validate $self->{'parts'}
my %x_negative = (4         => 1,
                  2         => 1,
                  1         => 0,
                  octant    => 0,
                  octant_up => 0,
                 );
sub x_negative {
  my ($self) = @_;
  return $x_negative{$self->{'parts'}};
}
sub y_negative {
  my ($self) = @_;
  return $self->{'parts'} eq '4';
}

sub x_negative_at_n {
  my ($self) = @_;
  return ($x_negative{$self->{'parts'}} ? $self->n_start + 3 : undef);
}
sub y_negative_at_n {
  my ($self) = @_;
  return ($self->{'parts'} eq '4' ? $self->n_start + 4 : undef);
}

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

{
  my %dir_maximum_dxdy = (4         => [1,-1],  # N=4  South-East
                          2         => [1,-1],  # N=44 South-East
                          1         => [2,-1],  # N=3  ESE
                          octant    => [10,-3], # N=51
                          octant_up => [2,-1],  # N=8  ESE
                         );
  sub dir_maximum_dxdy {
    my ($self) = @_;
    return @{$dir_maximum_dxdy{$self->{'parts'}}};
  }
}

{
  my %_UNDOCUMENTED__turn_any_right_at_n
    = (
       4         => 20,
       2         => 35,
       1         => 2,
       octant    => 4,
       octant_up => 2,
      );
  sub _UNDOCUMENTED__turn_any_right_at_n {
    my ($self) = @_;
    return $self->n_start
      + $_UNDOCUMENTED__turn_any_right_at_n{$self->{'parts'}};
  }
}

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

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

sub n_to_xy {
  my ($self, $n) = @_;
  ### UlamWarburton n_to_xy(): "$n  parts=$self->{'parts'}"

  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'};  # N=0 basis
  if ($n == 0) { return (0,0); }

  my $parts = $self->{'parts'};
  my ($depthsum, $factor, $nrem) = _n0_to_depthsum_factor_rem($n, $parts)
    or return $n;  # N=nan or +inf
  ### depthsum: join(',',@$depthsum)
  ### $factor
  ### n rem within row: $nrem

  if ($parts eq '4') {
    $factor /= 4;
  } elsif ($parts eq '2') {
    $factor /= 2;
    $nrem += ($factor-1)/2;
  } elsif ($parts eq 'octant_up') {
    $nrem += $factor;
  } else {
    $nrem += ($factor-1)/2;
  }
  (my $quad, $nrem) = _divrem ($nrem, $factor);

  ### factor modulus: $factor
  ### $quad
  ### n rem within quad: $nrem
  ### assert: $quad >= 0
  ### assert: $quad <= 3

  my $dhigh = shift @$depthsum;  # highest term
  my @ndigits = digit_split_lowtohigh($nrem,3);
  ### $dhigh
  ### ndigits low to high: join(',',@ndigits)

  my $x = 0;
  my $y = 0;
  foreach my $depthterm (reverse @$depthsum) { # depth terms low to high
    my $ndigit = shift @ndigits;              # N digits low to high
    ### $depthterm
    ### $ndigit

    $x += $depthterm;
    ### bit to x: "$x,$y"

    if ($ndigit) {
      if ($ndigit == 2) {
        ($x,$y) = (-$y,$x);   # rotate +90
      }
    } else {
      # $ndigit==0 (or undef when @ndigits shorter than @$depthsum)
      ($x,$y) = ($y,-$x);   # rotate -90
    }
    ### rotate to: "$x,$y"
  }
  $x += $dhigh;

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

  ### final: "$x,$y"
  return $x,$y;
}
# no Smart::Comments;

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

  $x = round_nearest ($x);
  $y = round_nearest ($y);
  if ($x == 0 && $y == 0) {
    return $self->{'n_start'};
  }

  my $parts = $self->{'parts'};
  if ($parts ne '4'
      && ($y < 0
          || ($parts ne '2' && $x < ($parts eq 'octant' ? $y : 0))
          || ($parts eq 'octant_up' && $x > $y))) {
    return undef;
  }

  my $quad;
  if ($y > $x) {
    ### quad above leading diagonal ...
    #        /
    # above /
    #      /
    if ($y > -$x) {
      ### quad above opposite diagonal, top quarter ...
      #  top
      # \  /
      #  \/
      $quad = 1;
      ($x,$y) = ($y,-$x);  # rotate -90
    } else  {
      ### quad below opposite diagonal, left quarter ...
      #      \
      # left  \
      #       /
      #      /
      $quad = 2;
      $x = -$x;  # rotate -180
      $y = -$y;
    }
  } else {
    ### quad below leading diagonal ...
    #   /
    #  / below
    # /
    if ($y > -$x) {
      ### quad above opposite diagonal, right quarter ...
      #   /
      #  / right
      #  \
      #   \
      $quad = 0;
    } else {
      ### quad below opposite diagonal, bottom quarter ...
      #  /\
      # /  \
      # bottom
      $quad = 3;
      ($x,$y) = (-$y,$x);  # rotate +90
    }
  }
  ### $quad
  ### quad rotated xy: "$x,$y"
  ### assert: ! ($y > $x)
  ### assert: ! ($y < -$x)

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


  my $depth =
    my $ndigits =
      my $n = ($x * 0 * $y);  # inherit bignum 0

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

    my $abs_y = abs($y);
    if ($x && $x == $abs_y) {
      return undef;
    }

    # right quarter diamond
    ### assert: $x >= 0
    ### assert: $x >= abs($y)
    ### assert: $x+abs($y) < 2*$len || $x==abs($y)

    if ($x + $abs_y >= $len) {
      # one of the three quarter diamonds away from the origin
      $x -= $len;
      ### shift to: "$x,$y"

      $depth += $len;
      if ($x || $y) {
        $n *= 3;
        $ndigits++;

        if ($y < -$x) {
          ### bottom, digit 0 ...
          ($x,$y) = (-$y,$x);  # rotate +90

        } elsif ($y > $x) {
          ### top, digit 2 ...
          ($x,$y) = ($y,-$x);  # rotate -90
          $n += 2;
        } else {
          ### right, digit 1 ...
          $n += 1;
        }
      }
    }

    $len /= 2;
  }

  ### $n
  ### $depth
  ### $ndigits
  ### npower: 3**$ndigits
  ### $quad
  ### quad powered: $quad*3**$ndigits

  my $npower = 3**$ndigits;
  if ($parts eq 'octant_up') {
     $n -= $npower;
  } elsif ($parts ne '4') {
     $n -= ($npower-1)/2;
  }

  return $n + $quad*$npower + $self->tree_depth_to_n($depth);
}

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

  my ($dlo, $dhi)
    = _rect_to_diamond_range (round_nearest($x1), round_nearest($y1),
                              round_nearest($x2), round_nearest($y2));
  ### $dlo
  ### $dhi

  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),
          $self->tree_depth_to_n(2*$dhi) - 1);
}

#     x1       |       x2
#     +--------|-------+ y2          xzero true, yzero false
#     |        |       |             diamond min is y1
#     +--------|-------+ y1
#              |
#    ----------O-------------
#
#     |   x1        x2
#     |    +--------+ y2          xzero false, yzero true
#     |    |        |             diamond min is x1
#    -O--------------------
#     |    |        |
#     |    +--------+ y1
#     |
#
sub _rect_to_diamond_range {
  my ($x1,$y1, $x2,$y2) = @_;

  my $xzero = ($x1 < 0) != ($x2 < 0);  # x range covers x=0
  my $yzero = ($y1 < 0) != ($y2 < 0);  # y range covers y=0

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

  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1) }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1) }

  return (($yzero ? 0 : $y1) + ($xzero ? 0 : $x1),
          $x2+$y2);
}


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

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

  if ($n < $self->{'n_start'}) {
    return;
  }
  my ($x,$y) = $self->n_to_xy($n);
  my @ret;
  my $dx = 1;
  my $dy = 0;
  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) = @_;
  ### UlamWarburton tree_n_parent(): $n

  if ($n <= $self->{'n_start'}) {
    return undef;
  }
  my ($x,$y) = $self->n_to_xy($n);
  my $dx = 1;
  my $dy = 0;
  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;
}
# sub tree_n_children {
#   my ($self, $n) = @_;
#   my ($power, $exp) = _round_down_pow (3*$n-2, 4);
#   $exp -= 1;
#   $power /= 4;
#
#   ### $power
#   ### $exp
#   ### pow base: 2 + 4*(4**$exp - 1)/3
#
#   $n -= ($power - 1)/3 * 4 + 2;
#   ### n less pow base: $n
#
#   my @$depthsum = (2**$exp);
#   $power = 3**$exp;
#
#   # find the cumulative levelpoints total <= $n, being the start of the
#   # level containing $n
#   #
#   my $factor = 4;
#   while (--$exp >= 0) {
#     $power /= 3;
#     my $sub = 4**$exp * $factor;
#     ### $sub
#     # $power*$factor;
#     my $rem = $n - $sub;
#
#     ### $n
#     ### $power
#     ### $factor
#     ### consider subtract: $sub
#     ### $rem
#
#     if ($rem >= 0) {
#       $n = $rem;
#       push @$depthsum, 2**$exp;
#       $factor *= 3;
#     }
#   }
#
#   $n += $factor;
#   if (1) {
#     return ($n,$n+1,$n+2);
#   } else {
#     return $n,$n+1,$n+2;
#   }
# }

# Converting quarter ...
# (N-start)*4+1+start = 4*N-4*start+1+start
#                     = 4*N-3*start+1
#
sub tree_depth_to_n {
  my ($self, $depth) = @_;
  ### UlamWarburton tree_depth_to_n(): $depth

  if ($depth == 0) {
    return $self->{'n_start'};
  }
  my $n = $self->Math::PlanePath::UlamWarburtonQuarter::tree_depth_to_n($depth-1);
  if (! defined $n) {
    return undef;
  }
  my $parts = $self->{'parts'};
  if ($parts eq '2') {
    return 2*$n - $self->{'n_start'} + $depth;
  }
  if ($parts eq '1') {
    return $n + $depth;
  }
  if ($parts eq 'octant' || $parts eq 'octant_up') {
    return ($n + 1);
  }
  ### assert: $parts eq '4'
  return 4*$n - 3*$self->{'n_start'} + 1;
}
# sub _NOTWORKING__tree_depth_to_n_range {
#   my ($self, $depth) = @_;
#   my ($nstart, $nend) = $self->Math::PlanePath::UlamWarburtonQuarter::tree_depth_to_n_range($self, $depth)
#     or return;
#   return (4*$nstart-3 + $self->{'n_start'}-1,
#           4*$nend-3 + $self->{'n_start'}-1);
# }


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

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


# 1+3+3+9=16
#
# 0 +1
# 1 +4        <- 0
# 5 +4        <- 1
# 9 +12
# 21 +4     <- 5 + 4+12 = 21 = 5 + 4*(1+3)
# 25 +12
# 37 +12
# 49 +36
# 85 +4     <- 21 + 4+12+12+36  = 21 + 4*(1+3+3+9)
# 89 +12      <- 8   +64
# 101 +12
# 113 +36
# 149
# 161
# 197
# 233
# 341
# 345         <- 16  +256
# 357
# 369

# 1+3 = 4  power 2
# 1+3+3+9 = 16    power 3
# 1+3+3+9+3+9+9+27 = 64    power 4
#
# 4*(1+4+...+4^(l-1)) = 4*(4^l-1)/3
#    l=1 total=4*(4-1)/3 = 4
#    l=2 total=4*(16-1)/3=4*5 = 20
#    l=3 total=4*(64-1)/3=4*63/3 = 4*21 = 84
#
# n = 2 + 4*(4^l-1)/3
# (n-2) = 4*(4^l-1)/3
# 3*(n-2) = 4*(4^l-1)
# 3n-6 = 4^(l+1)-4
# 3n-2 = 4^(l+1)
#
# 3^0+3^1+3^1+3^2 = 1+3+3+9=16
# x+3x+3x+9x = 16x = 256
# 4 quads is 4*16=64
#
# 1+1+3 = 5
# 1+1+3 +1+1+3 +3+3+9 = 25

# 1+4 = 5
# 1+4+4+12 = 21 = 1 + 4*(1+1+3)
# 2  +1
# 3  +3
# 6  +1
# 7  +1
# 10 +3
# 13


# parts=1
#   1+4+...+4^(l-1) + 2^l
#     = (4^l-1)/3 + 2^l
#     = (4^l-1 + 3*2^l)/3
#     = (2^l*(2^l + 3) - 1)/3
#   l=1 total= 3
#   l=2 total= 9
#   l=3 total= 29
#   l=4 total= 101
#
#   N = (4^l-1)/3 + 2^l
#   3*(N-2^l)+1 = 4^l
#   12*(N-2^l)+1 = 4 * 4^l
#
# parts=2
#   N = 2*(4^l-1)/3 + 2^l
#   3/2*(N-2^l)+1 = 4^l
#   6*(N-2^l)+1 = 4 * 4^l
#
# parts=4
#   N = (4^l-1)/3
#   3*N+1 = 4 * 4^l

# use Smart::Comments;

# Return ($aref, $factor, $remaining_n).
# sum(@$aref) = depth starting depth=1
#
sub _n0_to_depthsum_factor_rem {
  my ($n, $parts) = @_;
  ### _n0_to_depthsum_factor_rem(): "$n  parts=$parts"

  my $factor = ($parts eq '4' ? 4 : $parts eq '2' ? 2 : 1);
  if ($n == 0) {
    return ([], $factor, 0);
  }

  my $n3 = 3*$n + 1;
  my $ndepth = 0;
  my $power = $n3;
  my $exp;
  if ($parts eq '4') {
    $power /= 4;
  } elsif ($parts eq '2') {
    $power /= 2;
    $ndepth = -1;
  } elsif ($parts =~ /octant/) {
    $power *= 2;
    $ndepth = 2;
  }
  ($power, $exp) = round_down_pow ($power, 4);
  ### $n3
  ### $power
  ### $exp
  if (is_infinite($exp)) {
    return;
  }

  # ### pow base: ($power - 1)/3 * $factor + 1 + ($parts ne '4' && $exp)
  # $n -= ($power - 1)/3 * $factor + 1;
  # if ($parts ne '4') { $n -= $exp; }
  # ### n less pow base: $n

  my $twopow = 2**$exp;
  my @depthsum;

  for (;
       $exp-- >= 0;
       $power /= 4, $twopow /= 2) {
    ### at: "power=$power twopow=$twopow factor=$factor n3=$n3 ndepth=$ndepth depthsum=".join(',',@depthsum)

    my $nmore = $power * $factor;
    if ($parts ne '4') { $nmore += 3*$twopow; }
    if ($parts =~ /octant/) {
      ### assert: $nmore % 2 == 0
      $nmore = $nmore/2;
    }

    my $ncmp = $ndepth + $nmore;
    ### $nmore
    ### $ncmp

    if ($n3 >= $ncmp) {
      ### go to ncmp, remainder: $n3-$ncmp
      $factor *= 3;
      $ndepth = $ncmp;
      push @depthsum, $twopow;
    }
  }

  if ($parts eq '2') {
    $n3 += 1;
  }

  # ### assert: ($n3 - $ndepth)%3 == 0
  $n = ($n3 - $ndepth) / 3;
  $factor /= 3;

  ### $ndepth
  ### @depthsum
  ### remaining n: $n
  ### assert: $n >= 0
  ### assert: $n < $factor + ($parts ne '4')

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

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

  $n = int($n - $self->{'n_start'});  # N=0 basis
  if ($n < 0) {
    return undef;
  }
  my ($depthsum, $factor, $nrem) = _n0_to_depthsum_factor_rem($n, $self->{'parts'})
    or return $n;  # N=nan or +inf
  ### $depthsum
  ### $factor
  ### $nrem

  my $parts = $self->{'parts'};
  if ($parts eq '4') {
    $factor /= 4;
  } elsif ($parts eq '2') {
    $factor /= 2;
    $nrem += ($factor-1)/2;
  } elsif ($parts eq 'octant_up') {
  } else {
    $nrem += ($factor-1)/2;
  }
  (my $quad, $nrem) = _divrem ($nrem, $factor);

  my $sub = pop @$depthsum;
  while (_divrem_mutate($nrem,3) == 1) {  # low "1" ternary digits of Nrem
    $sub += pop @$depthsum;
  }
  if (@$depthsum) {
    return $depthsum->[-1] - 1 - $sub;
  } else {
    return undef;  # N all 1-digits, on central infinite spine
  }
}

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

sub level_to_n_range {
  my ($self, $level) = @_;
  return ($self->{'n_start'},
          $self->tree_depth_to_n_end(2**($level+1)-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_down_pow ($depth, 2);
  return $exp;
}

# parts=4
# Ndepth(2^a) = 2 + 4*(4^a-1)/3
# Nend(2^a-1) = 1 + 4*(4^a-1)/3 = (4^(a+1)-1)/3
# parts=2
#
# {
#   my %factor = (4         => 16,
#                 2         => 8,
#                 1         => 4,
#                 octant    => 2,
#                 octant_up => 2,
#                );
#   my %constant = (4         => -4,
#                   2         => -5,
#                   1         => -4,
#                   octant    => 0,
#                   octant_up => 0,
#                  );
#   my %spine = (4         => 0,
#                2         => 2,
#                1         => 2,
#                octant    => 1,
#                octant_up => 1,
#               );
#   sub level_to_n_range {
#     my ($self, $level) = @_;
#     my $parts = $self->{'parts'};
#     return ($self->{'n_start'},
#             $self->{'n_start'}
#             + (4**$level * $factor{$parts} + $constant{$parts}) / 3
#             + 2**$level * $spine{$parts});
#   }
# }

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