Math::Complex update

2 files changed, 565 insertions, 131 deletions

diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm
index fce53f7b81..30194eb9a7 100644
--- a/lib/Math/Complex.pm
+++ b/lib/Math/Complex.pm
@@ -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>>
diff --git a/t/lib/complex.t b/t/lib/complex.t
index 1ffd7d5447..5dcac6c6ce 100755
--- a/t/lib/complex.t
+++ b/t/lib/complex.t
@@ -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