File Coverage

blib/lib/Math/PlanePath/ArchimedeanChords.pm

Criterion	Covered	Total	%
statement	130	154	84.4
branch	14	24	58.3
condition			n/a
subroutine	31	32	96.8
pod	4	4	100.0
total	179	214	83.6

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


# Also possible would be circle involute spiral, unrolling string around
# centre of circumference 1, but is only very slightly different radius from
# an Archimedean spiral.


package Math::PlanePath::ArchimedeanChords;
use 5.004;
use strict;
use Math::Libm 'hypot', 'asinh';
use POSIX 'ceil';
#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');

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

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


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


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

use constant 1.02 _PI => 2*atan2(1,0);

# Starting at polar angle position t in radians,
#
#     r = t / 2pi
#
#     x = r * cos(t) = t * cos(t) / 2pi
#     y = r * sin(t) = t * sin(t) / 2pi
#
# Want a polar angle amount u to move by a chord distance of 1.  Hypot
# square distance from t to t+u is
#
#     dist(u) =   ( (t+u)/2pi*cos(t+u) - t/2pi*cos(t) )^2     # X
#               + ( (t+u)/2pi*sin(t+u) - t/2pi*sin(t) )^2     # Y
#           = [  (t+u)^2*cos^2(t+u) - 2*(t+u)*t*cos(t+u)*cos(t) + t^2*cos^2(t)
#              + (t+u)^2*sin^2(t+u) - 2*(t+u)*t*sin(t+u)*sin(t) + t^2*sin^2(t)
#             ] / (4*pi^2)
#
# and from sin^2 + cos^2 = 1
# and addition cosA*cosB + sinA*sinB = cos(A-B)
#
#           = [  (t+u)^2            - 2*(t+u)*t*cos((t+u)-t)    + t^2 ] /4pi^2
#           = [ (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) ] / (4*pi^2)
#
# then double angle cos(u) = 1 - 2*sin^2(u/2) to go to the sine since if u
# is small then cos(u) near 1.0 might lose accuracy
#
#     dist(u) = [(t+u)^2 + t^2 - 2*t*(t+u)*(1 - 2*sin^2(u/2))] / (4*pi^2)
#             = [(t+u)^2 + t^2 - 2*t*(t+u) + 2*t*(t+u)*2*sin^2(u/2)] / (4*pi^2)
#             = [((t+u) - t)^2 + 4*t*(t+u)*sin^2(u/2)] / (4*pi^2)
#             = [ u^2 + 4*t*(t+u)*sin^2(u/2) ] / (4*pi^2)
#
# Seek d(u) = 1 by letting f(u)=4*pi^2*(d(u)-1) and seeking f(u)=0
#
#     f(u) = u^2 + 4*t*(t+u)*sin^2(u/2) - 4*pi^2
#
# Derivative f'(u) for the slope, starting from the cosine form,
#
#     f(u)  = (t+u)^2 + t^2 - 2*t*(t+u)*cos(u) - 4*pi^2
#
#     f'(u) = 2*(t+u) - 2*t*[ cos(u) - (t+u)*sin(u) ]
#           = 2*(t+u) - 2*t*[ 1 - 2*sin^2(u/2) - (t+u)*sin(u) ]
#           = 2*t + 2*u - 2*t + 2*t*2*sin^2(u/2) + 2*t*(t+u)*sin(u)
#           = 2*[ u + 2*t*sin^2(u/2) + t*(t+u)*sin(u) ]
#           = 2*[ u + t * [2*sin^2(u/2) + (t+u)*sin(u) ] ]
#
# Newton's method
#                           */    <- f(x) high
#                          */|
#                        * / |
#                      *  /  |
#          ---------*------------------
#                        +---+  <- subtract
#
#      f(x) / sub = f'(x)
#      sub = f(x) / f'(x)
#
#
# _chord_angle_inc() takes $t is a polar angle around the Archimedean spiral.
# Returns an increment polar angle $u which may be added to $t to move around
# the spiral by a chord length 1 unit.
#
# The loop is Newton's method, $f=f(u), $slope=f'(u) so $u-$f/$slope is a
# better $u, ie. f($u) closer to 0.  Stop when the subtract becomes small,
# usually only about 3 iterations.
#
sub _chord_angle_inc {
  my ($t) = @_;
  # ### _chord_angle_inc(): $t

  my $u = 2*_PI/$t; # estimate

  foreach (0 .. 10) {
    my $shu = sin($u/2); $shu *= $shu;   # sin^2(u/2)
    my $tu = ($t+$u);

    my $f = $u*$u + 4*$t*$tu*$shu - 4*_PI*_PI;
    my $slope = 2 * ( $u + $t*(2*$shu + $tu*sin($u)));

    # unsimplified versions ...
    # $f = ($t+$u)**2 + $t**2 - 2*$t*($t+$u)*cos($u) - 4*_PI*_PI;
    # $slope = 2*($t+$u) - 2*$t*( cos($u) - ($t+$u)*sin($u) );

    my $sub = $f/$slope;
    $u -= $sub;

    # printf ("f=%.6f slope=%.6f sub=%.20f u=%.6f\n", $f, $slope, $sub, $u);
    last if (abs($sub) < 1e-15);
  }
  # printf ("return u=%.6f\n", $u);
  return $u;
}

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

my @save_n = (1);
my @save_t = (2*_PI);
my $next_save = $save_n[0] + _SAVE;

