view release on metacpan or  search on metacpan

lib/Math/PlanePath/ArchimedeanChords.pm  view on Meta::CPAN

#     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);