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>

diff --git a/src/math/cmplx/asin.go b/src/math/cmplx/asin.go
index 61880a257d49d340616f6427f5509ec1ad3b6faf..062f324ce2cec423e0bdcf77fe4bde0d2786decd 100644 (file)
--- a/src/math/cmplx/asin.go
+++ b/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 {

diff --git a/src/math/cmplx/cmath_test.go b/src/math/cmplx/cmath_test.go
index 7a5c485a3103bd565d0cad3a2a6fb80edc9e7178..8d705622fd4368809505882ee3b326569b8d0980 100644 (file)
--- a/src/math/cmplx/cmath_test.go
+++ b/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++ {

diff --git a/src/math/cmplx/sqrt.go b/src/math/cmplx/sqrt.go
index 72f81e907cf0444e14931e52820bbe3c0302700c..0fbdcdedd36cdd38584fe9883a9bd26a50323b4e 100644 (file)
--- a/src/math/cmplx/sqrt.go
+++ b/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 {