File Coverage

blib/lib/Math/PlanePath/KochSquareflakes.pm

Criterion	Covered	Total	%
statement	212	245	86.5
branch	68	90	75.5
condition	13	14	92.8
subroutine	23	25	92.0
pod	5	5	100.0
total	321	379	84.7

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 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::KochSquareflakes;
use 5.004;
use strict;

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

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

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


use constant n_frac_discontinuity => 0;

use constant parameter_info_array =>
  [ { name        => 'inward',
      display     => 'Inward',
      type        => 'boolean',
      default     => 0,
      description => 'Whether to direct the sides of the square inward, rather than outward.',
    } ];

use constant x_negative_at_n => 1;
use constant y_negative_at_n => 1;
use constant sumabsxy_minimum => 1;
use constant rsquared_minimum => 0.5; # minimum X=0.5,Y=0.5

# jump across rings is South-West, so
use constant dx_maximum => 1;
use constant dy_maximum => 1;
use constant dsumxy_maximum => 2; # diagonal NE
use constant ddiffxy_maximum => 2;
use constant ddiffxy_minimum => -2;
use constant dir_maximum_dxdy => (1,-1); # South-East

# N=1,2,3,4  gcd(1/2,1/2) = 1/2
use constant gcdxy_minimum => 1/2;

use constant turn_any_straight => 0; # never straight


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

# level 0 inner square
# sidelen = 4^level
# ring points 4*4^level
# Nend = 4 * [ 1 + ... + 4^level ]
#      = 4 * (4^(level+1) - 1) / 3
#      = (4^(level+2) - 4) / 3
# Nstart = Nend(level-1) + 1
#        = (4^(level+1) - 4) / 3 + 1
#        = (4^(level+1) - 4 + 3) / 3
#        = (4^(level+1) - 1) / 3
#
# level    Nstart             Nend
#    0     (4-1)/3=1          (16-4)/3=12/3=4
#    1     (16-1)/3=15/3=5    (64-4)/3=60/3=20
#    2     (64-1)/3=63/3=21   (256-4)/3=252/3=84
#    3     (256-1)/3=255/3=85
#

sub n_to_xy {
  my ($self, $n) = @_;
  ### KochSquareflakes n_to_xy(): $n
  if ($n < 1) { return; }

  my $frac;
  {
    my $int = int($n);
    $frac = $n - $int;
    $n = $int;       # BigFloat int() gives BigInt, use that
  }

  # (4^(level+1) - 1) / 3 = N
  # 4^(level+1) - 1 = 3*N
  # 4^(level+1) = 3*N+1
  #
  my ($pow,$level) = round_down_pow (3*$n + 1, 4);
  ### $level
  ### $pow
  if (is_infinite($level)) { return ($level,$level); }

  # Nstart = (4^(level+1)-1)/3 with $power=4^(level+1) here
  #
  $n -= ($pow-1)/3;

  ### base: ($pow-1)/3
  ### next base would be: (4*$pow-1)/3
  ### n remainder from base: $n

  my $sidelen = $pow/4;
  (my $rot, $n) = _divrem ($n, $sidelen);  # high part is rot
  ### $sidelen
  ### n remainder: $n
  ### $rot

  ### assert: $n>=0
  ### assert: $n < 4 ** $level

  my @horiz = (1);
  my @diag = (1);
  my $i = 0;
  while (--$level > 0) {
    $horiz[$i+1] = 2*$horiz[$i] + 2*$diag[$i];
    $diag[$i+1]  = $horiz[$i] + 2*$diag[$i];
    ### horiz: $horiz[$i+1]
    ### diag: $diag[$i+1]
    $i++;
  }

  ### horiz: join(', ',@horiz)
  ### $i
  my $x
    = my $y
      = ($n * 0) + $horiz[$i]/-2;  # inherit bignum
  if ($rot & 1) {
    ($x,$y) = (-$y,$x);
  }
  if ($rot & 2) {
    $x = -$x;
    $y = -$y;
  }
  $rot *= 2;

  my $inward = $self->{'inward'};
  my @digits = digit_split_lowtohigh($n,4);

  while ($i > 0) {
    $i--;
    my $digit = $digits[$i] || 0;

    my ($dx, $dy, $drot);
    if ($digit == 0) {
      $dx = 0;
      $dy = 0;
      $drot = 0;
    } elsif ($digit == 1) {
      if ($rot & 1) {
        $dx = $diag[$i];
        $dy = $diag[$i];
      } else {
        $dx = $horiz[$i];
        $dy = 0;
      }
      $drot = ($inward ? 1 : -1);
    } elsif ($digit == 2) {
      if ($rot & 1) {
        if ($inward) {
          $dx = $diag[$i];
          $dy = $diag[$i] + $horiz[$i];
        } else {
          $dx = $diag[$i] + $horiz[$i];
          $dy = $diag[$i];
        }
      } else {
        $dx = $horiz[$i] + $diag[$i];
        $dy = $diag[$i];
        unless ($inward) { $dy = -$dy; }
      }
      $drot = ($inward ? -1 : 1);
    } elsif ($digit == 3) {
      if ($rot & 1) {
        $dx = $dy = $diag[$i] + $horiz[$i];
      } else {
        $dx = $horiz[$i] + 2*$diag[$i];
        $dy = 0;
      }
      $drot = 0;
    }
    ### delta: "$dx,$dy   rot=$rot   drot=$drot"

    if ($rot & 2) {
      ($dx,$dy) = (-$dy,$dx);
    }
    if ($rot & 4) {
      $dx = -$dx;
      $dy = -$dy;
    }
    ### delta with rot: "$dx,$dy"

    $x += $dx;
    $y += $dy;
    $rot += $drot;
  }

  {
    my $dx = $frac;
    my $dy = ($rot & 1 ? $frac : 0);
    if ($rot & 2) {
      ($dx,$dy) = (-$dy,$dx);
    }
    if ($rot & 4) {
      $dx = -$dx;
      $dy = -$dy;
    }
    $x = $dx + $x;
    $y = $dy + $y;
  }

  return ($x,$y);
}

