File Coverage

blib/lib/Math/PlanePath/KochSnowflakes.pm

Criterion	Covered	Total	%
statement	125	150	83.3
branch	35	56	62.5
condition	11	14	78.5
subroutine	26	27	96.3
pod	7	7	100.0
total	204	254	80.3

# 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=KochSnowflakes --lines --scale=10
#
# area approaches sqrt(48)/10
#     *     height=sqrt(1-1/4)=sqrt(3)/2
#    /|\    halfbase=1/2
#   / | \   trianglearea = sqrt(3)/4
#  *-----*
# segments = 3*4^level = 3,12,48,192,...
#
# with initial triangle area=1
# add a new triangle onto each side
# x,y scale by 3* so 9*area
#
# area[level+1] = 9*area[level] + segments
#               = 9*area[level] + 3*4^level

# area[0] = 1
# area[1] = 9*area[0] + 3       = 9 + 3 = 12
# area[2] = 9*area[1] + 3*4
#         = 9*(9*1 + 3) + 3*4
#         = 9*9 + 3*9 + 3*4     = 120
# area[3] = 9*area[2] + 3*4
#         = 9*(9*9 + 3*9 + 3*4) + 3*4^2
#         = 9^3 + 3*9^2 + 3*0*4 + 3*4^2

# area[level+1]
#   = 9^(level+1) + (9^(level+1) - 4^(level+1)) * 3/5
#   = (5*9^(level+1) + 3*9^(level+1) - 3*4^(level+1)) / 5
#   = (8*9^(level+1) - 3*4^(level+1)) / 5
#
# area[level] = (8*9^level - 3*4^level) / 5
#             = 1,12,120,1128,10344,93864,847848
#
#         .
#        / \       area[0] = 1
#       .---.
#
#         .
#        / \       area=[1] = 12 = 9*area[0] + 3*4^0
#   .---.---.---.
#    \ / \ / \ /
#     .---.---.
#    / \ / \ / \
#   .---.---.---.
#        \ /
#         .
#
# area[level] / 9^level
#   = (8*9^level / 9^level - 3*4^level / 9^level) / 5
#   = (8 - 3*(4/0)^level)/5
#   -> 8/5   as level->infinity

# in integer coords
# initial triangle area
#         *       2/3        1*2 / 2 = 1 unit
#        / \
#       *---*     -1/3
#      -1   +1
#
# so area[level] / (sqrt(3)/2)



package Math::PlanePath::KochSnowflakes;
use 5.004;
use strict;
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use vars '$VERSION', '@ISA';
$VERSION = 129;
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;

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


use constant n_frac_discontinuity => 0;
use constant x_negative_at_n => 1;
use constant y_negative_at_n => 1;
use constant sumabsxy_minimum => 2/3; # minimum X=0,Y=2/3
use constant rsquared_minimum => 4/9; # minimum X=0,Y=2/3
# maybe: use constant radius_minimum => 2/3; # minimum X=0,Y=2/3

# jump across rings is WSW slope 2, so following maximums
use constant dx_maximum => 2;
use constant dy_maximum => 1;
use constant dsumxy_maximum => 2;
use constant ddiffxy_maximum => 2;

use constant absdx_minimum => 1; # never vertical
use constant dir_maximum_dxdy => (1,-1); # South-East

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

use constant turn_any_straight => 0; # never straight


#------------------------------------------------------------------------------
sub new {
  my $self = shift->SUPER::new (@_);
  $self->{'sides'} ||= 3; # default
  return $self;
}

# N=1 to 3      3 of, level=1
# N=4 to 15    12 of, level=2
# N=16 to ..   48 of, level=3
#
# each loop = 3*4^level
#
#     n_base = 1 + 3*4^0 + 3*4^1 + ... + 3*4^(level-1)
#            = 1 + 3*[ 4^0 + 4^1 + ... + 4^(level-1) ]
#            = 1 + 3*[ (4^level - 1)/3 ]
#            = 1 + (4^level - 1)
#            = 4^level
#
# each side = loop/3
#           = 3*4^level / 3
#           = 4^level
#
# 6 sides
# n_base = 1 + 2*3*4^0 + ...
#        = 2*4^level - 1
# level = log4 (n+1)/2

