Math::Complex update
Jarkko Hietaniemi [Thu, 2 Jan 1997 23:17:59 +0000 (11:17 +1200)]
lib/Math/Complex.pm
t/lib/complex.t

index fce53f7..30194eb 100644 (file)
@@ -1,34 +1,48 @@
 # $RCSFile$
 #
 # Complex numbers and associated mathematical functions
-# -- Raphael Manfredi, Sept 1996
+# -- Raphael Manfredi, September 1996
+# -- Jarkko Hietaniemi, March 1997
 
 require Exporter;
 package Math::Complex; @ISA = qw(Exporter);
 
+use strict;
+
+use vars qw(@EXPORT $package $display
+           $pi $i $ilog10 $logn %logn);
+
 @EXPORT = qw(
        pi i Re Im arg
+       sqrt exp log ln
        log10 logn cbrt root
-       tan cotan asin acos atan acotan
-       sinh cosh tanh cotanh asinh acosh atanh acotanh
+       tan
+       cosec csc sec cotan cot
+       asin acos atan
+       acosec acsc asec acotan acot
+       sinh cosh tanh
+       cosech csch sech cotanh coth
+       asinh acosh atanh
+       acosech acsch asech acotanh acoth
        cplx cplxe
 );
 
 use overload
-       '+'             => \&plus,
-       '-'             => \&minus,
-       '*'             => \&multiply,
-       '/'             => \&divide,
+       '+'     => \&plus,
+       '-'     => \&minus,
+       '*'     => \&multiply,
+       '/'     => \&divide,
        '**'    => \&power,
        '<=>'   => \&spaceship,
        'neg'   => \&negate,
-       '~'             => \&conjugate,
+       '~'     => \&conjugate,
        'abs'   => \&abs,
        'sqrt'  => \&sqrt,
        'exp'   => \&exp,
        'log'   => \&log,
        'sin'   => \&sin,
        'cos'   => \&cos,
+       'tan'   => \&tan,
        'atan2' => \&atan2,
        qw("" stringify);
 
@@ -87,7 +101,7 @@ sub new { &make }            # For backward compatibility only.
 #
 sub cplx {
        my ($re, $im) = @_;
-       return $package->make($re, $im);
+       return $package->make($re, defined $im ? $im : 0);
 }
 
 #
