From: Russ Cox Date: Tue, 4 Oct 2016 02:13:57 +0000 (-0400) Subject: math: fix Gamma(x) for x < -170.5 and other corner cases X-Git-Tag: go1.8beta1~1033 X-Git-Url: http://www.git.cypherpunks.su/?a=commitdiff_plain;h=a39920fdbbeca653a8b2ac52678378d260e2d396;p=gostls13.git math: fix Gamma(x) for x < -170.5 and other corner cases Fixes #11441. Test tables generated by package main import ( "bytes" "fmt" "log" "os/exec" "strconv" "strings" ) var inputs = []float64{ 0.5, 1.5, 2.5, 3.5, -0.5, -1.5, -2.5, -3.5, 0.1, 0.01, 1e-8, 1e-16, 1e-3, 1e-16, 1e-308, 5.6e-309, 5.5e-309, 1e-309, 1e-323, 5e-324, -0.1, -0.01, -1e-8, -1e-16, -1e-3, -1e-16, -1e-308, -5.6e-309, -5.5e-309, -1e-300 / 1e9, -1e-300 / 1e23, -5e-300 / 1e24, -0.9999999999999999, -1.0000000000000002, -1.9999999999999998, -2.0000000000000004, -100.00000000000001, -99.999999999999986, 17, 171, 171.6, 171.624, 171.625, 172, 2000, -100.5, -160.5, -170.5, -171.5, -176.5, -177.5, -178.5, -179.5, -201.0001, -202.9999, -1000.5, -1000000000.3, -4503599627370495.5, -63.349078729022985, -127.45117632943295, } func main() { var buf bytes.Buffer for _, v := range inputs { fmt.Fprintf(&buf, "gamma(%.1000g)\n", v) } cmd := exec.Command("gp", "-q") cmd.Stdin = &buf out, err := cmd.CombinedOutput() if err != nil { log.Fatalf("gp: %v", err) } f := strings.Split(string(out), "\n") if len(f) > 0 && f[len(f)-1] == "" { f = f[:len(f)-1] } if len(f) != len(inputs) { log.Fatalf("gp: wrong output count\n%s\n", out) } for i, g := range f { gf, err := strconv.ParseFloat(strings.Replace(g, " E", "e", -1), 64) if err != nil { if strings.Contains(err.Error(), "value out of range") { if strings.HasPrefix(g, "-") { fmt.Printf("\t{%g, Inf(-1)},\n", inputs[i]) } else { fmt.Printf("\t{%g, Inf(1)},\n", inputs[i]) } continue } log.Fatal(err) } if gf == 0 && strings.HasPrefix(g, "-") { fmt.Printf("\t{%g, Copysign(0, -1)},\n", inputs[i]) continue } fmt.Printf("\t{%g, %g},\n", inputs[i], gf) } } Change-Id: Ie98c7751d92b8ffb40e8313f5ea10df0890e2feb Reviewed-on: https://go-review.googlesource.com/30146 Run-TryBot: Russ Cox Reviewed-by: Quentin Smith --- diff --git a/src/math/all_test.go b/src/math/all_test.go index 882400527b..c77dcfea6e 100644 --- a/src/math/all_test.go +++ b/src/math/all_test.go @@ -1165,21 +1165,88 @@ var frexpSC = []fi{ {NaN(), 0}, } -var vfgammaSC = []float64{ - Inf(-1), - -3, - Copysign(0, -1), - 0, - Inf(1), - NaN(), -} -var gammaSC = []float64{ - NaN(), - NaN(), - Inf(-1), - Inf(1), - Inf(1), - NaN(), +var vfgamma = [][2]float64{ + {Inf(1), Inf(1)}, + {Inf(-1), NaN()}, + {0, Inf(1)}, + {Copysign(0, -1), Inf(-1)}, + {NaN(), NaN()}, + {-1, NaN()}, + {-2, NaN()}, + {-3, NaN()}, + {-1e16, NaN()}, + {-1e300, NaN()}, + {1.7e308, Inf(1)}, + + // Test inputs inspired by Python test suite. + // Outputs computed at high precision by PARI/GP. + // If recomputing table entries, be careful to use + // high-precision (%.1000g) formatting of the float64 inputs. + // For example, -2.0000000000000004 is the float64 with exact value + // -2.00000000000000044408920985626161695, and + // gamma(-2.0000000000000004) = -1249999999999999.5386078562728167651513, while + // gamma(-2.00000000000000044408920985626161695) = -1125899906826907.2044875028130093136826. + // Thus the table lists -1.1258999068426235e+15 as the answer. + {0.5, 1.772453850905516}, + {1.5, 0.886226925452758}, + {2.5, 1.329340388179137}, + {3.5, 3.3233509704478426}, + {-0.5, -3.544907701811032}, + {-1.5, 2.363271801207355}, + {-2.5, -0.9453087204829419}, + {-3.5, 0.2700882058522691}, + {0.1, 9.51350769866873}, + {0.01, 99.4325851191506}, + {1e-08, 9.999999942278434e+07}, + {1e-16, 1e+16}, + {0.001, 999.4237724845955}, + {1e-16, 1e+16}, + {1e-308, 1e+308}, + {5.6e-309, 1.7857142857142864e+308}, + {5.5e-309, Inf(1)}, + {1e-309, Inf(1)}, + {1e-323, Inf(1)}, + {5e-324, Inf(1)}, + {-0.1, -10.686287021193193}, + {-0.01, -100.58719796441078}, + {-1e-08, -1.0000000057721567e+08}, + {-1e-16, -1e+16}, + {-0.001, -1000.5782056293586}, + {-1e-16, -1e+16}, + {-1e-308, -1e+308}, + {-5.6e-309, -1.7857142857142864e+308}, + {-5.5e-309, Inf(-1)}, + {-1e-309, Inf(-1)}, + {-1e-323, Inf(-1)}, + {-5e-324, Inf(-1)}, + {-0.9999999999999999, -9.007199254740992e+15}, + {-1.0000000000000002, 4.5035996273704955e+15}, + {-1.9999999999999998, 2.2517998136852485e+15}, + {-2.0000000000000004, -1.1258999068426235e+15}, + {-100.00000000000001, -7.540083334883109e-145}, + {-99.99999999999999, 7.540083334884096e-145}, + {17, 2.0922789888e+13}, + {171, 7.257415615307999e+306}, + {171.6, 1.5858969096672565e+308}, + {171.624, 1.7942117599248104e+308}, + {171.625, Inf(1)}, + {172, Inf(1)}, + {2000, Inf(1)}, + {-100.5, -3.3536908198076787e-159}, + {-160.5, -5.255546447007829e-286}, + {-170.5, -3.3127395215386074e-308}, + {-171.5, 0}, // TODO: 1.9316265431712e-310 + {-176.5, Copysign(0, -1)}, // TODO: -1.196e-321 + {-177.5, 0}, // TODO: 5e-324 + {-178.5, Copysign(0, -1)}, + {-179.5, 0}, + {-201.0001, 0}, + {-202.9999, Copysign(0, -1)}, + {-1000.5, Copysign(0, -1)}, + {-1.0000000003e+09, Copysign(0, -1)}, + {-4.5035996273704955e+15, 0}, + {-63.349078729022985, 4.177797167776188e-88}, + {-127.45117632943295, 1.183111089623681e-214}, } var vfhypotSC = [][2]float64{ @@ -2147,9 +2214,18 @@ func TestGamma(t *testing.T) { t.Errorf("Gamma(%g) = %g, want %g", vf[i], f, gamma[i]) } } - for i := 0; i < len(vfgammaSC); i++ { - if f := Gamma(vfgammaSC[i]); !alike(gammaSC[i], f) { - t.Errorf("Gamma(%g) = %g, want %g", vfgammaSC[i], f, gammaSC[i]) + for _, g := range vfgamma { + f := Gamma(g[0]) + var ok bool + if IsNaN(g[1]) || IsInf(g[1], 0) || g[1] == 0 || f == 0 { + ok = alike(g[1], f) + } else if g[0] > -50 && g[0] <= 171 { + ok = veryclose(g[1], f) + } else { + ok = close(g[1], f) + } + if !ok { + t.Errorf("Gamma(%g) = %g, want %g", g[0], f, g[1]) } } } diff --git a/src/math/gamma.go b/src/math/gamma.go index 841ec114a0..c965f275e7 100644 --- a/src/math/gamma.go +++ b/src/math/gamma.go @@ -93,6 +93,9 @@ var _gamS = [...]float64{ // Gamma function computed by Stirling's formula. // The polynomial is valid for 33 <= x <= 172. func stirling(x float64) float64 { + if x > 171.625 { + return Inf(1) + } const ( SqrtTwoPi = 2.506628274631000502417 MaxStirling = 143.01608 @@ -130,8 +133,6 @@ func Gamma(x float64) float64 { return Inf(-1) } return Inf(1) - case x < -170.5674972726612 || x > 171.61447887182298: - return Inf(1) } q := Abs(x) p := Floor(q) @@ -139,8 +140,11 @@ func Gamma(x float64) float64 { if x >= 0 { return stirling(x) } + // Note: x is negative but (checked above) not a negative integer, + // so x must be small enough to be in range for conversion to int64. + // If |x| were >= 2⁶³ it would have to be an integer. signgam := 1 - if ip := int(p); ip&1 == 0 { + if ip := int64(p); ip&1 == 0 { signgam = -1 } z := q - p