00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
#if defined(LIBM_SCCS) && !defined(lint)
00014
static char rcsid[] =
"$NetBSD: e_lgamma_r.c,v 1.7 1995/05/10 20:45:42 jtc Exp $";
00015
#endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
#include "math.h"
00085
#include "mathP.h"
00086
00087
#ifdef __STDC__
00088
static const double
00089
#else
00090
static double
00091
#endif
00092 two52= 4.50359962737049600000e+15,
00093 half= 5.00000000000000000000e-01,
00094 one = 1.00000000000000000000e+00,
00095 pi = 3.14159265358979311600e+00,
00096 a0 = 7.72156649015328655494e-02,
00097 a1 = 3.22467033424113591611e-01,
00098 a2 = 6.73523010531292681824e-02,
00099 a3 = 2.05808084325167332806e-02,
00100 a4 = 7.38555086081402883957e-03,
00101 a5 = 2.89051383673415629091e-03,
00102 a6 = 1.19270763183362067845e-03,
00103 a7 = 5.10069792153511336608e-04,
00104 a8 = 2.20862790713908385557e-04,
00105 a9 = 1.08011567247583939954e-04,
00106 a10 = 2.52144565451257326939e-05,
00107 a11 = 4.48640949618915160150e-05,
00108 tc = 1.46163214496836224576e+00,
00109 tf = -1.21486290535849611461e-01,
00110
00111 tt = -3.63867699703950536541e-18,
00112 t0 = 4.83836122723810047042e-01,
00113 t1 = -1.47587722994593911752e-01,
00114 t2 = 6.46249402391333854778e-02,
00115 t3 = -3.27885410759859649565e-02,
00116 t4 = 1.79706750811820387126e-02,
00117 t5 = -1.03142241298341437450e-02,
00118 t6 = 6.10053870246291332635e-03,
00119 t7 = -3.68452016781138256760e-03,
00120 t8 = 2.25964780900612472250e-03,
00121 t9 = -1.40346469989232843813e-03,
00122 t10 = 8.81081882437654011382e-04,
00123 t11 = -5.38595305356740546715e-04,
00124 t12 = 3.15632070903625950361e-04,
00125 t13 = -3.12754168375120860518e-04,
00126 t14 = 3.35529192635519073543e-04,
00127 u0 = -7.72156649015328655494e-02,
00128 u1 = 6.32827064025093366517e-01,
00129 u2 = 1.45492250137234768737e+00,
00130 u3 = 9.77717527963372745603e-01,
00131 u4 = 2.28963728064692451092e-01,
00132 u5 = 1.33810918536787660377e-02,
00133 v1 = 2.45597793713041134822e+00,
00134 v2 = 2.12848976379893395361e+00,
00135 v3 = 7.69285150456672783825e-01,
00136 v4 = 1.04222645593369134254e-01,
00137 v5 = 3.21709242282423911810e-03,
00138 s0 = -7.72156649015328655494e-02,
00139 s1 = 2.14982415960608852501e-01,
00140 s2 = 3.25778796408930981787e-01,
00141 s3 = 1.46350472652464452805e-01,
00142 s4 = 2.66422703033638609560e-02,
00143 s5 = 1.84028451407337715652e-03,
00144 s6 = 3.19475326584100867617e-05,
00145 r1 = 1.39200533467621045958e+00,
00146 r2 = 7.21935547567138069525e-01,
00147 r3 = 1.71933865632803078993e-01,
00148 r4 = 1.86459191715652901344e-02,
00149 r5 = 7.77942496381893596434e-04,
00150 r6 = 7.32668430744625636189e-06,
00151 w0 = 4.18938533204672725052e-01,
00152 w1 = 8.33333333333329678849e-02,
00153 w2 = -2.77777777728775536470e-03,
00154 w3 = 7.93650558643019558500e-04,
00155 w4 = -5.95187557450339963135e-04,
00156 w5 = 8.36339918996282139126e-04,
00157 w6 = -1.63092934096575273989e-03;
00158
00159
#ifdef __STDC__
00160
static const double zero= 0.00000000000000000000e+00;
00161
#else
00162 static double zero= 0.00000000000000000000e+00;
00163
#endif
00164
00165
static
00166
#ifdef __GNUC__
00167
__inline__
00168
#endif
00169
#ifdef __STDC__
00170
double sin_pi(
double x)
00171 #
else
00172 double sin_pi(x)
00173 double x;
00174 #endif
00175 {
00176
double y,z;
00177
int n,ix;
00178
00179
GET_HIGH_WORD(ix,x);
00180 ix &= 0x7fffffff;
00181
00182
if(ix<0x3fd00000)
return __kernel_sin(pi*x,
zero,0);
00183 y = -x;
00184
00185
00186
00187
00188
00189 z =
floor(y);
00190
if(z!=y) {
00191 y *= 0.5;
00192 y = 2.0*(y -
floor(y));
00193 n = (
int) (y*4.0);
00194 }
else {
00195
if(ix>=0x43400000) {
00196 y =
zero; n = 0;
00197 }
else {
00198
if(ix<0x43300000) z = y+
two52;
00199
GET_LOW_WORD(n,z);
00200 n &= 1;
00201 y = n;
00202 n<<= 2;
00203 }
00204 }
00205
switch (n) {
00206
case 0: y =
__kernel_sin(pi*y,
zero,0);
break;
00207
case 1:
00208
case 2: y =
__kernel_cos(pi*(0.5-y),
zero);
break;
00209
case 3:
00210
case 4: y =
__kernel_sin(pi*(one-y),
zero,0);
break;
00211
case 5:
00212
case 6: y = -
__kernel_cos(pi*(y-1.5),
zero);
break;
00213
default: y =
__kernel_sin(pi*(y-2.0),
zero,0);
break;
00214 }
00215
return -y;
00216 }
00217
00218
00219
#ifdef __STDC__
00220
double __ieee754_lgamma_r(
double x,
int *signgamp)
00221 #
else
00222 double __ieee754_lgamma_r(x,signgamp)
00223 double x;
int *signgamp;
00224 #endif
00225 {
00226
double t,y,z,nadj=0.0,p,p1,p2,p3,q,r,w;
00227
int i,hx,lx,ix;
00228
00229
EXTRACT_WORDS(hx,lx,x);
00230
00231
00232 *signgamp = 1;
00233 ix = hx&0x7fffffff;
00234
if(ix>=0x7ff00000)
return x*x;
00235
if((ix|lx)==0)
return one/
zero;
00236
if(ix<0x3b900000) {
00237
if(hx<0) {
00238 *signgamp = -1;
00239
return -
__ieee754_log(-x);
00240 }
else return -
__ieee754_log(x);
00241 }
00242
if(hx<0) {
00243
if(ix>=0x43300000)
00244
return one/
zero;
00245 t =
sin_pi(x);
00246
if(t==
zero)
return one/
zero;
00247 nadj =
__ieee754_log(pi/
fabs(t*x));
00248
if(t<
zero) *signgamp = -1;
00249 x = -x;
00250 }
00251
00252
00253
if((((ix-0x3ff00000)|lx)==0)||(((ix-0x40000000)|lx)==0)) r = 0;
00254
00255
else if(ix<0x40000000) {
00256
if(ix<=0x3feccccc) {
00257 r = -
__ieee754_log(x);
00258
if(ix>=0x3FE76944) {y = one-x; i= 0;}
00259
else if(ix>=0x3FCDA661) {y= x-(tc-one); i=1;}
00260
else {y = x; i=2;}
00261 }
else {
00262 r =
zero;
00263
if(ix>=0x3FFBB4C3) {y=2.0-x;i=0;}
00264
else if(ix>=0x3FF3B4C4) {y=x-tc;i=1;}
00265
else {y=x-one;i=2;}
00266 }
00267
switch(i) {
00268
case 0:
00269 z = y*y;
00270 p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10))));
00271 p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11)))));
00272 p = y*p1+p2;
00273 r += (p-0.5*y);
break;
00274
case 1:
00275 z = y*y;
00276 w = z*y;
00277 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12)));
00278 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13)));
00279 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
00280 p = z*p1-(tt-w*(p2+y*p3));
00281 r += (tf + p);
break;
00282
case 2:
00283 p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
00284 p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
00285 r += (-0.5*y + p1/p2);
00286 }
00287 }
00288
else if(ix<0x40200000) {
00289 i = (
int)x;
00290 t =
zero;
00291 y = x-(
double)i;
00292 p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6))))));
00293 q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))));
00294 r = half*y+p/q;
00295 z = one;
00296
switch(i) {
00297
case 7: z *= (y+6.0);
00298
case 6: z *= (y+5.0);
00299
case 5: z *= (y+4.0);
00300
case 4: z *= (y+3.0);
00301
case 3: z *= (y+2.0);
00302 r +=
__ieee754_log(z);
break;
00303 }
00304
00305 }
else if (ix < 0x43900000) {
00306 t =
__ieee754_log(x);
00307 z = one/x;
00308 y = z*z;
00309 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
00310 r = (x-half)*(t-one)+w;
00311 }
else
00312
00313 r = x*(
__ieee754_log(x)-one);
00314
if(hx<0) r = nadj - r;
00315
return r;
00316 }