File Coverage

blib/lib/Math/PlanePath/KochPeaks.pm

Criterion	Covered	Total	%
statement	137	148	92.5
branch	26	34	76.4
condition	7	9	77.7
subroutine	26	27	96.3
pod	5	5	100.0
total	201	223	90.1

# 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::KochPeaks;
use 5.004;
use strict;
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

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

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_down_pow';
use Math::PlanePath::KochCurve;
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

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


use constant class_y_negative => 0;
use constant n_frac_discontinuity => .5;
use constant x_negative_at_n => 1;
use constant sumabsxy_minimum  => 1; # minimum X=1,Y=0
use constant absdiffxy_minimum => 1; # X=Y never occurs
use constant rsquared_minimum  => 1; # minimum X=1,Y=0

use constant dx_maximum => 2;
use constant dy_minimum => -1;
use constant dy_maximum => 1;
use constant absdx_minimum => 1; # never vertical
use constant dsumxy_maximum => 2; # diagonal NE
use constant ddiffxy_maximum => 2; # diagonal NW
use constant dir_maximum_dxdy => (1,-1); # South-East
use constant turn_any_straight => 0; # never straight


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

# N=1 to 3      3 of, level=0
# N=4 to 12     9 of, level=1
# N=13 to 45   33 of, level=2
#
# N=0.5 to 3.49     diff=3
# N=3.39 to 12.49   diff=9
# N=12.5 to 45.5    diff=33
#
# each length = 2*4^level + 1
#
#     Nstart = 1 + 2*4^0 + 1 + 2*4^1 + 1 + ... + 2*4^(level-1) + 1
#            = 1 + level + 2*[ 4^0 + 4^1 + ... + 4^(level-1) ]
#            = level+1 + 2*[ (4^level - 1)/3 ]
#            = level+1 + (2*4^level - 2)/3
#            = level + (2*4^level - 2 + 3)/3
#            = level + (2*4^level + 1)/3
#
#     3*n = 2*4^level + 1
#     3*n-1 = 2*4^level
#     (3*n-1)/2 = 4^level
#
# Nbase = 0.5 + 2*4^0 + 1 + 2*4^1 + 1 + ... + 2*4^(level-1) + 1
#       = level + (2*4^level + 1)/3 - 1/2
#       = level + 2/3*4^level + 1/3 - 1/2
#       = level + 2/3*4^level - 1/6
#       = level + 4/6*4^level - 1/6
#       = level + (4*4^level - 1)/6
#       = level + (4^(level+1) - 1)/6
#
#     6*N = 4^(level+1) - 1
#     6*N + 1 = 4^(level+1)
#     level+1 = log4(6*N + 1)
#     level = log4(6*N + 1) - 1
#
### loop 1: (2*4**1 + 1)/3
### loop 2: (2*4**2 + 1)/3
### loop 3: (2*4**3 + 1)/3

# sub _n_to_level {
#   my ($n) = @_;
#   my ($side, $level) = round_down_pow(6*$n + 1, 4);
#   my $base = $level + (2*$side + 1)/3 - .5;
#   ### $level
#   ### $base
#   if ($base > $n) {
#     $level--;
#     $side /= 4;
#     $base = $level + (2*$side + 1)/3 - .5;
#     ### $level
#     ### $base
#   }
#   return ($level, $base, $side + .5);
# }
# sub _level_to_base {
#   my ($level) = @_;
#   return $level + (2*$side + 1)/3 - .5;
# }

