File Coverage

blib/lib/Math/PlanePath/DragonCurve.pm

Criterion	Covered	Total	%
statement	97	225	43.1
branch	18	82	21.9
condition	11	16	68.7
subroutine	19	35	54.2
pod	11	11	100.0
total	156	369	42.2

# 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 .


# math-image --path=DragonCurve --lines --scale=20
# math-image --path=DragonCurve --all --scale=10
# math-image --path=DragonCurve --output=numbers_dash
#
# math-image --wx --path=DragonCurve,arms=4 --expression='i<16384?i:0'


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

use vars '$VERSION', '@ISA';
$VERSION = 127;
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',
  'bit_split_lowtohigh',
  'digit_split_lowtohigh';
*_divrem_mutate = \&Math::PlanePath::_divrem_mutate;

use Math::PlanePath::DragonMidpoint;

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



use constant n_start => 0;

use constant parameter_info_array => [ { name      => 'arms',
                                         share_key => 'arms_4',
                                         display   => 'Arms',
                                         type      => 'integer',
                                         minimum   => 1,
                                         maximum   => 4,
                                         default   => 1,
                                         width     => 1,
                                         description => 'Arms',
                                       } ];

{
  my @x_negative_at_n = (undef, 5,5,5,6);
  sub x_negative_at_n {
    my ($self) = @_;
    return $x_negative_at_n[$self->{'arms'}];
  }
}
{
  my @y_negative_at_n = (undef, 14,11,8,7);
  sub y_negative_at_n {
    my ($self) = @_;
    return $y_negative_at_n[$self->{'arms'}];
  }
}
{
  my @_UNDOCUMENTED__dxdy_list_at_n = (undef, 5, 5, 5, 3);
  sub _UNDOCUMENTED__dxdy_list_at_n {
    my ($self) = @_;
    return $_UNDOCUMENTED__dxdy_list_at_n[$self->{'arms'}];
  }
}

use constant turn_any_straight => 0; # never straight


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

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

{
  # sub state_string {
  #   my ($state) = @_;
  #   my $digit = $state & 3;  $state >>= 2;
  #   my $rot = $state & 3;  $state >>= 2;
  #   my $rev = $state & 1;  $state >>= 1;
  #   return "rot=$rot rev=$rev (digit=$digit)";
  # }

  # generated by tools/dragon-curve-table.pl
  # next_state length 32
  my @next_state = (12,16, 4,16,  0,20, 8,20,  4,24,12,24,  8,28, 0,28,
                    0,20, 0,28,  4,24, 4,16,  8,28, 8,20, 12,16,12,24);
  my @digit_to_x = ( 0, 0, 1, 1,  0, 1, 1, 0,  0, 0,-1,-1,  0,-1,-1, 0,
                     0, 1, 1, 2,  0, 0,-1,-1,  0,-1,-1,-2,  0, 0, 1, 1);
  my @digit_to_y = ( 0,-1,-1, 0,  0, 0, 1, 1,  0, 1, 1, 0,  0, 0,-1,-1,
                     0, 0, 1, 1,  0, 1, 1, 2,  0, 0,-1,-1,  0,-1,-1,-2);
  my @digit_to_dxdy = ( 1, 0,undef,undef,  0, 1,undef,undef, -1, 0,undef,undef,  0,-1,undef,undef,
                        1, 0,undef,undef,  0, 1,undef,undef, -1, 0,undef,undef,  0,-1);

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

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

    my $int = int($n);      # integer part
    $n -= $int;             # $n = fraction part
    my $zero = ($int * 0);  # inherit bignum 0

    my $arm = _divrem_mutate ($int, $self->{'arms'});
    my @digits = digit_split_lowtohigh($int,4);
    ### @digits

    # initial state from rotation by arm and number of digits
    my $state = ((scalar(@digits) + $arm) & 3) << 2;

    my $len = (2+$zero) ** $#digits;
    my $x = $zero;
    my $y = $zero;
    foreach my $digit (reverse @digits) {  # high to low
      ### at: "x=$x,y=$y  len=$len digit=$digit state=$state"
      # ### state is: state_string($state)

      $state += $digit;
      $x += $len * $digit_to_x[$state];
      $y += $len * $digit_to_y[$state];
      $state = $next_state[$state];
      $len /= 2;
    }

    ### final: "x=$x y=$y  state=$state"
    # ### state is: state_string($state)
    ### final: "frac dx=$digit_to_dxdy[$state], dy=$digit_to_dxdy[$state+1]"

    return ($n * $digit_to_dxdy[$state] + $x,
            $n * $digit_to_dxdy[$state+1] + $y);
  }
}


