GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: lib/libm/src/ld80/e_tgammal.c Lines: 0 74 0.0 %
Date: 2017-11-13 Branches: 0 42 0.0 %

Line Branch Exec Source
1
/*	$OpenBSD: e_tgammal.c,v 1.4 2013/11/12 20:35:19 martynas Exp $	*/
2
3
/*
4
 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5
 *
6
 * Permission to use, copy, modify, and distribute this software for any
7
 * purpose with or without fee is hereby granted, provided that the above
8
 * copyright notice and this permission notice appear in all copies.
9
 *
10
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17
 */
18
19
/*							tgammal.c
20
 *
21
 *	Gamma function
22
 *
23
 *
24
 *
25
 * SYNOPSIS:
26
 *
27
 * long double x, y, tgammal();
28
 * extern int signgam;
29
 *
30
 * y = tgammal( x );
31
 *
32
 *
33
 *
34
 * DESCRIPTION:
35
 *
36
 * Returns gamma function of the argument.  The result is
37
 * correctly signed, and the sign (+1 or -1) is also
38
 * returned in a global (extern) variable named signgam.
39
 * This variable is also filled in by the logarithmic gamma
40
 * function lgamma().
41
 *
42
 * Arguments |x| <= 13 are reduced by recurrence and the function
43
 * approximated by a rational function of degree 7/8 in the
44
 * interval (2,3).  Large arguments are handled by Stirling's
45
 * formula. Large negative arguments are made positive using
46
 * a reflection formula.
47
 *
48
 *
49
 * ACCURACY:
50
 *
51
 *                      Relative error:
52
 * arithmetic   domain     # trials      peak         rms
53
 *    IEEE     -40,+40      10000       3.6e-19     7.9e-20
54
 *    IEEE    -1755,+1755   10000       4.8e-18     6.5e-19
55
 *
56
 * Accuracy for large arguments is dominated by error in powl().
57
 *
58
 */
59
60
#include <float.h>
61
#include <math.h>
62
63
#include "math_private.h"
64
65
/*
66
tgamma(x+2)  = tgamma(x+2) P(x)/Q(x)
67
0 <= x <= 1
68
Relative error
69
n=7, d=8
70
Peak error =  1.83e-20
71
Relative error spread =  8.4e-23
72
*/
73
74
static long double P[8] = {
75
 4.212760487471622013093E-5L,
76
 4.542931960608009155600E-4L,
77
 4.092666828394035500949E-3L,
78
 2.385363243461108252554E-2L,
79
 1.113062816019361559013E-1L,
80
 3.629515436640239168939E-1L,
81
 8.378004301573126728826E-1L,
82
 1.000000000000000000009E0L,
83
};
84
static long double Q[9] = {
85
-1.397148517476170440917E-5L,
86
 2.346584059160635244282E-4L,
87
-1.237799246653152231188E-3L,
88
-7.955933682494738320586E-4L,
89
 2.773706565840072979165E-2L,
90
-4.633887671244534213831E-2L,
91
-2.243510905670329164562E-1L,
92
 4.150160950588455434583E-1L,
93
 9.999999999999999999908E-1L,
94
};
95
96
/*
97
static long double P[] = {
98
-3.01525602666895735709e0L,
99
-3.25157411956062339893e1L,
100
-2.92929976820724030353e2L,
101
-1.70730828800510297666e3L,
102
-7.96667499622741999770e3L,
103
-2.59780216007146401957e4L,
104
-5.99650230220855581642e4L,
105
-7.15743521530849602425e4L
106
};
107
static long double Q[] = {
108
 1.00000000000000000000e0L,
109
-1.67955233807178858919e1L,
110
 8.85946791747759881659e1L,
111
 5.69440799097468430177e1L,
112
-1.98526250512761318471e3L,
113
 3.31667508019495079814e3L,
114
 1.60577839621734713377e4L,
115
-2.97045081369399940529e4L,
116
-7.15743521530849602412e4L
117
};
118
*/
119
#define MAXGAML 1755.455L
120
/*static const long double LOGPI = 1.14472988584940017414L;*/
121
122
/* Stirling's formula for the gamma function
123
tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
124
z(x) = x
125
13 <= x <= 1024
126
Relative error
127
n=8, d=0
128
Peak error =  9.44e-21
129
Relative error spread =  8.8e-4
130
*/
131
132
static long double STIR[9] = {
133
 7.147391378143610789273E-4L,
134
-2.363848809501759061727E-5L,
135
-5.950237554056330156018E-4L,
136
 6.989332260623193171870E-5L,
137
 7.840334842744753003862E-4L,
138
-2.294719747873185405699E-4L,
139
-2.681327161876304418288E-3L,
140
 3.472222222230075327854E-3L,
141
 8.333333333333331800504E-2L,
142
};
143
144
#define MAXSTIR 1024.0L
145
static const long double SQTPI = 2.50662827463100050242E0L;
146
147
/* 1/tgamma(x) = z P(z)
148
 * z(x) = 1/x
149
 * 0 < x < 0.03125
150
 * Peak relative error 4.2e-23
151
 */
