From 26f0c83eb8ebed8e00ad7bdb4cd23672a9c9e40e Mon Sep 17 00:00:00 2001 From: "Charles L. Dorian" Date: Fri, 19 Mar 2010 15:29:22 -0700 Subject: [PATCH] math: add Gamma function R=rsc CC=golang-dev https://golang.org/cl/649041 --- src/pkg/math/Makefile | 1 + src/pkg/math/all_test.go | 46 ++++++++++ src/pkg/math/gamma.go | 188 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 235 insertions(+) create mode 100644 src/pkg/math/gamma.go diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile index 3b82a786b3e..a92a50e0c8b 100644 --- a/src/pkg/math/Makefile +++ b/src/pkg/math/Makefile @@ -53,6 +53,7 @@ ALLGOFILES=\ floor.go\ fmod.go\ frexp.go\ + gamma.go\ hypot.go\ hypot_port.go\ logb.go\ diff --git a/src/pkg/math/all_test.go b/src/pkg/math/all_test.go index 8cb575659af..b28a1f49a6c 100644 --- a/src/pkg/math/all_test.go +++ b/src/pkg/math/all_test.go @@ -286,6 +286,18 @@ var frexp = []fi{ fi{9.1265404584042750000e-01, 1}, fi{-5.4287029803597508250e-01, 4}, } +var gamma = []float64{ + 2.3254348370739963835386613898e+01, + 2.991153837155317076427529816e+03, + -4.561154336726758060575129109e+00, + 7.719403468842639065959210984e-01, + 1.6111876618855418534325755566e+05, + 1.8706575145216421164173224946e+00, + 3.4082787447257502836734201635e+01, + 1.579733951448952054898583387e+00, + 9.3834586598354592860187267089e-01, + -2.093995902923148389186189429e-05, +} var lgamma = []fi{ fi{3.146492141244545774319734e+00, 1}, fi{8.003414490659126375852113e+00, 1}, @@ -736,6 +748,21 @@ var frexpSC = []fi{ fi{NaN(), 0}, } +var vfgammaSC = []float64{ + Inf(-1), + -3, + 0, + Inf(1), + NaN(), +} +var gammaSC = []float64{ + Inf(-1), + Inf(1), + Inf(1), + Inf(1), + NaN(), +} + var vfhypotSC = [][2]float64{ [2]float64{Inf(-1), Inf(-1)}, [2]float64{Inf(-1), 0}, @@ -1278,6 +1305,19 @@ func TestFrexp(t *testing.T) { } } +func TestGamma(t *testing.T) { + for i := 0; i < len(vf); i++ { + if f := Gamma(vf[i]); !close(gamma[i], f) { + t.Errorf("Gamma(%g) = %g, want %g\n", 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\n", vfgammaSC[i], f, gammaSC[i]) + } + } +} + func TestHypot(t *testing.T) { for i := 0; i < len(vf); i++ { a := Fabs(1e200 * tanh[i] * Sqrt(2)) @@ -1748,6 +1788,12 @@ func BenchmarkFrexp(b *testing.B) { } } +func BenchmarkGamma(b *testing.B) { + for i := 0; i < b.N; i++ { + Gamma(2.5) + } +} + func BenchmarkHypot(b *testing.B) { for i := 0; i < b.N; i++ { Hypot(3, 4) diff --git a/src/pkg/math/gamma.go b/src/pkg/math/gamma.go new file mode 100644 index 00000000000..4c5b17d05c7 --- /dev/null +++ b/src/pkg/math/gamma.go @@ -0,0 +1,188 @@ +// 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. + +package math + +// The original C code, the long comment, and the constants +// below are from http://netlib.sandia.gov/cephes/cprob/gamma.c. +// The go code is a simplified version of the original C. +// +// tgamma.c +// +// Gamma function +// +// SYNOPSIS: +// +// double x, y, tgamma(); +// extern int signgam; +// +// y = tgamma( x ); +// +// DESCRIPTION: +// +// Returns gamma function of the argument. The result is +// correctly signed, and the sign (+1 or -1) is also +// returned in a global (extern) variable named signgam. +// This variable is also filled in by the logarithmic gamma +// function lgamma(). +// +// Arguments |x| <= 34 are reduced by recurrence and the function +// approximated by a rational function of degree 6/7 in the +// interval (2,3). Large arguments are handled by Stirling's +// formula. Large negative arguments are made positive using +// a reflection formula. +// +// ACCURACY: +// +// Relative error: +// arithmetic domain # trials peak rms +// DEC -34, 34 10000 1.3e-16 2.5e-17 +// IEEE -170,-33 20000 2.3e-15 3.3e-16 +// IEEE -33, 33 20000 9.4e-16 2.2e-16 +// IEEE 33, 171.6 20000 2.3e-15 3.2e-16 +// +// Error for arguments outside the test range will be larger +// owing to error amplification by the exponential function. +// +// Cephes Math Library Release 2.8: June, 2000 +// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier +// +// The readme file at http://netlib.sandia.gov/cephes/ says: +// Some software in this archive may be from the book _Methods and +// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster +// International, 1989) or from the Cephes Mathematical Library, a +// commercial product. In either event, it is copyrighted by the author. +// What you see here may be used freely but it comes with no support or +// guarantee. +// +// The two known misprints in the book are repaired here in the +// source listings for the gamma function and the incomplete beta +// integral. +// +// Stephen L. Moshier +// moshier@na-net.ornl.gov + +var _P = []float64{ + 1.60119522476751861407e-04, + 1.19135147006586384913e-03, + 1.04213797561761569935e-02, + 4.76367800457137231464e-02, + 2.07448227648435975150e-01, + 4.94214826801497100753e-01, + 9.99999999999999996796e-01, +} +var _Q = []float64{ + -2.31581873324120129819e-05, + 5.39605580493303397842e-04, + -4.45641913851797240494e-03, + 1.18139785222060435552e-02, + 3.58236398605498653373e-02, + -2.34591795718243348568e-01, + 7.14304917030273074085e-02, + 1.00000000000000000320e+00, +} +var _S = []float64{ + 7.87311395793093628397e-04, + -2.29549961613378126380e-04, + -2.68132617805781232825e-03, + 3.47222221605458667310e-03, + 8.33333333333482257126e-02, +} + +// Gamma function computed by Stirling's formula. +// The polynomial is valid for 33 <= x <= 172. +func stirling(x float64) float64 { + const ( + SqrtTwoPi = 2.506628274631000502417 + MaxStirling = 143.01608 + ) + w := 1 / x + w = 1 + w*((((_S[0]*w+_S[1])*w+_S[2])*w+_S[3])*w+_S[4]) + y := Exp(x) + if x > MaxStirling { // avoid Pow() overflow + v := Pow(x, 0.5*x-0.25) + y = v * (v / y) + } else { + y = Pow(x, x-0.5) / y + } + y = SqrtTwoPi * y * w + return y +} + +// Gamma(x) returns the Gamma function of x. +// +// Special cases are: +// Gamma(Inf) = Inf +// Gamma(-Inf) = -Inf +// Gamma(NaN) = NaN +// Large values overflow to +Inf. +// Negative integer values equal ±Inf. +func Gamma(x float64) float64 { + const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620 + // special cases + switch { + case x < -MaxFloat64 || x != x: // IsInf(x, -1) || IsNaN(x): + return x + case x < -170.5674972726612 || x > 171.61447887182298: + return Inf(1) + } + q := Fabs(x) + p := Floor(q) + if q > 33 { + if x >= 0 { + return stirling(x) + } + signgam := 1 + if ip := int(p); ip&1 == 0 { + signgam = -1 + } + z := q - p + if z > 0.5 { + p = p + 1 + z = q - p + } + z = q * Sin(Pi*z) + if z == 0 { + return Inf(signgam) + } + z = Pi / (Fabs(z) * stirling(q)) + return float64(signgam) * z + } + + // Reduce argument + z := float64(1) + for x >= 3 { + x = x - 1 + z = z * x + } + for x < 0 { + if x > -1e-09 { + goto small + } + z = z / x + x = x + 1 + } + for x < 2 { + if x < 1e-09 { + goto small + } + z = z / x + x = x + 1 + } + + if x == 2 { + return z + } + + x = x - 2 + p = (((((x*_P[0]+_P[1])*x+_P[2])*x+_P[3])*x+_P[4])*x+_P[5])*x + _P[6] + q = ((((((x*_Q[0]+_Q[1])*x+_Q[2])*x+_Q[3])*x+_Q[4])*x+_Q[5])*x+_Q[6])*x + _Q[7] + return z * p / q + +small: + if x == 0 { + return Inf(1) + } + return z / ((1 + Euler*x) * x) +}