{
  # generated by tools/dragon-curve-dxdy.pl
  # next_state length 32
  my @next_state = ( 0, 6,20, 2,  4,10,24, 6,  8,14,28,10, 12, 2,16,14,
                     0,22,20,18,  4,26,24,22,  8,30,28,26, 12,18,16,30);
  my @state_to_dxdy = ( 1, 0,-1, 1,  0, 1,-1,-1, -1, 0, 1,-1,  0,-1, 1, 1,
                        1, 0,-1,-1,  0, 1, 1,-1, -1, 0, 1, 1,  0,-1,-1, 1);

  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
    ### $int
    ### $n

    my $state = 4 * _divrem_mutate ($int, $self->{'arms'});
    ### arm as initial state: $state

    foreach my $bit (reverse bit_split_lowtohigh($int)) {  # high to low
      $state = $next_state[$state + $bit];
    }
    $state &= 0x1C;  # mask out "prevbit" from state, leaving state==0 mod 4

    ### final state: $state
    ### dx: $state_to_dxdy[$state]
    ### dy: $state_to_dxdy[$state+1],
    ### frac dx: $state_to_dxdy[$state+2],
    ### frac dy: $state_to_dxdy[$state+3],

    return ($state_to_dxdy[$state]   + $n * $state_to_dxdy[$state+2],
            $state_to_dxdy[$state+1] + $n * $state_to_dxdy[$state+3]);
  }
}

# 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
#
my @try_dx = (0,0,-1,-1);
my @try_dy = (0,1,1,0);

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

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

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

  if ($x == 0 && $y == 0) {
    return (0 .. $self->arms_count - 1);
  }

  my @n_list;
  my $xm = $x+$y;  # rotate -45 and mul sqrt(2)
  my $ym = $y-$x;
  foreach my $dx (0,-1) {
    foreach my $dy (0,1) {
      my $t = $self->Math::PlanePath::DragonMidpoint::xy_to_n
        ($xm+$dx, $ym+$dy);
      next unless defined $t;

      my ($tx,$ty) = $self->n_to_xy($t)
        or next;

      if ($tx == $x && $ty == $y) {
        ### found: $t
        if (@n_list && $t < $n_list[0]) {
          unshift @n_list, $t;
        } else {
          push @n_list, $t;
        }
        if (@n_list == 2) {
          return @n_list;
        }
      }
    }
  }
  return @n_list;
}

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

sub xy_is_visited {
  my ($self, $x, $y) = @_;

  my $arms_count = $self->{'arms'};
  if ($arms_count == 4) {
    # yes, whole plane visited
    return 1;
  }

  my $xm = $x+$y;
  my $ym = $y-$x;
  {
    my $arm = Math::PlanePath::DragonMidpoint::_xy_to_arm($xm,$ym);
    if ($arm < $arms_count) {
      # yes, segment $xm,$ym is on the desired arms
      return 1;
    }
    if ($arm == 2 && $arms_count == 1) {
      # no, segment $xm,$ym is on arm 2, which means its opposite is only on
      # arm 1,2,3 not arm 0 so arms_count==1 cannot be visited
      return 0;
    }
  }
  return (Math::PlanePath::DragonMidpoint::_xy_to_arm($xm-1,$ym+1)
          < $arms_count);
}


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

# 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
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### DragonCurve rect_to_n_range(): "$x1,$y1  $x2,$y2"
  my $xmax = int(max(abs($x1),abs($x2)));
  my $ymax = int(max(abs($y1),abs($y2)));
  return (0,
          $self->{'arms'} * ($xmax*$xmax + $ymax*$ymax + 1) * 7);
}