@@ -98,7 +112,7 @@ sub cplx {
 #
 sub cplxe {
        my ($rho, $theta) = @_;
-       return $package->emake($rho, $theta);
+       return $package->emake($rho, defined $theta ? $theta : 0);
 }
 
 #
@@ -129,8 +143,10 @@ sub i () {
 # Attribute access/set routines
 #
 
-sub cartesian {$_[0]->{c_dirty} ? $_[0]->update_cartesian : $_[0]->{'cartesian'}}
-sub polar     {$_[0]->{p_dirty} ? $_[0]->update_polar : $_[0]->{'polar'}}
+sub cartesian {$_[0]->{c_dirty} ?
+                  $_[0]->update_cartesian : $_[0]->{'cartesian'}}
+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] }
@@ -168,8 +184,9 @@ sub update_polar {
 #
 sub plus {
        my ($z1, $z2, $regular) = @_;
+       $z2 = cplx($z2, 0) unless ref $z2;
        my ($re1, $im1) = @{$z1->cartesian};
-       my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2);
+       my ($re2, $im2) = @{$z2->cartesian};
        unless (defined $regular) {
                $z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
                return $z1;
@@ -184,8 +201,9 @@ sub plus {
 #
 sub minus {
        my ($z1, $z2, $inverted) = @_;
+       $z2 = cplx($z2, 0) unless ref $z2;
        my ($re1, $im1) = @{$z1->cartesian};
-       my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2);
+       my ($re2, $im2) = @{$z2->cartesian};
        unless (defined $inverted) {
                $z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
                return $z1;
@@ -203,7 +221,8 @@ sub minus {
 sub multiply {
        my ($z1, $z2, $regular) = @_;
        my ($r1, $t1) = @{$z1->polar};
-       my ($r2, $t2) = ref $z2 ? @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
+       my ($r2, $t2) = ref $z2 ?
+           @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
        unless (defined $regular) {
                $z1->set_polar([$r1 * $r2, $t1 + $t2]);
                return $z1;
@@ -212,6 +231,20 @@ sub multiply {
 }
 
 #
+# divbyzero
+#
+# Die on division by zero.
+#
+sub divbyzero {
+    warn $package . '::' . "$_[0]: Division by zero.\n";
+    warn "(Because in the definition of $_[0], $_[1] is 0)\n"
+       if (defined $_[1]);
+    my @up = caller(1);
+    my $dmess = "Died at $up[1] line $up[2].\n";
+    die $dmess;
+}
+
+#
 # (divide)
 #
 # Computes z1/z2.
@@ -219,14 +252,20 @@ sub multiply {
 sub divide {
        my ($z1, $z2, $inverted) = @_;
        my ($r1, $t1) = @{$z1->polar};
-       my ($r2, $t2) = ref $z2 ? @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
+       my ($r2, $t2) = ref $z2 ?
+           @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
        unless (defined $inverted) {
+               divbyzero "$z1/0" if ($r2 == 0);
                $z1->set_polar([$r1 / $r2, $t1 - $t2]);
                return $z1;
        }
-       return $inverted ?
-               (ref $z1)->emake($r2 / $r1, $t2 - $t1) :
-               (ref $z1)->emake($r1 / $r2, $t1 - $t2);
+       if ($inverted) {
+               divbyzero "$z2/0" if ($r1 == 0);
+               return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
+       } else {
+               divbyzero "$z1/0" if ($r2 == 0);
+               return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
+       }
 }
 
 #
@@ -248,8 +287,9 @@ sub power {
 #
 sub spaceship {
        my ($z1, $z2, $inverted) = @_;
+       $z2 = cplx($z2, 0) unless ref $z2;
        my ($re1, $im1) = @{$z1->cartesian};
-       my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2);
+       my ($re2, $im2) = @{$z2->cartesian};
        my $sgn = $inverted ? -1 : 1;
        return $sgn * ($re1 <=> $re2) if $re1 != $re2;
        return $sgn * ($im1 <=> $im2);
@@ -292,6 +332,7 @@ sub conjugate {
 #
 sub abs {
        my ($z) = @_;
+       return abs($z) unless ref $z;
        my ($r, $t) = @{$z->polar};
        return abs($r);
 }
@@ -303,7 +344,7 @@ sub abs {
 #
 sub arg {
        my ($z) = @_;
-       return 0 unless ref $z;
+       return ($z < 0 ? pi : 0) unless ref $z;
        my ($r, $t) = @{$z->polar};
        return $t;
 }
@@ -311,10 +352,11 @@ sub arg {
 #
 # (sqrt)
 #
-# Compute sqrt(z) (positive only).
+# Compute sqrt(z).
 #
 sub sqrt {
        my ($z) = @_;
+       $z = cplx($z, 0) unless ref $z;
        my ($r, $t) = @{$z->polar};
        return (ref $z)->emake(sqrt($r), $t/2);
 }
@@ -322,11 +364,11 @@ sub sqrt {
 #
 # cbrt
 #
-# Compute cbrt(z) (cubic root, primary only).
+# Compute cbrt(z) (cubic root).
 #
 sub cbrt {
        my ($z) = @_;
-       return $z ** (1/3) unless ref $z;
+       return cplx($z, 0) ** (1/3) unless ref $z;
        my ($r, $t) = @{$z->polar};
        return (ref $z)->emake($r**(1/3), $t/3);
 }
@@ -389,6 +431,7 @@ sub Im {
 #
 sub exp {
        my ($z) = @_;
+       $z = cplx($z, 0) unless ref $z;
        my ($x, $y) = @{$z->cartesian};
        return (ref $z)->emake(exp($x), $y);
 }
@@ -400,21 +443,32 @@ sub exp {
 #
 sub log {
        my ($z) = @_;
+       $z = cplx($z, 0) unless ref $z;
        my ($r, $t) = @{$z->polar};
+       my ($x, $y) = @{$z->cartesian};
+       $t -= 2 * pi if ($t >  pi() and $x < 0);
+       $t += 2 * pi if ($t < -pi() and $x < 0);
        return (ref $z)->make(log($r), $t);
 }
 
 #
+# ln
+#
+# Alias for log().
+#
+sub ln { Math::Complex::log(@_) }
+
+#
 # log10
 #
 # Compute log10(z).
 #
 sub log10 {
        my ($z) = @_;
-       $log10 = log(10) unless defined $log10;
-       return log($z) / $log10 unless ref $z;
+       my $ilog10 = 1 / log(10) unless defined $ilog10;
+       return log(cplx($z, 0)) * $ilog10 unless ref $z;
        my ($r, $t) = @{$z->polar};
-       return (ref $z)->make(log($r) / $log10, $t / $log10);
+       return (ref $z)->make(log($r) * $ilog10, $t * $ilog10);
 }
 
 #
@@ -424,9 +478,10 @@ sub log10 {
 #
 sub logn {
        my ($z, $n) = @_;
+       $z = cplx($z, 0) unless ref $z;
        my $logn = $logn{$n};
        $logn = $logn{$n} = log($n) unless defined $logn;       # Cache log(n)
-       return log($z) / log($n);
+       return log($z) / $logn;
 }
 
 #
@@ -439,7 +494,8 @@ sub cos {
        my ($x, $y) = @{$z->cartesian};
        my $ey = exp($y);
        my $ey_1 = 1 / $ey;
-       return (ref $z)->make(cos($x) * ($ey + $ey_1)/2, sin($x) * ($ey_1 - $ey)/2);
+       return (ref $z)->make(cos($x) * ($ey + $ey_1)/2,
+                             sin($x) * ($ey_1 - $ey)/2);
 }
 
 #
@@ -452,7 +508,8 @@ sub sin {
        my ($x, $y) = @{$z->cartesian};
        my $ey = exp($y);
        my $ey_1 = 1 / $ey;
-       return (ref $z)->make(sin($x) * ($ey + $ey_1)/2, cos($x) * ($ey - $ey_1)/2);
+       return (ref $z)->make(sin($x) * ($ey + $ey_1)/2,
+                             cos($x) * ($ey - $ey_1)/2);
 }
 
 #
@@ -462,29 +519,70 @@ sub sin {
 #
 sub tan {
        my ($z) = @_;
-       return sin($z) / cos($z);
+       my $cz = cos($z);
+       divbyzero "tan($z)", "cos($z)" if ($cz == 0);
+       return sin($z) / $cz;
 }
 
 #
-# cotan
+# sec
+#
+# Computes the secant sec(z) = 1 / cos(z).
+#
+sub sec {
+       my ($z) = @_;
+       my $cz = cos($z);
+       divbyzero "sec($z)", "cos($z)" if ($cz == 0);
+       return 1 / $cz;
+}
+
+#
+# csc
+#
+# Computes the cosecant csc(z) = 1 / sin(z).
+#
+sub csc {
+       my ($z) = @_;
+       my $sz = sin($z);
+       divbyzero "csc($z)", "sin($z)" if ($sz == 0);
+       return 1 / $sz;
+}
+
 #
-# Computes cotan(z) = 1 / tan(z).
+# cosec
 #
-sub cotan {
+# Alias for csc().
+#
+sub cosec { Math::Complex::csc(@_) }
+
+#
+# cot
+#
+# Computes cot(z) = 1 / tan(z).
+#
+sub cot {
        my ($z) = @_;
-       return cos($z) / sin($z);
+       my $sz = sin($z);
+       divbyzero "cot($z)", "sin($z)" if ($sz == 0);
+       return cos($z) / $sz;
 }
 
 #
+# cotan
+#
+# Alias for cot().
+#
+sub cotan { Math::Complex::cot(@_) }
+
+#
 # acos
 #
 # Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
 #
 sub acos {
        my ($z) = @_;
-       my $cz = $z*$z - 1;
-       $cz = cplx($cz, 0) if !ref $cz && $cz < 0;      # Force complex if <0
-       return ~i * log($z + sqrt $cz);                         # ~i is -i
+       $z = cplx($z, 0) unless ref $z;
+       return ~i * log($z + (Re($z) * Im($z) > 0 ? 1 : -1) * sqrt($z*$z - 1));
 }
 
 #
@@ -494,43 +592,80 @@ sub acos {
 #
 sub asin {
        my ($z) = @_;
-       my $cz = 1 - $z*$z;
-       $cz = cplx($cz, 0) if !ref $cz && $cz < 0;      # Force complex if <0
-       return ~i * log(i * $z + sqrt $cz);                     # ~i is -i
+       $z = cplx($z, 0) unless ref $z;
+       return ~i * log(i * $z + sqrt(1 - $z*$z));
 }
 
 #
 # atan
 #
-# Computes the arc tagent atan(z) = i/2 log((i+z) / (i-z)).
+# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
 #
 sub atan {
        my ($z) = @_;
-       return i/2 * log((i + $z) / (i - $z));
+       divbyzero "atan($z)", "i - $z" if ($z == i);
+       return i/2*log((i + $z) / (i - $z));
 }
 
 #
-# acotan
+# asec
+#
+# Computes the arc secant asec(z) = acos(1 / z).
+#
+sub asec {
+       my ($z) = @_;
+       return acos(1 / $z);
+}
+
+#
+# acosec
+#
+# Computes the arc cosecant sec(z) = asin(1 / z).
+#
+sub acosec {
+       my ($z) = @_;
+       return asin(1 / $z);
+}
+
+#
+# acsc
 #
-# Computes the arc cotangent acotan(z) = -i/2 log((i+z) / (z-i))
+# Alias for acosec().
+#
+sub acsc { Math::Complex::acosec(@_) }
+
 #
-sub acotan {
+# acot
+#
+# Computes the arc cotangent acot(z) = -i/2 log((i+z) / (z-i))
+#
+sub acot {
        my ($z) = @_;
+       divbyzero "acot($z)", "$z - i" if ($z == i);
        return i/-2 * log((i + $z) / ($z - i));
 }
 
 #
+# acotan
+#
+# Alias for acot().
+#
+sub acotan { Math::Complex::acot(@_) }
+
+#
 # cosh
 #
 # Computes the hyperbolic cosine cosh(z) = (exp(z) + exp(-z))/2.
 #
 sub cosh {
        my ($z) = @_;
-       my ($x, $y) = ref $z ? @{$z->cartesian} : ($z);
+       $z = cplx($z, 0) unless ref $z;
+       my ($x, $y) = @{$z->cartesian};
        my $ex = exp($x);
        my $ex_1 = 1 / $ex;
        return ($ex + $ex_1)/2 unless ref $z;
-       return (ref $z)->make(cos($y) * ($ex + $ex_1)/2, sin($y) * ($ex - $ex_1)/2);
+       return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
+                             sin($y) * ($ex - $ex_1)/2);
 }
 
 #
@@ -540,11 +675,13 @@ sub cosh {
 #
 sub sinh {
        my ($z) = @_;
-       my ($x, $y) = ref $z ? @{$z->cartesian} : ($z);
+       $z = cplx($z, 0) unless ref $z;
+       my ($x, $y) = @{$z->cartesian};
        my $ex = exp($x);
        my $ex_1 = 1 / $ex;
        return ($ex - $ex_1)/2 unless ref $z;
-       return (ref $z)->make(cos($y) * ($ex - $ex_1)/2, sin($y) * ($ex + $ex_1)/2);
+       return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
+                             sin($y) * ($ex + $ex_1)/2);
 }
 
 #
@@ -554,29 +691,70 @@ sub sinh {
 #
 sub tanh {
        my ($z) = @_;
-       return sinh($z) / cosh($z);
+       my $cz = cosh($z);
+       divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
+       return sinh($z) / $cz;
 }
 
 #
-# cotanh
+# sech
+#
+# Computes the hyperbolic secant sech(z) = 1 / cosh(z).
+#
+sub sech {
+       my ($z) = @_;
+       my $cz = cosh($z);
+       divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
+       return 1 / $cz;
+}
+
+#
+# csch
+#
+# Computes the hyperbolic cosecant csch(z) = 1 / sinh(z).
 #
-# Comptutes the hyperbolic cotangent cotanh(z) = cosh(z) / sinh(z).
+sub csch {
+       my ($z) = @_;
+       my $sz = sinh($z);
+       divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
+       return 1 / $sz;
+}
+
+#
+# cosech
+#
+# Alias for csch().
+#
+sub cosech { Math::Complex::csch(@_) }
+
 #
-sub cotanh {
+# coth
+#
+# Computes the hyperbolic cotangent coth(z) = cosh(z) / sinh(z).
+#
+sub coth {
        my ($z) = @_;
-       return cosh($z) / sinh($z);
+       my $sz = sinh($z);
+       divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
+       return cosh($z) / $sz;
 }
 
 #
+# cotanh
+#
+# Alias for coth().
+#
+sub cotanh { Math::Complex::coth(@_) }
+
+#
 # acosh
 #
 # Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
 #
 sub acosh {
        my ($z) = @_;
-       my $cz = $z*$z - 1;
-       $cz = cplx($cz, 0) if !ref $cz && $cz < 0;      # Force complex if <0
-       return log($z + sqrt $cz);
+       $z = cplx($z, 0) unless ref $z; # asinh(-2)
+       return log($z + sqrt($z*$z - 1));
 }
 
 #
@@ -586,8 +764,8 @@ sub acosh {
 #
 sub asinh {
        my ($z) = @_;
-       my $cz = $z*$z + 1;                                                     # Already complex if <0
-       return log($z + sqrt $cz);
+       $z = cplx($z, 0) unless ref $z; # asinh(-2)
+       return log($z + sqrt($z*$z + 1));
 }
 
 #
@@ -597,24 +775,62 @@ sub asinh {
 #
 sub atanh {
        my ($z) = @_;
+       $z = cplx($z, 0) unless ref $z; # atanh(-2)
+       divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
        my $cz = (1 + $z) / (1 - $z);
-       $cz = cplx($cz, 0) if !ref $cz && $cz < 0;      # Force complex if <0
        return log($cz) / 2;
 }
 
 #
-# acotanh
+# asech
+#
+# Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
+#
+sub asech {
+       my ($z) = @_;
+       divbyzero 'asech(0)', $z if ($z == 0);
+       return acosh(1 / $z);
+}
+
+#
+# acsch
 #
-# Computes the arc hyperbolic cotangent acotanh(z) = 1/2 log((1+z) / (z-1)).
+# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
 #
-sub acotanh {
+sub acsch {
        my ($z) = @_;
+       divbyzero 'acsch(0)', $z if ($z == 0);
+       return asinh(1 / $z);
+}
+
+#
+# acosech
+#
+# Alias for acosh().
+#
+sub acosech { Math::Complex::acsch(@_) }
+
+#
+# acoth
+#
+# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
+#
+sub acoth {
+       my ($z) = @_;
+       $z = cplx($z, 0) unless ref $z; # acoth(-2)
+       divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
        my $cz = (1 + $z) / ($z - 1);
-       $cz = cplx($cz, 0) if !ref $cz && $cz < 0;      # Force complex if <0
        return log($cz) / 2;
 }
 
 #
+# acotanh
+#
+# Alias for acot().
+#
+sub acotanh { Math::Complex::acoth(@_) }
+
+#
 # (atan2)
 #
 # Compute atan(z1/z2).
@@ -622,7 +838,7 @@ sub acotanh {
 sub atan2 {
        my ($z1, $z2, $inverted) = @_;
        my ($re1, $im1) = @{$z1->cartesian};
-       my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2);
+       my ($re2, $im2) = @{$z2->cartesian};
        my $tan;
        if (defined $inverted && $inverted) {   # atan(z2/z1)
                return pi * ($re2 > 0 ? 1 : -1) if $re1 == 0 && $im1 == 0;
@@ -653,7 +869,7 @@ sub display_format {
 
        if (ref $self) {                        # Called as a method
                $format = shift;
-       } else {                                        # Regular procedure call
+       } else {                                # Regular procedure call
                $format = $self;
                undef $self;
        }
@@ -709,7 +925,7 @@ sub stringify_cartesian {
        elsif ($y == -1)                        { $im = '-i' }
        elsif (abs($y) >= 1e-14)        { $im = $y . "i" }
 
-       my $str;
+       my $str = '';
        $str = $re if defined $re;
        $str .= "+$im" if defined $im;
        $str =~ s/\+-/-/;
@@ -728,22 +944,24 @@ sub stringify_polar {
        my $z  = shift;
        my ($r, $t) = @{$z->polar};
        my $theta;
+       my $eps = 1e-14;
 
-       return '[0,0]' if $r <= 1e-14;
+       return '[0,0]' if $r <= $eps;
 
        my $tpi = 2 * pi;
        my $nt = $t / $tpi;
        $nt = ($nt - int($nt)) * $tpi;
        $nt += $tpi if $nt < 0;                 # Range [0, 2pi]
 
-       if (abs($nt) <= 1e-14)                  { $theta = 0 }
-       elsif (abs(pi-$nt) <= 1e-14)    { $theta = 'pi' }
+       if (abs($nt) <= $eps)           { $theta = 0 }
+       elsif (abs(pi-$nt) <= $eps)     { $theta = 'pi' }
 
        if (defined $theta) {
-               $r = int($r + ($r < 0 ? -1 : 1) * 1e-14)
-                       if int(abs($r)) != int(abs($r) + 1e-14);
-               $theta = int($theta + ($theta < 0 ? -1 : 1) * 1e-14)
-                       if int(abs($theta)) != int(abs($theta) + 1e-14);
+               $r = int($r + ($r < 0 ? -1 : 1) * $eps)
+                       if int(abs($r)) != int(abs($r) + $eps);
+               $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
+                       if ($theta ne 'pi' and
+                           int(abs($theta)) != int(abs($theta) + $eps));
                return "\[$r,$theta\]";
        }
 
@@ -756,18 +974,20 @@ sub stringify_polar {
        
        for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
                $n = int($kpi / $nt + ($nt > 0 ? 1 : -1) * 0.5);
-               if (abs($kpi/$n - $nt) <= 1e-14) {
-                       $theta = ($nt < 0 ? '-':'').($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
+               if (abs($kpi/$n - $nt) <= $eps) {
+                       $theta = ($nt < 0 ? '-':'').
+                                ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
                        last;
                }
        }
 
        $theta = $nt unless defined $theta;
 
-       $r = int($r + ($r < 0 ? -1 : 1) * 1e-14)
-               if int(abs($r)) != int(abs($r) + 1e-14);
-       $theta = int($theta + ($theta < 0 ? -1 : 1) * 1e-14)
-               if int(abs($theta)) != int(abs($theta) + 1e-14);
+       $r = int($r + ($r < 0 ? -1 : 1) * $eps)
+               if int(abs($r)) != int(abs($r) + $eps);
+       $theta = int($theta + ($theta < 0 ? -1 : 1) * $eps)
+               if ($theta !~ m(^-?\d*pi/\d+$) and
+                   int(abs($theta)) != int(abs($theta) + $eps));
 
        return "\[$r,$theta\]";
 }
@@ -974,24 +1194,41 @@ numbers:
        logn(z, n) = log(z) / log(n)
 
        tan(z) = sin(z) / cos(z)
-       cotan(z) = 1 / tan(z)
+
+        csc(z) = 1 / sin(z)
+        sec(z) = 1 / cos(z)
+       cot(z) = 1 / tan(z)
 
        asin(z) = -i * log(i*z + sqrt(1-z*z))
        acos(z) = -i * log(z + sqrt(z*z-1))
        atan(z) = i/2 * log((i+z) / (i-z))
-       acotan(z) = -i/2 * log((i+z) / (z-i))
+
+        acsc(z) = asin(1 / z)
+        asec(z) = acos(1 / z)
+       acot(z) = -i/2 * log((i+z) / (z-i))
 
        sinh(z) = 1/2 (exp(z) - exp(-z))
        cosh(z) = 1/2 (exp(z) + exp(-z))
-       tanh(z) = sinh(z) / cosh(z)
-       cotanh(z) = 1 / tanh(z)
+       tanh(z) = sinh(z) / cosh(z) = (exp(z) - exp(-z)) / (exp(z) + exp(-z))
+
+        csch(z) = 1 / sinh(z)
+        sech(z) = 1 / cosh(z)
+       coth(z) = 1 / tanh(z)
        
        asinh(z) = log(z + sqrt(z*z+1))
        acosh(z) = log(z + sqrt(z*z-1))
        atanh(z) = 1/2 * log((1+z) / (1-z))
-       acotanh(z) = 1/2 * log((1+z) / (z-1))
 
-The I<root> function is available to compute all the I<n>th
+        acsch(z) = asinh(1 / z)
+        asech(z) = acosh(1 / z)
+       acoth(z) = atanh(1 / z) = 1/2 * log((1+z) / (z-1))
+
+I<log>, I<csc>, I<cot>, I<acsc>, I<acot>, I<csch>, I<coth>,
+I<acosech>, I<acotanh>, have aliases I<ln>, I<cosec>, I<cotan>,
+I<acosec>, I<acotan>, I<cosech>, I<cotanh>, I<acosech>, I<acotanh>,
+respectively.
+
+The I<root> function is available to compute all the I<n>
 roots of some complex, where I<n> is a strictly positive integer.
 There are exactly I<n> such roots, returned as a list. Getting the
 number mathematicians call C<j> such that:
@@ -1006,10 +1243,11 @@ 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)
 
-The I<spaceshift> operation 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> comparison operator 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
 
@@ -1027,9 +1265,9 @@ if you like. To create a number using the trigonometric form, use either:
        $z = Math::Complex->emake(5, pi/3);
        $x = cplxe(5, pi/3);
 
-instead. The first argument is the modulus, the second is the angle (in radians).
-(Mnmemonic: C<e> is used as a notation for complex numbers in the trigonometric
-form).
+instead. The first argument is the modulus, the second is the angle
+(in radians, the full circle is 2*pi).  (Mnmemonic: C<e> is used as a
+notation for complex numbers in the trigonometric form).
 
 It is possible to write:
 
@@ -1104,6 +1342,7 @@ All routines expect to be given real or complex numbers. Don't attempt to
 use BigFloat, since Perl has currently no rule to disambiguate a '+'
 operation (for instance) between two overloaded entities.
 
-=head1 AUTHOR
+=head1 AUTHORS
 
-Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
+       Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>>
+       Jarkko Hietaniemi <F<jhi@iki.fi>>
index 1ffd7d5..5dcac6c 100755 (executable)
@@ -3,7 +3,8 @@
 # $RCSfile$
 #
 # Regression tests for the new Math::Complex pacakge
-# -- Raphael Manfredi, Sept 1996
+# -- Raphael Manfredi, Septemeber 1996
+# -- Jarkko Hietaniemi Manfredi, March 1997
 BEGIN {
     chdir 't' if -d 't';
     @INC = '../lib';
@@ -13,12 +14,13 @@ use Math::Complex;
 $test = 0;
 $| = 1;
 $script = '';
-$epsilon = 1e-10;
+my $eps = 1e-4; # for example root() is quite bad
 
 while (<DATA>) {
-       next if /^#/ || /^\s*$/;
-       chop;
-       $set_test = 0;                  # Assume not a test over a set of values
+       s/^\s+//;
+       next if $_ eq '' || /^\#/;
+       chomp;
+       $test_set = 0;          # Assume not a test over a set of values
        if (/^&(.*)/) {
                $op = $1;
                next;
@@ -28,14 +30,16 @@ while (<DATA>) {
                next;
        }
        elsif (s/^\|//) {
-               $set_test = 1;          # Requests we loop over the set...
+               $test_set = 1;  # Requests we loop over the set...
        }
        my @args = split(/:/);
-       if ($set_test) {
+       if ($test_set == 1) {
                my $i;
                for ($i = 0; $i < @set; $i++) {
-                       $target = $set[$i];             # complex number
-                       $zvalue = $val[$i];             # textual value as found in set definition
+                       # complex number
+                       $target = $set[$i];
+                       # textual value as found in set definition
+                       $zvalue = $val[$i];
                        test($zvalue, $target, @args);
                }
        } else {
@@ -57,7 +61,7 @@ sub test {
        }
        if (defined $z) {
                $args = "'$op'";                # Really the value
-               $try = "abs(\$z0 - \$z1) <= 1e-10 ? \$z1 : \$z0";
+               $try = "abs(\$z0 - \$z1) <= $eps ? \$z1 : \$z0";
                $script .= "\$res = $try; ";
                $script .= "check($test, $args[0], \$res, \$z$#args, $args);\n";
        } else {
@@ -114,7 +118,15 @@ sub value {
 
 sub check {
        my ($test, $try, $got, $expected, @z) = @_;
-       if ("$got" eq "$expected" || ($expected =~ /^-?\d/ && $got == $expected)) {
+
+#      print "# @_\n";
+
+       if ("$got" eq "$expected"
+           ||
+           ($expected =~ /^-?\d/ && $got == $expected)
+           ||
+           (abs($got - $expected) < $eps)
+           ) {
                print "ok $test\n";
        } else {
                print "not ok $test\n";
@@ -163,10 +175,24 @@ __END__
 [4,pi]:[2,pi/2]:[2,pi/2]
 [2,pi/2]:[4,pi]:[0.5,-(pi)/2]
 
+&Re
+(3,4):3
+(-3,4):-3
+[1,pi/2]:0
+
+&Im
+(3,4):4
+(3,-4):-4
+[1,pi/2]:1
+
 &abs
 (3,4):5
 (-3,4):5
 
+&arg
+[2,0]:0
+[-2,0]:pi
+
 &~
 (4,5):(4,-5)
 (-3,4):(-3,-4)
@@ -185,6 +211,7 @@ __END__
 (3,4):(3,4):1
 
 &sqrt
+-9:(0,3)
 (-100,0):(0,10)
 (16,-30):(5,-3)
 
@@ -209,46 +236,214 @@ __END__
 |'z - ~z':'2*i*Im(z)'
 |'z * ~z':'abs(z) * abs(z)'
 
-{ (4,3); [3,2]; (-3,4); (0,2); 3; 1; (-5, 0); [2,1] }
+{ (2,3); [3,2]; (-3,2); (0,2); 3; 1.2; -3; (-3, 0); (-2, -1); [2,1] }
 
-|'exp(z)':'exp(a) * exp(i * b)'
+|'(root(z, 4))[1] ** 4':'z'
+|'(root(z, 5))[3] ** 5':'z'
+|'(root(z, 8))[7] ** 8':'z'
 |'abs(z)':'r'
-|'sqrt(z) * sqrt(z)':'z'
-|'sqrt(z)':'sqrt(r) * exp(i * t/2)'
+|'acot(z)':'acotan(z)'
+|'acsc(z)':'acosec(z)'
+|'acsc(z)':'asin(1 / z)'
+|'asec(z)':'acos(1 / z)'
 |'cbrt(z)':'cbrt(r) * exp(i * t/3)'
-|'log(z)':'log(r) + i*t'
-|'sin(asin(z))':'z'
 |'cos(acos(z))':'z'
-|'tan(atan(z))':'z'
-|'cotan(acotan(z))':'z'
 |'cos(z) ** 2 + sin(z) ** 2':1
-|'cosh(z) ** 2 - sinh(z) ** 2':1
 |'cos(z)':'cosh(i*z)'
-|'cotan(z)':'1 / tan(z)'
-|'cotanh(z)':'1 / tanh(z)'
-|'i*sin(z)':'sinh(i*z)'
-|'z**z':'exp(z * log(z))'
-|'log(exp(z))':'z'
+|'cosh(z) ** 2 - sinh(z) ** 2':1
+|'cot(acot(z))':'z'
+|'cot(z)':'1 / tan(z)'
+|'cot(z)':'cotan(z)'
+|'csc(acsc(z))':'z'
+|'csc(z)':'1 / sin(z)'
+|'csc(z)':'cosec(z)'
 |'exp(log(z))':'z'
+|'exp(z)':'exp(a) * exp(i * b)'
+|'ln(z)':'log(z)'
+|'log(exp(z))':'z'
+|'log(z)':'log(r) + i*t'
 |'log10(z)':'log(z) / log(10)'
-|'logn(z, 3)':'log(z) / log(3)'
 |'logn(z, 2)':'log(z) / log(2)'
-|'(root(z, 4))[1] ** 4':'z'
-|'(root(z, 8))[7] ** 8':'z'
+|'logn(z, 3)':'log(z) / log(3)'
+|'sec(asec(z))':'z'
+|'sec(z)':'1 / cos(z)'
+|'sin(asin(z))':'z'
+|'sin(i * z)':'i * sinh(z)'
+|'sqrt(z) * sqrt(z)':'z'
+|'sqrt(z)':'sqrt(r) * exp(i * t/2)'
+|'tan(atan(z))':'z'
+|'z**z':'exp(z * log(z))'
 
-{ (1,1); [1,0.5]; (-2, -1); 2; (-1,0.5); (0,0.5); 0.5; (2, 0) }
+{ (1,1); [1,0.5]; (-2, -1); 2; -3; (-1,0.5); (0,0.5); 0.5; (2, 0); (-1, -2) }
 
-|'sinh(asinh(z))':'z'
 |'cosh(acosh(z))':'z'
+|'coth(acoth(z))':'z'
+|'coth(z)':'1 / tanh(z)'
+|'coth(z)':'cotanh(z)'
+|'csch(acsch(z))':'z'
+|'csch(z)':'1 / sinh(z)'
+|'csch(z)':'cosech(z)'
+|'sech(asech(z))':'z'
+|'sech(z)':'1 / cosh(z)'
+|'sinh(asinh(z))':'z'
 |'tanh(atanh(z))':'z'
-|'cotanh(acotanh(z))':'z'
 
-{ (0.2,-0.4); [1,0.5]; -1.2; (-1,0.5); (0,-0.5); 0.5; (1.1, 0) }
+{ (0.2,-0.4); [1,0.5]; -1.2; (-1,0.5); 0.5; (1.1, 0) }
 
-|'asin(sin(z))':'z'
 |'acos(cos(z)) ** 2':'z * z'
-|'atan(tan(z))':'z'
-|'asinh(sinh(z))':'z'
 |'acosh(cosh(z)) ** 2':'z * z'
+|'acoth(z)':'acotanh(z)'
+|'acoth(z)':'atanh(1 / z)'
+|'acsch(z)':'acosech(z)'
+|'acsch(z)':'asinh(1 / z)'
+|'asech(z)':'acosh(1 / z)'
+|'asin(sin(z))':'z'
+|'asinh(sinh(z))':'z'
+|'atan(tan(z))':'z'
 |'atanh(tanh(z))':'z'
 
+&sin
+( 2, 3):(  9.15449914691143, -4.16890695996656)
+(-2, 3):( -9.15449914691143, -4.16890695996656)
+(-2,-3):( -9.15449914691143,  4.16890695996656)
+( 2,-3):(  9.15449914691143,  4.16890695996656)
+
+&cos
+( 2, 3):( -4.18962569096881, -9.10922789375534)
+(-2, 3):( -4.18962569096881,  9.10922789375534)
+(-2,-3):( -4.18962569096881, -9.10922789375534)
+( 2,-3):( -4.18962569096881,  9.10922789375534)
+
+&tan
+( 2, 3):( -0.00376402564150,  1.00323862735361)
+(-2, 3):(  0.00376402564150,  1.00323862735361)
+(-2,-3):(  0.00376402564150, -1.00323862735361)
+( 2,-3):( -0.00376402564150, -1.00323862735361)
+
+&sec
+( 2, 3):( -0.04167496441114,  0.09061113719624)
+(-2, 3):( -0.04167496441114, -0.09061113719624)
+(-2,-3):( -0.04167496441114,  0.09061113719624)
+( 2,-3):( -0.04167496441114, -0.09061113719624)
+
+&csc
+( 2, 3):(  0.09047320975321,  0.04120098628857)
+(-2, 3):( -0.09047320975321,  0.04120098628857)
+(-2,-3):( -0.09047320975321, -0.04120098628857)
+( 2,-3):(  0.09047320975321, -0.04120098628857)
+
+&cot
+( 2, 3):( -0.00373971037634, -0.99675779656936)
+(-2, 3):(  0.00373971037634, -0.99675779656936)
+(-2,-3):(  0.00373971037634,  0.99675779656936)
+( 2,-3):( -0.00373971037634,  0.99675779656936)
+
+&asin
+( 2, 3):(  0.57065278432110,  1.98338702991654)
+(-2, 3):( -0.57065278432110,  1.98338702991654)
+(-2,-3):( -0.57065278432110, -1.98338702991654)
+( 2,-3):(  0.57065278432110, -1.98338702991654)
+
+&acos
+( 2, 3):(  1.00014354247380, -1.98338702991654)
+(-2, 3):(  2.14144911111600, -1.98338702991654)
+(-2,-3):(  2.14144911111600,  1.98338702991654)
+( 2,-3):(  1.00014354247380,  1.98338702991654)
+
+&atan
+( 2, 3):(  1.40992104959658,  0.22907268296854)
+(-2, 3):( -1.40992104959658,  0.22907268296854)
+(-2,-3):( -1.40992104959658, -0.22907268296854)
+( 2,-3):(  1.40992104959658, -0.22907268296854)
+
+&asec
+( 2, 3):(  1.42041072246703,  0.23133469857397)
+(-2, 3):(  1.72118193112276,  0.23133469857397)
+(-2,-3):(  1.72118193112276, -0.23133469857397)
+( 2,-3):(  1.42041072246703, -0.23133469857397)
+
+&acsc
+( 2, 3):(  0.15038560432786, -0.23133469857397)
+(-2, 3):( -0.15038560432786, -0.23133469857397)
+(-2,-3):( -0.15038560432786,  0.23133469857397)
+( 2,-3):(  0.15038560432786,  0.23133469857397)
+
+&acot
+( 2, 3):(  0.16087527719832, -0.22907268296854)
+(-2, 3):( -0.16087527719832, -0.22907268296854)
+(-2,-3):( -0.16087527719832,  0.22907268296854)
+( 2,-3):(  0.16087527719832,  0.22907268296854)
+
+&sinh
+( 2, 3):( -3.59056458998578,  0.53092108624852)
+(-2, 3):(  3.59056458998578,  0.53092108624852)
+(-2,-3):(  3.59056458998578, -0.53092108624852)
+( 2,-3):( -3.59056458998578, -0.53092108624852)
+
+&cosh
+( 2, 3):( -3.72454550491532,  0.51182256998738)
+(-2, 3):( -3.72454550491532, -0.51182256998738)
+(-2,-3):( -3.72454550491532,  0.51182256998738)
+( 2,-3):( -3.72454550491532, -0.51182256998738)
+
+&tanh
+( 2, 3):(  0.96538587902213, -0.00988437503832)
+(-2, 3):( -0.96538587902213, -0.00988437503832)
+(-2,-3):( -0.96538587902213,  0.00988437503832)
+( 2,-3):(  0.96538587902213,  0.00988437503832)
+
+&sech
+( 2, 3):( -0.26351297515839, -0.03621163655877)
+(-2, 3):( -0.26351297515839,  0.03621163655877)
+(-2,-3):( -0.26351297515839, -0.03621163655877)
+( 2,-3):( -0.26351297515839,  0.03621163655877)
+
+&csch
+( 2, 3):( -0.27254866146294, -0.04030057885689)
+(-2, 3):(  0.27254866146294, -0.04030057885689)
+(-2,-3):(  0.27254866146294,  0.04030057885689)
+( 2,-3):( -0.27254866146294,  0.04030057885689)
+
+&coth
+( 2, 3):(  1.03574663776500,  0.01060478347034)
+(-2, 3):( -1.03574663776500,  0.01060478347034)
+(-2,-3):( -1.03574663776500, -0.01060478347034)
+( 2,-3):(  1.03574663776500, -0.01060478347034)
+
+&asinh
+( 2, 3):(  1.96863792579310,  0.96465850440760)
+(-2, 3):( -1.96863792579310,  0.96465850440761)
+(-2,-3):( -1.96863792579310, -0.96465850440761)
+( 2,-3):(  1.96863792579310, -0.96465850440760)
+
+&acosh
+( 2, 3):(  1.98338702991654,  1.00014354247380)
+(-2, 3):( -1.98338702991653, -2.14144911111600)
+(-2,-3):( -1.98338702991653,  2.14144911111600)
+( 2,-3):(  1.98338702991654, -1.00014354247380)
+
+&atanh
+( 2, 3):(  0.14694666622553,  1.33897252229449)
+(-2, 3):( -0.14694666622553,  1.33897252229449)
+(-2,-3):( -0.14694666622553, -1.33897252229449)
+( 2,-3):(  0.14694666622553, -1.33897252229449)
+
+&asech
+( 2, 3):(  0.23133469857397, -1.42041072246703)
+(-2, 3):( -0.23133469857397,  1.72118193112276)
+(-2,-3):( -0.23133469857397, -1.72118193112276)
+( 2,-3):(  0.23133469857397,  1.42041072246703)
+
+&acsch
+( 2, 3):(  0.15735549884499, -0.22996290237721)
+(-2, 3):( -0.15735549884499, -0.22996290237721)
+(-2,-3):( -0.15735549884499,  0.22996290237721)
+( 2,-3):(  0.15735549884499,  0.22996290237721)
+
+&acoth
+( 2, 3):(  0.14694666622553, -0.23182380450040)
+(-2, 3):( -0.14694666622553, -0.23182380450040)
+(-2,-3):( -0.14694666622553,  0.23182380450040)
+( 2,-3):(  0.14694666622553,  0.23182380450040)
+
+# eof