Complex update (five patches)

2 files changed, 159 insertions, 59 deletions

diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm
index 7d5a014203..a501b99ddd 100644
--- a/lib/Math/Complex.pm
+++ b/lib/Math/Complex.pm
@@ -207,6 +207,7 @@ sub update_polar {
 sub plus {
 	my ($z1, $z2, $regular) = @_;
 	my ($re1, $im1) = @{$z1->cartesian};
+	$z2 = cplx($z2) unless ref $z2;
 	my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
 	unless (defined $regular) {
 		$z1->set_cartesian([$re1 + $re2, $im1 + $im2]);
@@ -223,7 +224,8 @@ sub plus {
 sub minus {
 	my ($z1, $z2, $inverted) = @_;
 	my ($re1, $im1) = @{$z1->cartesian};
-	my ($re2, $im2) = ref $z2 ? @{$z2->cartesian} : ($z2, 0);
+	$z2 = cplx($z2) unless ref $z2;
+	my ($re2, $im2) = @{$z2->cartesian};
 	unless (defined $inverted) {
 		$z1->set_cartesian([$re1 - $re2, $im1 - $im2]);
 		return $z1;
@@ -231,6 +233,7 @@ sub minus {
 	return $inverted ?
 		(ref $z1)->make($re2 - $re1, $im2 - $im1) :
 		(ref $z1)->make($re1 - $re2, $im1 - $im2);
+
 }
 
 #
@@ -241,8 +244,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);
+	$z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2;
+	my ($r2, $t2) = @{$z2->polar};
 	unless (defined $regular) {
 		$z1->set_polar([$r1 * $r2, $t1 + $t2]);
 		return $z1;
@@ -251,11 +254,11 @@ sub multiply {
 }
 
 #
-# divbyzero
+# _divbyzero
 #
 # Die on division by zero.
 #
-sub divbyzero {
+sub _divbyzero {
     my $mess = "$_[0]: Division by zero.\n";
 
     if (defined $_[1]) {
@@ -272,21 +275,6 @@ sub divbyzero {
 }
 
 #
-# zerotozero
-#
-# Die on zero raised to the zeroth.
-#
-sub zerotozero {
-    my $mess = "The zero raised to the zeroth power is not defined.\n";
-
-    my @up = caller(1);
-    
-    $mess .= "Died at $up[1] line $up[2].\n";
-
-    die $mess;
-}
-
-#
 # (divide)
 #
 # Computes z1/z2.
@@ -294,32 +282,55 @@ sub zerotozero {
 sub divide {
 	my ($z1, $z2, $inverted) = @_;
 	my ($r1, $t1) = @{$z1->polar};
-	my ($r2, $t2) = ref $z2 ?
-	    @{$z2->polar} : (abs($z2), $z2 >= 0 ? 0 : pi);
+	$z2 = cplxe(abs($z2), $z2 >= 0 ? 0 : pi) unless ref $z2;
+	my ($r2, $t2) = @{$z2->polar};
 	unless (defined $inverted) {
-		divbyzero "$z1/0" if ($r2 == 0);
+		_divbyzero "$z1/0" if ($r2 == 0);
 		$z1->set_polar([$r1 / $r2, $t1 - $t2]);
 		return $z1;
 	}
 	if ($inverted) {
-		divbyzero "$z2/0" if ($r1 == 0);
+		_divbyzero "$z2/0" if ($r1 == 0);
 		return (ref $z1)->emake($r2 / $r1, $t2 - $t1);
 	} else {
-		divbyzero "$z1/0" if ($r2 == 0);
+		_divbyzero "$z1/0" if ($r2 == 0);
 		return (ref $z1)->emake($r1 / $r2, $t1 - $t2);
 	}
 }
 
 #
+# _zerotozero
+#
+# Die on zero raised to the zeroth.
+#
+sub _zerotozero {
+    my $mess = "The zero raised to the zeroth power is not defined.\n";
+
+    my @up = caller(1);
+    
+    $mess .= "Died at $up[1] line $up[2].\n";
+
+    die $mess;
+}
+
+#
 # (power)
 #
 # Computes z1**z2 = exp(z2 * log z1)).
 #
 sub power {
 	my ($z1, $z2, $inverted) = @_;
-	zerotozero if ($z1 == 0 and $z2 == 0);
-	return exp($z1 * log $z2) if defined $inverted && $inverted;
-	return exp($z2 * log $z1);
+	_zerotozero if ($z1 == 0 and $z2 == 0);
+	return 1 if ($z1 == 1);
+	return 0 if ($z2 == 0);
+	$z2 = cplx($z2) unless ref $z2;
+	unless (defined $inverted) {
+		my $z3 = exp($z2 * log $z1);
+		$z1->set_cartesian([@{$z3->cartesian}]);
+		return $z1;
+	}
+	return exp($z2 * log $z1) unless $inverted;
+	return exp($z1 * log $z2);
 }
 
 #
@@ -416,6 +427,21 @@ sub cbrt {
 }
 
 #
+# _rootbad
+#
+# Die on bad root.
+#
+sub _rootbad {
+    my $mess = "Root $_[0] not defined, root must be positive integer.\n";
+
+    my @up = caller(1);
+    
+    $mess .= "Died at $up[1] line $up[2].\n";
+
+    die $mess;
+}
+
+#
 # root
 #
 # Computes all nth root for z, returning an array whose size is n.
@@ -427,8 +453,7 @@ sub cbrt {
 #
 sub root {
 	my ($z, $n) = @_;
-	$n = int($n + 0.5);
-	return undef unless $n > 0;
+	_rootbad($n) if ($n < 1 or int($n) != $n);
 	my ($r, $t) = ref $z ? @{$z->polar} : (abs($z), $z >= 0 ? 0 : pi);
 	my @root;
 	my $k;
@@ -565,7 +590,7 @@ sub sin {
 sub tan {
 	my ($z) = @_;
 	my $cz = cos($z);
-	divbyzero "tan($z)", "cos($z)" if ($cz == 0);
+	_divbyzero "tan($z)", "cos($z)" if ($cz == 0);
 	return sin($z) / $cz;
 }
 
@@ -577,7 +602,7 @@ sub tan {
 sub sec {
 	my ($z) = @_;
 	my $cz = cos($z);
-	divbyzero "sec($z)", "cos($z)" if ($cz == 0);
+	_divbyzero "sec($z)", "cos($z)" if ($cz == 0);
 	return 1 / $cz;
 }
 
@@ -589,7 +614,7 @@ sub sec {
 sub csc {
 	my ($z) = @_;
 	my $sz = sin($z);
-	divbyzero "csc($z)", "sin($z)" if ($sz == 0);
+	_divbyzero "csc($z)", "sin($z)" if ($sz == 0);
 	return 1 / $sz;
 }
 
@@ -608,7 +633,7 @@ sub cosec { Math::Complex::csc(@_) }
 sub cot {
 	my ($z) = @_;
 	my $sz = sin($z);
-	divbyzero "cot($z)", "sin($z)" if ($sz == 0);
+	_divbyzero "cot($z)", "sin($z)" if ($sz == 0);
 	return cos($z) / $sz;
 }
 
@@ -649,7 +674,7 @@ sub asin {
 sub atan {
 	my ($z) = @_;
 	$z = cplx($z, 0) unless ref $z;
-	divbyzero "atan($z)", "i - $z" if ($z == i);
+	_divbyzero "atan($z)", "i - $z" if ($z == i);
 	return i/2*log((i + $z) / (i - $z));
 }
 
@@ -660,7 +685,7 @@ sub atan {
 #
 sub asec {
 	my ($z) = @_;
-	divbyzero "asec($z)", $z if ($z == 0);
+	_divbyzero "asec($z)", $z if ($z == 0);
 	return acos(1 / $z);
 }
 
@@ -671,7 +696,7 @@ sub asec {
 #
 sub acsc {
 	my ($z) = @_;
-	divbyzero "acsc($z)", $z if ($z == 0);
+	_divbyzero "acsc($z)", $z if ($z == 0);
 	return asin(1 / $z);
 }
 
@@ -690,7 +715,7 @@ sub acosec { Math::Complex::acsc(@_) }
 sub acot {
 	my ($z) = @_;
 	$z = cplx($z, 0) unless ref $z;
-	divbyzero "acot($z)", "$z - i" if ($z == i);
+	_divbyzero "acot($z)", "$z - i" if ($z == i);
 	return i/-2 * log((i + $z) / ($z - i));
 }
 
@@ -708,10 +733,15 @@ sub acotan { Math::Complex::acot(@_) }
 #
 sub cosh {
 	my ($z) = @_;
-	my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
+	my $real;
+	unless (ref $z) {
+	    $z = cplx($z, 0);
+	    $real = 1;
+	}
+	my ($x, $y) = @{$z->cartesian};
 	my $ex = exp($x);
 	my $ex_1 = 1 / $ex;
-	return ($ex + $ex_1)/2 unless ref $z;
+	return ($ex + $ex_1)/2 if $real;
 	return (ref $z)->make(cos($y) * ($ex + $ex_1)/2,
 			      sin($y) * ($ex - $ex_1)/2);
 }
@@ -723,10 +753,15 @@ sub cosh {
 #
 sub sinh {
 	my ($z) = @_;
-	my ($x, $y) = ref $z ? @{$z->cartesian} : ($z, 0);
+	my $real;
+	unless (ref $z) {
+	    $z = cplx($z, 0);
+	    $real = 1;
+	}
+	my ($x, $y) = @{$z->cartesian};
 	my $ex = exp($x);
 	my $ex_1 = 1 / $ex;
-	return ($ex - $ex_1)/2 unless ref $z;
+	return ($ex - $ex_1)/2 if $real;
 	return (ref $z)->make(cos($y) * ($ex - $ex_1)/2,
 			      sin($y) * ($ex + $ex_1)/2);
 }
@@ -739,7 +774,7 @@ sub sinh {
 sub tanh {
 	my ($z) = @_;
 	my $cz = cosh($z);
-	divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
+	_divbyzero "tanh($z)", "cosh($z)" if ($cz == 0);
 	return sinh($z) / $cz;
 }
 
@@ -751,7 +786,7 @@ sub tanh {
 sub sech {
 	my ($z) = @_;
 	my $cz = cosh($z);
-	divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
+	_divbyzero "sech($z)", "cosh($z)" if ($cz == 0);
 	return 1 / $cz;
 }
 
@@ -763,7 +798,7 @@ sub sech {
 sub csch {
 	my ($z) = @_;
 	my $sz = sinh($z);
-	divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
+	_divbyzero "csch($z)", "sinh($z)" if ($sz == 0);
 	return 1 / $sz;
 }
 
@@ -782,7 +817,7 @@ sub cosech { Math::Complex::csch(@_) }
 sub coth {
 	my ($z) = @_;
 	my $sz = sinh($z);
-	divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
+	_divbyzero "coth($z)", "sinh($z)" if ($sz == 0);
 	return cosh($z) / $sz;
 }
 
@@ -822,7 +857,7 @@ sub asinh {
 #
 sub atanh {
 	my ($z) = @_;
-	divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
+	_divbyzero 'atanh(1)', "1 - $z" if ($z == 1);
 	$z = cplx($z, 0) unless ref $z;
 	my $cz = (1 + $z) / (1 - $z);
 	return log($cz) / 2;
@@ -835,7 +870,7 @@ sub atanh {
 #
 sub asech {
 	my ($z) = @_;
-	divbyzero 'asech(0)', $z if ($z == 0);
+	_divbyzero 'asech(0)', $z if ($z == 0);
 	return acosh(1 / $z);
 }
 
@@ -846,7 +881,7 @@ sub asech {
 #
 sub acsch {
 	my ($z) = @_;
-	divbyzero 'acsch(0)', $z if ($z == 0);
+	_divbyzero 'acsch(0)', $z if ($z == 0);
 	return asinh(1 / $z);
 }
 
@@ -864,7 +899,7 @@ sub acosech { Math::Complex::acsch(@_) }
 #
 sub acoth {
 	my ($z) = @_;
-	divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
+	_divbyzero 'acoth(1)', "$z - 1" if ($z == 1);
 	$z = cplx($z, 0) unless ref $z;
 	my $cz = (1 + $z) / ($z - 1);
 	return log($cz) / 2;
diff --git a/t/lib/complex.t b/t/lib/complex.t
index 310e6f5e3c..dfc6cb7eca 100755
--- a/t/lib/complex.t
+++ b/t/lib/complex.t
@@ -23,11 +23,11 @@ while (<DATA>) {
 	next if $_ eq '' || /^\#/;
 	chomp;
 	$test_set = 0;		# Assume not a test over a set of values
-	if (/^&(.*)/) {
+	if (/^&(.+)/) {
 		$op = $1;
 		next;
 	}
-	elsif (/^\{(.*)\}/) {
+	elsif (/^\{(.+)\}/) {
 		set($1, \@set, \@val);
 		next;
 	}
@@ -51,6 +51,17 @@ while (<DATA>) {
 
 # test the divbyzeros
 
+sub test_dbz {
+    for my $op (@_) {
+	$test++;
+
+#	push(@script, qq(print "# '$op'\n";));
+	push(@script, qq(eval '$op';));
+	push(@script, qq(print 'not ' unless (\$@ =~ /Division by zero/);));
+	push(@script, qq(print "ok $test\n";));
+    }
+}
+
 test_dbz(
 	 'i/0',
 #	 'tan(pi/2)',	# may succeed thanks to floating point inaccuracies
@@ -71,24 +82,50 @@ test_dbz(
 	 'acoth(1)'
 	);
 
-sub test_dbz {
+# test the 0**0
+
+sub test_ztz {
+	$test++;
+
+#	push(@script, qq(print "# 0**0\n";));
+	push(@script, qq(eval 'cplx(0)**cplx(0)';));
+	push(@script, qq(print 'not ' unless (\$@ =~ /zero raised to the/);));
+	push(@script, qq(print "ok $test\n";));
+}
+
+test_ztz;
+
+# test the bad roots
+
+sub test_broot {
     for my $op (@_) {
 	$test++;
 
-	push(@script, qq(eval '$op';));
-	push(@script, qq(print 'not ' unless (\$@ =~ /Division by zero/);));
+#	push(@script, qq(print "# root(2, $op)\n";));
+	push(@script, qq(eval 'root(2, $op)';));
+	push(@script, qq(print 'not ' unless (\$@ =~ /root must be/);));
 	push(@script, qq(print "ok $test\n";));
     }
 }
 
+test_broot(qw(-3 -2.1 0 0.99));
+
 print "1..$test\n";
 eval join '', @script;
 die $@ if $@;
 
+sub abop {
+	my ($op) = @_;
+
+	push(@script, qq(print "# $op=\n";));
+}
+
 sub test {
 	my ($op, $z, @args) = @_;
+	my ($baop) = 0;
 	$test++;
 	my $i;
+	$baop = 1 if ($op =~ s/;=$//);
 	for ($i = 0; $i < @args; $i++) {
 		$val = value($args[$i]);
 		push @script, "\$z$i = $val;\n";
@@ -109,6 +146,28 @@ sub test {
 		}
 		push @script, "\$res = $try; ";
 		push @script, "check($test, '$try', \$res, \$z$#args, $args);\n";
+		if (@args > 2 and $baop) { # binary assignment ops
+			$test++;
+			# check the op= works
+			push @script, <<EOB;
+{
+        my \$za = cplx(ref \$z0 ? \@{\$z0->cartesian} : (\$z0, 0));
+
+	my (\$z1r, \$z1i) = ref \$z1 ? \@{\$z1->cartesian} : (\$z1, 0);
+
+        my \$zb = cplx(\$z1r, \$z1i);
+
+	\$za $op= \$zb;
+	my (\$zbr, \$zbi) = \@{\$zb->cartesian};
+
+	check($test, '\$z0 $op= \$z1', \$za, \$z$#args, $args);
+EOB
+			$test++;
+			# check that the rhs has not changed
+			push @script, qq(print "not " unless (\$zbr == \$z1r and \$zbi == \$z1i););
+			push @script, qq(print "ok $test\n";);
+			push @script, "}\n";
+		}
 	}
 }
 
@@ -169,7 +228,7 @@ sub check {
 	}
 }
 __END__
-&+
+&+;=
 (3,4):(3,4):(6,8)
 (-3,4):(3,-4):(0,0)
 (3,4):-3:(0,4)
@@ -179,19 +238,20 @@ __END__
 &++
 (2,1):(3,1)
 
-&-
+&-;=
 (2,3):(-2,-3)
 [2,pi/2]:[2,-(pi)/2]
 2:[2,0]:(0,0)
 [3,0]:2:(1,0)
 3:(4,5):(-1,-5)
 (4,5):3:(1,5)
+(2,1):(3,5):(-1,-4)
 
 &--
 (1,2):(0,2)
 [2,pi]:[3,pi]
 
-&*
+&*;=
 (0,1):(0,1):(-1,0)
 (4,5):(1,0):(4,5)
 [2,2*pi/3]:(1,0):[2,2*pi/3]
@@ -200,7 +260,7 @@ __END__
 (0,1):(4,1):(-1,4)
 (2,1):(4,-1):(9,2)
 
-&/
+&/;=
 (3,4):(3,4):(1,0)
 (4,-5):1:(4,-5)
 1:(0,1):(0,-1)
@@ -209,6 +269,11 @@ __END__
 [4,pi]:[2,pi/2]:[2,pi/2]
 [2,pi/2]:[4,pi]:[0.5,-(pi)/2]
 
+&**;=
+(2,0):(3,0):(8,0)
+(3,0):(2,0):(9,0)
+(2,3):(4,0):(-119,-120)
+
 &Re
 (3,4):3
 (-3,4):-3