# Not quite right yet ...
#
# sub rect_to_n_range {
#   my ($self, $x1,$y1, $x2,$y2) = @_;
#   ### DragonCurve rect_to_n_range(): "$x1,$y1  $x2,$y2"
#
#
#    my ($length, $level_limit) = round_down_pow
#      ((max(abs($x1),abs($x2))**2 + max(abs($y1),abs($y2))**2 + 1) * 7,
#       2);
#    $level_limit += 2;
#    ### $level_limit
#
#    if (is_infinite($level_limit)) {
#      return ($level_limit,$level_limit);
#    }
#
#    $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;
#    ### sorted range: "$x1,$y1  $x2,$y2"
#
#
#    my @xend = (0, 1);
#    my @yend = (0, 0);
#    my @xmin = (0, 0);
#    my @xmax = (0, 1);
#    my @ymin = (0, 0);
#    my @ymax = (0, 0);
#    my @sidemax = (0, 1);
#    my $extend = sub {
#      my ($i) = @_;
#      ### extend(): $i
#      while ($i >= $#xend) {
#        ### extend from: $#xend
#        my $xend = $xend[-1];
#        my $yend = $yend[-1];
#        ($xend,$yend) = ($xend-$yend,  # rotate +45
#                         $xend+$yend);
#        push @xend, $xend;
#        push @yend, $yend;
#        my $xmax = $xmax[-1];
#        my $xmin = $xmin[-1];
#        my $ymax = $ymax[-1];
#        my $ymin = $ymin[-1];
#        ### assert: $xmax >= $xmin
#        ### assert: $ymax >= $ymin
#
#        #    ### at: "end=$xend,$yend   $xmin..$xmax  $ymin..$ymax"
#        push @xmax, max($xmax, $xend + $ymax);
#        push @xmin, min($xmin, $xend + $ymin);
#
#        push @ymax, max($ymax, $yend - $xmin);
#        push @ymin, min($ymin, $yend - $xmax);
#
#        push @sidemax, max ($xmax[-1], -$xmin[-1],
#                             $ymax[-1], -$ymin[-1],
#                             abs($xend),
#                             abs($yend));
#      }
#      ### @sidemax
#    };
#
#    my $rect_dist = sub {
#      my ($x,$y) = @_;
#      my $xd = ($x < $x1 ? $x1 - $x
#                : $x > $x2 ? $x - $x2
#                : 0);
#      my $yd = ($y < $y1 ? $y1 - $y
#                : $y > $y2 ? $y - $y2
#                : 0);
#      return max($xd,$yd);
#    };
#
#    my $arms = $self->{'arms'};
#    ### $arms
#    my $n_lo;
#    {
#      my $top = 0;
#      for (;;) {
#      ARM_LO: foreach my $arm (0 .. $arms-1) {
#          my $i = 0;
#          my @digits;
#          if ($top > 0) {
#            @digits = ((0)x($top-1), 1);
#          } else {
#            @digits = (0);
#          }
#
#          for (;;) {
#            my $n = 0;
#            foreach my $digit (reverse @digits) { # high to low
#              $n = 2*$n + $digit;
#            }
#            $n = $n*$arms + $arm;
#            my ($nx,$ny) = $self->n_to_xy($n);
#            my $nh = &$rect_dist ($nx,$ny);
#
#            ### lo consider: "i=$i  digits=".join(',',reverse @digits)."  is n=$n xy=$nx,$ny nh=$nh"
#
#            if ($i == 0 && $nh == 0) {
#              ### lo found inside: $n
#              if (! defined $n_lo || $n < $n_lo) {
#                $n_lo = $n;
#              }
#              next ARM_LO;
#            }
#
#            if ($i == 0 || $nh > $sidemax[$i+2]) {
#              ### too far away: "nxy=$nx,$ny   nh=$nh vs ".$sidemax[$i+2]." at i=$i"
#
#              while (++$digits[$i] > 1) {
#                $digits[$i] = 0;
#                if (++$i <= $top) {
#                  ### backtrack up ...
#                } else {
#                  ### not found within this top and arm, next arm ...
#                  next ARM_LO;
#                }
#              }
#            } else {
#              ### lo descend ...
#              ### assert: $i > 0
#              $i--;
#              $digits[$i] = 0;
#            }
#          }
#        }
#
#        # if an $n_lo was found on any arm within this $top then done
#        if (defined $n_lo) {
#          last;
#        }
#
#        ### lo extend top ...
#        if (++$top > $level_limit) {
#          ### nothing below level limit ...
#          return (1,0);
#        }
#        &$extend($top+3);
#      }
#    }
#
#    my $n_hi = 0;
#   ARM_HI: foreach my $arm (reverse 0 .. $arms-1) {
#      &$extend($level_limit+2);
#      my @digits = ((1) x $level_limit);
#      my $i = $#digits;
#      for (;;) {
#        my $n = 0;
#        foreach my $digit (reverse @digits) { # high to low
#          $n = 2*$n + $digit;
#        }
#
#        $n = $n*$arms + $arm;
#        my ($nx,$ny) = $self->n_to_xy($n);
#        my $nh = &$rect_dist ($nx,$ny);
#
#        ### hi consider: "arm=$arm  i=$i  digits=".join(',',reverse @digits)."  is n=$n xy=$nx,$ny nh=$nh"
#
#        if ($i == 0 && $nh == 0) {
#          ### hi found inside: $n
#          if ($n > $n_hi) {
#            $n_hi = $n;
#            next ARM_HI;
#          }
#        }
#
#        if ($i == 0 || $nh > $sidemax[$i+2]) {
#          ### too far away: "$nx,$ny   nh=$nh vs ".$sidemax[$i+2]." at i=$i"
#
#          while (--$digits[$i] < 0) {
#            $digits[$i] = 1;
#            if (++$i < $level_limit) {
#              ### hi backtrack up ...
#            } else {
#              ### hi nothing within level limit for this arm ...
#              next ARM_HI;
#            }
#          }
#
#        } else {
#          ### hi descend
#          ### assert: $i > 0
#          $i--;
#          $digits[$i] = 1;
#        }
#      }
#    }
#
#    if ($n_hi == 0) {
#      ### oops, lo found but hi not found
#      $n_hi = $n_lo;
#    }
#
#    return ($n_lo, $n_hi);
# }


