File Coverage

blib/lib/Math/PlanePath/SacksSpiral.pm

Criterion	Covered	Total	%
statement	94	103	91.2
branch	19	30	63.3
condition	0	3	0.0
subroutine	28	28	100.0
pod	4	4	100.0
total	145	168	86.3

# Copyright 2010, 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 .


# could loop by more or less, eg. 4*n^2 each time like a square spiral
# (Kevin Vicklund at the_surprises_never_eend_the_u.php)


package Math::PlanePath::SacksSpiral;
use 5.004;
use strict;
use Math::Libm 'hypot';
use POSIX 'floor';
#use List::Util 'max';
*max = \&Math::PlanePath::_max;

use Math::PlanePath;
use Math::PlanePath::MultipleRings;

use vars '$VERSION', '@ISA';
$VERSION = 128;
@ISA = ('Math::PlanePath');


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


use constant n_start => 0;
use constant figure => 'circle';
use constant x_negative_at_n => 2;
use constant y_negative_at_n => 3;

use constant 1.02; # for leading underscore
use constant _TWO_PI => 4*atan2(1,0);

# at N=k^2 polygon of 2k+1 sides R=k
# dX -> sin(2pi/(2k+1))*k
#    -> 2pi/(2k+1) * k
#    -> pi
use constant dx_minimum => - 2*atan2(1,0);  # -pi
use constant dx_maximum =>   2*atan2(1,0);  # +pi
use constant dy_minimum => - 2*atan2(1,0);
use constant dy_maximum =>   2*atan2(1,0);

use constant turn_any_right    => 0; # left always
use constant turn_any_straight => 0; # left always


#------------------------------------------------------------------------------
# sub _as_float {
#   my ($x) = @_;
#   if (ref $x) {
#     if ($x->isa('Math::BigInt')) {
#       return Math::BigFloat->new($x);
#     }
#     if ($x->isa('Math::BigRat')) {
#       return $x->as_float;
#     }
#   }
#   return $x;
# }

# Note: this is "use Math::BigFloat" not "require Math::BigFloat" because
# BigFloat 1.997 does some setups in its import() needed to tie-in to the
# BigInt back-end, or something.
use constant::defer _bigfloat => sub {
  eval "use Math::BigFloat; 1" or die $@;
  return "Math::BigFloat";
};

sub n_to_xy {
  my ($self, $n) = @_;
  if ($n < 0) {
    return;
  }
  my $two_pi = _TWO_PI();

  if (ref $n) {
    if ($n->isa('Math::BigInt')) {
      $n = _bigfloat()->new($n);
    }
    if ($n->isa('Math::BigRat')) {
      $n = $n->as_float;
    }
    if ($n->isa('Math::BigFloat')) {
      $two_pi = 2 * Math::BigFloat->bpi ($n->accuracy
                                         || $n->precision
                                         || $n->div_scale);
    }
  }

  my $r = sqrt($n);
  my $theta = $two_pi * ($r - int($r));  # 0 <= $theta < 2*pi
  return ($r * cos($theta),
          $r * sin($theta));

}

sub n_to_rsquared {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  return $n;  # exactly RSquared=$n
}

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

  my $theta_frac = Math::PlanePath::MultipleRings::_xy_to_angle_frac($x,$y);
  ### assert: 0 <= $theta_frac && $theta_frac < 1

  # the nearest arc, integer
  my $s = floor (hypot($x,$y) - $theta_frac + 0.5);

  # the nearest N on the arc
  my $n = floor ($s*$s + $theta_frac * (2*$s + 1) + 0.5);

  # check within 0.5 radius
  my ($nx, $ny) = $self->n_to_xy($n);

  ### $theta_frac
  ### raw hypot: hypot($x,$y)
  ### $s
  ### $n
  ### hypot: hypot($nx-$x, $ny-$y)
  if (hypot($nx-$x,$ny-$y) <= 0.5) {
    return $n;
  } else {
    return undef;
  }
}

