File Coverage

blib/lib/Math/PlanePath/SierpinskiCurve.pm

Criterion	Covered	Total	%
statement	159	236	67.3
branch	56	102	54.9
condition	27	35	77.1
subroutine	23	38	60.5
pod	18	18	100.0
total	283	429	65.9

# Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 .


# Sierpinski, "Sur Une Nouvelle Courbe Qui Remplit Toute Une Aire Plaine",
# Bull. Acad. Sci. Cracovie Series A, 1912, pages 462-478.


package Math::PlanePath::SierpinskiCurve;
use 5.004;
use strict;
use List::Util 'sum','first';
#use List::Util 'min','max';
*min = \&Math::PlanePath::_min;
*max = \&Math::PlanePath::_max;

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_up_pow',
  'round_down_pow',
  'digit_split_lowtohigh';

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


use constant n_start => 0;

sub x_negative {
  my ($self) = @_;
  return ($self->{'arms'} >= 3);
}
sub y_negative {
  my ($self) = @_;
  return ($self->{'arms'} >= 5);
}

use constant parameter_info_array =>
  [
   { name        => 'arms',
     share_key   => 'arms_8',
     display     => 'Arms',
     type        => 'integer',
     minimum     => 1,
     maximum     => 8,
     default     => 1,
     width       => 1,
     description => 'Arms',
   },

   { name        => 'straight_spacing',
     display     => 'Straight Spacing',
     type        => 'integer',
     minimum     => 1,
     default     => 1,
     width       => 1,
     description => 'Spacing of the straight line points.',
   },
   { name        => 'diagonal_spacing',
     display     => 'Diagonal Spacing',
     type        => 'integer',
     minimum     => 1,
     default     => 1,
     width       => 1,
     description => 'Spacing of the diagonal points.',
   },
  ];

# Ntop = (4^level)/2 - 1
# Xtop = 3*2^(level-1) - 1
# fill = Ntop / (Xtop*(Xtop-1)/2)
#      -> 2 * ((4^level)/2 - 1) / (3*2^(level-1) - 1)^2
#      -> 2 * ((4^level)/2) / (3*2^(level-1))^2
#      =  4^level / (9*4^(level-1)
#      =  4/9 = 0.444

sub x_negative_at_n {
  my ($self) = @_;
  return $self->arms_count >= 3 ? 2 : undef;
}
sub y_negative_at_n {
  my ($self) = @_;
  return $self->arms_count >= 5 ? 4 : undef;
}

{
  # Note: shared by Math::PlanePath::SierpinskiCurveStair
  my @x_minimum = (undef,
                   1,  # 1 arm
                   0,  # 2 arms
                  );   # more than 2 arm, X goes negative
  sub x_minimum {
    my ($self) = @_;
    return $x_minimum[$self->arms_count];
  }
}
{
  # Note: shared by Math::PlanePath::SierpinskiCurveStair
  my @sumxy_minimum = (undef,
                       1,  # 1 arm, octant and X>=1 so X+Y>=1
                       1,  # 2 arms, X>=1 or Y>=1 so X+Y>=1
                       0,  # 3 arms, Y>=1 and X>=Y, so X+Y>=0
                      );   # more than 3 arm, Sum goes negative so undef
  sub sumxy_minimum {
    my ($self) = @_;
    return $sumxy_minimum[$self->arms_count];
  }
}
use constant sumabsxy_minimum => 1;

# Note: shared by Math::PlanePath::SierpinskiCurveStair
#                 Math::PlanePath::AlternatePaper
#                 Math::PlanePath::AlternatePaperMidpoint
sub diffxy_minimum {
  my ($self) = @_;
  return ($self->arms_count == 1
          ? 1       # octant Y<=X-1 so X-Y>=1
          : undef); # more than 1 arm, DiffXY goes negative
}
use constant absdiffxy_minimum => 1; # X=Y never occurs
use constant rsquared_minimum => 1; # minimum X=1,Y=0

sub dx_minimum {
  my ($self) = @_;
  return - max($self->{'straight_spacing'},
               $self->{'diagonal_spacing'});
}
*dy_minimum = \&dx_minimum;