#------------------------------------------------------------------------------
# level ranges

#           arms=1       arms=2      arms=4
# level 0  0..1 = 2    0..3  = 4    0..7  = 8
# level 1  0..2 = 3    0..5  = 6    0..11 = 12
# level 2  0..4 = 5    0..9  = 10   0..19 = 20
# level 3  0..8 = 9    0..17 = 18   0..35 = 36
#            2^k       2*2^k+1      4*2^k+3
#
sub level_to_n_range {
  my ($self, $level) = @_;
  return (0, (2**$level + 1) * $self->{'arms'} - 1);
}
# 0 .. 2^level
# -1 .. 2^level-1
# level = round_up_pow(N)
# eg N=13 -> 2^4=16 level=4
#
sub n_to_level {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  if (is_infinite($n)) { return $n; }
  $n = round_nearest($n);
  _divrem_mutate ($n, $self->{'arms'});
  my ($pow, $exp) = round_up_pow ($n, 2);
  return $exp;
}

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

{
  my @_UNDOCUMENTED_level_to_left_line_boundary = (1,2,4);
  sub _UNDOCUMENTED_level_to_left_line_boundary {
    my ($self, $level) = @_;
    if ($level < 0) { return undef; }
    if ($level <= 2) { return $_UNDOCUMENTED_level_to_left_line_boundary[$level]; }
    if (is_infinite($level)) { return $level; }

    my $l0 = 2;
    my $l1 = 4;
    my $l2 = 8;
    foreach (4 .. $level) {
      ($l2,$l1,$l0) = ($l2 + 2*$l0, $l2, $l1);
    }
    return $l2;
  }
}

{
  my @level_to_right_line_boundary = (1,2,4,8,undef);
  sub _UNDOCUMENTED_level_to_right_line_boundary {
    my ($self, $level) = @_;
    if ($level < 0) { return undef; }
    if ($level <= 3) { return $level_to_right_line_boundary[$level]; }
    if (is_infinite($level)) { return $level; }

    my $r0 =  2;
    my $r1 =  4;
    my $r2 =  8;
    my $r3 = 16;
    foreach (5 .. $level) {
      ($r3,$r2,$r1,$r0) = (2*$r3 - $r2 + 2*$r1 - 2*$r0,  $r3, $r2, $r1);
    }
    return $r3;
  }
}
sub _UNDOCUMENTED_level_to_line_boundary {
  my ($self, $level) = @_;
  if ($level < 0) { return undef; }
  return $self->_UNDOCUMENTED_level_to_right_line_boundary($level+1);
}

