File Coverage

blib/lib/Math/PlanePath/CCurve.pm

Criterion	Covered	Total	%
statement	140	192	72.9
branch	35	70	50.0
condition	20	22	90.9
subroutine	20	25	80.0
pod	8	8	100.0
total	223	317	70.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 .


package Math::PlanePath::CCurve;
use 5.004;
use strict;
use List::Util 'min','max','sum';

use vars '$VERSION', '@ISA';
$VERSION = 129;
use Math::PlanePath;
use Math::PlanePath::Base::NSEW;
@ISA = ('Math::PlanePath::Base::NSEW',
        'Math::PlanePath');

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';
use Math::PlanePath::Base::Digits
  'round_up_pow',
  'round_down_pow',
  'bit_split_lowtohigh',
  'digit_split_lowtohigh',
  'digit_join_lowtohigh';
*_divrem = \&Math::PlanePath::_divrem;
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

use Math::PlanePath::KochCurve;
*_digit_join_hightolow = \&Math::PlanePath::KochCurve::_digit_join_hightolow;

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


# Not sure about this yet ... 2 or 4?  With mirror images too 8 arms would
# fill the plane everywhere 4-visited points double-traversed segments.
# use constant parameter_info_array => [ { name      => 'arms',
#                                          share_key => 'arms_2',
#                                          display   => 'Arms',
#                                          type      => 'integer',
#                                          minimum   => 1,
#                                          maximum   => 2,
#                                          default   => 1,
#                                          width     => 1,
#                                          description => 'Arms',
#                                        } ];

use constant n_start => 0;
use constant x_negative_at_n => 6;
use constant y_negative_at_n => 22;
use constant 1.02 _UNDOCUMENTED__dxdy_list_at_n => 7;


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

sub new {
  my $self = shift->SUPER::new(@_);
  $self->{'arms'} = max(1, min(2, $self->{'arms'} || 1));
  return $self;
}


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

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

  my $zero = ($n * 0);  # inherit bignum 0
  my $x = $zero;
  my $y = $zero;
  {
    my $int = int($n);
    $x = $n - $int;  # inherit possible BigFloat
    $n = $int;        # BigFloat int() gives BigInt, use that
  }

  # initial rotation from arm number $n mod $arms
  my $rot = _divrem_mutate ($n, $self->{'arms'});

  my $len = $zero+1;
  foreach my $digit (digit_split_lowtohigh($n,4)) {
    ### $digit

    if ($digit == 0) {
      ($x,$y) = ($y,-$x);    # rotate -90
    } elsif ($digit == 1) {
      $y -= $len;            # at Y=-len
    } elsif ($digit == 2) {
      $x += $len;            # at X=len,Y=-len
      $y -= $len;
    } else {
      ### assert: $digit == 3
      ($x,$y) = (2*$len - $y,  # at X=2len,Y=-len and rotate +90
                 $x-$len);
    }
    $rot++; # to keep initial direction
    $len *= 2;
  }

  if ($rot & 2) {
    $x = -$x;
    $y = -$y;
  }
  if ($rot & 1) {
    ($x,$y) = (-$y,$x);
  }

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

# point N=2^(2k) at XorY=+/-2^k  radius 2^k
#       N=2^(2k-1) at X=Y=+/-2^(k-1) radius sqrt(2)*2^(k-1)
# radius = sqrt(2^level)
# R(l)-R(l-1) = sqrt(2^level) - sqrt(2^(level-1))
#             = sqrt(2^level) * (1 - 1/sqrt(2))
# about 0.29289

# len=1 extent of lower level 0
# len=4 extent of lower level 2
# len=8 extent of lower level 4+1 = 5
# len=16 extent of lower level 8+3
# len/2 + len/4-1

my @digit_to_rot = (-1, 1, 0, 1);
my @dir4_to_dsdd = ([1,-1],[1,1],[-1,1],[-1,-1]);