sub new {
  ### ArchimedeanChords new() ...
  return shift->SUPER::new (i => $save_n[0],
                            t => $save_t[0],
                            @_);
}

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

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

  if ($n <= 1) {
    return ($n, 0);  # exactly Y=0
  }

  {
    # ENHANCE-ME: look at the N+1 position for the frac directly, without
    # the full call for N+1
    my $int = int($n);
    if ($n != $int) {
      my $frac = $n - $int;  # inherit possible BigFloat/BigRat
      my ($x1,$y1) = $self->n_to_xy($int);
      my ($x2,$y2) = $self->n_to_xy($int+1);
      my $dx = $x2-$x1;
      my $dy = $y2-$y1;
      return ($frac*$dx + $x1, $frac*$dy + $y1);
    }
  }

  my $i = $self->{'i'};
  my $t = $self->{'t'};

  if ($i > $n) {
    my $pos = min ($#save_n, int (($n - $save_n[0]) / _SAVE));
    $i = $save_n[$pos];
    $t = $save_t[$pos];
    ### resume: "$i  t=$t"
  }

  while ($i < $n) {
    $t += _chord_angle_inc($t);
    if (++$i == $next_save) {
      push @save_n, $i;
      push @save_t, $t;
      $next_save += _SAVE;
    }
  }

  $self->{'i'} = $i;
  $self->{'t'} = $t;

  my $r = $t * (1 / (2*_PI));
  return ($r*cos($t),
          $r*sin($t));
}

sub _xy_to_nearest_r {
  my ($x, $y) = @_;
  my $frac = Math::PlanePath::MultipleRings::_xy_to_angle_frac($x,$y);
  ### assert: 0 <= $frac && $frac < 1

  # if $frac > 0.5 then 0.5-$frac is negative and int() rounds towards zero
  # giving $r==$frac
  return int(hypot($x,$y) + 0.5 - $frac) + $frac;
}

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

  my $r = _xy_to_nearest_r($x,$y);
  my $r_limit = 1.001 * $r;

  ### hypot: hypot($x,$y)
  ### $r
  ### $r_limit
  ### save_t: "end index=$#save_t  save_t[0]=".($save_t[0]//'undef')

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

  my $theta = 0.999 * 2*_PI*$r;
  my $n_lo = 0;
  foreach my $i (1 .. $#save_t) {
    if ($save_t[$i] > $theta) {
      $n_lo = $save_n[$i-1];
      if ($n_lo == 1) { $n_lo = 0; } # for finding X=0,Y=0
      last;
    }
  }
  ### $n_lo

  # loop with for(;;) since $n_lo..$n_hi limited to IV range
  for (my $n = $n_lo; ; $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) {
      ### hypot in range ...
      return $n;
    }
    if (hypot($nx,$ny) >= $r_limit) {
      last;
    }
  }

  ### n not found ...
  return undef;
}

  # int (max (0, int(_PI*$r2) - 4*$r));
  #
  #  my $r2 = $r * $r;
  #  my $n_lo =  int (max (0, int(_PI*$r2) - 4*$r));
  #  my $n_hi = $n_lo + 7*$r + 2;
  # ### $r2
  # $n_lo == $n_lo-1 ||




# x,y has radius hypot(x,y), then want the next higher spiral arc which is r
# >= hypot(x,y)+0.5, with the 0.5 being the width of the circle figure on
# the spiral.
#
# The polar angle of x,y is a=atan2(y,x) and frac=a/2pi is the extra away
# from an integer radius for the spiral.  So seek integer k with k+a/2pi >=
# h with h=hypot(x,y)+0.5.
#
#     k + a/2pi >= h
#     k >= h-a/2pi
#     k = ceil(h-a/2pi)
#       = ceil(hypot(x,y) + 0.5 - atan2(y,x)/2pi)
#
#
# circle radius i has circumference 2*pi*i and at most that many N on it
# rectangle corner at radius Rcorner = hypot(x,y)
#
# sum i=1 to i=Rlimit of 2*pi*i = 2*pi/2 * Rlimit*(Rlimit+1)
#                               = pi * Rlimit*(Rlimit+1)
# is an upper bound, though a fairly slack one
#
#
# cf. arc length along the spiral r=a*theta with a=1/2pi
#     arclength = (1/2) * a * (theta*sqrt(1+theta^2) + asinh(theta))
#               = (1/4*pi) * (theta*sqrt(1+theta^2) + asinh(theta))
# and theta = 2*pi*r
#               = (1/4*pi) * (4*pi^2*r^2 * sqrt(1+1/theta^2) + asinh(theta))
#               = pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
#
# and to compare to the circles formula
#
#               = pi * (r*(r+1) * r/(r+1) * sqrt(1+1/r^2)
#                       + asinh(theta)/(4*pi^2))
#
# so it's smaller hence better upper bound.  Only a little smaller than the
# squaring once get to 100 loops or so.
#
#
# not exact
sub rect_to_n_range {
  my ($self, $x1,$y1, $x2,$y2) = @_;
  ### rect_to_n_range() ...

  my $rhi = 0;
  foreach my $x ($x1, $x2) {
    foreach my $y ($y1, $y2) {
      my $frac = atan2($y,$x) / (2*_PI);  # -0.5 <= $frac <= 0.5
      $frac += ($frac < 0);                #    0 <= $frac < 1
      $rhi = max ($rhi, ceil(hypot($x,$y)+0.5 - $frac) + $frac);
    }
  }
  ### $rhi

  # arc length pi * (r^2 * sqrt(1+1/r^2) + asinh(theta)/(4*pi^2))
  #          = pi*r^2*sqrt(1+1/r^2) + asinh(theta)/4pi
  my $rhi2 = $rhi*$rhi;
  return (0,
          ceil (_PI * $rhi2 * sqrt(1+1/$rhi2)
                + asinh(2*_PI*$rhi) / (4*_PI)));


  # # each loop begins at N = pi*k^2 - 2 or thereabouts
  # return (0,
  #         int(_PI*$rhi*($rhi+1) + 1));
}

1;
__END__