sub _UNDOCUMENTED_level_to_u_left_line_boundary {
  my ($self, $level) = @_;
  if ($level < 0) { return undef; }
  return ($level == 0 ? 3
          : $self->_UNDOCUMENTED_level_to_right_line_boundary($level) + 4);
}
sub _UNDOCUMENTED_level_to_u_right_line_boundary {
  my ($self, $level) = @_;
  if ($level < 0) { return undef; }
  return ($self->_UNDOCUMENTED_level_to_right_line_boundary($level)
          + $self->_UNDOCUMENTED_level_to_right_line_boundary($level+1));
}
sub _UNDOCUMENTED_level_to_u_line_boundary {
  my ($self, $level) = @_;
  if ($level < 0) { return undef; }
  return ($self->_UNDOCUMENTED_level_to_u_left_line_boundary($level)
          + $self->_UNDOCUMENTED_level_to_u_right_line_boundary($level));
}

sub _UNDOCUMENTED_level_to_enclosed_area {
  my ($self, $level) = @_;
  # A[k] = 2^(k-1) - B[k]/4
  if ($level < 0) { return undef; }
  if ($level == 0) { return 0; } # avoid 2**(-1)
  return 2**($level-1) - $self->_UNDOCUMENTED_level_to_line_boundary($level) / 4;
}
*_UNDOCUMENTED_level_to_doubled_points = \&_UNDOCUMENTED_level_to_enclosed_area;

{
  my @_UNDOCUMENTED_level_to_single_points = (2,3,5);
  sub _UNDOCUMENTED_level_to_single_points {
    my ($self, $level) = @_;
    if ($level < 0) { return undef; }
    if ($level <= 2) { return $_UNDOCUMENTED_level_to_single_points[$level]; }
    if (is_infinite($level)) { return $level; }

    my $l0 = 3;
    my $l1 = 5;
    my $l2 = 9;
    foreach (4 .. $level) {
      ($l2,$l1,$l0) = ($l2 + 2*$l0, $l2, $l1);
    }
    return $l2;
  }
}

{
  my @_UNDOCUMENTED_level_to_enclosed_area_join = (0,0,0,1);
  sub _UNDOCUMENTED_level_to_enclosed_area_join {
    my ($self, $level) = @_;
    if ($level < 0) { return undef; }
    if ($level <= 3) { return $_UNDOCUMENTED_level_to_enclosed_area_join[$level]; }
    if (is_infinite($level)) { return $level; }

    my ($j0,$j1,$j2,$j3) = @_UNDOCUMENTED_level_to_enclosed_area_join;
    $j3 += $level*0;
    foreach (4 .. $level) {
      ($j3,$j2,$j1,$j0) = (2*$j3 - $j2 + 2*$j1 - 2*$j0,  $j3, $j2, $j1);
    }
    return $j3;
  }
}

#------------------------------------------------------------------------------
# points visited

{
  my @_UNDOCUMENTED_level_to_visited = (2, 3, 5, 9, 16);
  sub _UNDOCUMENTED_level_to_visited {
    my ($self, $level) = @_;

    if ($level < 0) { return undef; }
    if ($level <= $#_UNDOCUMENTED_level_to_visited) { return $_UNDOCUMENTED_level_to_visited[$level]; }
    if (is_infinite($level)) { return $level; }

    my ($p0,$p1,$p2,$p3,$p4) = @_UNDOCUMENTED_level_to_visited;
    foreach (5 .. $level) {
      ($p4,$p3,$p2,$p1,$p0) = (4*$p4 - 5*$p3 + 4*$p2 - 6*$p1 + 4*$p0,  $p4, $p3, $p2, $p1);
    }
    return $p4;
  }
}

#------------------------------------------------------------------------------
{
  my @_UNDOCUMENTED__n_segment_is_right_boundary
    #         R M A B C D F G H
    #         1 2 3 4 5 6 7 8 9
    = ([undef,1,3,1,6,7,9,3     ],
       [undef,2,4,5,4,8,5,0,0,4 ]);

  sub _UNDOCUMENTED__n_segment_is_right_boundary {
    my ($self, $n) = @_;
    if (is_infinite($n)) { return 0; }
    unless ($n >= 0) { return 0; }
    $n = int($n);

    my $state = 1;
    foreach my $bit (reverse bit_split_lowtohigh($n)) { # high to low
      $state = $_UNDOCUMENTED__n_segment_is_right_boundary[$bit][$state]
        || return 0;
    }
    return 1;
  }
}

1;
__END__