runtime: make complex division c99 compatible

- changes tests to check that the real and imaginary part of the go complex
  division result is equal to the result gcc produces for c99
- changes complex division code to satisfy new complex division test
- adds float functions isNan, isFinite, isInf, abs and copysign
  in the runtime package

Fixes #14644.

name                   old time/op  new time/op  delta
Complex128DivNormal-4  21.8ns ± 6%  13.9ns ± 6%  -36.37%  (p=0.000 n=20+20)
Complex128DivNisNaN-4  14.1ns ± 1%  15.0ns ± 1%   +5.86%  (p=0.000 n=20+19)
Complex128DivDisNaN-4  12.5ns ± 1%  16.7ns ± 1%  +33.79%  (p=0.000 n=19+20)
Complex128DivNisInf-4  10.1ns ± 1%  13.0ns ± 1%  +28.25%  (p=0.000 n=20+19)
Complex128DivDisInf-4  11.0ns ± 1%  20.9ns ± 1%  +90.69%  (p=0.000 n=16+19)
ComplexAlgMap-4        86.7ns ± 1%  86.8ns ± 2%     ~     (p=0.804 n=20+20)

Change-Id: I261f3b4a81f6cc858bc7ff48f6fd1b39c300abf0
Reviewed-on: https://go-review.googlesource.com/37441
Reviewed-by: Robert Griesemer <gri@golang.org>

1 files changed, 49 insertions, 59 deletions

diff --git a/src/runtime/complex.go b/src/runtime/complex.go
index 73f1161a50..07c596fc0b 100644
--- a/src/runtime/complex.go
+++ b/src/runtime/complex.go
@@ -4,68 +4,58 @@
 
 package runtime
 
-func isposinf(f float64) bool { return f > maxFloat64 }
-func isneginf(f float64) bool { return f < -maxFloat64 }
-func isnan(f float64) bool    { return f != f }
-
-func nan() float64 {
-	var f float64 = 0
-	return f / f
-}
-
-func posinf() float64 {
-	var f float64 = maxFloat64
-	return f * f
-}
-
-func neginf() float64 {
-	var f float64 = maxFloat64
-	return -f * f
+// inf2one returns a signed 1 if f is an infinity and a signed 0 otherwise.
+// The sign of the result is the sign of f.
+func inf2one(f float64) float64 {
+	g := 0.0
+	if isInf(f) {
+		g = 1.0
+	}
+	return copysign(g, f)
 }
 
-func complex128div(n complex128, d complex128) complex128 {
-	// Special cases as in C99.
-	ninf := isposinf(real(n)) || isneginf(real(n)) ||
-		isposinf(imag(n)) || isneginf(imag(n))
-	dinf := isposinf(real(d)) || isneginf(real(d)) ||
-		isposinf(imag(d)) || isneginf(imag(d))
-
-	nnan := !ninf && (isnan(real(n)) || isnan(imag(n)))
-	dnan := !dinf && (isnan(real(d)) || isnan(imag(d)))
+func complex128div(n complex128, m complex128) complex128 {
+	var e, f float64 // complex(e, f) = n/m
+
+	// Algorithm for robust complex division as described in
+	// Robert L. Smith: Algorithm 116: Complex division. Commun. ACM 5(8): 435 (1962).
+	if abs(real(m)) >= abs(imag(m)) {
+		ratio := imag(m) / real(m)
+		denom := real(m) + ratio*imag(m)
+		e = (real(n) + imag(n)*ratio) / denom
+		f = (imag(n) - real(n)*ratio) / denom
+	} else {
+		ratio := real(m) / imag(m)
+		denom := imag(m) + ratio*real(m)
+		e = (real(n)*ratio + imag(n)) / denom
+		f = (imag(n)*ratio - real(n)) / denom
+	}
 
-	switch {
-	case nnan || dnan:
-		return complex(nan(), nan())
-	case ninf && !dinf:
-		return complex(posinf(), posinf())
-	case !ninf && dinf:
-		return complex(0, 0)
-	case real(d) == 0 && imag(d) == 0:
-		if real(n) == 0 && imag(n) == 0 {
-			return complex(nan(), nan())
-		} else {
-			return complex(posinf(), posinf())
-		}
-	default:
-		// Standard complex arithmetic, factored to avoid unnecessary overflow.
-		a := real(d)
-		if a < 0 {
-			a = -a
-		}
-		b := imag(d)
-		if b < 0 {
-			b = -b
-		}
-		if a <= b {
-			ratio := real(d) / imag(d)
-			denom := real(d)*ratio + imag(d)
-			return complex((real(n)*ratio+imag(n))/denom,
-				(imag(n)*ratio-real(n))/denom)
-		} else {
-			ratio := imag(d) / real(d)
-			denom := imag(d)*ratio + real(d)
-			return complex((imag(n)*ratio+real(n))/denom,
-				(imag(n)-real(n)*ratio)/denom)
+	if isNaN(e) && isNaN(f) {
+		// Correct final result to infinities and zeros if applicable.
+		// Matches C99: ISO/IEC 9899:1999 - G.5.1  Multiplicative operators.
+
+		a, b := real(n), imag(n)
+		c, d := real(m), imag(m)
+
+		switch {
+		case m == 0 && (!isNaN(a) || !isNaN(b)):
+			e = copysign(inf, c) * a
+			f = copysign(inf, c) * b
+
+		case (isInf(a) || isInf(b)) && isFinite(c) && isFinite(d):
+			a = inf2one(a)
+			b = inf2one(b)
+			e = inf * (a*c + b*d)
+			f = inf * (b*c - a*d)
+
+		case (isInf(c) || isInf(d)) && isFinite(a) && isFinite(b):
+			c = inf2one(c)
+			d = inf2one(d)
+			e = 0 * (a*c + b*d)
+			f = 0 * (b*c - a*d)
 		}
 	}
+
+	return complex(e, f)
 }