sub dx_maximum {
  my ($self) = @_;
  return max($self->{'straight_spacing'},
             $self->{'diagonal_spacing'});
}
*dy_maximum = \&dx_maximum;

sub _UNDOCUMENTED__dxdy_list {
  my ($self) = @_;
  my $s = $self->{'straight_spacing'};
  my $d = $self->{'diagonal_spacing'};
  return ($s,0,                    # E     eight scaled
          ($d ? ( $d, $d) : ()),   # NE    except s=0
          ($s ? (  0, $s) : ()),   # N     or d=0 skips
          ($d ? (-$d, $d) : ()),   # NW
          ($s ? (-$s,  0) : ()),   # W
          ($d ? (-$d,-$d) : ()),   # SW
          ($s ? (  0,-$s) : ()),   # S
          ($d ? ( $d,-$d) : ()));  # SE

}
{
  my @_UNDOCUMENTED__dxdy_list_at_n = (undef,
                                       21, 20, 27, 36,
                                       29, 12, 12, 13);
  sub _UNDOCUMENTED__dxdy_list_at_n {
    my ($self) = @_;
    return $_UNDOCUMENTED__dxdy_list_at_n[$self->{'arms'}];
  }
}
sub dsumxy_minimum {
  my ($self) = @_;
  return - max($self->{'straight_spacing'},
               2*$self->{'diagonal_spacing'});
}
sub dsumxy_maximum {
  my ($self) = @_;
  return max($self->{'straight_spacing'},
             2*$self->{'diagonal_spacing'});
}
*ddiffxy_minimum = \&dsumxy_minimum;
*ddiffxy_maximum = \&dsumxy_maximum;

use constant dir_maximum_dxdy => (1,-1); # South-East
use constant turn_any_straight => 0; # never straight


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

sub new {
  my $self = shift->SUPER::new(@_);

  $self->{'arms'} = max(1, min(8, $self->{'arms'} || 1));
  $self->{'straight_spacing'} ||= 1;
  $self->{'diagonal_spacing'} ||= 1;
  return $self;
}

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

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

  my $int = int($n); # BigFloat int() gives BigInt, use that
  $n -= $int;   # preserve possible BigFloat
  ### $int
  ### $n

  my $arm = _divrem_mutate ($int, $self->{'arms'});

  my $s = $self->{'straight_spacing'};
  my $d = $self->{'diagonal_spacing'};
  my $base = 2*$d+$s;
  my $x = my $y = ($int * 0);  # inherit big 0
  my $len = $x + $base;      # inherit big

  foreach my $digit (digit_split_lowtohigh($int,4)) {
    ### at: "$x,$y  digit=$digit"

    if ($digit == 0) {
      $x = $n*$d + $x;
      $y = $n*$d + $y;
      $n = 0;

    } elsif ($digit == 1) {
      ($x,$y) = ($n*$s - $y + $len-$d-$s,   # rotate +90
                 $x + $d);
      $n = 0;

    } elsif ($digit == 2) {
      # rotate -90
      ($x,$y) = ($n*$d + $y  + $len-$d,
                 -$n*$d - $x + $len-$d-$s);
      $n = 0;

    } else { # digit==3
      $x += $len;
    }
    $len *= 2;
  }

  # n=0 or n=33..33
  $x = $n*$d + $x;
  $y = $n*$d + $y;

  $x += 1;
  if ($arm & 1) {
    ($x,$y) = ($y,$x);   # mirror 45
  }
  if ($arm & 2) {
    ($x,$y) = (-1-$y,$x);   # rotate +90
  }
  if ($arm & 4) {
    $x = -1-$x;   # rotate 180
    $y = -1-$y;
  }

  # use POSIX 'floor';
  # $x += floor($x/3);
  # $y += floor($y/3);

  # $x += floor(($x-1)/3) + floor(($x-2)/3);
  # $y += floor(($y-1)/3) + floor(($y-2)/3);


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

my @digit_to_dir = (0, -2, 2, 0);
my @dir8_to_dx = (1, 1, 0,-1, -1, -1,  0, 1);
my @dir8_to_dy = (0, 1, 1, 1,  0, -1, -1,-1);
my @digit_to_nextturn = (-1,   # after digit=0
                         2,    #       digit=1
                         -1);  #       digit=2