152
153
static long double S[9] = {
154
-1.193945051381510095614E-3L,
155
 7.220599478036909672331E-3L,
156
-9.622023360406271645744E-3L,
157
-4.219773360705915470089E-2L,
158
 1.665386113720805206758E-1L,
159
-4.200263503403344054473E-2L,
160
-6.558780715202540684668E-1L,
161
 5.772156649015328608253E-1L,
162
 1.000000000000000000000E0L,
163
};
164
165
/* 1/tgamma(-x) = z P(z)
166
 * z(x) = 1/x
167
 * 0 < x < 0.03125
168
 * Peak relative error 5.16e-23
169
 * Relative error spread =  2.5e-24
170
 */
171
172
static long double SN[9] = {
173
 1.133374167243894382010E-3L,
174
 7.220837261893170325704E-3L,
175
 9.621911155035976733706E-3L,
176
-4.219773343731191721664E-2L,
177
-1.665386113944413519335E-1L,
178
-4.200263503402112910504E-2L,
179
 6.558780715202536547116E-1L,
180
 5.772156649015328608727E-1L,
181
-1.000000000000000000000E0L,
182
};
183
184
static const long double PIL = 3.1415926535897932384626L;
185
186
static long double stirf ( long double );
187
188
/* Gamma function computed by Stirling's formula.
189
 */
190
static long double stirf(long double x)
191
{
192
long double y, w, v;
193
194
w = 1.0L/x;
195
/* For large x, use rational coefficients from the analytical expansion.  */
196
if( x > 1024.0L )
197
	w = (((((6.97281375836585777429E-5L * w
198
		+ 7.84039221720066627474E-4L) * w
199
		- 2.29472093621399176955E-4L) * w
200
		- 2.68132716049382716049E-3L) * w
201
		+ 3.47222222222222222222E-3L) * w
202
		+ 8.33333333333333333333E-2L) * w
203
		+ 1.0L;
204
else
205
	w = 1.0L + w * __polevll( w, STIR, 8 );
206
y = expl(x);
207
if( x > MAXSTIR )
208
	{ /* Avoid overflow in pow() */
209
	v = powl( x, 0.5L * x - 0.25L );
210
	y = v * (v / y);
211
	}
212
else
213
	{
214
	y = powl( x, x - 0.5L ) / y;
215
	}
216
y = SQTPI * y * w;
217
return( y );
218
}
219
220
long double
221
tgammal(long double x)
222
{
223
long double p, q, z;
224
int i;
225
226
signgam = 1;
227
if( isnan(x) )
228
	return(NAN);
229
if(x == INFINITY)
230
	return(INFINITY);
231
if(x == -INFINITY)
232
	return(x - x);
233
if( x == 0.0L )
234
	return( 1.0L / x );
235
q = fabsl(x);
236
237
if( q > 13.0L )
238
	{
239
	if( q > MAXGAML )
240
		goto goverf;
241
	if( x < 0.0L )
242
		{
243
		p = floorl(q);
244
		if( p == q )
245
			return (x - x) / (x - x);
246
		i = p;
247
		if( (i & 1) == 0 )
248
			signgam = -1;
249
		z = q - p;
250
		if( z > 0.5L )
251
			{
252
			p += 1.0L;
253
			z = q - p;
254
			}
255
		z = q * sinl( PIL * z );
256
		z = fabsl(z) * stirf(q);
257
		if( z <= PIL/LDBL_MAX )
258
			{
259
goverf:
260
			return( signgam * INFINITY);
261
			}
262
		z = PIL/z;
263
		}
264
	else
265
		{
266
		z = stirf(x);
267
		}
268
	return( signgam * z );
269
	}
270
271
z = 1.0L;
272
while( x >= 3.0L )
273
	{
274
	x -= 1.0L;
275
	z *= x;
276
	}
277
278
while( x < -0.03125L )
279
	{
280
	z /= x;
281
	x += 1.0L;
282
	}
283
284
if( x <= 0.03125L )
285
	goto small;
286
287
while( x < 2.0L )
288
	{
289
	z /= x;
290
	x += 1.0L;
291
	}
292
293
if( x == 2.0L )
294
	return(z);
295
296
x -= 2.0L;
297
p = __polevll( x, P, 7 );
298
q = __polevll( x, Q, 8 );
299
z = z * p / q;
300
if( z < 0 )
301
	signgam = -1;
302
return z;
303
304
small:
305
if( x == 0.0L )
306
	return (x - x) / (x - x);
307
else
308
	{
309
	if( x < 0.0L )
310
		{
311
		x = -x;
312
		q = z / (x * __polevll( x, SN, 8 ));
313
		signgam = -1;
314
		}
315
	else
316
		q = z / (x * __polevll( x, S, 8 ));
317
	}
318
return q;
319
}