OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgo/] [go/] [math/] [gamma.go] - Blame information for rev 867

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 747 jeremybenn
// 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 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
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.