File Coverage

blib/lib/Math/PlanePath/KochCurve.pm

Criterion	Covered	Total	%
statement	198	241	82.1
branch	52	80	65.0
condition	35	52	67.3
subroutine	30	37	81.0
pod	6	6	100.0
total	321	416	77.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 .


# math-image --path=KochCurve --lines --scale=10
# math-image --path=KochCurve --all --scale=10

# continuous but nowhere differentiable
#
# Sur une Courbe Continue sans Tangente, Obtenue par une Construction
# Géométrique Élémentaire
#
# http://www.nku.edu/~curtin/grenouille.html
# http://www.nku.edu/~curtin/koch_171.jpg
#
# Cesàro, "Remarques sur la Courbe de von Koch." Atti della
# R. Accad. della Scienze Fisiche e Matem. Napoli 12, No. 15, 1-12,
# 1905. Reprinted as §228 in Opere scelte, a cura dell'Unione matematica
# italiana e col contributo del Consiglio nazionale delle ricerche, Vol. 2:
# Geometria, Analisi, Fisica Matematica. Rome: Edizioni Cremonese,
# pp. 464-479, 1964.
#
# Thue-Morse count 1s mod 2 is net direction
# Toeplitz first diffs is turn sequence +1 or -1
#
# J. Ma and J.A. Holdener. When Thue-Morse Meets Koch. In Fractals:
# Complex Geometry, Patterns, and Scaling in Nature and Society, volume 13,
# pages 191-206, 2005.
# http://personal.kenyon.edu/holdenerj/StudentResearch/WhenThueMorsemeetsKochJan222005.pdf
#
# F.M. Dekking. On the Distribution of Digits In Arithmetic Sequences.
# In Seminaire de Theorie des Nombres de Bordeaux, volume 12, 1983, pages
# 3201-3212,
#



package Math::PlanePath::KochCurve;
use 5.004;
use strict;
use List::Util 'sum','first';

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',
  'round_up_pow',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';

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


use constant n_start => 0;
use constant class_x_negative => 0;
use constant class_y_negative => 0;
use constant diffxy_minimum => 0;  # X>=Y octant so X-Y>=0
use constant dx_minimum => -2;
use constant dx_maximum => 2;
use constant dy_minimum => -1;
use constant dy_maximum => 1;
*_UNDOCUMENTED__dxdy_list = \&Math::PlanePath::_UNDOCUMENTED__dxdy_list_six;
use constant absdx_minimum => 1; # never vertical
use constant dsumxy_minimum => -2; # diagonals
use constant dsumxy_maximum => 2;
use constant ddiffxy_minimum => -2;
use constant ddiffxy_maximum => 2;
use constant dir_maximum_dxdy => (1,-1); # South-East
use constant turn_any_straight => 0; # never straight


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

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

  # secret negatives to -.5
  if (2*$n < -1) { return; }
  if (is_infinite($n)) { return ($n,$n); }

  my $x;
  my $y;
  {
    my $int = int($n);
    $x = 2 * ($n - $int);  # usually positive, but n=-0.5 gives x=-0.5
    $y = $x * 0;           # inherit possible bigrat 0
    $n = $int;             # BigFloat int() gives BigInt, use that
  }

  my $len = $y+1;  # inherit bignum 1
  foreach my $digit (digit_split_lowtohigh($n,4)) {
    ### at: "$x,$y  digit=$digit"

    if ($digit == 0) {

    } elsif ($digit == 1) {
      ($x,$y) = (($x-3*$y)/2 + 2*$len,     # rotate +60
                 ($x+$y)/2);

    } elsif ($digit == 2) {
      ($x,$y) = (($x+3*$y)/2 + 3*$len,    # rotate -60
                 ($y-$x)/2   + $len);

    } else {
      ### assert: $digit==3
      $x += 4*$len;
    }
    $len *= 3;
  }

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

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 < 0 || (($x ^ $y) & 1)) {
    ### neg y or parity different ...
    return undef;
  }
  my ($len,$level) = round_down_pow(($x/2)||1, 3);
  ### $level
  ### $len
  if (is_infinite($level)) {
    return $level;
  }

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

