File Coverage

blib/lib/Math/PlanePath/TheodorusSpiral.pm

Criterion	Covered	Total	%
statement	111	125	88.8
branch	13	20	65.0
condition	0	3	0.0
subroutine	27	28	96.4
pod	4	4	100.0
total	155	180	86.1

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


# Hlawka, angles of point N is
# phi(n) = sum k=1 to n of arcsin 1/sqrt(k+1)
# is equidistributed mod 2pi


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

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

use Math::PlanePath::Base::Generic
  'is_infinite',
  'round_nearest';

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


use constant n_start => 0;
use constant figure => 'circle';
use constant x_negative_at_n => 4;
use constant y_negative_at_n => 7;
use constant gcdxy_maximum => 1;
use constant dx_minimum => -1; # supremum when straight
use constant dx_maximum => 1;  # at N=0
use constant dy_minimum => -1;
use constant dy_maximum => 1;  # at N=1
use constant dsumxy_minimum => -sqrt(2); # supremum diagonal
use constant dsumxy_maximum => sqrt(2);
use constant ddiffxy_minimum => -sqrt(2); # supremum diagonal
use constant ddiffxy_maximum => sqrt(2);
use constant turn_any_right    => 0; # left always
use constant turn_any_straight => 0; # left always


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

# This adding up of unit steps isn't very good.  The last x,y,n is kept
# anticipating successively higher n, not necessarily consecutive, plus past
# x,y,n at _SAVE intervals for going backwards.
#
# The simplest formulas for the polar angle, possibly with the analytic
# continuation version don't seem much better, but theta approaches
# 2*sqrt(N) + const, or 2*sqrt(N) + 1/(6*sqrt(N+1)) + const + O(n^(3/2)), so
# more terms of that might have tolerably rapid convergence.
#
# The arctan sums for the polar angle end up as the generalized Riemann
# zeta, or the generalized minus the plain.  Is there a good formula for
# that which would converge quickly?

use constant 1.02; # for leading underscore
use constant _SAVE => 1000;

my @save_n = (1);
my @save_x = (1);
my @save_y = (0);
my $next_save = _SAVE;

sub new {
  return shift->SUPER::new (i => 1,
                            x => 1,
                            y => 0,
                            @_);
}

# r = sqrt(int)
# (frac r)^2
#   = hypot(r, frac)^2   frac at right angle to radial
#   = r^2 + $frac^2
#   = sqrt(int)^2 + $frac^2
#   = $int + $frac^2
#
sub n_to_rsquared {
  my ($self, $n) = @_;
  if ($n < 0) { return undef; }
  my $int = int($n);
  $n -= $int;  # fractional part
  return $n*$n + $int;
}

# r = sqrt(i)
# x,y angle
# r*x/hypot, r*y/hypot
#
# newx = x - y/r
# newy = y + x/r
# (x-y/r)^2 + (y+x/r)^2
#   =   x^2 - 2y/r + y^2/r^2
#     + y^2 + 2x/r + x^2/r^2

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

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

  if ($n < 1) {
    return ($n, 0);
  }
  my $frac = $n;
  $n = int($n);
  $frac -= $n;

  my $i = $self->{'i'};
  my $x = $self->{'x'};
  my $y = $self->{'y'};
  #### n_to_xy(): "$n from state $i $x,$y"

  if ($i > $n) {
    for (my $pos = $#save_n; $pos >= 0; $pos--) {
      if ($save_n[$pos] <= $n) {
        $i = $save_n[$pos];
        $x = $save_x[$pos];
        $y = $save_y[$pos];
        last;
      }
    }
    ### resume: "$i  $x,$y"
  }

  while ($i < $n) {
    my $r = sqrt($i);
    ($x,$y) = ($x - $y/$r, $y + $x/$r);
    $i++;

    if ($i == $next_save) {
      push @save_n, $i;
      push @save_x, $x;
      push @save_y, $y;
      $next_save += _SAVE;

      ### save: $i
      ### @save_n
      ### @save_x
      ### @save_y
    }
  }

  $self->{'i'} = $i;
  $self->{'x'} = $x;
  $self->{'y'} = $y;

  if ($frac) {
    my $r = sqrt($n);
    return ($x - $frac*$y/$r,
            $y + $frac*$x/$r);
  } else {
    #### integer return: "$i  $x,$y"
    return ($x,$y);
  }
}

sub xy_to_n {
  my ($self, $x, $y) = @_;
  ### TheodorusSpiral xy_to_n(): "$x, $y"
  my $r = hypot ($x,$y);
  my $n_lo = int (max (0, $r - .51) ** 2);
  my $n_hi = int (($r + .51) ** 2);
  ### $n_lo
  ### $n_hi

  if (is_infinite($n_lo) || is_infinite($n_hi)) {
    ### infinite range, r inf or too big ...
    return undef;
  }

  # for(;;) loop since $n_lo..$n_hi limited to IV range
  for (my $n = $n_lo; $n <= $n_hi; $n += 1) {
    my ($nx,$ny) = $self->n_to_xy($n);
    #### $n
    #### $nx
    #### $ny
    #### hypot: hypot ($x-$nx,$y-$ny)
    if (hypot ($x-$nx,$y-$ny) <= 0.5) {
      return $n;
    }
  }
  return undef;
}

use Math::PlanePath::SacksSpiral;
# not exact
*rect_to_n_range = \&Math::PlanePath::SacksSpiral::rect_to_n_range;

1;
__END__