# r^2 = x^2 + y^2
# (r+1)^2 = r^2 + 2r + 1
# r < x+y
# (r+1)^2 < x^2+y^2 + x + y + 1
#         < (x+.5)^2 + (y+.5)^2 + 1
# (x+1)^2 + (y+1)^2 = x^2+y^2 + 2x+2y+2
#
# (x+1)^2 + (y+1)^2 - (r+1)^2
#   = x^2+y^2 + 2x+2y+2 - (r^2 + 2r + 1)
#   = x^2+y^2 + 2x+2y+2 - x^2-y^2 - 2*sqrt(x^2+y^2) - 1
#   = 2x+2y+1 - 2*sqrt(x^2+y^2)
#   >= 2x+2y+1 - 2*(x+y)
#   = 1
#
# (x+e)^2 + (y+e)^2 - (r+e)^2
#   = x^2+y^2 + 2xe+2ye + 2e^2 - (r^2 + 2re + e^2)
#   = x^2+y^2 + 2xe+2ye + 2e^2 - x^2-y^2 - 2*e*sqrt(x^2+y^2) - e^2
#   = 2xe+2ye + e^2 - 2*e*sqrt(x^2+y^2)
#   >= 2xe+2ye + e^2 - 2*e*(x+y)
#   = e^2
#
# x+1,y+1 increases the radius by at least 1 thus pushing it to the outside
# of a ring.  Actually it's more, as much as sqrt(2)=1.4142 on the leading
# diagonal X=Y.  But the over-estimate is close enough for now.
#

# r = hypot(xmin,ymin)
# Nlo = (r-1/2)^2
#     = r^2 - r + 1/4
#     >= x^2+y^2 - (x+y)    because x+y >= r
#     = x(x-1) + y(y-1)
#     >= floorx(floorx-1) + floory(floory-1)
# in integers if round down to x=0 then x*(x-1)=0 too, so not negative
#
# r = hypot(xmax,ymax)
# Nhi = (r+1/2)^2
#     = r^2 + r + 1/4
#     <= x^2+y^2 + (x+y) + 1
#     = x(x+1) + y(y+1) + 1
#     <= ceilx(ceilx+1) + ceily(ceily+1) + 1

# Note: this code shared by TheodorusSpiral.  If start using the polar angle
# for more accuracy here then unshare it first.
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ($x1,$y1, $x2,$y2) = _rect_to_radius_corners ($x1,$y1, $x2,$y2);

  ### $x_min
  ### $y_min
  ### N min: $x_min*($x_min-1) + $y_min*($y_min-1)

  ### $x_max
  ### $y_max
  ### N max: $x_max*($x_max+1) + $y_max*($y_max+1) + 1

  return ($x1*($x1-1) + $y1*($y1-1),
          $x2*($x2+1) + $y2*($y2+1) + 1);
}

#------------------------------------------------------------------------------
# generic

# $x1,$y1, $x2,$y2 is a rectangle.
# Return ($xmin,$ymin, $xmax,$ymax).
#
# The two points are respectively minimum and maximum radius from the
# origin, rounded down or up to integers.
#
# If the rectangle is entirely one quadrant then the points are two opposing
# corners.  But if an axis is crossed then the minimum is on that axis and
# if the origin is covered then the minimum is 0,0.
#
# Currently the return is abs() absolute values of the places.  Could change
# that if there was any significance to the quadrant containing the min/max
# corners.
#
sub _rect_to_radius_corners {
  my ($x1,$y1, $x2,$y2) = @_;

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

  return (int($x2 < 0 ? -$x2
              : $x1 > 0 ? $x1
              : 0),
          int($y2 < 0 ? -$y2
              : $y1 > 0 ? $y1
              : 0),

          max(_ceil(abs($x1)), _ceil(abs($x2))),
          max(_ceil(abs($y1)), _ceil(abs($y2))));
}

sub _ceil {
  my ($x) = @_;
  my $int = int($x);
  return ($x > $int ? $int+1 : $int);
}

# FIXME: prefer to stay in integers if possible
# return ($rlo,$rhi) which is the radial distance range found in the rectangle
sub _rect_to_radius_range {
  my ($x1,$y1, $x2,$y2) = @_;

  ($x1,$y1, $x2,$y2) = _rect_to_radius_corners ($x1,$y1, $x2,$y2);
  return (hypot($x1,$y1),
          hypot($x2,$y2));
}

1;
__END__