Complex update (five patches)

diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm
index 7d5a014..a501b99 100644 (file)
--- 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 310e6f5..dfc6cb7 100755 (executable)
--- 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