Upgrade to Math-Complex-1.42
authorSteve Peters <steve@fisharerojo.org>
Wed, 9 Jan 2008 02:36:33 +0000 (02:36 +0000)
committerSteve Peters <steve@fisharerojo.org>
Wed, 9 Jan 2008 02:36:33 +0000 (02:36 +0000)
p4raw-id: //depot/perl@32908

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

index 4006b6f..315bea4 100644 (file)
@@ -9,26 +9,37 @@ package Math::Complex;
 
 use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $Inf);
 
-$VERSION = 1.38;
+$VERSION = 1.42;
 
 BEGIN {
-    unless ($^O eq 'unicosmk') {
+    my $IEEE_DBL_MAX = eval "1.7976931348623157e+308";
+    if ($^O eq 'unicosmk') {
+       $Inf = $IEEE_DBL_MAX;
+    } else {
         local $!;
        # We do want an arithmetic overflow, Inf INF inf Infinity:.
-        undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
-         local $SIG{FPE} = sub {die};
-         my $t = CORE::exp 30;
-         $Inf = CORE::exp $t;
-EOE
-       if (!defined $Inf) {            # Try a different method
-         undef $Inf unless eval <<'EOE' and $Inf =~ /^inf(?:inity)?$/i;
-           local $SIG{FPE} = sub {die};
-           my $t = 1;
-           $Inf = $t + "1e99999999999999999999999999999999";
-EOE
+       for my $t (
+           'exp(999)',
+           '9**9**9',
+           '1e999',
+           'inf',
+           'Inf',
+           'INF',
+           'infinity',
+           'Infinity',
+           'INFINITY',
+           ) {
+           local $SIG{FPE} = { };
+           local $^W = 0;
+           my $i = eval "$t+1.0";
+           if ($i =~ /inf/i && $i > 1e+99) {
+               $Inf = $i;
+               last;
+           }
        }
+       $Inf = $IEEE_DBL_MAX unless defined $Inf;  # Oh well, close enough.
+       die "Could not get Infinity" unless $Inf > 1e99;
     }
-    $Inf = "Inf" if !defined $Inf || !($Inf > 0); # Desperation.
 }
 
 use strict;
@@ -65,7 +76,7 @@ my @trig = qw(
             ),
           @trig);
 
-my @pi = qw(pi pi2 pi4 pip2 pip4);
+my @pi = qw(pi pi2 pi4 pip2 pip4 Inf);
 
 @EXPORT_OK = @pi;
 
@@ -109,8 +120,6 @@ my $eps            = 1e-14;         # Epsilon
 #      c_dirty         cartesian form not up-to-date
 #      p_dirty         polar form not up-to-date
 #      display         display format (package's global when not set)
-#      bn_cartesian
-#       bnc_dirty
 #
 
 # Die on bad *make() arguments.
@@ -858,7 +867,7 @@ sub cos {
        my $ey = CORE::exp($y);
        my $sx = CORE::sin($x);
        my $cx = CORE::cos($x);
-       my $ey_1 = $ey ? 1 / $ey : $Inf;
+       my $ey_1 = $ey ? 1 / $ey : Inf();
        return (ref $z)->make($cx * ($ey + $ey_1)/2,
                              $sx * ($ey_1 - $ey)/2);
 }
@@ -875,7 +884,7 @@ sub sin {
        my $ey = CORE::exp($y);
        my $sx = CORE::sin($x);
        my $cx = CORE::cos($x);
-       my $ey_1 = $ey ? 1 / $ey : $Inf;
+       my $ey_1 = $ey ? 1 / $ey : Inf();
        return (ref $z)->make($sx * ($ey + $ey_1)/2,
                              $cx * ($ey - $ey_1)/2);
 }
@@ -1069,11 +1078,11 @@ sub cosh {
        my $ex;
        unless (ref $z) {
            $ex = CORE::exp($z);
-           return $ex ? ($ex + 1/$ex)/2 : $Inf;
+           return $ex ? ($ex + 1/$ex)/2 : Inf();
        }
        my ($x, $y) = @{$z->_cartesian};
        $ex = CORE::exp($x);
-       my $ex_1 = $ex ? 1 / $ex : $Inf;
+       my $ex_1 = $ex ? 1 / $ex : Inf();
        return (ref $z)->make(CORE::cos($y) * ($ex + $ex_1)/2,
                              CORE::sin($y) * ($ex - $ex_1)/2);
 }
@@ -1089,13 +1098,13 @@ sub sinh {
        unless (ref $z) {
            return 0 if $z == 0;
            $ex = CORE::exp($z);
-           return $ex ? ($ex - 1/$ex)/2 : "-$Inf";
+           return $ex ? ($ex - 1/$ex)/2 : -Inf();
        }
        my ($x, $y) = @{$z->_cartesian};
        my $cy = CORE::cos($y);
        my $sy = CORE::sin($y);
        $ex = CORE::exp($x);
-       my $ex_1 = $ex ? 1 / $ex : $Inf;
+       my $ex_1 = $ex ? 1 / $ex : Inf();
        return (ref $z)->make(CORE::cos($y) * ($ex - $ex_1)/2,
                              CORE::sin($y) * ($ex + $ex_1)/2);
 }
@@ -1109,7 +1118,10 @@ sub tanh {
        my ($z) = @_;
        my $cz = cosh($z);
        _divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
-       return sinh($z) / $cz;
+       my $sz = sinh($z);
+       return  1 if $cz ==  $sz;
+       return -1 if $cz == -$sz;
+       return $sz / $cz;
 }
 
 #
@@ -1152,7 +1164,10 @@ sub coth {
        my ($z) = @_;
        my $sz = sinh($z);
        _divbyzero "coth($z)", "sinh($z)" if $sz == 0;
-       return cosh($z) / $sz;
+       my $cz = cosh($z);
+       return  1 if $cz ==  $sz;
+       return -1 if $cz == -$sz;
+       return $cz / $sz;
 }
 
 #
@@ -1488,6 +1503,10 @@ sub _stringify_polar {
        return "[$r,$theta]";
 }
 
+sub Inf {
+    return $Inf;
+}
+
 1;
 __END__
 
@@ -1755,11 +1774,11 @@ The I<k>th root for C<z = [r,t]> is given by:
 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
-number first, and imaginary parts are compared only when the real
-parts match.
+The I<spaceship> numeric 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 number first, and imaginary parts are compared only when
+the real parts match.
 
 =head1 CREATION
 
@@ -1916,6 +1935,25 @@ exported:
     use Math::Complex ':pi'; 
     $third_of_circle = pi2 / 3;
 
+=head2 Inf
+
+The floating point infinity can be exported as a subroutine Inf():
+
+    use Math::Complex qw(Inf sinh);
+    my $AlsoInf = Inf() + 42;
+    my $AnotherInf = sinh(1e42);
+    print "$AlsoInf is $AnotherInf\n" if $AlsoInf == $AnotherInf;
+
+Note that the stringified form of infinity varies between platforms:
+it can be for example any of
+
+   inf
+   infinity
+   INF
+   1.#INF
+
+or it can be something else. 
+
 =head1 ERRORS DUE TO DIVISION BY ZERO OR LOGARITHM OF ZERO
 
 The division (/) and the following functions
@@ -1986,6 +2024,11 @@ Daniel S. Lewart <F<lewart!at!uiuc.edu>>
 Jarkko Hietaniemi <F<jhi!at!iki.fi>>
 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>
 
+=head1 LICENSE
+
+This library is free software; you can redistribute it and/or modify
+it under the same terms as Perl itself. 
+
 =cut
 
 1;
index 4c5ab37..25c068d 100755 (executable)
@@ -13,7 +13,7 @@ BEGIN {
     }
 }
 
