complex divide: match C99 implementation

R=iant, ken2, r, r2, ken3
CC=golang-dev
https://golang.org/cl/1686044

Author: rsc

src/pkg/runtime/complex.c

@@ -4,33 +4,57 @@
 
 #include "runtime.h"
 
-// complex128div(num, den complex128) (quo complex128)
+typedef struct Complex128 Complex128;
+
 void
-·complex128div(float64 numreal, float64 numimag,
-	float64 denreal, float64 denimag,
-	float64 quoreal, float64 quoimag)
+·complex128div(Complex128 n, Complex128 d, Complex128 q)
 {
+	int32 ninf, dinf, nnan, dnan;
 	float64 a, b, ratio, denom;
 
-	a = denreal;
-	if(a < 0)
-		a = -a;
-	b = denimag;
-	if(b < 0)
-		b = -b;
-	if(a <= b) {
-		if(b == 0)
-			panicstring("complex divide by zero");
-		ratio = denreal/denimag;
-		denom = denreal*ratio + denimag;
-		quoreal = (numreal*ratio + numimag) / denom;
-		quoimag = (numimag*ratio - numreal) / denom;
+	// Special cases as in C99.
+	ninf = isInf(n.real, 0) || isInf(n.imag, 0);
+	dinf = isInf(d.real, 0) || isInf(d.imag, 0);
+
+	nnan = !ninf && (isNaN(n.real) || isNaN(n.imag));
+	dnan = !dinf && (isNaN(d.real) || isNaN(d.imag));
+
+	if(nnan || dnan) {
+		q.real = NaN();
+		q.imag = NaN();
+	} else if(ninf && !dinf && !dnan) {
+		q.real = Inf(0);
+		q.imag = Inf(0);
+	} else if(!ninf && !nnan && dinf) {
+		q.real = 0;
+		q.imag = 0;
+	} else if(d.real == 0 && d.imag == 0) {
+		if(n.real == 0 && n.imag == 0) {
+			q.real = NaN();
+			q.imag = NaN();
+		} else {
+			q.real = Inf(0);
+			q.imag = Inf(0);
+		}
 	} else {
-		ratio = denimag/denreal;
-		denom = denimag*ratio + denreal;
-		quoreal = (numimag*ratio + numreal) / denom;
-		quoimag = (numimag - numreal*ratio) / denom;
+		// Standard complex arithmetic, factored to avoid unnecessary overflow.
+		a = d.real;
+		if(a < 0)
+			a = -a;
+		b = d.imag;
+		if(b < 0)
+			b = -b;
+		if(a <= b) {
+			ratio = d.real/d.imag;
+			denom = d.real*ratio + d.imag;
+			q.real = (n.real*ratio + n.imag) / denom;
+			q.imag = (n.imag*ratio - n.real) / denom;
+		} else {
+			ratio = d.imag/d.real;
+			denom = d.imag*ratio + d.real;
+			q.real = (n.imag*ratio + n.real) / denom;
+			q.imag = (n.imag - n.real*ratio) / denom;
+		}
 	}
-	FLUSH(&quoreal);
-	FLUSH(&quoimag);
+	FLUSH(&q);
 }

test/cmplxdivide.c

@@ -0,0 +1,61 @@
+// Copyright 2010 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.
+
+// gcc '-std=c99' cmplxdivide.c && a.out >cmplxdivide1.go
+
+#include <complex.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+
+#define nelem(x) (sizeof(x)/sizeof((x)[0]))
+
+double f[] = {
+	0,
+	1,
+	-1,
+	2,
+	NAN,
+	INFINITY,
+	-INFINITY,
+};
+
+char*
+fmt(double g)
+{
+	static char buf[10][30];
+	static int n;
+	char *p;
+	
+	p = buf[n++];
+	if(n == 10)
+		n = 0;
+	sprintf(p, "%g", g);
+	if(strcmp(p, "-0") == 0)
+		strcpy(p, "negzero");
+	return p;
+}	
+
+int
+main(void)
+{
+	int i, j, k, l;
+	double complex n, d, q;
+	
+	printf("// # generated by cmplxdivide.c\n");
+	printf("\n");
+	printf("package main\n");
+	printf("var tests = []Test{\n");
+	for(i=0; i<nelem(f); i++)
+	for(j=0; j<nelem(f); j++)
+	for(k=0; k<nelem(f); k++)
+	for(l=0; l<nelem(f); l++) {
+		n = f[i] + f[j]*I;
+		d = f[k] + f[l]*I;
+		q = n/d;
+		printf("\tTest{cmplx(%s, %s), cmplx(%s, %s), cmplx(%s, %s)},\n", fmt(creal(n)), fmt(cimag(n)), fmt(creal(d)), fmt(cimag(d)), fmt(creal(q)), fmt(cimag(q)));
+	}
+	printf("}\n");
+	return 0;
+}

test/cmplxdivide.go

@@ -0,0 +1,43 @@
+// $G $D/$F.go $D/cmplxdivide1.go && $L $D/$F.$A && ./$A.out
+
+// Copyright 2010 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.
+
+// Driver for complex division table defined in cmplxdivide1.go
+
+package main
+
+import (
+	"cmath"
+	"fmt"
+	"math"
+)
+
+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 cmath.IsInf(a) && cmath.IsInf(b):
+		return true
+	case cmath.IsNaN(a) && cmath.IsNaN(b):
+		return true
+	}
+	return a == b
+}
+
+func main() {
+	for _, t := range tests {
+		x := t.f/t.g
+		if !calike(x, t.out) {
+			fmt.Printf("%v/%v: expected %v error; got %v\n", t.f, t.g, t.out, x)
+		}
+	}
+}