sub _n_to_side_level_base {
  my ($n) = @_;
  my ($side, $level) = round_down_pow((3*$n-1)/2, 4);
  my $base = $level + (2*$side + 1)/3;
  ### $level
  ### $base
  if (2*$n+1 < 2*$base) {
    $level--;
    $side /= 4;
    $base = $level + (2*$side + 1)/3;
    ### $level
    ### $base
  }
  return ($side, $level, $base);
}

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

  # $n<0.5 no good for Math::BigInt circa Perl 5.12, compare in integers
  return if 2*$n < 1;

  if (is_infinite($n)) { return ($n,$n); }

  my ($side, $level, $base) = _n_to_side_level_base($n);

  my $rem = $n - $base;
  my $frac;
  if ($rem < 0) {
    ### neg frac
    $frac = $rem;
    $rem = 0;
  } elsif ($rem > 2*$side) {
    ### excess frac
    $frac = $rem - 2*$side;
    $rem -= $frac;
  } else {
    ### no frac
    $frac = 0;
  }
  ### $frac
  ### $rem
  ### $n

  ### next base would be: ($level+1) + (2*4**($level+1) + 1)/3
  ### assert: $n-$frac >= $base
  ### assert: $n-$frac < ($level+1) + (2*4**($level+1) + 1)/3

  ### assert: $rem>=0
  ### assert: $rem < 2 * 4 ** $level + 1
  ### assert: $rem <= 2*$side+1

  my $pos = 3**$level;
  if ($rem < $side) {
    my ($x, $y) = Math::PlanePath::KochCurve->n_to_xy($rem);
    ### left side: $rem
    ### flat: "$x,$y"
    $x += 2*$frac;
    return (($x-3*$y)/2 - $pos,    # rotate +60
            ($x+$y)/2);
  } else {
    my ($x, $y) = Math::PlanePath::KochCurve->n_to_xy($rem-$side);
    ### right side: $rem-$side
    ### flat: "$x,$y"
    $x += 2*$frac;
    return (($x+3*$y)/2,           # rotate -60
            ($y-$x)/2 + $pos);
  }
}

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

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

  if ($y < 0 || ! (($x ^ $y) & 1)) {
    ### neg y or parity...
    return undef;
  }
  my ($len,$level) = round_down_pow ($y+abs($x), 3);
  ### $level
  ### $len
  if (is_infinite($level)) {
    return $level;
  }

  my $n;
  if ($x < 0) {
    $x += $len;
    ($x,$y) = (($x+3*$y)/2,   # rotate -60
               ($y-$x)/2);
    $n = 0;
    ### left rotate -60 to: "x=$x,y=$y   n=$n"
  } else {
    $y -= $len;
    ($x,$y) = (($x-3*$y)/2,   # rotate +60
               ($x+$y)/2);
    $n = 1;
    ### right rotate +60 to: "x=$x,y=$y   n=$n"
  }

  foreach (1 .. $level) {
    $n *= 4;
    ### at: "level=$level len=$len   x=$x,y=$y  n=$n"
    if ($x < $len) {
      $len /= 3;
      my $rel = 2*$len;
      if ($x < $rel) {
        ### digit 0
      } else {
        ### digit 1 sub: "$rel to x=".($x-$rel)
        $x -= $rel;
        ($x,$y) = (($x+3*$y)/2,   # rotate -60
                   ($y-$x)/2);
        $n += 1;
      }
    } else {
      $len /= 3;
      $x -= 4*$len;
      if ($x < $y) {   # before diagonal
        ### digit 2...
        ($x,$y) = (($x-3*$y)/2 + 2*$len,     # rotate +60
                   ($x+$y)/2);
        $n += 2;
      } else {
        #### digit 3...
        $n += 3;
      }
    }
  }
  ### end at: "x=$x,y=$y   n=$n"
  if ($x) {
    ### endmost point
    $n += 1;
    $x -= 2;
  }
  if ($x != 0 || $y != 0) {
    return undef;
  }
  return $n + $level + (2*4**$level + 1)/3 + ($x == 2);
}


# level extends to x= +/- 3^level
#                  y= 0 to 3^level
#
# diagonal X=Y or Y=-X is lowest in a level, so round down abs(X)+Y to pow 3
#
# end of level is 1 before base of level+1
#     basenext = (level+1) + (2*4^(level+1) + 1)/3
#     basenext-1 = level + (2*4^(level+1) + 1)/3
#                = level + (8*4^level + 1)/3

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

  $x1 = round_nearest ($x1);
  $y1 = round_nearest ($y1);
  $x2 = round_nearest ($x2);
  $y2 = round_nearest ($y2);
  ### rounded: "$x1,$y1  $x2,$y2"

  if ($y1 < 0 && $y2 < 0) {
    return (1,0);
  }

  # can't make use of the len=3**$level returned by round_down_pow()
  my ($len, $level) = round_down_pow (max(abs($x1),abs($x2))
                                      + max($y1, $y2),
                                      3);
  ### $level
  return (1, $level + (8 * 4**$level + 1)/3);
}

# peak Y is at N = Nstart + (count-1)/2
#                = level + (2*4^level + 1)/3 + (2*4^level + 1 - 1)/2
#                = level + (2*4^level + 1)/3 + (2*4^level)/2
#                = level + (2*4^level + 1)/3 + 4^level
#                = level + (2*4^level + 1 + 3*4^level)/3
#                = level + (5*4^level + 1)/3

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

sub level_to_n_range {
  my ($self, $level) = @_;
  my $pow = 4**$level;
  return ((2*$pow + 1)/3 + $level,
          (8*$pow + 1)/3 + $level);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 1) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  my ($side, $level, $base) = _n_to_side_level_base($n);
  return $level;
}

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