sub n_to_dxdy {
  my ($self, $n) = @_;
  ### n_to_dxdy(): $n

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

  my $int = int($n);
  $n -= $int;

  my $arm = _divrem_mutate($int,$self->{'arms'});
  my $lowbit = _divrem_mutate($int,2);
  ### $lowbit
  ### $int

  if (is_infinite($int)) {
    return ($int,$int);
  }
  my @ndigits = digit_split_lowtohigh($int,4);
  ### @ndigits

  my $dir8 = sum(0, map {$digit_to_dir[$_]} @ndigits);
  if ($arm & 1) {
    $dir8 = - $dir8;  # mirrored on second,fourth,etc arm
  }
  $dir8 += ($arm|1);  # NE,NW,SW, or SE

  my $turn;
  if ($n || $lowbit) {
    # next turn

    # lowest non-3 digit, or zero if all 3s (implicit 0 above high digit)
    $turn = $digit_to_nextturn[ first {$_!=3} @ndigits, 0 ];
    if ($arm & 1) {
      $turn = - $turn;  # mirrored on second,fourth,etc arm
    }
  }

  if ($lowbit) {
    $dir8 += $turn;
  }

  my $s = $self->{'straight_spacing'};
  my $d = $self->{'diagonal_spacing'};

  $dir8 &= 7;
  my $spacing = ($dir8 & 1 ? $d : $s);
  my $dx = $spacing * $dir8_to_dx[$dir8];
  my $dy = $spacing * $dir8_to_dy[$dir8];

  if ($n) {
    $dir8 += $turn;
    $dir8 &= 7;
    $spacing = ($dir8 & 1 ? $d : $s);
    $dx += $n*($spacing * $dir8_to_dx[$dir8]
               - $dx);
    $dy += $n*($spacing * $dir8_to_dy[$dir8]
               - $dy);
  }

  return ($dx, $dy);
}

# 2| . 3 .
# 1| 1 . 2
# 0| . 0 .
#  +------
#    0 1 2
#
# 4| . . . 3 .          # diagonal_spacing == 3
# 3| . . . . 2 4        # mod=2*3+1=7
# 2| . . . . . . .
# 1| 1 . . . . . . .
# 0| . 0 . . . . . . 6
#  +------------------
#    0 1 2 3 4 5 6 7 8
#
sub _NOTWORKING__xy_is_visited {
  my ($self, $x, $y) = @_;
  $x = round_nearest($x);
  $y = round_nearest($y);
  my $mod = 2*$self->{'diagonal_spacing'} + $self->{'straight_spacing'};
  return (_rect_within_arms($x,$y, $x,$y, $self->{'arms'})
          && ((($x%$mod)+($y%$mod)) & 1));
}