-use Math::Complex 1.38;
+use Math::Complex 1.42;
 
 use vars qw($VERSION);
 
index b7bfbb9..5bd85c5 100644 (file)
@@ -10,14 +10,14 @@ package Math::Trig;
 use 5.005;
 use strict;
 
-use Math::Complex 1.37;
+use Math::Complex 1.42;
 use Math::Complex qw(:trig :pi);
 
 use vars qw($VERSION $PACKAGE @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS);
 
 @ISA = qw(Exporter);
 
-$VERSION = 1.05;
+$VERSION = 1.07;
 
 my @angcnv = qw(rad2deg rad2grad
                deg2rad deg2grad
@@ -44,7 +44,7 @@ my @greatcircle = qw(
 
 my @pi = qw(pi pi2 pi4 pip2 pip4);
 
-@EXPORT_OK = (@rdlcnv, @greatcircle, @pi);
+@EXPORT_OK = (@rdlcnv, @greatcircle, @pi, 'Inf');
 
 # See e.g. the following pages:
 # http://www.movable-type.co.uk/scripts/LatLong.html
@@ -674,6 +674,11 @@ Do not attempt navigation using these formulas.
 Jarkko Hietaniemi <F<jhi!at!iki.fi>> and 
 Raphael Manfredi <F<Raphael_Manfredi!at!pobox.com>>.
 
+=head1 LICENSE
+
+This library is free software; you can redistribute it and/or modify
+it under the same terms as Perl itself. 
+
 =cut
 
 # eof
index eab2ed7..93012b4 100755 (executable)
@@ -26,9 +26,10 @@ BEGIN {
     }
 }
 
-plan(tests => 69);
+plan(tests => 135);
 
-use Math::Trig 1.05;
+use Math::Trig 1.07;
+use Math::Trig 1.07 qw(Inf);
 
 my $pip2 = pi / 2;
 
@@ -49,6 +50,42 @@ sub near ($$;$) {
     $_[1] ? ($d < $e) : abs($_[0]) < $e;
 }
 
+print "# Sanity checks\n";
+
+ok(near(sin(1), 0.841470984807897));
+ok(near(cos(1), 0.54030230586814));
+ok(near(tan(1), 1.5574077246549));
+
+ok(near(sec(1), 1.85081571768093));
+ok(near(csc(1), 1.18839510577812));
+ok(near(cot(1), 0.642092615934331));
+
+ok(near(asin(1), 1.5707963267949));
+ok(near(acos(1), 0));
+ok(near(atan(1), 0.785398163397448));
+
+ok(near(asec(1), 0));
+ok(near(acsc(1), 1.5707963267949));
+ok(near(acot(1), 0.785398163397448));
+
+ok(near(sinh(1), 1.1752011936438));
+ok(near(cosh(1), 1.54308063481524));
+ok(near(tanh(1), 0.761594155955765));
+
+ok(near(sech(1), 0.648054273663885));
+ok(near(csch(1), 0.850918128239322));
+ok(near(coth(1), 1.31303528549933));
+
+ok(near(asinh(1), 0.881373587019543));
+ok(near(acosh(1), 0));
+ok(near(atanh(0.9), 1.47221948958322)); # atanh(1.0) would be an error.
+
+ok(near(asech(0.9), 0.467145308103262));
+ok(near(acsch(2), 0.481211825059603));
+ok(near(acoth(2), 0.549306144334055));
+
+print "# Basics\n";
+
 $x = 0.9;
 ok(near(tan($x), sin($x) / cos($x)));
 
@@ -270,4 +307,58 @@ use Math::Trig ':radial';
     ok(near($dst1, $dst2));
 }
 