### loop 1: 3* 4**1
### loop 2: 3* 4**2
### loop 3: 3* 4**3

# sub _level_to_base {
#   my ($level) = @_;
#   return -3*$level + 4**($level+1) - 2;
# }
# ### level_to_base(1): _level_to_base(1)
# ### level_to_base(2): _level_to_base(2)
# ### level_to_base(3): _level_to_base(3)

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

  my $sides = $self->{'sides'};
  my ($sidelen, $level) = round_down_pow (($sides == 6 ? ($n+1)/2 : $n),
                                          4);
  my $base = ($sides == 6 ? 2*$sidelen - 1 : $sidelen);
  my $rem = $n - $base;

  ### $level
  ### $base
  ### $sidelen
  ### $rem
  ### assert: $n >= $base
  ### assert: $rem >= 0
  ### assert: $rem < $sidelen * $sides

  my $side = int($rem / $sidelen);
  ### $side
  ### $rem
  ### subtract: $side*$sidelen
  $rem -= $side*$sidelen;

  ### assert: $side >= 0 && $side < $sides

  my ($x, $y) = Math::PlanePath::KochCurve->n_to_xy ($rem);
  ### $x
  ### $y

  if ($sides == 3) {
    my $len = 3**($level-1);
    if ($side < 1) {
      ### horizontal rightwards
      return ($x - 3*$len,
              -$y - $len);
    } elsif ($side < 2) {
      ### right slope upwards
      return (($x-3*$y)/-2 + 3*$len,  # flip vert and rotate +120
              ($x+$y)/2 - $len);
    } else {
      ### left slope downwards
      return ((-3*$y-$x)/2,  # flip vert and rotate -120
              ($y-$x)/2 + 2*$len);
    }
  } else {

    #          3
    #     5-----4
    #  4 /       \
    #   /         \ 2
    #  6     o     3
    # 5 \   . .   /
    #    \ .   . / 1
    #     1-----2
    #      0
    #  7
    #
    my $len = 3**$level;
    $x -= $len;    # -y flip vert and offset
    $y = -$y - $len;
    if ($side >= 3) {
      ### rotate 180 ...
      $x = -$x;   # rotate 180
      $y = -$y;
      $side -= 3;
    }
    if ($side >= 2) {
      ### upper right slope upwards ...
      return (($x+3*$y)/-2,  # rotate +120
              ($x-$y)/2);
    }
    if ($side >= 1) {
      ### lower right slope upwards ...
      return (($x-3*$y)/2,  # rotate +60
              ($x+$y)/2);
    }
    ### horizontal ...
    return ($x,$y);
  }
}


