From: Jarkko Hietaniemi <jhi@iki.fi>
Date: Wed, 24 Jun 1998 15:19:05 +0000 (+0300)
Subject: Complex.pm update
X-Git-Url: http://git.shadowcat.co.uk/gitweb/gitweb.cgi?a=commitdiff_plain;h=d09ae4e65fd5bf88945016a3fc05dfaedbf59acc;p=p5sagit%2Fp5-mst-13.2.git

Complex.pm update
	Message-Id: <199806241219.PAA04061@alpha.hut.fi>
	Subject: [PATCH] 5.004_68: Complex.pm, complex.t

p4raw-id: //depot/perl@1235
---

diff --git a/lib/Math/Complex.pm b/lib/Math/Complex.pm
index 36ca060..aca85c6 100644
--- a/lib/Math/Complex.pm
+++ b/lib/Math/Complex.pm
@@ -196,6 +196,14 @@ use constant pit2 => 2 * pi;
 use constant pip2 => pi / 2;
 
 #
+# deg1
+#
+# One degree in radians, used in stringify_polar.
+#
+
+use constant deg1 => pi / 180;
+
+#
 # uplog10
 #
 # Used in log10().
@@ -424,7 +432,11 @@ sub power {
 	    return 0 if ($z1z);
 	    return 1 if ($z2z or $z1 == 1);
 	}
-	return $inverted ? exp($z1 * log $z2) : exp($z2 * log $z1);
+	my $w = $inverted ? exp($z1 * log $z2) : exp($z2 * log $z1);
+	# If both arguments cartesian, return cartesian, else polar.
+	return $z1->{c_dirty} == 0 &&
+	       (not ref $z2 or $z2->{c_dirty} == 0) ?
+	       cplx(@{$w->cartesian}) : $w;
 }
 
 #
@@ -590,9 +602,11 @@ sub root {
 	my $theta_inc = pit2 / $n;
 	my $rho = $r ** (1/$n);
 	my $theta;
-	my $complex = ref($z) || $package;
+	my $cartesian = ref $z && $z->{c_dirty} == 0;
 	for ($k = 0, $theta = $t / $n; $k < $n; $k++, $theta += $theta_inc) {
-		push(@root, $complex->emake($rho, $theta));
+	    my $w = cplxe($rho, $theta);
+	    # Yes, $cartesian is loop invariant.
+	    push @root, $cartesian ? cplx(@{$w->cartesian}) : $w;
 	}
 	return @root;
 }
@@ -1232,11 +1246,42 @@ sub stringify_cartesian {
 	$str .= "+$im" if defined $im;
 	$str =~ s/\+-/-/;
 	$str =~ s/^\+//;
+	$str =~ s/([-+])1i/$1i/; # Not redundant with the above 1/-1 tests.
 	$str = '0' unless $str;
 
 	return $str;
 }
 
+
+# Helper for stringify_polar, a Greatest Common Divisor with a memory.
+
+sub _gcd {
+    my ($a, $b) = @_;
+
+    use integer;
+
+    # Loops forever if given negative inputs.
+
+    if    ($b and $a > $b) { return gcd($a % $b, $b) }
+    elsif ($a and $b > $a) { return gcd($b % $a, $a) }
+    else                   { return $a ? $a : $b     }
+}
+
+my %gcd;
+
+sub gcd {
+    my ($a, $b) = @_;
+
+    my $id = "$a $b";
+    
+    unless (exists $gcd{$id}) {
+	$gcd{$id} = _gcd($a, $b);
+	$gcd{"$b $a"} = $gcd{$id};
+    }
+
+    return $gcd{$id};
+}
+
 #
 # ->stringify_polar
 #
@@ -1270,15 +1315,26 @@ sub stringify_polar {
 	#
 
 	$nt -= pit2 if $nt > pi;
-	my ($n, $k, $kpi);
 
-	for ($k = 1, $kpi = pi; $k < 10; $k++, $kpi += pi) {
+	if (abs($nt) >= deg1) {
+	    my ($n, $k, $kpi);
+
+	    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) <= $eps) {
-			$theta = ($nt < 0 ? '-':'').
-				 ($k == 1 ? 'pi':"${k}pi").'/'.abs($n);
-			last;
+		    $n = abs $n;
+		    my $gcd = gcd($k, $n);
+		    if ($gcd > 1) {
+			$k /= $gcd;
+			$n /= $gcd;
+		    }
+		    next if $n > 360;
+		    $theta = ($nt < 0 ? '-':'').
+			     ($k == 1 ? 'pi':"${k}pi");
+		    $theta .= '/'.$n if $n > 1;
+		    last;
 		}
+	    }
 	}
 
 	$theta = $nt unless defined $theta;
@@ -1700,6 +1756,11 @@ 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.
 
+In Cray UNICOS there is some strange numerical instability that results
+in root(), cos(), sin(), cosh(), sinh(), losing accuracy fast.  Beware.
+The bug may be in UNICOS math libs, in UNICOS C compiler, in Math::Complex.
+Whatever it is, it does not manifest itself anywhere else where Perl runs.
+
 =head1 AUTHORS
 
 Raphael Manfredi <F<Raphael_Manfredi@grenoble.hp.com>> and
diff --git a/t/lib/complex.t b/t/lib/complex.t
index 2783e42..2bb14f0 100755
--- a/t/lib/complex.t
+++ b/t/lib/complex.t
@@ -26,6 +26,11 @@ my @script = (
 );
 my $eps = 1e-13;
 
+if ($^O eq 'unicos') { 	# For some reason root() produces very inaccurate
+    $eps = 1e-11;	# results in Cray UNICOS, and occasionally also
+}			# cos(), sin(), cosh(), sinh().  The division
+			# of doubles is the current suspect.
+
 while (<DATA>) {
 	s/^\s+//;
 	next if $_ eq '' || /^\#/;