# level extends to x= 2*3^level
#                  level = log3(x/2)
#
# exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### KochCurve rect_to_n_range(): "$x1,$y1  $x2,$y2"

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

  if ($x2 < 0 || $y2 < 0
      || 3*$y1 > $x2 ) {   # above line Y=X/3
    return (1,0);
  }

  #        \
  #          \
  #       *    \
  #      / \     \
  # o-+-*   *-+-e  \
  # 0     3     6
  #
  # 3*Y+X/2 - (Y!=0)
  #
  #                  /
  #             *-+-*
  #              \
  #       *       *
  #      / \     /
  # o-+-*   *-+-*
  # 0     3     6   X/2
  #
  my ($len, $level) = round_down_pow ($x2/2, 3);
  return _rect_to_n_range_rot ($len, $level, 0, $x1,$y1, $x2,$y2);



  # (undef, my $level) = round_down_pow ($x2/2, 3);
  # ### $level
  # return (0, 4**($level+1)-1);
}


my @dir6_to_dx = (2, 1,-1,-2, -1, 1);
my @dir6_to_dy = (0, 1, 1, 0, -1,-1);
my @max_digit_to_rot = (1, -2, 1, 0);
my @min_digit_to_rot = (0, 1, -2, 1);
my @max_digit_to_offset = (-1, -1, -1, 2);