sub xy_to_n {
  return scalar((shift->xy_to_n_list(@_))[0]);
}
sub xy_to_n_list {
  my ($self, $x, $y) = @_;
  ### CCurve xy_to_n(): "$x, $y"

  $x = round_nearest($x);
  $y = round_nearest($y);
  my $zero = $x*0*$y;

  ($x,$y) = ($x + $y, $y - $x);  # sum and diff
  if (is_infinite($x)) { return $x; }
  if (is_infinite($y)) { return $y; }

  my @n_list;
  foreach my $dsdd (@dir4_to_dsdd) {
    my ($ds,$dd) = @$dsdd;
    ### attempt: "ds=$ds  dd=$dd"
    my $s = $x;  # sum X+Y
    my $d = $y;  # diff Y-X
    my @nbits;

    until ($s >= -1 && $s <= 1 && $d >= -1 && $d <= 1) {
      ### at: "s=$s, d=$d   nbits=".join('',reverse @nbits)
      my $bit = $s % 2;
      push @nbits, $bit;
      if ($bit) {
        $s -= $ds;
        $d -= $dd;
        ($ds,$dd) = ($dd,-$ds); # rotate -90
      }

      # divide 1/(1+i) = (1-i)/(1^2 - i^2)
      #                = (1-i)/2
      # so multiply (s + i*d) * (1-i)/2
      #   s = (s + d)/2
      #   d = (d - s)/2
      #
      ### assert: (($s+$d)%2)==0

      # this form avoids overflow near DBL_MAX
      my $odd = $s % 2;
      $s -= $odd;
      $d -= $odd;
      $s /= 2;
      $d /= 2;
      ($s,$d) = ($s+$d+$odd, $d-$s);
    }

    # five final positions
    #      .   0,1   .       ds,dd
    #           |
    #    -1,0--0,0--1,0
    #           |
    #      .   0,-1  .
    #
    ### end: "s=$s d=$d  ds=$ds dd=$dd"

    # last step must be East dx=1,dy=0
    unless ($ds == 1 && $dd == -1) { next; }

    if ($s == $ds && $d == $dd) {
      push @nbits, 1;
    } elsif ($s != 0 || $d != 0) {
      next;
    }
    # ended s=0,d=0 or s=ds,d=dd, found an N
    push @n_list, digit_join_lowtohigh(\@nbits, 2, $zero);
    ### found N: "$n_list[-1]"
  }
  ### @n_list
  return sort {$a<=>$b} @n_list;
}

# f = (1 - 1/sqrt(2) = .292
# 1/f = 3.41
# N = 2^level
# Rend = sqrt(2)^level
# Rmin = Rend / 2  maybe
# Rmin^2 = (2^level)/4
# N = 4 * Rmin^2
#
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### CCurve rect_to_n_range(): "$x1,$y1  $x2,$y2"

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

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

  my ($len,$level) = _rect_to_k ($x1,$y1, $x2,$y2);
  if (is_infinite($level)) {
    return (0, $level);
  }
  return (0, 4*$len*$len*$self->{'arms'} - 1);
}

# N=16 is Y=4 away   k=2
# N=64 is Y=-8+1=-7 away  k=3
# N=256=4^4 is X=2^4=16-3=-7 away  k=4
# dist = 2^k - (2^(k-2)-1)
#      = 2^k - 2^(k-2) + 1
#      = 4*2^(k-2) - 2^(k-2) + 1
#      = 3*2^(k-2) + 1
#   k=2 3*2^(2-2)+1=4   len=4^2=16
#   k=3 3*2^(3-2)+1=7   len=4^3=64
#   k=4 3*2^(4-2)+1=13
# 2^(k-2) = (dist-1)/3
# 2^k = (dist-1)*4/3
#
# up = 3*2^(k-2+1) + 1
# 2^(k+1) = (dist-1)*4/3
# 2^k = (dist-1)*2/3
#
# left = 3*2^(k-2+1) + 1
# 2^(k+1) = (dist-1)*4/3
# 2^k = (dist-1)*2/3
#
# down = 3*2^(k-2+1) + 1
# 2^(k+1) = (dist-1)*4/3
# 2^k = (dist-1)*2/3
#
# m=2 4*(2-1)/3=4/3=1
# m=4 4*(4-1)/3=4
sub _rect_to_k {
  my ($x1,$y1, $x2,$y2) = @_;
  ### _rect_to_k(): $x1,$y1

  {
    my $m = max(abs($x1),abs($y1),abs($x2),abs($y2));
    if ($m < 2) {
      return (2, 1);
    }
    if ($m < 4) {
      return (4, 2);
    }
    ### round_down: 4*($m-1)/3
    my ($len, $k) = round_down_pow (4*($m-1)/3, 2);
    return ($len, $k);
  }

  my $len;
  my $k = 0;

  my $offset = -1;
  foreach my $m ($x2, $y2, -$x1, -$y1) {
    $offset++;
    ### $offset
    ### $m
    next if $m < 0;

    my ($len1, $k1);
    # if ($m < 2) {
    #   $len1 = 1;
    #   $k1 = 0;
    # } else {
    # }

    ($len1, $k1) = round_down_pow (($m-1)/3, 2);
    next if $k1 < $offset;
    my $sub = ($offset-$k1) % 4;
    $k1 -= $sub;  # round down to k1 == offset mod 4

    if ($k1 > $k) {
      $k = $k1;
      $len = $len1 / 2**$sub;
    }
  }

  ### result: "k=$k  len=$len"
  return ($len, 2*$k);
}



