math/cmplx: use signed zero to correct branch cuts

Branch cuts for the elementary complex functions along real or imaginary axes
should be resolved in floating point calculations by one-sided continuity with
signed zero as described in:

"Branch Cuts for Complex Elementary Functions or Much Ado About Nothing's Sign Bit"
W. Kahan

Available at: https://people.freebsd.org/~das/kahan86branch.pdf

And as described in the C99 standard which is claimed as the original cephes source.

Sqrt did not return the correct branch when imag(x) == 0. The branch is now
determined by sign(imag(x)).  This incorrect branch choice was affecting the behavior
of the Trigonometric/Hyperbolic functions that use Sqrt in intermediate calculations.

Asin, Asinh and Atan had spurious domain checks, whereas the functions should be valid
over the whole complex plane with appropriate branch cuts.

Fixes #6888

Change-Id: I9b1278af54f54bfb4208276ae345bbd3ddf3ec83
Reviewed-on: https://go-review.googlesource.com/46492
Run-TryBot: Brad Fitzpatrick <bradfitz@golang.org>
TryBot-Result: Gobot Gobot <gobot@golang.org>
Reviewed-by: Russ Cox <rsc@golang.org>

src/math/cmplx/asin.go

@@ -49,11 +49,8 @@ import "math"
 
 // Asin returns the inverse sine of x.
 func Asin(x complex128) complex128 {
-	if imag(x) == 0 {
-		if math.Abs(real(x)) > 1 {
-			return complex(math.Pi/2, 0) // DOMAIN error
-		}
-		return complex(math.Asin(real(x)), 0)
+	if imag(x) == 0 && math.Abs(real(x)) <= 1 {
+		return complex(math.Asin(real(x)), imag(x))
 	}
 	ct := complex(-imag(x), real(x)) // i * x
 	xx := x * x
@@ -65,12 +62,8 @@ func Asin(x complex128) complex128 {
 
 // Asinh returns the inverse hyperbolic sine of x.
 func Asinh(x complex128) complex128 {
-	// TODO check range
-	if imag(x) == 0 {
-		if math.Abs(real(x)) > 1 {
-			return complex(math.Pi/2, 0) // DOMAIN error
-		}
-		return complex(math.Asinh(real(x)), 0)
+	if imag(x) == 0 && math.Abs(real(x)) <= 1 {
+		return complex(math.Asinh(real(x)), imag(x))
 	}
 	xx := x * x
 	x1 := complex(1+real(xx), imag(xx)) // 1 + x*x
@@ -140,10 +133,6 @@ func Acosh(x complex128) complex128 {
 
 // Atan returns the inverse tangent of x.
 func Atan(x complex128) complex128 {
-	if real(x) == 0 && imag(x) > 1 {
-		return NaN()
-	}
-
 	x2 := real(x) * real(x)
 	a := 1 - x2 - imag(x)*imag(x)
 	if a == 0 {

src/math/cmplx/cmath_test.go

@@ -435,6 +435,24 @@ var tanhSC = []complex128{
 	NaN(),
 }
 
+// branch cut continuity checks
+// points on each axis at |z| > 1 are checked for one-sided continuity from both the positive and negative side
+// all possible branch cuts for the elementary functions are at one of these points
+
+var zero = 0.0
+var eps = 1.0 / (1 << 53)
+
+var branchPoints = [][2]complex128{
+	{complex(2.0, zero), complex(2.0, eps)},
+	{complex(2.0, -zero), complex(2.0, -eps)},
+	{complex(-2.0, zero), complex(-2.0, eps)},
+	{complex(-2.0, -zero), complex(-2.0, -eps)},
+	{complex(zero, 2.0), complex(eps, 2.0)},
+	{complex(-zero, 2.0), complex(-eps, 2.0)},
+	{complex(zero, -2.0), complex(eps, -2.0)},
+	{complex(-zero, -2.0), complex(-eps, -2.0)},
+}
+
 // functions borrowed from pkg/math/all_test.go
 func tolerance(a, b, e float64) bool {
 	d := a - b
@@ -508,6 +526,11 @@ func TestAcos(t *testing.T) {
 			t.Errorf("Acos(%g) = %g, want %g", vcAcosSC[i], f, acosSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Acos(pt[0]), Acos(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Acos(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestAcosh(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -520,6 +543,11 @@ func TestAcosh(t *testing.T) {
 			t.Errorf("Acosh(%g) = %g, want %g", vcAcoshSC[i], f, acoshSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Acosh(pt[0]), Acosh(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Acosh(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestAsin(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -532,6 +560,11 @@ func TestAsin(t *testing.T) {
 			t.Errorf("Asin(%g) = %g, want %g", vcAsinSC[i], f, asinSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Asin(pt[0]), Asin(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Asin(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestAsinh(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -544,6 +577,11 @@ func TestAsinh(t *testing.T) {
 			t.Errorf("Asinh(%g) = %g, want %g", vcAsinhSC[i], f, asinhSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Asinh(pt[0]), Asinh(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Asinh(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestAtan(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -556,6 +594,11 @@ func TestAtan(t *testing.T) {
 			t.Errorf("Atan(%g) = %g, want %g", vcAtanSC[i], f, atanSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Atan(pt[0]), Atan(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Atan(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestAtanh(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -568,6 +611,11 @@ func TestAtanh(t *testing.T) {
 			t.Errorf("Atanh(%g) = %g, want %g", vcAtanhSC[i], f, atanhSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Atanh(pt[0]), Atanh(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Atanh(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestConj(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -635,6 +683,11 @@ func TestLog(t *testing.T) {
 			t.Errorf("Log(%g) = %g, want %g", vcLogSC[i], f, logSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Log(pt[0]), Log(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Log(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestLog10(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -685,6 +738,11 @@ func TestPow(t *testing.T) {
 			t.Errorf("Pow(%g, %g) = %g, want %g", vcPowSC[i][0], vcPowSC[i][0], f, powSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Pow(pt[0], 0.1), Pow(pt[1], 0.1); !cVeryclose(f0, f1) {
+			t.Errorf("Pow(%g, 0.1) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestRect(t *testing.T) {
 	for i := 0; i < len(vc); i++ {
@@ -733,6 +791,11 @@ func TestSqrt(t *testing.T) {
 			t.Errorf("Sqrt(%g) = %g, want %g", vcSqrtSC[i], f, sqrtSC[i])
 		}
 	}
+	for _, pt := range branchPoints {
+		if f0, f1 := Sqrt(pt[0]), Sqrt(pt[1]); !cVeryclose(f0, f1) {
+			t.Errorf("Sqrt(%g) not continuous, got %g want %g", pt[0], f0, f1)
+		}
+	}
 }
 func TestTan(t *testing.T) {
 	for i := 0; i < len(vc); i++ {

src/math/cmplx/sqrt.go

@@ -57,13 +57,14 @@ import "math"
 // The result r is chosen so that real(r) ≥ 0 and imag(r) has the same sign as imag(x).
 func Sqrt(x complex128) complex128 {
 	if imag(x) == 0 {
+		// Ensure that imag(r) has the same sign as imag(x) for imag(x) == signed zero.
 		if real(x) == 0 {
-			return complex(0, 0)
+			return complex(0, imag(x))
 		}
 		if real(x) < 0 {
-			return complex(0, math.Sqrt(-real(x)))
+			return complex(0, math.Copysign(math.Sqrt(-real(x)), imag(x)))
 		}
-		return complex(math.Sqrt(real(x)), 0)
+		return complex(math.Sqrt(real(x)), imag(x))
 	}
 	if real(x) == 0 {
 		if imag(x) < 0 {