math/cmplx: handle special cases

Implement special case handling and testing to ensure
conformance with the C99 standard annex G.6 Complex arithmetic.

Fixes #29320

Change-Id: Ieb0527191dd7fdea5b1aecb42b9e23aae3f74260
Reviewed-on: https://go-review.googlesource.com/c/go/+/169501
Run-TryBot: Brian Kessler <brian.m.kessler@gmail.com>
TryBot-Result: Gobot Gobot <gobot@golang.org>
Reviewed-by: Robert Griesemer <gri@golang.org>

src/math/cmplx/abs.go

@@ -3,7 +3,8 @@
 // license that can be found in the LICENSE file.
 
 // Package cmplx provides basic constants and mathematical functions for
-// complex numbers.
+// complex numbers. Special case handling conforms to the C99 standard
+// Annex G IEC 60559-compatible complex arithmetic.
 package cmplx
 
 import "math"

src/math/cmplx/asin.go

@@ -49,8 +49,31 @@ import "math"
 
 // Asin returns the inverse sine of x.
 func Asin(x complex128) complex128 {
-	if imag(x) == 0 && math.Abs(real(x)) <= 1 {
-		return complex(math.Asin(real(x)), imag(x))
+	switch re, im := real(x), imag(x); {
+	case im == 0 && math.Abs(re) <= 1:
+		return complex(math.Asin(re), im)
+	case re == 0 && math.Abs(im) <= 1:
+		return complex(re, math.Asinh(im))
+	case math.IsNaN(im):
+		switch {
+		case re == 0:
+			return complex(re, math.NaN())
+		case math.IsInf(re, 0):
+			return complex(math.NaN(), re)
+		default:
+			return NaN()
+		}
+	case math.IsInf(im, 0):
+		switch {
+		case math.IsNaN(re):
+			return x
+		case math.IsInf(re, 0):
+			return complex(math.Copysign(math.Pi/4, re), im)
+		default:
+			return complex(math.Copysign(0, re), im)
+		}
+	case math.IsInf(re, 0):
+		return complex(math.Copysign(math.Pi/2, re), math.Copysign(re, im))
 	}
 	ct := complex(-imag(x), real(x)) // i * x
 	xx := x * x
@@ -62,8 +85,31 @@ func Asin(x complex128) complex128 {
 
 // Asinh returns the inverse hyperbolic sine of x.
 func Asinh(x complex128) complex128 {
-	if imag(x) == 0 && math.Abs(real(x)) <= 1 {
-		return complex(math.Asinh(real(x)), imag(x))
+	switch re, im := real(x), imag(x); {
+	case im == 0 && math.Abs(re) <= 1:
+		return complex(math.Asinh(re), im)
+	case re == 0 && math.Abs(im) <= 1:
+		return complex(re, math.Asin(im))
+	case math.IsInf(re, 0):
+		switch {
+		case math.IsInf(im, 0):
+			return complex(re, math.Copysign(math.Pi/4, im))
+		case math.IsNaN(im):
+			return x
+		default:
+			return complex(re, math.Copysign(0.0, im))
+		}
+	case math.IsNaN(re):
+		switch {
+		case im == 0:
+			return x
+		case math.IsInf(im, 0):
+			return complex(im, re)
+		default:
+			return NaN()
+		}
+	case math.IsInf(im, 0):
+		return complex(math.Copysign(im, re), math.Copysign(math.Pi/2, im))
 	}
 	xx := x * x
 	x1 := complex(1+real(xx), imag(xx)) // 1 + x*x
@@ -91,6 +137,9 @@ func Acos(x complex128) complex128 {
 
 // Acosh returns the inverse hyperbolic cosine of x.
 func Acosh(x complex128) complex128 {
+	if x == 0 {
+		return complex(0, math.Copysign(math.Pi/2, imag(x)))
+	}
 	w := Acos(x)
 	if imag(w) <= 0 {
 		return complex(-imag(w), real(w)) // i * w
@@ -133,6 +182,17 @@ func Acosh(x complex128) complex128 {
 
 // Atan returns the inverse tangent of x.
 func Atan(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case im == 0:
+		return complex(math.Atan(re), im)
+	case re == 0 && math.Abs(im) <= 1:
+		return complex(re, math.Atanh(im))
+	case math.IsInf(im, 0) || math.IsInf(re, 0):
+		if math.IsNaN(re) {
+			return complex(math.NaN(), math.Copysign(0, im))
+		}
+		return complex(math.Copysign(math.Pi/2, re), math.Copysign(0, im))
+	}
 	x2 := real(x) * real(x)
 	a := 1 - x2 - imag(x)*imag(x)
 	if a == 0 {

src/math/cmplx/exp.go

@@ -49,6 +49,23 @@ import "math"
 
 // Exp returns e**x, the base-e exponential of x.
 func Exp(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case math.IsInf(re, 0):
+		switch {
+		case re > 0 && im == 0:
+			return x
+		case math.IsInf(im, 0) || math.IsNaN(im):
+			if re < 0 {
+				return complex(0, math.Copysign(0, im))
+			} else {
+				return complex(math.Inf(1.0), math.NaN())
+			}
+		}
+	case math.IsNaN(re):
+		if im == 0 {
+			return complex(math.NaN(), im)
+		}
+	}
 	r := math.Exp(real(x))
 	s, c := math.Sincos(imag(x))
 	return complex(r*c, r*s)

src/math/cmplx/log.go

@@ -60,5 +60,6 @@ func Log(x complex128) complex128 {
 
 // Log10 returns the decimal logarithm of x.
 func Log10(x complex128) complex128 {
-	return math.Log10E * Log(x)
+	z := Log(x)
+	return complex(math.Log10E*real(z), math.Log10E*imag(z))
 }

src/math/cmplx/sin.go

@@ -51,6 +51,19 @@ import "math"
 
 // Sin returns the sine of x.
 func Sin(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case im == 0 && (math.IsInf(re, 0) || math.IsNaN(re)):
+		return complex(math.NaN(), im)
+	case math.IsInf(im, 0):
+		switch {
+		case re == 0:
+			return x
+		case math.IsInf(re, 0) || math.IsNaN(re):
+			return complex(math.NaN(), im)
+		}
+	case re == 0 && math.IsNaN(im):
+		return x
+	}
 	s, c := math.Sincos(real(x))
 	sh, ch := sinhcosh(imag(x))
 	return complex(s*ch, c*sh)
@@ -71,6 +84,19 @@ func Sin(x complex128) complex128 {
 
 // Sinh returns the hyperbolic sine of x.
 func Sinh(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case re == 0 && (math.IsInf(im, 0) || math.IsNaN(im)):
+		return complex(re, math.NaN())
+	case math.IsInf(re, 0):
+		switch {
+		case im == 0:
+			return complex(re, im)
+		case math.IsInf(im, 0) || math.IsNaN(im):
+			return complex(re, math.NaN())
+		}
+	case im == 0 && math.IsNaN(re):
+		return complex(math.NaN(), im)
+	}
 	s, c := math.Sincos(imag(x))
 	sh, ch := sinhcosh(real(x))
 	return complex(c*sh, s*ch)
@@ -96,6 +122,19 @@ func Sinh(x complex128) complex128 {
 
 // Cos returns the cosine of x.
 func Cos(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case im == 0 && (math.IsInf(re, 0) || math.IsNaN(re)):
+		return complex(math.NaN(), -im*math.Copysign(0, re))
+	case math.IsInf(im, 0):
+		switch {
+		case re == 0:
+			return complex(math.Inf(1), -re*math.Copysign(0, im))
+		case math.IsInf(re, 0) || math.IsNaN(re):
+			return complex(math.Inf(1), math.NaN())
+		}
+	case re == 0 && math.IsNaN(im):
+		return complex(math.NaN(), 0)
+	}
 	s, c := math.Sincos(real(x))
 	sh, ch := sinhcosh(imag(x))
 	return complex(c*ch, -s*sh)
@@ -115,6 +154,19 @@ func Cos(x complex128) complex128 {
 
 // Cosh returns the hyperbolic cosine of x.
 func Cosh(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case re == 0 && (math.IsInf(im, 0) || math.IsNaN(im)):
+		return complex(math.NaN(), re*math.Copysign(0, im))
+	case math.IsInf(re, 0):
+		switch {
+		case im == 0:
+			return complex(math.Inf(1), im*math.Copysign(0, re))
+		case math.IsInf(im, 0) || math.IsNaN(im):
+			return complex(math.Inf(1), math.NaN())
+		}
+	case im == 0 && math.IsNaN(re):
+		return complex(math.NaN(), im)
+	}
 	s, c := math.Sincos(imag(x))
 	sh, ch := sinhcosh(real(x))
 	return complex(c*ch, s*sh)

src/math/cmplx/sqrt.go

@@ -65,6 +65,8 @@ func Sqrt(x complex128) complex128 {
 			return complex(0, math.Copysign(math.Sqrt(-real(x)), imag(x)))
 		}
 		return complex(math.Sqrt(real(x)), imag(x))
+	} else if math.IsInf(imag(x), 0) {
+		return complex(math.Inf(1.0), imag(x))
 	}
 	if real(x) == 0 {
 		if imag(x) < 0 {

src/math/cmplx/tan.go

@@ -57,6 +57,16 @@ import "math"
 
 // Tan returns the tangent of x.
 func Tan(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case math.IsInf(im, 0):
+		switch {
+		case math.IsInf(re, 0) || math.IsNaN(re):
+			return complex(math.Copysign(0, re), math.Copysign(1, im))
+		}
+		return complex(math.Copysign(0, math.Sin(2*re)), math.Copysign(1, im))
+	case re == 0 && math.IsNaN(im):
+		return x
+	}
 	d := math.Cos(2*real(x)) + math.Cosh(2*imag(x))
 	if math.Abs(d) < 0.25 {
 		d = tanSeries(x)
@@ -81,6 +91,16 @@ func Tan(x complex128) complex128 {
 
 // Tanh returns the hyperbolic tangent of x.
 func Tanh(x complex128) complex128 {
+	switch re, im := real(x), imag(x); {
+	case math.IsInf(re, 0):
+		switch {
+		case math.IsInf(im, 0) || math.IsNaN(im):
+			return complex(math.Copysign(1, re), math.Copysign(0, im))
+		}
+		return complex(math.Copysign(1, re), math.Copysign(0, math.Sin(2*im)))
+	case im == 0 && math.IsNaN(re):
+		return x
+	}
 	d := math.Cosh(2*real(x)) + math.Cos(2*imag(x))
 	if d == 0 {
 		return Inf()