my @inner_to_n = (1,2,4,3);

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

  # +/- 0.75
  if (4*$x < 3 && 4*$y < 3 && 4*$x >= -3 && 4*$y >= -3) {
    return $inner_to_n[($x >= 0) + 2*($y >= 0)];
  }

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

  # quarter curve segment and high digit
  my $n;
  {
    my $negx = -$x;
    if (($y > 0 ? $x > $y : $x >= $y)) {
      ### below leading diagonal ...
      if ($negx > $y) {
        ### bottom quarter ...
        $n = 1;
      } else {
        ### right quarter ...
        $n = 2;
        ($x,$y) = ($y, $negx);   # rotate -90
      }
    } else {
      ### above leading diagonal
      if ($y > $negx) {
        ### top quarter ...
        $n = 3;
        $x = $negx;   # rotate 180
        $y = -$y;
      } else {
        ### right quarter ...
        $n = 4;
        ($x,$y) = (-$y, $x);   # rotate +90
      }
    }
  }
  $y = -$y;
  ### rotate to: "$x,$y   n=$n"

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

  my @horiz;
  my @diag;
  my $horiz = 1;
  my $diag = 1;
  for (;;) {
    push @horiz, $horiz;
    push @diag, $diag;
    my $offset = $horiz+$diag;
    my $nextdiag = $offset + $diag;  # horiz + 2*diag
    ### $horiz
    ### $diag
    ### $offset
    ### $nextdiag

    if ($y <= $nextdiag) {
      ### found level at: "top=$nextdiag vs y=$y"
      $y -= $offset;
      $x += $offset;
      last;
    }
    $horiz = 2*$offset;  # 2*horiz+2*diag
    $diag = $nextdiag;
  }
  ### base subtract to: "$x,$y"


  if ($self->{'inward'}) {
    $y = -$y;
    ### inward invert to: "$x,$y"
  }

  ### origin based side: "$x,$y   horiz=$horiz diag=$diag  with levels ".scalar(@horiz)

  # loop 4*1, 4*4, 4*4^2 etc, extra +1 on the digits to include that in the sum
  #
  my $slope;
  while (@horiz) {
    ### at: "$x,$y slope=".($slope||0)." n=$n"
    $horiz = pop @horiz;
    $diag = pop @diag;
    $n *= 4;

    if ($slope) {
      if ($y < $diag) {
        ### slope digit 0 ...
        $n += 1;
      } else {
        $x -= $diag;
        $y -= $diag;
        ### slope not digit 0, move to: "$x,$y"

        if ($y < $horiz) {
          ### digit 1 ...
          $n += 2;
          ($x,$y) = ($y, -$x);   # rotate -90
          $slope = 0;
        } else {
          $y -= $horiz;
          ### slope not digit 1, move to: "$x,$y"

          if ($x < $horiz) {
            ### digit 2 ...
            $n += 3;
            $slope = 0;

          } else {
            ### digit 3 ...
            $n += 4;
            $x -= $horiz;
          }
        }
      }

    } else {
      if ($x < $horiz) {
        ### digit 0 ...
        $n += 1;
      } else {
        $x -= $horiz;
        ### not digit 0, move to: "$x,$y"

        if ($x < $diag) {
          ### digit 1 ...
          $n += 2;
          $slope = 1;
        } else {
          $x -= $diag;
          ### not digit 1, move to: "$x,$y"

          if ($x < $diag) {
            ### digit 2 ...
            $n += 3;
            $slope = 1;
            ($x,$y) = ($diag-$y, $x);   # offset and rotate +90

          } else {
            ### digit 3 ...
            $n += 4;
            $x -= $diag;
          }
        }
      }
    }
  }
  ### final: "$x,$y n=$n"

  if ($x == 0 && $y == 0) {
    return $n;
  } else {
    return undef;
  }
}

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

  foreach ($x1,$y1, $x2,$y2) {
    if (is_infinite($_)) {
      return (0, $_);
    }
    $_ = abs(round_nearest($_));
  }
  if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  my $max = ($x2 > $y2 ? $x2 : $y2);

  # Nend = 4 * [ 1 + ... + 4^level ]
  #      = 4 + 16 + ... + 4^(level+1)
  #
  my $horiz = 4;
  my $diag = 3;
  my $nhi = 4;
  for (;;) {
    $nhi += 1;
    $nhi *= 4;
    my $nextdiag = $horiz + 2*$diag;
    if (($self->{'inward'} ? $horiz : $nextdiag) >= 2*$max) {
      return (1, $nhi);
    }
    $horiz = $nextdiag + $horiz;   # 2*$horiz + 2*$diag;
    $diag = $nextdiag;
  }
}


#------------------------------------------------------------------------------
# Nstart = (4^(k+1) - 1)/3
# Nend = Nstart(k+1) - 1
#      = (4*4^(k+1) - 1)/3 - 1
#      = (4*4^(k+1) - 1 - 3)/3
#      = (4*4^(k+1) - 4)/3
#      = 4*(4^(k+1) - 1)/3
#      = 4*Nstart(k)

sub level_to_n_range {
  my ($self, $level) = @_;
  my $n_lo = (4**($level+1) - 1)/3;
  return ($n_lo, 4*$n_lo);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 1) { return undef; }
  if (is_infinite($n)) { return $n; }
  my ($pow,$exp) = round_down_pow (3*$n + 1, 4);
  return $exp-1;
}

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