e_lgammal_r.c revision 755318548d787c7ca36901be064c69efba41c6cf 
1/* @(#)e_lgamma_r.c 1.3 95/01/18 */
2/*
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5 *
6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
9 * is preserved.
10 * ====================================================
11 */
12
13#include <sys/cdefs.h>
14__FBSDID("$FreeBSD$");
15
16/*
17 * See e_lgamma_r.c for complete comments.
18 *
19 * Converted to long double by Steven G. Kargl.
20 */
21
22#include "fpmath.h"
23#include "math.h"
24#include "math_private.h"
25
26static const volatile double vzero = 0;
27
28static const double
29zero=  0,
30half=  0.5,
31one =  1;
32
33static const long double
34pi  =  3.14159265358979323846264338327950288e+00L;
35/*
36 * Domain y in [0x1p-119, 0.28], range ~[-1.4065e-36, 1.4065e-36]:
37 * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-119.1
38 */
39static const long double
40a0  =  7.72156649015328606065120900824024296e-02L,
41a1  =  3.22467033424113218236207583323018498e-01L,
42a2  =  6.73523010531980951332460538330282217e-02L,
43a3  =  2.05808084277845478790009252803463129e-02L,
44a4  =  7.38555102867398526627292839296001626e-03L,
45a5  =  2.89051033074152328576829509522483468e-03L,
46a6  =  1.19275391170326097618357349881842913e-03L,
47a7  =  5.09669524743042462515256340206203019e-04L,
48a8  =  2.23154758453578096143609255559576017e-04L,
49a9  =  9.94575127818397632126978731542755129e-05L,
50a10 =  4.49262367375420471287545895027098145e-05L,
51a11 =  2.05072127845117995426519671481628849e-05L,
52a12 =  9.43948816959096748454087141447939513e-06L,
53a13 =  4.37486780697359330303852050718287419e-06L,
54a14 =  2.03920783892362558276037363847651809e-06L,
55a15 =  9.55191070057967287877923073200324649e-07L,
56a16 =  4.48993286185740853170657139487620560e-07L,
57a17 =  2.13107543597620911675316728179563522e-07L,
58a18 =  9.70745379855304499867546549551023473e-08L,
59a19 =  5.61889970390290257926487734695402075e-08L,
60a20 =  6.42739653024130071866684358960960951e-09L,
61a21 =  3.34491062143649291746195612991870119e-08L,
62a22 = -1.57068547394315223934653011440641472e-08L,
63a23 =  1.30812825422415841213733487745200632e-08L;
64/*
65 * Domain x in [tc-0.24, tc+0.28], range ~[-6.3201e-37, 6.3201e-37]:
66 * |(lgamma(x) - tf) - t(x - tc)| < 2**-120.3.
67 */
68static const long double
69tc  =  1.46163214496836234126265954232572133e+00L,
70tf  = -1.21486290535849608095514557177691584e-01L,
71tt  =  1.57061739945077675484237837992951704e-36L,
72t0  = -1.99238329499314692728655623767019240e-36L,
73t1  = -6.08453430711711404116887457663281416e-35L,
74t2  =  4.83836122723810585213722380854828904e-01L,
75t3  = -1.47587722994530702030955093950668275e-01L,
76t4  =  6.46249402389127526561003464202671923e-02L,
77t5  = -3.27885410884813055008502586863748063e-02L,
78t6  =  1.79706751152103942928638276067164935e-02L,
79t7  = -1.03142230366363872751602029672767978e-02L,
80t8  =  6.10053602051788840313573150785080958e-03L,
81t9  = -3.68456960831637325470641021892968954e-03L,
82t10 =  2.25976482322181046611440855340968560e-03L,
83t11 = -1.40225144590445082933490395950664961e-03L,
84t12 =  8.78232634717681264035014878172485575e-04L,
85t13 = -5.54194952796682301220684760591403899e-04L,
86t14 =  3.51912956837848209220421213975000298e-04L,
87t15 = -2.24653443695947456542669289367055542e-04L,
88t16 =  1.44070395420840737695611929680511823e-04L,
89t17 = -9.27609865550394140067059487518862512e-05L,
90t18 =  5.99347334438437081412945428365433073e-05L,
91t19 = -3.88458388854572825603964274134801009e-05L,
92t20 =  2.52476631610328129217896436186551043e-05L,
93t21 = -1.64508584981658692556994212457518536e-05L,
94t22 =  1.07434583475987007495523340296173839e-05L,
95t23 = -7.03070407519397260929482550448878399e-06L,
96t24 =  4.60968590693753579648385629003100469e-06L,
97t25 = -3.02765473778832036018438676945512661e-06L,
98t26 =  1.99238771545503819972741288511303401e-06L,
99t27 = -1.31281299822614084861868817951788579e-06L,
100t28 =  8.60844432267399655055574642052370223e-07L,
101t29 = -5.64535486432397413273248363550536374e-07L,
102t30 =  3.99357783676275660934903139592727737e-07L,
103t31 = -2.95849029193433121795495215869311610e-07L,
104t32 =  1.37790144435073124976696250804940384e-07L;
105/*
106 * Domain y in [-0.1, 0.232], range ~[-1.4046e-37, 1.4181e-37]:
107 * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-122.8
108 */
109static const long double
110u0  = -7.72156649015328606065120900824024311e-02L,
111u1  =  4.24082772271938167430983113242482656e-01L,
112u2  =  2.96194003481457101058321977413332171e+00L,
113u3  =  6.49503267711258043997790983071543710e+00L,
114u4  =  7.40090051288150177152835698948644483e+00L,
115u5  =  4.94698036296756044610805900340723464e+00L,
116u6  =  2.00194224610796294762469550684947768e+00L,
117u7  =  4.82073087750608895996915051568834949e-01L,
118u8  =  6.46694052280506568192333848437585427e-02L,
119u9  =  4.17685526755100259316625348933108810e-03L,
120u10 =  9.06361003550314327144119307810053410e-05L,
121v1  =  5.15937098592887275994320496999951947e+00L,
122v2  =  1.14068418766251486777604403304717558e+01L,
123v3  =  1.41164839437524744055723871839748489e+01L,
124v4  =  1.07170702656179582805791063277960532e+01L,
125v5  =  5.14448694179047879915042998453632434e+00L,
126v6  =  1.55210088094585540637493826431170289e+00L,
127v7  =  2.82975732849424562719893657416365673e-01L,
128v8  =  2.86424622754753198010525786005443539e-02L,
129v9  =  1.35364253570403771005922441442688978e-03L,
130v10 =  1.91514173702398375346658943749580666e-05L,
131v11 = -3.25364686890242327944584691466034268e-08L;
132/*
133 * Domain x in (2, 3], range ~[-1.3341e-36, 1.3536e-36]:
134 * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-120.1
135 * with y = x - 2.
136 */
137static const long double
138s0  = -7.72156649015328606065120900824024297e-02L,
139s1  =  1.23221687850916448903914170805852253e-01L,
140s2  =  5.43673188699937239808255378293820020e-01L,
141s3  =  6.31998137119005233383666791176301800e-01L,
142s4  =  3.75885340179479850993811501596213763e-01L,
143s5  =  1.31572908743275052623410195011261575e-01L,
144s6  =  2.82528453299138685507186287149699749e-02L,
145s7  =  3.70262021550340817867688714880797019e-03L,
146s8  =  2.83374000312371199625774129290973648e-04L,
147s9  =  1.15091830239148290758883505582343691e-05L,
148s10 =  2.04203474281493971326506384646692446e-07L,
149s11 =  9.79544198078992058548607407635645763e-10L,
150r1  =  2.58037466655605285937112832039537492e+00L,
151r2  =  2.86289413392776399262513849911531180e+00L,
152r3  =  1.78691044735267497452847829579514367e+00L,
153r4  =  6.89400381446725342846854215600008055e-01L,
154r5  =  1.70135865462567955867134197595365343e-01L,
155r6  =  2.68794816183964420375498986152766763e-02L,
156r7  =  2.64617234244861832870088893332006679e-03L,
157r8  =  1.52881761239180800640068128681725702e-04L,
158r9  =  4.63264813762296029824851351257638558e-06L,
159r10 =  5.89461519146957343083848967333671142e-08L,
160r11 =  1.79027678176582527798327441636552968e-10L;
161/*
162 * Domain z in [8, 0x1p70], range ~[-9.8214e-35, 9.8214e-35]:
163 * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-113.0
164 */
165static const long double
166w0  =  4.18938533204672741780329736405617738e-01L,
167w1  =  8.33333333333333333333333333332852026e-02L,
168w2  = -2.77777777777777777777777727810123528e-03L,
169w3  =  7.93650793650793650791708939493907380e-04L,
170w4  = -5.95238095238095234390450004444370959e-04L,
171w5  =  8.41750841750837633887817658848845695e-04L,
172w6  = -1.91752691752396849943172337347259743e-03L,
173w7  =  6.41025640880333069429106541459015557e-03L,
174w8  = -2.95506530801732133437990433080327074e-02L,
175w9  =  1.79644237328444101596766586979576927e-01L,
176w10 = -1.39240539108367641920172649259736394e+00L,
177w11 =  1.33987701479007233325288857758641761e+01L,
178w12 = -1.56363596431084279780966590116006255e+02L,
179w13 =  2.14830978044410267201172332952040777e+03L,
180w14 = -3.28636067474227378352761516589092334e+04L,
181w15 =  5.06201257747865138432663574251462485e+05L,
182w16 = -6.79720123352023636706247599728048344e+06L,
183w17 =  6.57556601705472106989497289465949255e+07L,
184w18 = -3.26229058141181783534257632389415580e+08L;
185
186static long double
187sin_pil(long double x)
188{
189	volatile long double vz;
190	long double y,z;
191	uint64_t lx, n;
192	uint16_t hx;
193
194	y = -x;
195
196	vz = y+0x1.p112;
197	z = vz-0x1.p112;
198	if (z == y)
199	    return zero;
200
201	vz = y+0x1.p110;
202	EXTRACT_LDBL128_WORDS(hx,lx,n,vz);
203	z = vz-0x1.p110;
204	if (z > y) {
205	    z -= 0.25;
206	    n--;
207	}
208	n &= 7;
209	y = y - z + n * 0.25L;
210
211	switch (n) {
212	    case 0:   y =  __kernel_sinl(pi*y,zero,0); break;
213	    case 1:
214	    case 2:   y =  __kernel_cosl(pi*(0.5-y),zero); break;
215	    case 3:
216	    case 4:   y =  __kernel_sinl(pi*(one-y),zero,0); break;
217	    case 5:
218	    case 6:   y = -__kernel_cosl(pi*(y-1.5),zero); break;
219	    default:  y =  __kernel_sinl(pi*(y-2.0),zero,0); break;
220	    }
221	return -y;
222}
223
224
225long double
226lgammal_r(long double x, int *signgamp)
227{
228	long double nadj,p,p1,p2,p3,q,r,t,w,y,z;
229	uint64_t llx,lx;
230	int i;
231	uint16_t hx;
232
233	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
234
235	if((hx & 0x7fff) == 0x7fff) {	/* erfl(nan)=nan */
236		i = (hx>>15)<<1;
237		return (1-i)+one/x;	/* erfl(+-inf)=+-1 */
238	}
239
240    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
241	*signgamp = 1;
242	if((hx & 0x7fff) == 0x7fff)	/* x is +-Inf or NaN */
243		return x*x;
244	if((hx==0||hx==0x8000)&&lx==0) {
245	    if (hx&0x8000)
246		*signgamp = -1;
247	    return one/vzero;
248	}
249
250   /* purge off tiny and negative arguments */
251	if(fabsl(x)<0x1p-119L) {
252	    if(hx&0x8000) {
253	        *signgamp = -1;
254	        return -logl(-x);
255	    } else return -logl(x);
256	}
257	if(hx&0x8000) {
258	    if(fabsl(x)>=0x1p112)
259		return one/vzero;
260	    t = sin_pil(x);
261	    if(t==zero) return one/vzero;
262	    nadj = logl(pi/fabsl(t*x));
263	    if(t<zero) *signgamp = -1;
264	    x = -x;
265	}
266
267	if(x == 1 || x ==2) r = 0;
268	else if(x<2) {
269	    if(x<=0.8999996185302734) {
270		r = -logl(x);
271		if(x>=0.7315998077392578) {y = 1-x; i= 0;}
272		else if(x>=0.2316399812698364) {y= x-(tc-1); i=1;}
273	  	else {y = x; i=2;}
274	    } else {
275	  	r = 0;
276	        if(x>=1.7316312789916992) {y=2-x;i=0;}
277	        else if(x>=1.2316322326660156) {y=x-tc;i=1;}
278		else {y=x-1;i=2;}
279	    }
280	    switch(i) {
281	      case 0:
282		z = y*y;
283		p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*(a12+z*(a14+z*(a16+
284		    z*(a18+z*(a20+z*a22))))))))));
285		p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*(a13+z*(a15+
286		    z*(a17+z*(a19+z*(a21+z*a23)))))))))));
287		p  = y*p1+p2;
288		r  += (p-y/2); break;
289	      case 1:
290		p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+
291		    y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+
292		    y*(t17+y*(t18+y*(t19+y*(t20+y*(t21+y*(t22+y*(t23+
293		    y*(t24+y*(t25+y*(t26+y*(t27+y*(t28+y*(t29+y*(t30+
294		    y*(t31+y*t32))))))))))))))))))))))))))))));
295		r += (tf + p); break;
296	      case 2:
297		p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*(u7+
298		    y*(u8+y*(u9+y*u10))))))))));
299		p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7+
300		    y*(v8+y*(v9+y*(v10+y*v11))))))))));
301		r += (-y/2 + p1/p2);
302	    }
303	}
304	else if(x<8) {
305	    i = x;
306	    y = x-i;
307	    p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*(s6+y*(s7+y*(s8+
308		y*(s9+y*(s10+y*s11)))))))))));
309	    q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*(r7+y*(r8+
310		y*(r9+y*(r10+y*r11))))))))));
311	    r = y/2+p/q;
312	    z = 1;	/* lgamma(1+s) = log(s) + lgamma(s) */
313	    switch(i) {
314	    case 7: z *= (y+6);		/* FALLTHRU */
315	    case 6: z *= (y+5);		/* FALLTHRU */
316	    case 5: z *= (y+4);		/* FALLTHRU */
317	    case 4: z *= (y+3);		/* FALLTHRU */
318	    case 3: z *= (y+2);		/* FALLTHRU */
319		    r += logl(z); break;
320	    }
321	} else if (x < 0x1p119L) {
322	    t = logl(x);
323	    z = one/x;
324	    y = z*z;
325	    w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*(w8+
326		y*(w9+y*(w10+y*(w11+y*(w12+y*(w13+y*(w14+y*(w15+y*(w16+
327		y*(w17+y*w18)))))))))))))))));
328	    r = (x-half)*(t-one)+w;
329	} else
330	    r =  x*(logl(x)-1);
331	if(hx&0x8000) r = nadj - r;
332	return r;
333}
334