libgo/go/math/gamma.go

   1 // Copyright 2010 The Go Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package math
   6
   7 // The original C code, the long comment, and the constants
   8 // below are from http://netlib.sandia.gov/cephes/cprob/gamma.c.
   9 // The go code is a simplified version of the original C.
  10 //
  11 //      tgamma.c
  12 //
  13 //      Gamma function
  14 //
  15 // SYNOPSIS:
  16 //
  17 // double x, y, tgamma();
  18 // extern int signgam;
  19 //
  20 // y = tgamma( x );
  21 //
  22 // DESCRIPTION:
  23 //
  24 // Returns gamma function of the argument. The result is
  25 // correctly signed, and the sign (+1 or -1) is also
  26 // returned in a global (extern) variable named signgam.
  27 // This variable is also filled in by the logarithmic gamma
  28 // function lgamma().
  29 //
  30 // Arguments |x| <= 34 are reduced by recurrence and the function
  31 // approximated by a rational function of degree 6/7 in the
  32 // interval (2,3).  Large arguments are handled by Stirling's
  33 // formula. Large negative arguments are made positive using
  34 // a reflection formula.
  35 //
  36 // ACCURACY:
  37 //
  38 //                      Relative error:
  39 // arithmetic   domain     # trials      peak         rms
  40 //    DEC      -34, 34      10000       1.3e-16     2.5e-17
  41 //    IEEE    -170,-33      20000       2.3e-15     3.3e-16
  42 //    IEEE     -33,  33     20000       9.4e-16     2.2e-16
  43 //    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
  44 //
  45 // Error for arguments outside the test range will be larger
  46 // owing to error amplification by the exponential function.
  47 //
  48 // Cephes Math Library Release 2.8:  June, 2000
  49 // Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
  50 //
  51 // The readme file at http://netlib.sandia.gov/cephes/ says:
  52 //    Some software in this archive may be from the book _Methods and
  53 // Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
  54 // International, 1989) or from the Cephes Mathematical Library, a
  55 // commercial product. In either event, it is copyrighted by the author.
  56 // What you see here may be used freely but it comes with no support or
  57 // guarantee.
  58 //
  59 //   The two known misprints in the book are repaired here in the
  60 // source listings for the gamma function and the incomplete beta
  61 // integral.
  62 //
  63 //   Stephen L. Moshier
  64 //   moshier@na-net.ornl.gov
  65
  66 var _gamP = [...]float64{
  67         1.60119522476751861407e-04,
  68         1.19135147006586384913e-03,
  69         1.04213797561761569935e-02,
  70         4.76367800457137231464e-02,
  71         2.07448227648435975150e-01,
  72         4.94214826801497100753e-01,
  73         9.99999999999999996796e-01,
  74 }
  75 var _gamQ = [...]float64{
  76         -2.31581873324120129819e-05,
  77         5.39605580493303397842e-04,
  78         -4.45641913851797240494e-03,
  79         1.18139785222060435552e-02,
  80         3.58236398605498653373e-02,
  81         -2.34591795718243348568e-01,
  82         7.14304917030273074085e-02,
  83         1.00000000000000000320e+00,
  84 }
  85 var _gamS = [...]float64{
  86         7.87311395793093628397e-04,
  87         -2.29549961613378126380e-04,
  88         -2.68132617805781232825e-03,
  89         3.47222221605458667310e-03,
  90         8.33333333333482257126e-02,
  91 }
  92
  93 // Gamma function computed by Stirling's formula.
  94 // The pair of results must be multiplied together to get the actual answer.
  95 // The multiplication is left to the caller so that, if careful, the caller can avoid
  96 // infinity for 172 <= x <= 180.
  97 // The polynomial is valid for 33 <= x <= 172; larger values are only used
  98 // in reciprocal and produce denormalized floats. The lower precision there
  99 // masks any imprecision in the polynomial.
 100 func stirling(x float64) (float64, float64) {
 101         if x > 200 {
 102                 return Inf(1), 1
 103         }
 104         const (
 105                 SqrtTwoPi   = 2.506628274631000502417
 106                 MaxStirling = 143.01608
 107         )
 108         w := 1 / x
 109         w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
 110         y1 := Exp(x)
 111         y2 := 1.0
 112         if x > MaxStirling { // avoid Pow() overflow
 113                 v := Pow(x, 0.5*x-0.25)
 114                 y1, y2 = v, v/y1
 115         } else {
 116                 y1 = Pow(x, x-0.5) / y1
 117         }
 118         return y1, SqrtTwoPi * w * y2
 119 }
 120
 121 // Gamma returns the Gamma function of x.
 122 //
 123 // Special cases are:
 124 //      Gamma(+Inf) = +Inf
 125 //      Gamma(+0) = +Inf
 126 //      Gamma(-0) = -Inf
 127 //      Gamma(x) = NaN for integer x < 0
 128 //      Gamma(-Inf) = NaN
 129 //      Gamma(NaN) = NaN
 130 func Gamma(x float64) float64 {
 131         const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
 132         // special cases
 133         switch {
 134         case isNegInt(x) || IsInf(x, -1) || IsNaN(x):
 135                 return NaN()
 136         case IsInf(x, 1):
 137                 return Inf(1)
 138         case x == 0:
 139                 if Signbit(x) {
 140                         return Inf(-1)
 141                 }
 142                 return Inf(1)
 143         }
 144         q := Abs(x)
 145         p := Floor(q)
 146         if q > 33 {
 147                 if x >= 0 {
 148                         y1, y2 := stirling(x)
 149                         return y1 * y2
 150                 }
 151                 // Note: x is negative but (checked above) not a negative integer,
 152                 // so x must be small enough to be in range for conversion to int64.
 153                 // If |x| were >= 2⁶³ it would have to be an integer.
 154                 signgam := 1
 155                 if ip := int64(p); ip&1 == 0 {
 156                         signgam = -1
 157                 }
 158                 z := q - p
 159                 if z > 0.5 {
 160                         p = p + 1
 161                         z = q - p
 162                 }
 163                 z = q * Sin(Pi*z)
 164                 if z == 0 {
 165                         return Inf(signgam)
 166                 }
 167                 sq1, sq2 := stirling(q)
 168                 absz := Abs(z)
 169                 d := absz * sq1 * sq2
 170                 if IsInf(d, 0) {
 171                         z = Pi / absz / sq1 / sq2
 172                 } else {
 173                         z = Pi / d
 174                 }
 175                 return float64(signgam) * z
 176         }
 177
 178         // Reduce argument
 179         z := 1.0
 180         for x >= 3 {
 181                 x = x - 1
 182                 z = z * x
 183         }
 184         for x < 0 {
 185                 if x > -1e-09 {
 186                         goto small
 187                 }
 188                 z = z / x
 189                 x = x + 1
 190         }
 191         for x < 2 {
 192                 if x < 1e-09 {
 193                         goto small
 194                 }
 195                 z = z / x
 196                 x = x + 1
 197         }
 198
 199         if x == 2 {
 200                 return z
 201         }
 202
 203         x = x - 2
 204         p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
 205         q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
 206         return z * p / q
 207
 208 small:
 209         if x == 0 {
 210                 return Inf(1)
 211         }
 212         return z / ((1 + Euler*x) * x)
 213 }
 214
 215 func isNegInt(x float64) bool {
 216         if x < 0 {
 217                 _, xf := Modf(x)
 218                 return xf == 0
 219         }
 220         return false
 221 }