1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5 * Permission to use, copy, modify, and distribute this software for any
6 * purpose with or without fee is hereby granted, provided that the above
7 * copyright notice and this permission notice appear in all copies.
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
11 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
12 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
13 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
14 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
15 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
23 * long double x, y, tgammal();
31 * Returns gamma function of the argument. The result is
32 * correctly signed, and the sign (+1 or -1) is also
33 * returned in a global (extern) variable named signgam.
34 * This variable is also filled in by the logarithmic gamma
37 * Arguments |x| <= 13 are reduced by recurrence and the function
38 * approximated by a rational function of degree 7/8 in the
39 * interval (2,3). Large arguments are handled by Stirling's
40 * formula. Large negative arguments are made positive using
41 * a reflection formula.
47 * arithmetic domain # trials peak rms
48 * IEEE -40,+40 10000 3.6e-19 7.9e-20
49 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19
51 * Accuracy for large arguments is dominated by error in powl().
57 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
58 long double tgammal(long double x)
62 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
64 tgamma(x+2) = tgamma(x+2) P(x)/Q(x)
69 Relative error spread = 8.4e-23
71 static const long double P[8] = {
72 4.212760487471622013093E-5L,
73 4.542931960608009155600E-4L,
74 4.092666828394035500949E-3L,
75 2.385363243461108252554E-2L,
76 1.113062816019361559013E-1L,
77 3.629515436640239168939E-1L,
78 8.378004301573126728826E-1L,
79 1.000000000000000000009E0L,
81 static const long double Q[9] = {
82 -1.397148517476170440917E-5L,
83 2.346584059160635244282E-4L,
84 -1.237799246653152231188E-3L,
85 -7.955933682494738320586E-4L,
86 2.773706565840072979165E-2L,
87 -4.633887671244534213831E-2L,
88 -2.243510905670329164562E-1L,
89 4.150160950588455434583E-1L,
90 9.999999999999999999908E-1L,
94 static const long double P[] = {
95 -3.01525602666895735709e0L,
96 -3.25157411956062339893e1L,
97 -2.92929976820724030353e2L,
98 -1.70730828800510297666e3L,
99 -7.96667499622741999770e3L,
100 -2.59780216007146401957e4L,
101 -5.99650230220855581642e4L,
102 -7.15743521530849602425e4L
104 static const long double Q[] = {
105 1.00000000000000000000e0L,
106 -1.67955233807178858919e1L,
107 8.85946791747759881659e1L,
108 5.69440799097468430177e1L,
109 -1.98526250512761318471e3L,
110 3.31667508019495079814e3L,
111 1.60577839621734713377e4L,
112 -2.97045081369399940529e4L,
113 -7.15743521530849602412e4L
116 #define MAXGAML 1755.455L
117 /*static const long double LOGPI = 1.14472988584940017414L;*/
119 /* Stirling's formula for the gamma function
120 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
125 Peak error = 9.44e-21
126 Relative error spread = 8.8e-4
128 static const long double STIR[9] = {
129 7.147391378143610789273E-4L,
130 -2.363848809501759061727E-5L,
131 -5.950237554056330156018E-4L,
132 6.989332260623193171870E-5L,
133 7.840334842744753003862E-4L,
134 -2.294719747873185405699E-4L,
135 -2.681327161876304418288E-3L,
136 3.472222222230075327854E-3L,
137 8.333333333333331800504E-2L,
140 #define MAXSTIR 1024.0L
141 static const long double SQTPI = 2.50662827463100050242E0L;
143 /* 1/tgamma(x) = z P(z)
146 * Peak relative error 4.2e-23
148 static const long double S[9] = {
149 -1.193945051381510095614E-3L,
150 7.220599478036909672331E-3L,
151 -9.622023360406271645744E-3L,
152 -4.219773360705915470089E-2L,
153 1.665386113720805206758E-1L,
154 -4.200263503403344054473E-2L,
155 -6.558780715202540684668E-1L,
156 5.772156649015328608253E-1L,
157 1.000000000000000000000E0L,
160 /* 1/tgamma(-x) = z P(z)
163 * Peak relative error 5.16e-23
164 * Relative error spread = 2.5e-24
166 static const long double SN[9] = {
167 1.133374167243894382010E-3L,
168 7.220837261893170325704E-3L,
169 9.621911155035976733706E-3L,
170 -4.219773343731191721664E-2L,
171 -1.665386113944413519335E-1L,
172 -4.200263503402112910504E-2L,
173 6.558780715202536547116E-1L,
174 5.772156649015328608727E-1L,
175 -1.000000000000000000000E0L,
178 static const long double PIL = 3.1415926535897932384626L;
180 /* Gamma function computed by Stirling's formula.
182 static long double stirf(long double x)
187 /* For large x, use rational coefficients from the analytical expansion. */
189 w = (((((6.97281375836585777429E-5L * w
190 + 7.84039221720066627474E-4L) * w
191 - 2.29472093621399176955E-4L) * w
192 - 2.68132716049382716049E-3L) * w
193 + 3.47222222222222222222E-3L) * w
194 + 8.33333333333333333333E-2L) * w
197 w = 1.0L + w * __polevll(w, STIR, 8);
199 if (x > MAXSTIR) { /* Avoid overflow in pow() */
200 v = powl(x, 0.5L * x - 0.25L);
203 y = powl(x, x - 0.5L) / y;
209 long double tgammal(long double x)
228 return (x - x) / (x - x);
237 z = q * sinl(PIL * z);
238 z = fabsl(z) * stirf(q);
239 if (z <= PIL/LDBL_MAX) {
241 return signgam * INFINITY;
255 while (x < -0.03125L) {
269 p = __polevll(x, P, 7);
270 q = __polevll(x, Q, 8);
278 return (x - x) / (x - x);
281 q = z / (x * __polevll(x, SN, 8));
284 q = z / (x * __polevll(x, S, 8));