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 polynomial is valid for 33 <= x <= 172.
  95 func stirling(x float64) float64 {
  96         const (
  97                 SqrtTwoPi   = 2.506628274631000502417
  98                 MaxStirling = 143.01608
  99         )
 100         w := 1 / x
 101         w = 1 + w*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[4])
 102         y := Exp(x)
 103         if x > MaxStirling { // avoid Pow() overflow
 104                 v := Pow(x, 0.5*x-0.25)
 105                 y = v * (v / y)
 106         } else {
 107                 y = Pow(x, x-0.5) / y
 108         }
 109         y = SqrtTwoPi * y * w
 110         return y
 111 }
 112
 113 // Gamma(x) returns the Gamma function of x.
 114 //
 115 // Special cases are:
 116 //      Gamma(±Inf) = ±Inf
 117 //      Gamma(NaN) = NaN
 118 // Large values overflow to +Inf.
 119 // Negative integer values equal ±Inf.
 120 func Gamma(x float64) float64 {
 121         const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
 122         // special cases
 123         switch {
 124         case x < -MaxFloat64 || x != x: // IsInf(x, -1) || IsNaN(x):
 125                 return x
 126         case x < -170.5674972726612 || x > 171.61447887182298:
 127                 return Inf(1)
 128         }
 129         q := Abs(x)
 130         p := Floor(q)
 131         if q > 33 {
 132                 if x >= 0 {
 133                         return stirling(x)
 134                 }
 135                 signgam := 1
 136                 if ip := int(p); ip&1 == 0 {
 137                         signgam = -1
 138                 }
 139                 z := q - p
 140                 if z > 0.5 {
 141                         p = p + 1
 142                         z = q - p
 143                 }
 144                 z = q * Sin(Pi*z)
 145                 if z == 0 {
 146                         return Inf(signgam)
 147                 }
 148                 z = Pi / (Abs(z) * stirling(q))
 149                 return float64(signgam) * z
 150         }
 151
 152         // Reduce argument
 153         z := 1.0
 154         for x >= 3 {
 155                 x = x - 1
 156                 z = z * x
 157         }
 158         for x < 0 {
 159                 if x > -1e-09 {
 160                         goto small
 161                 }
 162                 z = z / x
 163                 x = x + 1
 164         }
 165         for x < 2 {
 166                 if x < 1e-09 {
 167                         goto small
 168                 }
 169                 z = z / x
 170                 x = x + 1
 171         }
 172
 173         if x == 2 {
 174                 return z
 175         }
 176
 177         x = x - 2
 178         p = (((((x*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]
 179         q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]
 180         return z * p / q
 181
 182 small:
 183         if x == 0 {
 184                 return Inf(1)
 185         }
 186         return z / ((1 + Euler*x) * x)
 187 }