+print "# Infinity\n";
+
+my $BigDouble = 1e40;
+
+ok(Inf() > $BigDouble);
+ok(Inf() + $BigDouble > $BigDouble);
+ok(Inf() + $BigDouble == Inf());
+ok(Inf() - $BigDouble > $BigDouble);
+ok(Inf() - $BigDouble == Inf());
+ok(Inf() * $BigDouble > $BigDouble);
+ok(Inf() * $BigDouble == Inf());
+ok(Inf() / $BigDouble > $BigDouble);
+ok(Inf() / $BigDouble == Inf());
+
+ok(-Inf() < -$BigDouble);
+ok(-Inf() + $BigDouble < $BigDouble);
+ok(-Inf() + $BigDouble == -Inf());
+ok(-Inf() - $BigDouble < -$BigDouble);
+ok(-Inf() - $BigDouble == -Inf());
+ok(-Inf() * $BigDouble < -$BigDouble);
+ok(-Inf() * $BigDouble == -Inf());
+ok(-Inf() / $BigDouble < -$BigDouble);
+ok(-Inf() / $BigDouble == -Inf());
+
+print "# sinh/sech/cosh/csch/tanh/coth unto infinity\n";
+
+ok(near(sinh(100), 1.3441e+43, 1e-3));
+ok(near(sech(100), 7.4402e-44, 1e-3));
+ok(near(cosh(100), 1.3441e+43, 1e-3));
+ok(near(csch(100), 7.4402e-44, 1e-3));
+ok(near(tanh(100), 1));
+ok(near(coth(100), 1));
+
+ok(near(sinh(-100), -1.3441e+43, 1e-3));
+ok(near(sech(-100),  7.4402e-44, 1e-3));
+ok(near(cosh(-100),  1.3441e+43, 1e-3));
+ok(near(csch(-100), -7.4402e-44, 1e-3));
+ok(near(tanh(-100), -1));
+ok(near(coth(-100), -1));
+
+ok(sinh(1e4) == Inf());
+ok(sech(1e4) == 0);
+ok(cosh(1e4) == Inf());
+ok(csch(1e4) == 0);
+ok(tanh(1e4) == 1);
+ok(coth(1e4) == 1);
+
+ok(sinh(-1e4) == -Inf());
+ok(sech(-1e4) == 0);
+ok(cosh(-1e4) == Inf());
+ok(csch(-1e4) == 0);
+ok(tanh(-1e4) == -1);
+ok(coth(-1e4) == -1);
+
 # eof