# N=1 overlaps N=5
# N=2 overlaps N=7
#      +---------+         +---------+   Y=1.5
#      |         |         |         |
#      |         +---------+         |   Y=7/6 = 1.166
#      |         |         |         |
#      |    * 13 |         |    * 11 |   Y=1
#      |         |         |         |
#      |         |    * 3  |         |   Y=2/3 = 0.666
#      |         |         |         |
#      +---------+         +---------+   Y=0.5
#                |         |
#      +---------+---------+---------+   Y=1/6 = 0.166
#      |         |    O    |         | --Y=0
#      |         |         |         |
#      |         |         |         |
#      |    * 1  |         |    * 2  |   Y=-1/3 = -0.333
#      |         |         |         |
#      +---------+         +---------+   Y=-3/6 = -0.5
#      |         |         |         |
#      +---------+         +---------+   Y=-5/6 = -0.833
#      |         |         |         |
#      |    * 5  |         |    * 7  |   Y=-1
#      |         |         |         |
#      |         |         |         |
#      +---------+         +---------+   Y=-1.5
#
sub xy_to_n {
  return scalar((shift->xy_to_n_list(@_))[0]);
}
sub xy_to_n_list {
  my ($self, $x, $y) = @_;
  ### KochSnowflakes xy_to_n(): "$x, $y"

  $x = round_nearest ($x);
  if (abs($x) <= 1) {
    if ($x == 0) {
      my $y6 = 6*$y;
      if ($y6 >= 1 && $y6 < 7) {
        # Y = 2/3-1/2=1/6 to 2/3+1/2=7/6
        return 3;
      }
    } else {
      my $y6 = 6*$y;
      if ($y6 >= -5 && $y6 < 1) {
        # Y = -1/3-1/2=-5/6 to -1/3+1/2=+1/6
        return (1 + ($x > 0),
                ($y6 < -3 ? (5+2*($x>0)) : ()));   # 5 or 7 up to Y<-1/2
      }
    }
  }

  $y = round_nearest ($y);
  if (($x % 2) != ($y % 2)) {
    ### diff parity...
    return;
  }

  my $high;
  if ($x > 0 && $x >= -3*$y) {
    ### right upper third n=2 ...
    ($x,$y) = ((3*$y-$x)/2,   # rotate -120 and flip vert
               ($x+$y)/2);
    $high = 2;
  } elsif ($x <= 0 && 3*$y > $x) {
    ### left upper third n=3 ...
    ($x,$y) = (($x+3*$y)/-2,             # rotate +120 and flip vert
               ($y-$x)/2);
    $high = 3;
  } else {
    ### lower third n=1 ...
    $y = -$y;  # flip vert
    $high = 1;
  }
  ### rotate/flip is: "$x,$y"

  if ($y <= 0) {
    return;
  }

  my ($len,$level) = round_down_pow($y, 3);
  $level += 1;
  ### $level
  ### $len
  if (is_infinite($level)) {
    return $level;
  }


  $y -= $len;  # shift to Y=0 basis
  $len *= 3;

  ### compare for end: ($x+$y)." >= 3*len=".$len
  if ($x + $y >= $len) {
    ### past end of this level, no points ...
    return;
  }
  $x += $len;  # shift to X=0 basis

  my $n = Math::PlanePath::KochCurve->xy_to_n($x, $y);

  ### plain curve on: ($x+3*$len).",".($y-$len)."  n=".(defined $n && $n)
  ### $high
  ### high: (4**$level)*$high

  if (defined $n) {
    return (4**$level)*$high + $n;
  } else {
    return;
  }
}

# level extends to x= +/- 3^level
#                  y= +/- 2*3^(level-1)
#                   =     2/3 * 3^level
#                  1.5*y = 3^level
#
# ENHANCE-ME: share KochCurve segment checker to find actual min/max
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### KochSnowflakes 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;

  #
  #          |
  # +------  .   -----+
  # |x1,y2  /|\  x2,y2|
  #        / | \
  #       /  |  \
  # -----/---m---\-----
  #     /    |    \
  #    .-----------.
  #          |
  #           y1
  #        -------
  #
  # -y1 bottom horizontal
  # (x2+y2)/2 right side
  # (-x1+y2)/2 left side
  # each giving a power of 3 of the level
  #
  ### right: ($x2+$y2)/2
  ### left: (-$x1+$y2)/2
  ### bottom: -$y1

  my $sides = $self->{'sides'};
  my ($pow, $level) = round_down_pow (max ($sides == 6
                                           ? ($x1/-2,
                                              $x2/2,
                                              -$y1,
                                              $y2)
                                           : (int(($x2+$y2)/2),
                                              int((-$x1+$y2)/2),
                                              -$y1)),
                                      3);
  ### $level
  # end of $level is 1 before base of $level+1
  return (1, 4**($level+2) - 1);
}

#------------------------------------------------------------------------------
# Nstart = 4^k
# length = 3*4^k many points
# Nend = Nstart + length-1
#      = 4^k + 3*4^k - 1
#      = 4*4^k - 1
#      = Nstart(k+1) - 1
sub level_to_n_range {
  my ($self, $level) = @_;
  my $pow = 4**$level;
  return ($pow, 4*$pow-1);
}
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 1) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  my ($sidelen, $level) = round_down_pow (($self->{'sides'} == 6 ? ($n+1)/2 : $n),
                                          4);
  return $level;
}

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