my @dir4_to_dx = (1,0,-1,0);
my @dir4_to_dy = (0,1,0,-1);

# Not yet supporting "arms" parameter ...
sub n_to_dxdy {
  my ($self, $n) = @_;
  ### n_to_dxdy(): $n

  if ($n < 0) { return; }
  if (is_infinite($n)) { return ($n,$n); }
  my $int = int($n);
  $n -= $int;  # $n fraction part

  my @digits = bit_split_lowtohigh($int);
  my $dir = (sum(@digits)||0) & 3;  # count of 1-bits
  my $dx = $dir4_to_dx[$dir];
  my $dy = $dir4_to_dy[$dir];

  if ($n) {
    # apply fraction part $n

    # count low 1-bits is right turn of N+1, apply as dir-(turn-1) so decr $dir
    while (shift @digits) {
      $dir--;
    }

    # this with turn=count-1 turn which is dir++ worked into swap and negate
    # of dir4_to_dy parts
    $dir &= 3;
    $dx -= $n*($dir4_to_dy[$dir] + $dx);  # with rot-90 instead of $dir+1
    $dy += $n*($dir4_to_dx[$dir] - $dy);

    # this the equivalent with explicit dir++ for turn=count-1
    # $dir++;
    # $dir &= 3;
    # $dx += $n*($dir4_to_dx[$dir] - $dx);
    # $dy += $n*($dir4_to_dy[$dir] - $dy);
  }

  ### result: "$dx, $dy"
  return ($dx,$dy);
}

#------------------------------------------------------------------------------
# k even
#         S[h]
#       ---------
#      /          \  Z[h-1]
#     /            \
#    |              |  S[h-1]
#     \            / Z[h-2]
#      --        --
# Hb[k] = S[h] + 2*S[h-1] + S[h] + 2*(Z[h-1]/2 - Z[h-2]/2)
#          + sqrt(2)*(2*Z[h-1]/2 + 2*Z[h-2]/2)
#       = 2*S[h] + 2*S[h-1] + Z[h-1]-Z[h-2]   + sqrt(2) * (Z[h-1] + Z[h-2])
#       = 2*2^h + 2*2^(h-1) + 2*2^(h-1)-2 - (2*2^(h-2)-2)   + sqrt(2) * (2*2^(h-1)-2 + 2*2^(h-2)-2)
#       = 3*2^h             + 2*2^(h-1)-2 - 2*2^(h-2) + 2   + sqrt(2) * (3*2^(h-1) - 4)
#       = 3*2^h             +   2^(h-1)                     + sqrt(2) * (3*2^(h-1) - 4)
#       = 7*2^(h-1) + sqrt(2) * (3*2^(h-1) - 4)
#       = 7*sqrt(2)^(2h-2) + sqrt(2) * (3*sqrt(2)^(2h-2) - 4)
#       = 7*sqrt(2)^(k-2) + sqrt(2) * (3*sqrt(2)^(k-2) - 4)
#       = 7*sqrt(2)^(k-2) + sqrt(2)*3*sqrt(2)^(k-2) - 4*sqrt(2)
#       = 7*sqrt(2)^(k-2) + 3*sqrt(2)*sqrt(2)^(k-2) - 4*sqrt(2)
#       = (7 + 3*sqrt(2))*sqrt(2)^(k-2) - 4*sqrt(2)
#
#             S[2]=4
#       11--10--7,9--6---5  Z[1]=2         k=4 h=2
#        |       |       |
#   13--12       8       4---3             4 + 2*2 + 4+(2-0) = 14
#    |                       |  S[1]=2     (2+0) = 2
#   14                       2
#    |                       |
#   15---16              0---1  Z[0] = 0
#

