This is a live mirror of the Perl 5 development currently hosted at https://github.com/perl/perl5
Math::Complex and Math::Trig updates (Re: [perl #37117] Math::Complex atan2 bug)
authorJarkko Hietaniemi <jhi@iki.fi>
Wed, 14 Sep 2005 09:26:11 +0000 (12:26 +0300)
committerRafael Garcia-Suarez <rgarciasuarez@gmail.com>
Wed, 14 Sep 2005 12:49:58 +0000 (12:49 +0000)
Message-ID: <4327C283.80706@gmail.com>

p4raw-id: //depot/perl@25414

lib/Math/Complex.pm
lib/Math/Complex.t
lib/Math/Trig.pm
lib/Math/Trig.t
pod/perlfunc.pod

index 400366c..e4b980b 100644 (file)
@@ -7,9 +7,9 @@
 
 package Math::Complex;
 
-our($VERSION, @ISA, @EXPORT, %EXPORT_TAGS, $Inf);
+use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf);
 
-$VERSION = 1.34;
+$VERSION = 1.35;
 
 BEGIN {
     unless ($^O eq 'unicosmk') {
@@ -38,7 +38,8 @@ my $i;
 my %LOGN;
 
 # Regular expression for floating point numbers.
-my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?)))';
+# These days we could use Scalar::Util::lln(), I guess.
+my $gre = qr'\s*([\+\-]?(?:(?:(?:\d+(?:_\d+)*(?:\.\d*(?:_\d+)*)?|\.\d+(?:_\d+)*)(?:[eE][\+\-]?\d+(?:_\d+)*)?))|inf)'i;
 
 require Exporter;
 
@@ -61,9 +62,12 @@ my @trig = qw(
             sqrt log ln
             log10 logn cbrt root
             cplx cplxe
+            atan2
             ),
           @trig);
 
+@EXPORT_OK = qw(decplx);
+
 %EXPORT_TAGS = (
     'trig' => [@trig],
 );
@@ -108,27 +112,53 @@ my $eps            = 1e-14;               # Epsilon
 # Die on bad *make() arguments.
 
 sub _cannot_make {
-    die "@{[(caller(1))[3]]}: Cannot take $_[0] of $_[1].\n";
+    die "@{[(caller(1))[3]]}: Cannot take $_[0] of '$_[1]'.\n";
 }
 
