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 _P = [...]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 _Q = [...]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 _S = [...]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*((((_S[0]*w+_S[1])*w+_S[2])*w+_S[3])*w+_S[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(-Inf) = -Inf
 118 //      Gamma(NaN) = NaN
 119 // Large values overflow to +Inf.
 120 // Negative integer values equal ±Inf.
 121 func Gamma(x float64) float64 {
 122         const Euler = 0.57721566490153286060651209008240243104215933593992 // A001620
 123         // special cases
 124         switch {
 125         case x < -MaxFloat64 || x != x: // IsInf(x, -1) || IsNaN(x):
 126                 return x
 127         case x < -170.5674972726612 || x > 171.61447887182298:
 128                 return Inf(1)
 129         }
 130         q := Abs(x)
 131         p := Floor(q)
 132         if q > 33 {
 133                 if x >= 0 {
 134                         return stirling(x)
 135                 }
 136                 signgam := 1
 137                 if ip := int(p); ip&1 == 0 {
 138                         signgam = -1
 139                 }
 140                 z := q - p
 141                 if z > 0.5 {
 142                         p = p + 1
 143                         z = q - p
 144                 }
 145                 z = q * Sin(Pi*z)
 146                 if z == 0 {
 147                         return Inf(signgam)
 148                 }
 149                 z = Pi / (Abs(z) * stirling(q))
 150                 return float64(signgam) * z
 151         }
 152
 153         // Reduce argument
 154         z := 1.0
 155         for x >= 3 {
 156                 x = x - 1
 157                 z = z * x
 158         }
 159         for x < 0 {
 160                 if x > -1e-09 {
 161                         goto small
 162                 }
 163                 z = z / x
 164                 x = x + 1
 165         }
 166         for x < 2 {
 167                 if x < 1e-09 {
 168                         goto small
 169                 }
 170                 z = z / x
 171                 x = x + 1
 172         }
 173
 174         if x == 2 {
 175                 return z
 176         }
 177
 178         x = x - 2
 179         p = (((((x*_P[0]+_P[1])*x+_P[2])*x+_P[3])*x+_P[4])*x+_P[5])*x + _P[6]
 180         q = ((((((x*_Q[0]+_Q[1])*x+_Q[2])*x+_Q[3])*x+_Q[4])*x+_Q[5])*x+_Q[6])*x + _Q[7]
 181         return z * p / q
 182
 183 small:
 184         if x == 0 {
 185                 return Inf(1)
 186         }
 187         return z / ((1 + Euler*x) * x)
 188 }