#   x1    *  x2 *
#    +-----*-+y2*
#    |      *|  *
#    |       *  *
#    |       |* *
#    |       | **
#    +-------+y1*
#   ----------------
#
# arms=5 x1,y2 after X=Y-1 line, so x1 > y2-1, x1 >= y2
# ************
#      x1   *   x2
#      +---*----+y2
#      |  *     |
#      | *      |
#      |*       |
#      *        |
#     *+--------+y1
#    *
#
# arms=7 x1,y1 after X=-2-Y line, so x1 > -2-y1
# ************
# ** +------+
# * *|      |
# *  *      |
# *  |*     |
# *  | *    |
# *y1+--*---+
# * x1   *
#
# _rect_within_arms() returns true if rectangle x1,y1,x2,y2 has some part
# within the extent of the $arms set of octants.
#
sub _rect_within_arms {
  my ($x1,$y1, $x2,$y2, $arms) = @_;
  return ($arms <= 4
          ? ($y2 >= 0  # y2 top edge must be positive
             && ($arms <= 2
                 ? ($arms == 1 ? $x2 > $y1   # arms==1  bottom right
                    :            $x2 >= 0)   # arms==2  right edge
                 : ($arms == 4               # arms==4  anything
                    || $x2 >= -$y2)))        # arms==3  top right

          # arms >= 5
          : ($y2 >= 0  # y2 top edge positive is good, otherwise check
             || ($arms <= 6
                 ? ($arms == 5 ? $x1 < $y2   # arms==5  top left
                    :            $x1 < 0)    # arms==6  left edge
                 : ($arms == 8               # arms==8  anything
                    || $x1 <= -2-$y1))));    # arms==7  bottom left
}

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

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

  my $arm = 0;
  if ($y < 0) {
    $arm = 4;
    $x = -1-$x;  # rotate -180
    $y = -1-$y;
  }
  if ($x < 0) {
    $arm += 2;
    ($x,$y) = ($y, -1-$x);  # rotate -90
  }
  if ($y > $x) {       # second octant
    $arm++;
    ($x,$y) = ($y,$x); # mirror 45
  }

  my $arms = $self->{'arms'};
  if ($arm >= $arms) {
    return undef;
  }

  $x -= 1;
  if ($x < 0 || $x < $y) {
    return undef;
  }
  ### x adjust to zero: "$x,$y"
  ### assert: $x >= 0
  ### assert: $y >= 0

  my $s = $self->{'straight_spacing'};
  my $d = $self->{'diagonal_spacing'};
  my $base = (2*$d+$s);
  my ($len,$level) = round_down_pow (($x+$y)/$base || 1,  2);
  ### $level
  ### $len
  if (is_infinite($level)) {
    return $level;
  }

  # Xtop = 3*2^(level-1)-1
  #
  $len *= 2*$base;
  ### initial len: $len

  my $n = 0;
  foreach (0 .. $level) {
    $n *= 4;
    ### at: "loop=$_ len=$len   x=$x,y=$y  n=$n"
    ### assert: $x >= 0
    ### assert: $y >= 0

    my $len_sub_d = $len - $d;
    if ($x < $len_sub_d) {
      ### digit 0 or 1...
      if ($x+$y+$s < $len) {
        ### digit 0 ...
      } else {
        ### digit 1 ...
        ($x,$y) = ($y-$d, $len-$s-$d-$x);   # shift then rotate -90
        $n += 1;
      }
    } else {
      $x -= $len_sub_d;
      ### digit 2 or 3 to: "x=$x y=$y"
      if ($x < $y) {   # before diagonal
        ### digit 2...
        ($x,$y) = ($len-$d-$s-$y, $x);     # shift y-len then rotate +90
        $n += 2;
      } else {
        #### digit 3...
        $x -= $d;
        $n += 3;
      }
      if ($x < 0) {
        return undef;
      }
    }
    $len /= 2;
  }

  ### end at: "x=$x,y=$y   n=$n"
  ### assert: $x >= 0
  ### assert: $y >= 0

  $n *= 4;
  if ($y == 0 && $x == 0) {
    ### final digit 0 ...
  } elsif ($x == $d && $y == $d) {
    ### final digit 1 ...
    $n += 1;
  } elsif ($x == $d+$s && $y == $d) {
    ### final digit 2 ...
    $n += 2;
  } elsif ($x == $base && $y == 0) {
    ### final digit 3 ...
    $n += 3;
  } else {
    return undef;
  }

  return $n*$arms + $arm;
}

# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### SierpinskiCurve 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 $arms = $self->{'arms'};
  unless (_rect_within_arms($x1,$y1, $x2,$y2, $arms)) {
    ### rect outside octants, for arms: $arms
    return (1,0);
  }

  my $max = ($x2 + $y2);
  if ($arms >= 3) {
    _apply_max ($max, -1-$x1 + $y2);

    if ($arms >= 5) {
      _apply_max ($max, -1-$x1 - $y1-1);

      if ($arms >= 7) {
        _apply_max ($max, $x2 - $y1-1);
      }
    }
  }

  # base=2d+s
  # level begins at
  #   base*(2^level-1)-s = X+Y     ... maybe
  #   base*2^level = X+base
  #   2^level = (X+base)/base
  #   level = log2((X+base)/base)
  # then
  #   Nlevel = 4^level-1

  my $base = 2 * $self->{'diagonal_spacing'} + $self->{'straight_spacing'};
  my ($power) = round_down_pow (int(($max+$base-2)/$base),
                                2);
  return (0, 4*$power*$power * $arms - 1);
}

sub _apply_max {
  ### _apply_max(): "$_[0] cf $_[1]"
  unless ($_[0] > $_[1]) {
    $_[0] = $_[1];
  }
}

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

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

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