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>

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
+// 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 posinf() float64 {
-	var f float64 = maxFloat64
-	return f * f
-}
+func complex128div(n complex128, m complex128) complex128 {
+	var e, f float64 // complex(e, f) = n/m
 
-func neginf() float64 {
-	var f float64 = maxFloat64
-	return -f * f
-}
+	// 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
+	}
 
-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))
+	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.
 
-	nnan := !ninf && (isnan(real(n)) || isnan(imag(n)))
-	dnan := !dinf && (isnan(real(d)) || isnan(imag(d)))
+		a, b := real(n), imag(n)
+		c, d := real(m), imag(m)
 
-	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)
+		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)
 }

src/runtime/fastlog2.go

@@ -4,8 +4,6 @@
 
 package runtime
 
-import "unsafe"
-
 // fastlog2 implements a fast approximation to the base 2 log of a
 // float64. This is used to compute a geometric distribution for heap
 // sampling, without introducing dependencies into package math. This
@@ -27,7 +25,3 @@ func fastlog2(x float64) float64 {
 	low, high := fastlog2Table[xManIndex], fastlog2Table[xManIndex+1]
 	return float64(xExp) + low + (high-low)*float64(xManScale)*fastlogScaleRatio
 }
-
-// float64bits returns the IEEE 754 binary representation of f.
-// Taken from math.Float64bits to avoid dependencies into package math.
-func float64bits(f float64) uint64 { return *(*uint64)(unsafe.Pointer(&f)) }

src/runtime/float.go

@@ -0,0 +1,53 @@
+// Copyright 2017 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package runtime
+
+import "unsafe"
+
+var inf = float64frombits(0x7FF0000000000000)
+
+// isNaN reports whether f is an IEEE 754 ``not-a-number'' value.
+func isNaN(f float64) (is bool) {
+	// IEEE 754 says that only NaNs satisfy f != f.
+	return f != f
+}
+
+// isFinite reports whether f is neither NaN nor an infinity.
+func isFinite(f float64) bool {
+	return !isNaN(f - f)
+}
+
+// isInf reports whether f is an infinity.
+func isInf(f float64) bool {
+	return !isNaN(f) && !isFinite(f)
+}
+
+// Abs returns the absolute value of x.
+//
+// Special cases are:
+//	Abs(±Inf) = +Inf
+//	Abs(NaN) = NaN
+func abs(x float64) float64 {
+	const sign = 1 << 63
+	return float64frombits(float64bits(x) &^ sign)
+}
+
+// copysign returns a value with the magnitude
+// of x and the sign of y.
+func copysign(x, y float64) float64 {
+	const sign = 1 << 63
+	return float64frombits(float64bits(x)&^sign | float64bits(y)&sign)
+}
+
+// Float64bits returns the IEEE 754 binary representation of f.
+func float64bits(f float64) uint64 {
+	return *(*uint64)(unsafe.Pointer(&f))
+}
+
+// Float64frombits returns the floating point number corresponding
+// the IEEE 754 binary representation b.
+func float64frombits(b uint64) float64 {
+	return *(*float64)(unsafe.Pointer(&b))
+}

test/cmplxdivide.c

@@ -23,50 +23,63 @@
 #define nelem(x) (sizeof(x)/sizeof((x)[0]))
 
 double f[] = {
-	0,
-	1,
-	-1,
-	2,
+	0.0,
+	-0.0,
+	1.0,
+	-1.0,
+	2.0,
 	NAN,
 	INFINITY,
 	-INFINITY,
 };
 
-char*
-fmt(double g)
-{
+char* fmt(double g) {
 	static char buf[10][30];
 	static int n;
 	char *p;
-	
+
 	p = buf[n++];
-	if(n == 10)
+	if(n == 10) {
 		n = 0;
+	}
+
 	sprintf(p, "%g", g);
-	if(strcmp(p, "-0") == 0)
-		strcpy(p, "negzero");
+
+	if(strcmp(p, "0") == 0) {
+		strcpy(p, "zero");
+		return p;
+	}
+
+	if(strcmp(p, "-0") == 0) {
+		strcpy(p, "-zero");
+		return p;
+	}
+
 	return p;
 }
 
-int
-iscnan(double complex d)
-{
-	return !isinf(creal(d)) && !isinf(cimag(d)) && (isnan(creal(d)) || isnan(cimag(d)));
-}
-
-double complex zero;	// attempt to hide zero division from gcc
-
-int
-main(void)
-{
+int main(void) {
 	int i, j, k, l;
 	double complex n, d, q;
-	
+
 	printf("// skip\n");
 	printf("// # generated by cmplxdivide.c\n");
 	printf("\n");
 	printf("package main\n");
-	printf("var tests = []Test{\n");
+	printf("\n");
+	printf("import \"math\"\n");
+	printf("\n");
+	printf("var (\n");
+	printf("\tnan     = math.NaN()\n");
+	printf("\tinf     = math.Inf(1)\n");
+	printf("\tzero    = 0.0\n");
+	printf(")\n");
+	printf("\n");
+	printf("var tests = []struct {\n");
+	printf("\tf, g complex128\n");
+	printf("\tout  complex128\n");
+	printf("}{\n");
+
 	for(i=0; i<nelem(f); i++)
 	for(j=0; j<nelem(f); j++)
 	for(k=0; k<nelem(f); k++)
@@ -74,17 +87,8 @@ main(void)
 		n = f[i] + f[j]*I;
 		d = f[k] + f[l]*I;
 		q = n/d;
-		
-		// BUG FIX.
-		// Gcc gets the wrong answer for NaN/0 unless both sides are NaN.
-		// That is, it treats (NaN+NaN*I)/0 = NaN+NaN*I (a complex NaN)
-		// but it then computes (1+NaN*I)/0 = Inf+NaN*I (a complex infinity).
-		// Since both numerators are complex NaNs, it seems that the
-		// results should agree in kind.  Override the gcc computation in this case.
-		if(iscnan(n) && d == 0)
-			q = (NAN+NAN*I) / zero;
 
-		printf("\tTest{complex(%s, %s), complex(%s, %s), complex(%s, %s)},\n",
+		printf("\t{complex(%s, %s), complex(%s, %s), complex(%s, %s)},\n",
 			fmt(creal(n)), fmt(cimag(n)),
 			fmt(creal(d)), fmt(cimag(d)),
 			fmt(creal(q)), fmt(cimag(q)));

test/cmplxdivide.go

@@ -5,33 +5,25 @@
 // license that can be found in the LICENSE file.
 
 // Driver for complex division table defined in cmplxdivide1.go
-// For details, see the comment at the top of in cmplxdivide.c.
+// For details, see the comment at the top of cmplxdivide.c.
 
 package main
 
 import (
 	"fmt"
 	"math"
-	"math/cmplx"
 )
 
-type Test struct {
-	f, g complex128
-	out  complex128
-}
-
-var nan = math.NaN()
-var inf = math.Inf(1)
-var negzero = math.Copysign(0, -1)
-
 func calike(a, b complex128) bool {
-	switch {
-	case cmplx.IsInf(a) && cmplx.IsInf(b):
-		return true
-	case cmplx.IsNaN(a) && cmplx.IsNaN(b):
-		return true
+	if imag(a) != imag(b) && !(math.IsNaN(imag(a)) && math.IsNaN(imag(b))) {
+		return false
 	}
-	return a == b
+
+	if real(a) != real(b) && !(math.IsNaN(real(a)) && math.IsNaN(real(b))) {
+		return false
+	}
+
+	return true
 }
 
 func main() {