# k odd
#            S[h]
#            ----
#   Z[h-1] /     \    middle Z[h]
# S[h-1]  |       \
#          \       \
#                   |  S[h]
#                   |
#             \    /  Z[h-1]
#               --
#              S[h-1]
#
# Hb[k] = 2*S[h] + 2*S[h-1]  + sqrt(2)*( Z[h]/2 + Z[h-1] + Z[h]/2 + S[h]-S[h-1] )
#       = 2*S[h] + 2*S[h-1]  + sqrt(2)*( Z[h]   + Z[h-1]          + S[h]-S[h-1] )
#       = 2*2^h  + 2*2^(h-1) + sqrt(2)*( 2*2^h-2 + 2*2^(h-1)-2 + 2^h - 2^(h-1) )
#       = 3*2^h              + sqrt(2)*( 3*2^h                 + 2^(h-1)       - 4 )
#       = 3*2^h + sqrt(2)*( 7*2^(h-1) - 4 )

sub _UNDOCUMENTED_level_to_hull_boundary {
  my ($self, $level) = @_;
  my ($a, $b) = $self->_UNDOCUMENTED_level_to_hull_boundary_sqrt2($level)
    or return undef;
  return $a + $b*sqrt(2);
}
sub _UNDOCUMENTED_level_to_hull_boundary_sqrt2 {
  my ($self, $level) = @_;
  if ($level <= 2) {
    if ($level < 0) { return; }
    if ($level == 2) { return (6,0); }
    return (2, ($level == 0 ? 0 : 1));
  }

  my ($h, $rem) = _divrem($level, 2);
  my $pow = 2**($h-1);

  if ($rem) {
    return (6*$pow, 7*$pow-4);

    # return (2*S_formula($h) + 2*S_formula($h-1),
    #         Z_formula($h)/2 + Z_formula($h-1)
    #         + Z_formula($h)/2 + (S_formula($h)-S_formula($h-1)) );

  } else {
    return (7*$pow, 3*$pow-4);

    # return (S_formula($h) + 2*S_formula($h-1) + S_formula($h)+(Z_formula($h-1)-Z_formula($h-2)),
    #         (Z_formula($h-1) + Z_formula($h-2)));
  }
}

#------------------------------------------------------------------------------
{
  my @_UNDOCUMENTED_level_to_hull_area = (0, 1/2, 2);

  sub _UNDOCUMENTED_level_to_hull_area {
    my ($self, $level) = @_;

    if ($level < 3) {
      if ($level < 0) { return undef; }
      return $_UNDOCUMENTED_level_to_hull_area[$level];
    }
    my ($h, $rem) = _divrem($level, 2);
    return 35*2**($level-4) - ($rem ? 13 : 10)*2**($h-1) + 2;

    #   if ($rem) {
    #     return 35*2**($level-4) - 13*$pow + 2;
    #
    #     my $width = S_formula($h) + Z_formula($h)/2 + Z_formula($h-1)/2;
    #     my $ul = Z_formula($h-1)/2;
    #     my $ur = Z_formula($h)/2;
    #     my $bl = $width - Z_formula($h-1)/2 - S_formula($h-1);
    #     my $br = Z_formula($h-1)/2;
    #     return $width**2 - $ul**2/2 - $ur**2/2 - $bl**2/2 - $br**2/2;
    #
    #   } else {
    #     return 35*2**($level-4) - 10*$pow + 2;
    #     return 0;
    #     return 35*2**($level-4) - 5*2**$h + 2;
    #
    #     # my $width = S_formula($h) + Z_formula($h-1);
    #     # my $upper = Z_formula($h-1)/2;
    #     # my $lower = Z_formula($h-2)/2;
    #     # my $height = S_formula($h-1) + $upper + $lower;
    #     # return $width; # * $height - $upper*$upper - $lower*$lower;
    #   }
    # }
  }
}

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

sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, 2**$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, 2);
  return $exp;
}

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