sub _rect_to_n_range_rot {
  my ($initial_len, $level_max, $initial_rot, $x1,$y1, $x2,$y2) = @_;
  ### KochCurve _rect_to_n_range_rot(): "$x1,$y1  $x2,$y2  len=$initial_len level=$level_max rot=$initial_rot"

  my ($rot, $len, $x, $y);
  my $overlap = sub {
    ### overlap: "$x,$y len=$len rot=$rot"

    if ($len == 1) {
      return ($x >= $x1 && $x <= $x2
              && $y >= $y1 && $y <= $y2);
    }
    my $len = $len / 3;

    if ($rot < 3) {
      if ($rot == 0) {
        #       *
        #      / \
        # o-+-*   *-+-.
        return ($y <= $y2               # bottom before end
                && $y+$len >= $y1
                && $x <= $x2
                && $x+6*$len > $x1);    # right before end, exclusive
      } elsif ($rot == 1) {
        #       .
        #      /
        # *-+-*
        #  \
        #   *  +-----
        #  /   |x1,y2
        # o
        return ($x <= $x2              # left before end
                && $y+3*$len > $y1     # top after start, exclusive
                && $y-$x <= $y2-$x1);  # diag before corner
      } else {
        # .    |x1,y1
        #  \   +-----
        #   *
        #  /
        # *-+-*
        #      \
        #       o
        return ($y <= $y2              # bottom before end
                && $x-3*$len <=$x2     # left before end
                && $y+$x >= $y1+$x1);  # diag after corner
      }
    } else {
      if ($rot == 3) {
        # .-+-*   *-+-o
        #      \ /
        #       *
        return ($y >= $y1              # top after start
                && $y-$len <= $y2      # bottom before end
                && $x >= $x1           # right after start
                && $x-6*$len < $x2);   # left before end, exclusive
      } elsif ($rot == 4) {
        # x2,y1|    o
        # -----+   /
        #         *
        #          \
        #       *-+-*
        #      /
        #     .
        return ($x >= $x1              # right after start
                && $y-3*$len < $y2     # bottom before end, exclusive
                && $y-$x >= $y1-$x2);  # diag after corner
      } else {
        #    o
        #     \
        #      *-+-*
        #         /
        #        *
        # -----+  \
        # x2,y2|   .
        return ($y >= $y1              # top after start
                && $x+3*$len >= $x1    # right after start
                && $y+$x <= $y2+$x2);  # diag before corner
      }
    }
  };

  my $zero = 0*$x1*$x2*$y1*$y2;
  my @lens = ($initial_len);
  my $n_hi;
  $rot = $initial_rot;
  $len = $initial_len;
  $x = $zero;
  $y = $zero;
  my @digits = (4);

  for (;;) {
    my $digit = --$digits[-1];
    ### max at: "digits=".join(',',@digits)."  xy=$x,$y   len=$len"

    if ($digit < 0) {
      pop @digits;
      if (! @digits) {
        ### nothing found to level_max ...
        return (1, 0);
      }
      ### end of digits, backtrack ...
      $len = $lens[$#digits];
      next;
    }

    my $offset = $max_digit_to_offset[$digit];
    $rot = ($rot - $max_digit_to_rot[$digit]) % 6;
    $x += $dir6_to_dx[$rot] * $offset * $len;
    $y += $dir6_to_dy[$rot] * $offset * $len;

    ### $offset
    ### $rot

    if (&$overlap()) {
      if ($#digits >= $level_max) {
        ### yes overlap, found n_hi ...
        ### digits: join(',',@digits)
        ### n_hi: _digit_join_hightolow (\@digits, 4, $zero)
        $n_hi = _digit_join_hightolow (\@digits, 4, $zero);
        last;
      }
      ### yes overlap, descend ...
      push @digits, 4;
      $len = ($lens[$#digits] ||= $len/3);
    } else {
      ### no overlap, next digit ...
    }
  }

  $rot = $initial_rot;
  $x = $zero;
  $y = $zero;
  $len = $initial_len;
  @digits = (-1);

  for (;;) {
    my $digit = ++$digits[-1];
    ### min at: "digits=".join(',',@digits)."  xy=$x,$y   len=$len"

    if ($digit > 3) {
      pop @digits;
      if (! @digits) {
        ### oops, n_lo not found to level_max ...
        return (1, 0);
      }
      ### end of digits, backtrack ...
      $len = $lens[$#digits];
      next;
    }

    ### $digit
    ### rot increment: $min_digit_to_rot[$digit]
    $rot = ($rot + $min_digit_to_rot[$digit]) % 6;

    if (&$overlap()) {
      if ($#digits >= $level_max) {
        ### yes overlap, found n_lo ...
        ### digits: join(',',@digits)
        ### n_lo: _digit_join_hightolow (\@digits, 4, $zero)
        return (_digit_join_hightolow (\@digits, 4, $zero),
                $n_hi);
      }
      ### yes overlap, descend ...
      push @digits, -1;
      $len = ($lens[$#digits] ||= $len/3);

    } else {
      ### no overlap, next digit ...
      $x += $dir6_to_dx[$rot] * $len;
      $y += $dir6_to_dy[$rot] * $len;
    }
  }
}

# $aref->[0] high digit
sub _digit_join_hightolow {
  my ($aref, $radix, $zero) = @_;
  my @lowtohigh = reverse @$aref;
  return digit_join_lowtohigh(\@lowtohigh, $radix, $zero);
}


my @digit_to_dir = (0, 1, -1, 0);
my @digit_to_nextturn = (1,  # digit=1 (with +1 for "next" N)
                         -2, # digit=2
                         1); # digit=3
sub n_to_dxdy {
  my ($self, $n) = @_;
  ### n_to_dxdy(): $n

  if ($n < 0) {
    return;  # first direction at N=0
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }

  my $int = int($n);
  $n -= $int;
  my @ndigits = digit_split_lowtohigh($int,4);

  my $dir6 = sum(0, map {$digit_to_dir[$_]} @ndigits) % 6;
  my $dx = $dir6_to_dx[$dir6];
  my $dy = $dir6_to_dy[$dir6];

  if ($n) {
    # fraction part

    # lowest non-3 digit, or zero if all 3s (0 above high digit)
    $dir6 += $digit_to_nextturn[ first {$_!=3} @ndigits, 0 ];
    $dir6 %= 6;
    $dx += $n*($dir6_to_dx[$dir6] - $dx);
    $dy += $n*($dir6_to_dy[$dir6] - $dy);
  }
  return ($dx, $dy);
}

sub _UNTESTED__n_to_dir6 {
  my ($self, $n) = @_;
  if ($n < 0) {
    return undef;  # first direction at N=0
  }
  if (is_infinite($n)) {
    return ($n,$n);
  }
  return (sum (map {$digit_to_dir[$_]} digit_split_lowtohigh($n,4))
          || 0) # if empty
    % 6;
}

my @n_to_turn6 = (undef,
                      1,  # +60 degrees
                      -2, # -120 degrees
                      1); # +60 degrees
sub _UNTESTED__n_to_turn6 {
  my ($self, $n) = @_;
  if (is_infinite($n)) {
    return undef;
  }
  while ($n) {
    my $digit = _divrem_mutate($n,4);
    if ($digit) {
      # lowest non-zero digit
      return $n_to_turn6[$digit];
    }
  }
  return 0;
}
sub _UNTESTED__n_to_turn_LSR {
  my ($self, $n) = @_;
  my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
  return ($turn6 > 0 ? 1 : -1);
}
sub _UNTESTED__n_to_turn_left {
  my ($self, $n) = @_;
  my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
  return ($turn6 > 0 ? 1 : 0);
}
sub _UNTESTED__n_to_turn_right {
  my ($self, $n) = @_;
  my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
  return ($turn6 < 0 ? 1 : 0);
}

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

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


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