Subversion Repositories openrisc

// 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 _gamP = [...]float64{

        1.60119522476751861407e-04,

        1.19135147006586384913e-03,

        1.04213797561761569935e-02,

        4.76367800457137231464e-02,

        2.07448227648435975150e-01,

        4.94214826801497100753e-01,

        9.99999999999999996796e-01,

}

var _gamQ = [...]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 _gamS = [...]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*((((_gamS[0]*w+_gamS[1])*w+_gamS[2])*w+_gamS[3])*w+_gamS[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(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 IsInf(x, -1) || IsNaN(x):

                return x

        case x < -170.5674972726612 || x > 171.61447887182298:

                return Inf(1)

        }

        q := Abs(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 / (Abs(z) * stirling(q))

                return float64(signgam) * z

        }

        // Reduce argument

        z := 1.0

        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*_gamP[0]+_gamP[1])*x+_gamP[2])*x+_gamP[3])*x+_gamP[4])*x+_gamP[5])*x + _gamP[6]

        q = ((((((x*_gamQ[0]+_gamQ[1])*x+_gamQ[2])*x+_gamQ[3])*x+_gamQ[4])*x+_gamQ[5])*x+_gamQ[6])*x + _gamQ[7]

        return z * p / q

small:

        if x == 0 {

                return Inf(1)

        }

        return z / ((1 + Euler*x) * x)

}