-sub _remake {
+sub _make {
     my $arg = shift;
-    my ($made, $p, $q);
+    my ($p, $q);
 
-    if ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
+    if ($arg =~ /^$gre$/) {
+       ($p, $q) = ($1, 0);
+    } elsif ($arg =~ /^(?:$gre)?$gre\s*i\s*$/) {
        ($p, $q) = ($1 || 0, $2);
-       $made = 'cart';
-    } elsif ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
+    } elsif ($arg =~ /^\s*\(\s*$gre\s*(?:,\s*$gre\s*)?\)\s*$/) {
        ($p, $q) = ($1, $2 || 0);
-       $made = 'exp';
     }
 
-    if ($made) {
+    if (defined $p) {
        $p =~ s/^\+//;
+       $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
        $q =~ s/^\+//;
+       $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
     }
 
-    return ($made, $p, $q);
+    return ($p, $q);
+}
+
+sub _emake {
+    my $arg = shift;
+    my ($p, $q);
+
+    if ($arg =~ /^\s*\[\s*$gre\s*(?:,\s*$gre\s*)?\]\s*$/) {
+       ($p, $q) = ($1, $2 || 0);
+    } elsif ($arg =~ m!^\s*\[\s*$gre\s*(?:,\s*([-+]?\d*\s*)?pi(?:/\s*(\d+))?\s*)?\]\s*$!) {
+       ($p, $q) = ($1, ($2 eq '-' ? -1 : ($2 || 1)) * pi() / ($3 || 1));
+    } elsif ($arg =~ /^\s*\[\s*$gre\s*\]\s*$/) {
+       ($p, $q) = ($1, 0);
+    } elsif ($arg =~ /^\s*$gre\s*$/) {
+       ($p, $q) = ($1, 0);
+    }
+
+    if (defined $p) {
+       $p =~ s/^\+//;
+       $q =~ s/^\+//;
+       $p =~ s/^(-?)inf$/"${1}9**9**9"/e;
+       $q =~ s/^(-?)inf$/"${1}9**9**9"/e;
+    }
+
+    return ($p, $q);
 }
 
 #
@@ -137,42 +167,26 @@ sub _remake {
 # Create a new complex number (cartesian form)
 #
 sub make {
-       my $self = bless {}, shift;
-       my ($re, $im) = @_;
-       if (@_ == 1) {
-           my ($remade, $p, $q) = _remake($re);
-           if ($remade) {
-               if ($remade eq 'cart') {
-                   ($re, $im) = ($p, $q);
-               } else {
-                   return (ref $self)->emake($p, $q);
-               }
-           }
-       }
-       my $rre = ref $re;
-       if ( $rre ) {
-           if ( $rre eq ref $self ) {
-               $re = Re($re);
-           } else {
-               _cannot_make("real part", $rre);
-           }
-       }
-       my $rim = ref $im;
-       if ( $rim ) {
-           if ( $rim eq ref $self ) {
-               $im = Im($im);
-           } else {
-               _cannot_make("imaginary part", $rim);
-           }
-       }
+    my $self = bless {}, shift;
+    my ($re, $im);
+    if (@_ == 0) {
+       ($re, $im) = (0, 0);
+    } elsif (@_ == 1) {
+       return (ref $self)->emake($_[0])
+           if ($_[0] =~ /^\s*\[/);
+       ($re, $im) = _make($_[0]);
+    } elsif (@_ == 2) {
+       ($re, $im) = @_;
+    }
+    if (defined $re) {
        _cannot_make("real part",      $re) unless $re =~ /^$gre$/;
-       $im ||= 0;
-       _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
-       $self->{'cartesian'} = [ $re, $im ];
-       $self->{c_dirty} = 0;
-       $self->{p_dirty} = 1;
-       $self->display_format('cartesian');
-       return $self;
+    }
+    $im ||= 0;
+    _cannot_make("imaginary part", $im) unless $im =~ /^$gre$/;
+    $self->set_cartesian([$re, $im ]);
+    $self->display_format('cartesian');
+
+    return $self;
 }
 
 #
@@ -181,46 +195,32 @@ sub make {
 # Create a new complex number (exponential form)
 #
 sub emake {
-       my $self = bless {}, shift;
-       my ($rho, $theta) = @_;
-       if (@_ == 1) {
-           my ($remade, $p, $q) = _remake($rho);
-           if ($remade) {
-               if ($remade eq 'exp') {
-                   ($rho, $theta) = ($p, $q);
-               } else {
-                   return (ref $self)->make($p, $q);
-               }
-           }
-       }
-       my $rrh = ref $rho;
-       if ( $rrh ) {
-           if ( $rrh eq ref $self ) {
-               $rho = rho($rho);
-           } else {
-               _cannot_make("rho", $rrh);
-           }
-       }
-       my $rth = ref $theta;
-       if ( $rth ) {
-           if ( $rth eq ref $self ) {
-               $theta = theta($theta);
-           } else {
-               _cannot_make("theta", $rth);
-           }
-       }
+    my $self = bless {}, shift;
+    my ($rho, $theta);
+    if (@_ == 0) {
+       ($rho, $theta) = (0, 0);
+    } elsif (@_ == 1) {
+       return (ref $self)->make($_[0])
+           if ($_[0] =~ /^\s*\(/ || $_[0] =~ /i\s*$/);
+       ($rho, $theta) = _emake($_[0]);
+    } elsif (@_ == 2) {
+       ($rho, $theta) = @_;
+    }
+    if (defined $rho && defined $theta) {
        if ($rho < 0) {
            $rho   = -$rho;
            $theta = ($theta <= 0) ? $theta + pi() : $theta - pi();
        }
+    }
+    if (defined $rho) {
        _cannot_make("rho",   $rho)   unless $rho   =~ /^$gre$/;
-       $theta ||= 0;
-       _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
-       $self->{'polar'} = [$rho, $theta];
-       $self->{p_dirty} = 0;
-       $self->{c_dirty} = 1;
-       $self->display_format('polar');
-       return $self;
+    }
+    $theta ||= 0;
+    _cannot_make("theta", $theta) unless $theta =~ /^$gre$/;
+    $self->set_polar([$rho, $theta]);
+    $self->display_format('polar');
+
+    return $self;
 }
 
 sub new { &make }              # For backward compatibility only.
@@ -312,8 +312,10 @@ sub cartesian {$_[0]->{c_dirty} ?
 sub polar     {$_[0]->{p_dirty} ?
                   $_[0]->update_polar : $_[0]->{'polar'}}
 
-sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{'cartesian'} = $_[1] }
-sub set_polar     { $_[0]->{c_dirty}++; $_[0]->{'polar'} = $_[1] }
+sub set_cartesian { $_[0]->{p_dirty}++; $_[0]->{c_dirty} = 0;
+                   $_[0]->{'cartesian'} = $_[1] }
+sub set_polar     { $_[0]->{c_dirty}++; $_[0]->{p_dirty} = 0;
+                   $_[0]->{'polar'} = $_[1] }
 
 #
 # ->update_cartesian
@@ -659,7 +661,7 @@ sub cbrt {
 # Die on bad root.
 #
 sub _rootbad {
-    my $mess = "Root $_[0] illegal, root rank must be positive integer.\n";
+    my $mess = "Root '$_[0]' illegal, root rank must be positive integer.\n";
 
     my @up = caller(1);
 
@@ -679,22 +681,27 @@ sub _rootbad {
 # z^(1/n) = r^(1/n) (cos ((t+2 k pi)/n) + i sin ((t+2 k pi)/n))
 #
 sub root {
-       my ($z, $n) = @_;
+       my ($z, $n, $k) = @_;
        _rootbad($n) if ($n < 1 or int($n) != $n);
        my ($r, $t) = ref $z ?
            @{$z->polar} : (CORE::abs($z), $z >= 0 ? 0 : pi);
-       my @root;
-       my $k;
        my $theta_inc = pit2 / $n;
        my $rho = $r ** (1/$n);
-       my $theta;
        my $cartesian = ref $z && $z->{c_dirty} == 0;
-       for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
-           my $w = cplxe($rho, $theta);
-           # Yes, $cartesian is loop invariant.
-           push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
+       if (@_ == 2) {
+           my @root;
+           for (my $i = 0, my $theta = $t / $n;
+                $i < $n;
+                $i++, $theta += $theta_inc) {
+               my $w = cplxe($rho, $theta);
+               # Yes, $cartesian is loop invariant.
+               push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
+           }
+           return @root;
+       } elsif (@_ == 3) {
+           my $w = cplxe($rho, $t / $n + $k * $theta_inc);
+           return $cartesian ? cplx(@{$w->cartesian}) : $w;
        }
-       return @root;
 }
 
 #
@@ -1265,27 +1272,30 @@ sub acotanh { Math::Complex::acoth(@_) }
 #
 # (atan2)
 #
-# Compute atan(z1/z2).
+# Compute atan(z1/z2), minding the right quadrant.
 #
 sub atan2 {
        my ($z1, $z2, $inverted) = @_;
        my ($re1, $im1, $re2, $im2);
        if ($inverted) {
            ($re1, $im1) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
-           ($re2, $im2) = @{$z1->cartesian};
+           ($re2, $im2) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
        } else {
-           ($re1, $im1) = @{$z1->cartesian};
+           ($re1, $im1) = ref $z1 ? @{$z1->cartesian} : ($z1, 0);
            ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
        }
-       if ($im2 == 0) {
-           return CORE::atan2($re1, $re2) if $im1 == 0;
-           return ($im1<=>0) * pip2 if $re2 == 0;
+       if ($im1 || $im2) {
+           # In MATLAB the imaginary parts are ignored.
+           # warn "atan2: Imaginary parts ignored";
+           # http://documents.wolfram.com/mathematica/functions/ArcTan
+           # NOTE: Mathematica ArcTan[x,y] while atan2(y,x)
+           my $s = $z1 * $z1 + $z2 * $z2;
+           _divbyzero("atan2") if $s == 0;
+           my $i = &i;
+           my $r = $z2 + $z1 * $i;
+           return -$i * &log($r / &sqrt( $s ));
        }
-       my $w = atan($z1/$z2);
-       my ($u, $v) = ref $w ? @{$w->cartesian} : ($w, 0);
-       $u += pi   if $re2 < 0;
-       $u -= pit2 if $u > pi;
-       return cplx($u, $v);
+       return CORE::atan2($re1, $re2);
 }
 
 #
@@ -1659,7 +1669,11 @@ the following (overloaded) operations are supported on complex numbers:
        log(z) = log(r1) + i*t
        sin(z) = 1/2i (exp(i * z1) - exp(-i * z))
        cos(z) = 1/2 (exp(i * z1) + exp(-i * z))
-       atan2(z1, z2) = atan(z1/z2)
+       atan2(y, x) = atan(y / x) # Minding the right quadrant, note the order.
+
+The definition used for complex arguments of atan2() is
+
+       -i log((x + iy)/sqrt(x*x+y*y))
 
 The following extra operations are supported on both real and complex
 numbers:
@@ -1726,6 +1740,9 @@ The I<k>th root for C<z = [r,t]> is given by:
 
        (root(z, n))[k] = r**(1/n) * exp(i * (t + 2*k*pi)/n)
 
+You can return the I<k>th root directly by C<root(z, n, k)>,
+indexing starting from I<zero> and ending at I<n - 1>.
+
 The I<spaceship> comparison operator, E<lt>=E<gt>, is also defined. In
 order to ensure its restriction to real numbers is conform to what you
 would expect, the comparison is run on the real part of the complex
@@ -1773,17 +1790,22 @@ understand a single (string) argument of the forms
        2-3i
        -3i
        [2,3]
+       [2,-3pi/4]
        [2]
 
 in which case the appropriate cartesian and exponential components
 will be parsed from the string and used to create new complex numbers.
 The imaginary component and the theta, respectively, will default to zero.
 
-=head1 STRINGIFICATION
+The C<new>, C<make>, C<emake>, C<cplx>, and C<cplxe> will also
+understand the case of no arguments: this means plain zero or (0, 0).
+
+=head1 DISPLAYING
 
 When printed, a complex number is usually shown under its cartesian
 style I<a+bi>, but there are legitimate cases where the polar style
-I<[r,t]> is more appropriate.
+I<[r,t]> is more appropriate.  The process of converting the complex
+number into a string that can be displayed is known as I<stringification>.
 
 By calling the class method C<Math::Complex::display_format> and
 supplying either C<"polar"> or C<"cartesian"> as an argument, you
@@ -1809,6 +1831,8 @@ The polar style attempts to emphasize arguments like I<k*pi/n>
 (where I<n> is a positive integer and I<k> an integer within [-9, +9]),
 this is called I<polar pretty-printing>.
 
+For the reverse of stringifying, see the C<make> and C<emake>.
+
 =head2 CHANGED IN PERL 5.6
 
 The C<display_format> class method and the corresponding
@@ -1902,7 +1926,8 @@ C<-1> (minus one).  For the C<atan>, C<acot>, the argument cannot be
 C<i> (the imaginary unit).  For the C<atan>, C<acoth>, the argument
 cannot be C<-i> (the negative imaginary unit).  For the C<tan>,
 C<sec>, C<tanh>, the argument cannot be I<pi/2 + k * pi>, where I<k>
-is any integer.
+is any integer.  atan2(0, 0) is undefined, and if the complex arguments
+are used for atan2(), a division by zero will happen if z1**2+z2**2 == 0.
 
 Note that because we are operating on approximations of real numbers,
 these errors can happen when merely `too close' to the singularities
@@ -1922,7 +1947,7 @@ messages like the following
 =head1 BUGS
 
 Saying C<use Math::Complex;> exports many mathematical routines in the
-caller environment and even overrides some (C<sqrt>, C<log>).
+caller environment and even overrides some (C<sqrt>, C<log>, C<atan2>).
 This is construed as a feature by the Authors, actually... ;-)
 
 All routines expect to be given real or complex numbers. Don't attempt to
index 555d5b5..b82350f 100755 (executable)
@@ -1,6 +1,5 @@
 #!./perl
 
-# $RCSfile: complex.t,v $
 #
 # Regression tests for the Math::Complex pacakge
 # -- Raphael Manfredi  since Sep 1996
@@ -8,8 +7,10 @@
 # -- Daniel S. Lewart  since Sep 1997
 
 BEGIN {
-    chdir 't' if -d 't';
-    @INC = '../lib';
+    if ($ENV{PERL_CORE}) {
+       chdir 't' if -d 't';
+       @INC = '../lib';
+    }
 }
 
 use Math::Complex;
@@ -23,7 +24,7 @@ my ($args, $op, $target, $test, $test_set, $try, $val, $zvalue, @set, @val);
 $test = 0;
 $| = 1;
 my @script = (
-    'my ($res, $s0,$s1,$s2,$s3,$s4,$s5,$s6,$s7,$s8,$s9,$s10, $z0,$z1,$z2);' .
+    'my ($res, $s0,$s1,$s2,$s3,$s4,$s5,$s6,$s7,$s8,$s9,$s10,$z0,$z1,$z2);' .
        "\n\n"
 );
 my $eps = 1e-13;
@@ -114,7 +115,9 @@ my $i    = cplx(0,  1);
 my $pi   = cplx(pi, 0);
 my $pii  = cplx(0, pi);
 my $pip2 = cplx(pi/2, 0);
+my $pip4 = cplx(pi/4, 0);
 my $zero = cplx(0, 0);
+my $inf  = 9**9**9;
 ';
 
 push(@script, $constants);
@@ -168,6 +171,7 @@ test_dbz(
         'coth(0)',
         'csc(0)',
         'csch(0)',
+        'atan(cplx(0, 1), cplx(1, 0))',
        );
 
 test_loz(
@@ -307,39 +311,152 @@ sub test_remake {
     $test++;
     push @script, <<EOS;
     print "# remake 2+3i\n";
-    my \$z = cplx('2+3i');
+    \$z = cplx('2+3i');
     print "not " unless \$z == Math::Complex->make(2,3);
     print "ok $test\n";
 EOS
 
     $test++;
     push @script, <<EOS;
-    print "# remake 3i\n";
-    my \$z = Math::Complex->make('3i');
+    print "# make 3i\n";
+    \$z = Math::Complex->make('3i');
     print "not " unless \$z == cplx(0,3);
     print "ok $test\n";
 EOS
 
     $test++;
     push @script, <<EOS;
-    print "# remake [2,3]\n";
-    my \$z = cplxe('[2,3]');
-    print "not " unless \$z == Math::Complex->emake(2,3);
+    print "# emake [2,3]\n";
+    \$z = Math::Complex->emake('[2,3]');
+    print "not " unless \$z == cplxe(2,3);
+    print "ok $test\n";
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    print "# make (2,3)\n";
+    \$z = Math::Complex->make('(2,3)');
+    print "not " unless \$z == cplx(2,3);
+    print "ok $test\n";
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    print "# emake [2,3pi/8]\n";
+    \$z = Math::Complex->emake('[2,3pi/8]');
+    print "not " unless \$z == cplxe(2,3*\$pi/8);
     print "ok $test\n";
 EOS
 
     $test++;
     push @script, <<EOS;
-    print "# remake [2]\n";
-    my \$z = Math::Complex->emake('[2]');
+    print "# emake [2]\n";
+    \$z = Math::Complex->emake('[2]');
     print "not " unless \$z == cplxe(2);
     print "ok $test\n";
 EOS
 }
 
+sub test_no_args {
+    push @script, <<'EOS';
+{
+    print "# cplx, cplxe, make, emake without arguments\n";
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    my \$z0 = cplx();
+    print ((\$z0->Re()  == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    print ((\$z0->Im()  == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    my \$z1 = cplxe();
+    print ((\$z1->rho()   == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    print ((\$z1->theta() == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    my \$z2 = Math::Complex->make();
+    print ((\$z2->Re()  == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    print ((\$z2->Im()  == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    my \$z3 = Math::Complex->emake();
+    print ((\$z3->rho()   == 0) ? "ok $test\n" : "not ok $test\n");
+EOS
+
+    $test++;
+    push @script, <<EOS;
+    print ((\$z3->theta() == 0) ? "ok $test\n" : "not ok $test\n");
+}
+EOS
+}
+
+sub test_atan2 {
+    push @script, <<'EOS';
+print "# atan2() with some real arguments\n";
+EOS
+    my @real = (-1, 0, 1);
+    for my $x (@real) {
+       for my $y (@real) {
+           next if $x == 0 && $y == 0;
+           $test++;
+           push @script, <<EOS;
+print ((Math::Complex::atan2($y, $x) == CORE::atan2($y, $x)) ? "ok $test\n" : "not ok $test\n");
+EOS
+        }
+    }
+    push @script, <<'EOS';
+    print "# atan2() with some complex arguments\n";
+EOS
+    $test++;
+    push @script, <<EOS;
+    print (abs(atan2(0, cplx(0, 1))) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+    $test++;
+    push @script, <<EOS;
+    print (abs(atan2(cplx(0, 1), 0) - \$pip2) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+    $test++;
+    push @script, <<EOS;
+    print (abs(atan2(cplx(0, 1), cplx(0, 1)) - \$pip4) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+    $test++;
+    push @script, <<EOS;
+    print (abs(atan2(cplx(0, 1), cplx(1, 1)) - cplx(0.553574358897045, 0.402359478108525)) < $eps ? "ok $test\n" : "not ok $test\n");
+EOS
+}
+
+sub test_decplx {
+}
+
 test_remake();
 
+test_no_args();
+
+test_atan2();
+
+test_decplx();
+
 print "1..$test\n";
+#print @script, "\n";
 eval join '', @script;
 die $@ if $@;
 
@@ -450,6 +567,8 @@ sub check {
            ||
            ($expected =~ /^-?\d/ && $got == $expected)
            ||
+           (abs(Math::Complex->make($got) - Math::Complex->make($expected)) < $eps)
+           ||
            (abs($got - $expected) < $eps)
            ) {
                print "ok $test\n";
@@ -591,6 +710,8 @@ __END__
 |'(root(z, 4))[1] ** 4':'z'
 |'(root(z, 5))[3] ** 5':'z'
 |'(root(z, 8))[7] ** 8':'z'
+|'(root(z, 8, 0)) ** 8':'z'
+|'(root(z, 8, 7)) ** 8':'z'
 |'abs(z)':'r'
 |'acot(z)':'acotan(z)'
 |'acsc(z)':'acosec(z)'
index 7560df5..9ff016f 100644 (file)
@@ -10,13 +10,14 @@ package Math::Trig;
 use 5.006;
 use strict;
 
+use Math::Complex 1.35;
 use Math::Complex qw(:trig);
 
 our($VERSION, $PACKAGE, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
 
 @ISA = qw(Exporter);
 
-$VERSION = 1.02;
+$VERSION = 1.03;
 
 my @angcnv = qw(rad2deg rad2grad
                deg2rad deg2grad
@@ -32,12 +33,30 @@ my @rdlcnv = qw(cartesian_to_cylindrical
                spherical_to_cartesian
                spherical_to_cylindrical);
 
-@EXPORT_OK = (@rdlcnv, 'great_circle_distance', 'great_circle_direction');
+my @greatcircle = qw(
+                    great_circle_distance
+                    great_circle_direction
+                    great_circle_bearing
+                    great_circle_waypoint
+                    great_circle_midpoint
+                    great_circle_destination
+                   );
 
-%EXPORT_TAGS = ('radial' => [ @rdlcnv ]);
+my @pi = qw(pi2 pip2 pip4);
+
+@EXPORT_OK = (@rdlcnv, @greatcircle, @pi);
+
+# See e.g. the following pages:
+# http://www.movable-type.co.uk/scripts/LatLong.html
+# http://williams.best.vwh.net/avform.htm
+
+%EXPORT_TAGS = ('radial' => [ @rdlcnv ],
+               'great_circle' => [ @greatcircle ],
+               'pi'     => [ @pi ]);
 
 sub pi2  () { 2 * pi }
 sub pip2 () { pi / 2 }
+sub pip4 () { pi / 4 }
 
 sub DR  () { pi2/360 }
 sub RD  () { 360/pi2 }
@@ -148,6 +167,57 @@ sub great_circle_direction {
     return rad2rad($direction);
 }
 
+*great_circle_bearing = \&great_circle_direction;
+
+sub great_circle_waypoint {
+    my ( $theta0, $phi0, $theta1, $phi1, $point ) = @_;
+
+    $point = 0.5 unless defined $point;
+
+    my $d = great_circle_distance( $theta0, $phi0, $theta1, $phi1 );
+
+    return undef if $d == pi;
+
+    my $sd = sin($d);
+
+    return ($theta0, $phi0) if $sd == 0;
+
+    my $A = sin((1 - $point) * $d) / $sd;
+    my $B = sin(     $point  * $d) / $sd;
+
+    my $lat0 = pip2 - $phi0;
+    my $lat1 = pip2 - $phi1;
+
+    my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1);
+    my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1);
+    my $z = $A * sin($lat0)                + $B * sin($lat1);
+
+    my $theta = atan2($y, $x);
+    my $phi   = atan2($z, sqrt($x*$x + $y*$y));
+    
+    return ($theta, $phi);
+}
+
+sub great_circle_midpoint {
+    great_circle_waypoint(@_[0..3], 0.5);
+}
+
+sub great_circle_destination {
+    my ( $theta0, $phi0, $dir0, $dst ) = @_;
+
+    my $lat0 = pip2 - $phi0;
+
+    my $phi1   = asin(sin($lat0)*cos($dst)+cos($lat0)*sin($dst)*cos($dir0));
+    my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
+                                cos($dst)-sin($lat0)*sin($phi1));
+
+    my $dir1 = great_circle_bearing($theta1, $phi1, $theta0, $phi0) + pi;
+
+    $dir1 -= pi2 if $dir1 > pi2;
+
+    return ($theta1, $phi1, $dir1);
+}
+
 1;
 
 __END__
@@ -169,12 +239,21 @@ Math::Trig - trigonometric functions
 
        $rad = deg2rad(120);
 
+        # Import constants pi2, pip2, pip4 (2*pi, pi/2, pi/4).
+       use Math::Trig ':pi';
+
+        # Import the conversions between cartesian/spherical/cylindrical.
+       use Math::Trig ':radial';
+
+        # Import the great circle formulas.
+       use Math::Trig ':great_circle';
+
 =head1 DESCRIPTION
 
 C<Math::Trig> defines many trigonometric functions not defined by the
 core Perl which defines only the C<sin()> and C<cos()>.  The constant
 B<pi> is also defined as are a few convenience functions for angle
-conversions.
+conversions, and I<great circle formulas> for spherical movement.
 
 =head1 TRIGONOMETRIC FUNCTIONS
 
@@ -265,7 +344,7 @@ C<asech>, C<acsch>, the argument cannot be C<0> (zero).  For the
 C<atanh>, C<acoth>, the argument cannot be C<1> (one).  For the
 C<atanh>, C<acoth>, the argument cannot be C<-1> (minus one).  For the
 C<tan>, C<sec>, C<tanh>, C<sech>, the argument cannot be I<pi/2 + k *
-pi>, where I<k> is any integer.
+pi>, where I<k> is any integer.  atan2(0, 0) is undefined.
 
 =head2 SIMPLE (REAL) ARGUMENTS, COMPLEX RESULTS
 
@@ -338,8 +417,7 @@ B<All angles are in radians>.
 
 =head2 COORDINATE SYSTEMS
 
-B<Cartesian> coordinates are the usual rectangular I<(x, y,
-z)>-coordinates.
+B<Cartesian> coordinates are the usual rectangular I<(x, y, z)>-coordinates.
 
 Spherical coordinates, I<(rho, theta, pi)>, are three-dimensional
 coordinates which define a point in three-dimensional space.  They are
@@ -430,16 +508,56 @@ degrees).
   $distance = great_circle_distance($lon0, pi/2 - $lat0,
                                     $lon1, pi/2 - $lat1, $rho);
 
-The direction you must follow the great circle can be computed by the
-great_circle_direction() function:
+The direction you must follow the great circle (also known as I<bearing>)
+can be computed by the great_circle_direction() function:
 
   use Math::Trig 'great_circle_direction';
 
   $direction = great_circle_direction($theta0, $phi0, $theta1, $phi1);
 
+(Alias 'great_circle_bearing' is also available.)
 The result is in radians, zero indicating straight north, pi or -pi
 straight south, pi/2 straight west, and -pi/2 straight east.
 
+You can inversely compute the destination if you know the
+starting point, direction, and distance:
+
+  use Math::Trig 'great_circle_destination';
+
+  # thetad and phid are the destination coordinates,
+  # dird is the final direction at the destination.
+
+  ($thetad, $phid, $dird) =
+    great_circle_destination($theta, $phi, $direction, $distance);
+
+or the midpoint if you know the end points:
+
+  use Math::Trig 'great_circle_midpoint';
+
+  ($thetam, $phim) =
+    great_circle_midpoint($theta0, $phi0, $theta1, $phi1);
+
+The great_circle_midpoint() is just a special case of
+
+  use Math::Trig 'great_circle_waypoint';
+
+  ($thetai, $phii) =
+    great_circle_waypoint($theta0, $phi0, $theta1, $phi1, $way);
+
+Where the $way is a value from zero ($theta0, $phi0) to one ($theta1,
+$phi1).  Note that antipodal points (where their distance is I<pi>
+radians) do not have waypoints between them (they would have an an
+"equator" between them), and therefore C<undef> is returned for
+antipodal points.  If the points are the same and the distance
+therefore zero and all waypoints therefore identical, the first point
+(either point) is returned.
+
+The thetas, phis, direction, and distance in the above are all in radians.
+
+You can import all the great circle formulas by
+
+  use Math::Trig ':great_circle';
+
 Notice that the resulting directions might be somewhat surprising if
 you are looking at a flat worldmap: in such map projections the great
 circles quite often do not look like the shortest routes-- but for
@@ -454,31 +572,31 @@ To calculate the distance between London (51.3N 0.5W) and Tokyo
         use Math::Trig qw(great_circle_distance deg2rad);
 
         # Notice the 90 - latitude: phi zero is at the North Pole.
-       @L = (deg2rad(-0.5), deg2rad(90 - 51.3));
-        @T = (deg2rad(139.8),deg2rad(90 - 35.7));
-
-        $km = great_circle_distance(@L, @T, 6378);
+       sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
+       my @L = NESW( -0.5, 51.3);
+        my @T = NESW(139.8, 35.7);
+        my $km = great_circle_distance(@L, @T, 6378); # About 9600 km.
 
-The direction you would have to go from London to Tokyo
+The direction you would have to go from London to Tokyo (in radians,
+straight north being zero, straight east being pi/2).
 
         use Math::Trig qw(great_circle_direction);
 
-        $rad = great_circle_direction(@L, @T);
+        my $rad = great_circle_direction(@L, @T); # About 0.547 or 0.174 pi.
 
-=head2 CAVEAT FOR GREAT CIRCLE FORMULAS
+The midpoint between London and Tokyo being
 
-The answers may be off by few percentages because of the irregular
-(slightly aspherical) form of the Earth.  The formula used for
-grear circle distances
+        use Math::Trig qw(great_circle_midpoint);
+
+        my @M = great_circle_midpoint(@L, @T);
+
+or about 68.11N 24.74E, in the Finnish Lapland.
 
-       lat0 = 90 degrees - phi0
-       lat1 = 90 degrees - phi1
-       d = R * arccos(cos(lat0) * cos(lat1) * cos(lon1 - lon01) +
-                       sin(lat0) * sin(lat1))
+=head2 CAVEAT FOR GREAT CIRCLE FORMULAS
 
-is also somewhat unreliable for small distances (for locations
-separated less than about five degrees) because it uses arc cosine
-which is rather ill-conditioned for values close to zero.
+The answers may be off by few percentages because of the irregular
+(slightly aspherical) form of the Earth.  The errors are at worst
+about 0.55%, but generally below 0.3%.
 
 =head1 BUGS
 
@@ -492,6 +610,8 @@ the computations even when the arguments are not. This, however,
 cannot be completely avoided if we want things like C<asin(2)> to give
 an answer instead of giving a fatal runtime error.
 
+Do not attempt navigation using these formulas.
+
 =head1 AUTHORS
 
 Jarkko Hietaniemi <F<jhi@iki.fi>> and 
index 5b1f12e..4777240 100755 (executable)
@@ -3,17 +3,21 @@
 #
 # Regression tests for the Math::Trig package
 #
-# The tests are quite modest as the Math::Complex tests exercise
-# these quite vigorously.
+# The tests here are quite modest as the Math::Complex tests exercise
+# these interfaces quite vigorously.
 # 
 # -- Jarkko Hietaniemi, April 1997
 
 BEGIN {
-    chdir 't' if -d 't';
-    @INC = '../lib';
+    if ($ENV{PERL_CORE}) {
+       chdir 't' if -d 't';
+       @INC = '../lib';
+    }
 }
 
-use Math::Trig;
+use Math::Trig 1.03;
+
+my $pip2 = pi / 2;
 
 use strict;
 
@@ -31,7 +35,7 @@ sub near ($$;$) {
     $_[1] ? (abs($_[0]/$_[1] - 1) < $e) : abs($_[0]) < $e;
 }
 
-print "1..29\n";
+print "1..49\n";
 
 $x = 0.9;
 print 'not ' unless (near(tan($x), sin($x) / cos($x)));
@@ -184,7 +188,7 @@ use Math::Trig ':radial';
        unless (near(great_circle_direction(0, 0, 0, pi/2), pi));
     print "ok 24\n";
 
-# Retired test: Relies on atan(0, 0), which is not portable.
+# Retired test: Relies on atan2(0, 0), which is not portable.
 #    print 'not '
 #      unless (near(great_circle_direction(0, 0, pi, pi), -pi()/2));
     print "ok 25\n";
@@ -197,8 +201,8 @@ use Math::Trig ':radial';
     print 'not '
        unless (near(rad2deg(great_circle_direction(@London, @Tokyo)),
                     31.791945393073));
-    
     print "ok 26\n";
+
     print 'not '
        unless (near(rad2deg(great_circle_direction(@Tokyo, @London)),
                     336.069766430326));
@@ -209,10 +213,107 @@ use Math::Trig ':radial';
                     246.800348034667));
     
     print "ok 28\n";
+
     print 'not '
        unless (near(rad2deg(great_circle_direction(@Paris, @Berlin)),
                     58.2079877553156));
     print "ok 29\n";
+
+    use Math::Trig 'great_circle_bearing';
+
+    print 'not '
+       unless (near(rad2deg(great_circle_bearing(@Paris, @Berlin)),
+                    58.2079877553156));
+    print "ok 30\n";
+
+    use Math::Trig 'great_circle_waypoint';
+    use Math::Trig 'great_circle_midpoint';
+
+    my ($lon, $lat);
+
+    ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.0);
+
+    print 'not ' unless (near($lon, $London[0]));
+    print "ok 31\n";
+
+    print 'not ' unless (near($lat, $pip2 - $London[1]));
+    print "ok 32\n";
+
+    ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 1.0);
+
+    print 'not ' unless (near($lon, $Tokyo[0]));
+    print "ok 33\n";
+
+    print 'not ' unless (near($lat, $pip2 - $Tokyo[1]));
+    print "ok 34\n";
+
+    ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.5);
+
+    print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E
+    print "ok 35\n";
+
+    print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N
+    print "ok 36\n";
+
+    ($lon, $lat) = great_circle_midpoint(@London, @Tokyo);
+
+    print 'not ' unless (near($lon, 1.55609593577679)); # 89.1577 E
+    print "ok 37\n";
+
+    print 'not ' unless (near($lat, 1.20296099733328)); # 68.9246 N
+    print "ok 38\n";
+
+    ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.25);
+
+    print 'not ' unless (near($lon, 0.516073562850837)); # 29.5688 E
+    print "ok 39\n";
+
+    print 'not ' unless (near($lat, 1.170565013391510)); # 67.0684 N
+    print "ok 40\n";
+    ($lon, $lat) = great_circle_waypoint(@London, @Tokyo, 0.75);
+
+    print 'not ' unless (near($lon, 2.17494903805952)); # 124.6154 E
+    print "ok 41\n";
+
+    print 'not ' unless (near($lat, 0.952987032741305)); # 54.6021 N
+    print "ok 42\n";
+
+    use Math::Trig 'great_circle_destination';
+
+    my $dir1 = great_circle_direction(@London, @Tokyo);
+    my $dst1 = great_circle_distance(@London,  @Tokyo);
+
+    ($lon, $lat) = great_circle_destination(@London, $dir1, $dst1);
+
+    print 'not ' unless (near($lon, $Tokyo[0]));
+    print "ok 43\n";
+
+    print 'not ' unless (near($lat, $pip2 - $Tokyo[1]));
+    print "ok 44\n";
+
+    my $dir2 = great_circle_direction(@Tokyo, @London);
+    my $dst2 = great_circle_distance(@Tokyo,  @London);
+
+    ($lon, $lat) = great_circle_destination(@Tokyo, $dir2, $dst2);
+
+    print 'not ' unless (near($lon, $London[0]));
+    print "ok 45\n";
+
+    print 'not ' unless (near($lat, $pip2 - $London[1]));
+    print "ok 46\n";
+
+    my $dir3 = (great_circle_destination(@London, $dir1, $dst1))[2];
+
+    print 'not ' unless (near($dir3, 2.69379263839118)); # about 154.343 deg
+    print "ok 47\n";
+
+    my $dir4 = (great_circle_destination(@Tokyo,  $dir2, $dst2))[2];
+
+    print 'not ' unless (near($dir4, 3.6993902625701)); # about 211.959 deg
+    print "ok 48\n";
+
+    print 'not ' unless (near($dst1, $dst2));
+    print "ok 49\n";
 }
 
 # eof
index f687728..2d6c40c 100644 (file)
@@ -446,6 +446,8 @@ function, or use the familiar relation:
 
     sub tan { sin($_[0]) / cos($_[0])  }
 
+Note that atan2(0, 0) is not well-defined.
+
 =item bind SOCKET,NAME
 
 Binds a network address to a socket, just as the bind system call