From 21ff75bc0efc34675775168ccae118cd286f904c Mon Sep 17 00:00:00 2001 From: Russ Cox Date: Fri, 18 Jun 2010 15:46:00 -0700 Subject: [PATCH] complex divide: match C99 implementation R=iant, ken2, r, r2, ken3 CC=golang-dev https://golang.org/cl/1686044 --- src/pkg/runtime/complex.c | 70 +- test/cmplxdivide.c | 61 + test/cmplxdivide.go | 43 + test/cmplxdivide1.go | 2406 +++++++++++++++++++++++++++++++++++++ test/golden.out | 8 - test/recover2.go | 8 +- test/zerodivide.go | 14 +- 7 files changed, 2569 insertions(+), 41 deletions(-) create mode 100644 test/cmplxdivide.c create mode 100644 test/cmplxdivide.go create mode 100644 test/cmplxdivide1.go diff --git a/src/pkg/runtime/complex.c b/src/pkg/runtime/complex.c index ca6ed79ba3..2240d9fb85 100644 --- a/src/pkg/runtime/complex.c +++ b/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); } diff --git a/test/cmplxdivide.c b/test/cmplxdivide.c new file mode 100644 index 0000000000..63473ba6cb --- /dev/null +++ b/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 +#include +#include +#include + +#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