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/* 15 * See e_lgamma_r.c for complete comments. 16 * 17 * Converted to long double by Steven G. Kargl. 18 */ 19 20#include "math.h" 21#include "math_private.h" 22 23static const volatile double vzero = 0; 24 25static const double 26zero= 0, 27half= 0.5, 28one = 1; 29 30static const union ieee_ext_u 31piu = LD80C(0xc90fdaa22168c235, 1, 3.14159265358979323851e+00L); 32#define pi (piu.extu_ld) 33/* 34 * Domain y in [0x1p-70, 0.27], range ~[-4.5264e-22, 4.5264e-22]: 35 * |(lgamma(2 - y) + y / 2) / y - a(y)| < 2**-70.9 36 */ 37static const union ieee_ext_u 38a0u = LD80C(0x9e233f1bed863d26, -4, 7.72156649015328606028e-02L), 39a1u = LD80C(0xa51a6625307d3249, -2, 3.22467033424113218889e-01L), 40a2u = LD80C(0x89f000d2abafda8c, -4, 6.73523010531979398946e-02L), 41a3u = LD80C(0xa8991563eca75f26, -6, 2.05808084277991211934e-02L), 42a4u = LD80C(0xf2027e10634ce6b6, -8, 7.38555102796070454026e-03L), 43a5u = LD80C(0xbd6eb76dd22187f4, -9, 2.89051035162703932972e-03L), 44a6u = LD80C(0x9c562ab05e0458ed, -10, 1.19275351624639999297e-03L), 45a7u = LD80C(0x859baed93ee48e46, -11, 5.09674593842117925320e-04L), 46a8u = LD80C(0xe9f28a4432949af2, -13, 2.23109648015769155122e-04L), 47a9u = LD80C(0xd12ad0d9b93c6bb0, -14, 9.97387167479808509830e-05L), 48a10u= LD80C(0xb7522643c78a219b, -15, 4.37071076331030136818e-05L), 49a11u= LD80C(0xca024dcdece2cb79, -16, 2.40813493372040143061e-05L), 50a12u= LD80C(0xbb90fb6968ebdbf9, -19, 2.79495621083634031729e-06L), 51a13u= LD80C(0xba1c9ffeeae07b37, -17, 1.10931287015513924136e-05L); 52#define a0 (a0u.extu_ld) 53#define a1 (a1u.extu_ld) 54#define a2 (a2u.extu_ld) 55#define a3 (a3u.extu_ld) 56#define a4 (a4u.extu_ld) 57#define a5 (a5u.extu_ld) 58#define a6 (a6u.extu_ld) 59#define a7 (a7u.extu_ld) 60#define a8 (a8u.extu_ld) 61#define a9 (a9u.extu_ld) 62#define a10 (a10u.extu_ld) 63#define a11 (a11u.extu_ld) 64#define a12 (a12u.extu_ld) 65#define a13 (a13u.extu_ld) 66/* 67 * Domain x in [tc-0.24, tc+0.28], range ~[-6.1205e-22, 6.1205e-22]: 68 * |(lgamma(x) - tf) - t(x - tc)| < 2**-70.5 69 */ 70static const union ieee_ext_u 71tcu = LD80C(0xbb16c31ab5f1fb71, 0, 1.46163214496836234128e+00L), 72tfu = LD80C(0xf8cdcde61c520e0f, -4, -1.21486290535849608093e-01L), 73ttu = LD80C(0xd46ee54b27d4de99, -69, -2.81152980996018785880e-21L), 74t0u = LD80C(0x80b9406556a62a6b, -68, 3.40728634996055147231e-21L), 75t1u = LD80C(0xc7e9c6f6df3f8c39, -67, -1.05833162742737073665e-20L), 76t2u = LD80C(0xf7b95e4771c55d51, -2, 4.83836122723810583532e-01L), 77t3u = LD80C(0x97213c6e35e119ff, -3, -1.47587722994530691476e-01L), 78t4u = LD80C(0x845a14a6a81dc94b, -4, 6.46249402389135358063e-02L), 79t5u = LD80C(0x864d46fa89997796, -5, -3.27885410884846056084e-02L), 80t6u = LD80C(0x93373cbd00297438, -6, 1.79706751150707171293e-02L), 81t7u = LD80C(0xa8fcfca7eddc8d1d, -7, -1.03142230361450732547e-02L), 82t8u = LD80C(0xc7e7015ff4bc45af, -8, 6.10053603296546099193e-03L), 83t9u = LD80C(0xf178d2247adc5093, -9, -3.68456964904901200152e-03L), 84t10u = LD80C(0x94188d58f12e5e9f, -9, 2.25976420273774583089e-03L), 85t11u = LD80C(0xb7cbaef14e1406f1, -10, -1.40224943666225639823e-03L), 86t12u = LD80C(0xe63a671e6704ea4d, -11, 8.78250640744776944887e-04L), 87t13u = LD80C(0x914b6c9cae61783e, -11, -5.54255012657716808811e-04L), 88t14u = LD80C(0xb858f5bdb79276fe, -12, 3.51614951536825927370e-04L), 89t15u = LD80C(0xea73e744c34b9591, -13, -2.23591563824520112236e-04L), 90t16u = LD80C(0x99aeabb0d67ba835, -13, 1.46562869351659194136e-04L), 91t17u = LD80C(0xd7c6938325db2024, -14, -1.02889866046435680588e-04L), 92t18u = LD80C(0xe24cb1e3b0474775, -15, 5.39540265505221957652e-05L); 93#define tc (tcu.extu_ld) 94#define tf (tfu.extu_ld) 95#define tt (ttu.extu_ld) 96#define t0 (t0u.extu_ld) 97#define t1 (t1u.extu_ld) 98#define t2 (t2u.extu_ld) 99#define t3 (t3u.extu_ld) 100#define t4 (t4u.extu_ld) 101#define t5 (t5u.extu_ld) 102#define t6 (t6u.extu_ld) 103#define t7 (t7u.extu_ld) 104#define t8 (t8u.extu_ld) 105#define t9 (t9u.extu_ld) 106#define t10 (t10u.extu_ld) 107#define t11 (t11u.extu_ld) 108#define t12 (t12u.extu_ld) 109#define t13 (t13u.extu_ld) 110#define t14 (t14u.extu_ld) 111#define t15 (t15u.extu_ld) 112#define t16 (t16u.extu_ld) 113#define t17 (t17u.extu_ld) 114#define t18 (t18u.extu_ld) 115/* 116 * Domain y in [-0.1, 0.232], range ~[-8.1938e-22, 8.3815e-22]: 117 * |(lgamma(1 + y) + 0.5 * y) / y - u(y) / v(y)| < 2**-71.2 118 */ 119static const union ieee_ext_u 120u0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L), 121u1u = LD80C(0x98280ee45e4ddd3d, -1, 5.94361239198682739769e-01L), 122u2u = LD80C(0xe330c8ead4130733, 0, 1.77492629495841234275e+00L), 123u3u = LD80C(0xd4a213f1a002ec52, 0, 1.66119622514818078064e+00L), 124u4u = LD80C(0xa5a9ca6f5bc62163, -1, 6.47122051417476492989e-01L), 125u5u = LD80C(0xc980e49cd5b019e6, -4, 9.83903751718671509455e-02L), 126u6u = LD80C(0xff636a8bdce7025b, -9, 3.89691687802305743450e-03L), 127v1u = LD80C(0xbd109c533a19fbf5, 1, 2.95413883330948556544e+00L), 128v2u = LD80C(0xd295cbf96f31f099, 1, 3.29039286955665403176e+00L), 129v3u = LD80C(0xdab8bcfee40496cb, 0, 1.70876276441416471410e+00L), 130v4u = LD80C(0xd2f2dc3638567e9f, -2, 4.12009126299534668571e-01L), 131v5u = LD80C(0xa07d9b0851070f41, -5, 3.91822868305682491442e-02L), 132v6u = LD80C(0xe3cd8318f7adb2c4, -11, 8.68998648222144351114e-04L); 133#define u0 (u0u.extu_ld) 134#define u1 (u1u.extu_ld) 135#define u2 (u2u.extu_ld) 136#define u3 (u3u.extu_ld) 137#define u4 (u4u.extu_ld) 138#define u5 (u5u.extu_ld) 139#define u6 (u6u.extu_ld) 140#define v1 (v1u.extu_ld) 141#define v2 (v2u.extu_ld) 142#define v3 (v3u.extu_ld) 143#define v4 (v4u.extu_ld) 144#define v5 (v5u.extu_ld) 145#define v6 (v6u.extu_ld) 146/* 147 * Domain x in (2, 3], range ~[-3.3648e-22, 3.4416e-22]: 148 * |(lgamma(y+2) - 0.5 * y) / y - s(y)/r(y)| < 2**-72.3 149 * with y = x - 2. 150 */ 151static const union ieee_ext_u 152s0u = LD80C(0x9e233f1bed863d27, -4, -7.72156649015328606095e-02L), 153s1u = LD80C(0xd3ff0dcc7fa91f94, -3, 2.07027640921219389860e-01L), 154s2u = LD80C(0xb2bb62782478ef31, -2, 3.49085881391362090549e-01L), 155s3u = LD80C(0xb49f7438c4611a74, -3, 1.76389518704213357954e-01L), 156s4u = LD80C(0x9a957008fa27ecf9, -5, 3.77401710862930008071e-02L), 157s5u = LD80C(0xda9b389a6ca7a7ac, -9, 3.33566791452943399399e-03L), 158s6u = LD80C(0xbc7a2263faf59c14, -14, 8.98728786745638844395e-05L), 159r1u = LD80C(0xbf5cff5b11477d4d, 0, 1.49502555796294337722e+00L), 160r2u = LD80C(0xd9aec89de08e3da6, -1, 8.50323236984473285866e-01L), 161r3u = LD80C(0xeab7ae5057c443f9, -3, 2.29216312078225806131e-01L), 162r4u = LD80C(0xf29707d9bd2b1e37, -6, 2.96130326586640089145e-02L), 163r5u = LD80C(0xd376c2f09736c5a3, -10, 1.61334161411590662495e-03L), 164r6u = LD80C(0xc985983d0cd34e3d, -16, 2.40232770710953450636e-05L), 165r7u = LD80C(0xe5c7a4f7fc2ef13d, -25, -5.34997929289167573510e-08L); 166#define s0 (s0u.extu_ld) 167#define s1 (s1u.extu_ld) 168#define s2 (s2u.extu_ld) 169#define s3 (s3u.extu_ld) 170#define s4 (s4u.extu_ld) 171#define s5 (s5u.extu_ld) 172#define s6 (s6u.extu_ld) 173#define r1 (r1u.extu_ld) 174#define r2 (r2u.extu_ld) 175#define r3 (r3u.extu_ld) 176#define r4 (r4u.extu_ld) 177#define r5 (r5u.extu_ld) 178#define r6 (r6u.extu_ld) 179#define r7 (r7u.extu_ld) 180/* 181 * Domain z in [8, 0x1p70], range ~[-3.0235e-22, 3.0563e-22]: 182 * |lgamma(x) - (x - 0.5) * (log(x) - 1) - w(1/x)| < 2**-71.7 183 */ 184static const union ieee_ext_u 185w0u = LD80C(0xd67f1c864beb4a69, -2, 4.18938533204672741776e-01L), 186w1u = LD80C(0xaaaaaaaaaaaaaaa1, -4, 8.33333333333333332678e-02L), 187w2u = LD80C(0xb60b60b60b5491c9, -9, -2.77777777777760927870e-03L), 188w3u = LD80C(0xd00d00cf58aede4c, -11, 7.93650793490637233668e-04L), 189w4u = LD80C(0x9c09bf626783d4a5, -11, -5.95238023926039051268e-04L), 190w5u = LD80C(0xdca7cadc5baa517b, -11, 8.41733700408000822962e-04L), 191w6u = LD80C(0xfb060e361e1ffd07, -10, -1.91515849570245136604e-03L), 192w7u = LD80C(0xcbd5101bb58d1f2b, -8, 6.22046743903262649294e-03L), 193w8u = LD80C(0xad27a668d32c821b, -6, -2.11370706734662081843e-02L); 194#define w0 (w0u.extu_ld) 195#define w1 (w1u.extu_ld) 196#define w2 (w2u.extu_ld) 197#define w3 (w3u.extu_ld) 198#define w4 (w4u.extu_ld) 199#define w5 (w5u.extu_ld) 200#define w6 (w6u.extu_ld) 201#define w7 (w7u.extu_ld) 202#define w8 (w8u.extu_ld) 203 204static long double 205sin_pil(long double x) 206{ 207 volatile long double vz; 208 long double y,z; 209 uint64_t n; 210 uint16_t hx __unused; 211 212 y = -x; 213 214 vz = y+0x1p63; 215 z = vz-0x1p63; 216 if (z == y) 217 return zero; 218 219 vz = y+0x1p61; 220 EXTRACT_LDBL80_WORDS(hx,n,vz); 221 z = vz-0x1p61; 222 if (z > y) { 223 z -= 0.25; /* adjust to round down */ 224 n--; 225 } 226 n &= 7; /* octant of y mod 2 */ 227 y = y - z + n * 0.25; /* y mod 2 */ 228 229 switch (n) { 230 case 0: y = __kernel_sinl(pi*y,zero,0); break; 231 case 1: 232 case 2: y = __kernel_cosl(pi*(0.5-y),zero); break; 233 case 3: 234 case 4: y = __kernel_sinl(pi*(one-y),zero,0); break; 235 case 5: 236 case 6: y = -__kernel_cosl(pi*(y-1.5),zero); break; 237 default: y = __kernel_sinl(pi*(y-2.0),zero,0); break; 238 } 239 return -y; 240} 241 242long double 243lgammal_r(long double x, int *signgamp) 244{ 245 long double nadj,p,p1,p2,q,r,t,w,y,z; 246 uint64_t lx; 247 int i; 248 uint16_t hx,ix; 249 250 EXTRACT_LDBL80_WORDS(hx,lx,x); 251 252 /* purge +-Inf and NaNs */ 253 *signgamp = 1; 254 ix = hx&0x7fff; 255 if(ix==0x7fff) return x*x; 256 257 ENTERI(); 258 259 /* purge +-0 and tiny arguments */ 260 *signgamp = 1-2*(hx>>15); 261 if(ix<0x3fff-67) { /* |x|<2**-(p+3), return -log(|x|) */ 262 if((ix|lx)==0) 263 RETURNI(one/vzero); 264 RETURNI(-logl(fabsl(x))); 265 } 266 267 /* purge negative integers and start evaluation for other x < 0 */ 268 if(hx&0x8000) { 269 *signgamp = 1; 270 if(ix>=0x3fff+63) /* |x|>=2**(p-1), must be -integer */ 271 RETURNI(one/vzero); 272 t = sin_pil(x); 273 if(t==zero) RETURNI(one/vzero); /* -integer */ 274 nadj = logl(pi/fabsl(t*x)); 275 if(t<zero) *signgamp = -1; 276 x = -x; 277 } 278 279 /* purge 1 and 2 */ 280 if((ix==0x3fff || ix==0x4000) && lx==0x8000000000000000ULL) r = 0; 281 /* for x < 2.0 */ 282 else if(ix<0x4000) { 283 /* 284 * XXX Supposedly, one can use the following information to replace the 285 * XXX FP rational expressions. A similar approach is appropriate 286 * XXX for ld128, but one (may need?) needs to consider llx, too. 287 * 288 * 8.9999961853027344e-01 3ffe e666600000000000 289 * 7.3159980773925781e-01 3ffe bb4a200000000000 290 * 2.3163998126983643e-01 3ffc ed33080000000000 291 * 1.7316312789916992e+00 3fff dda6180000000000 292 * 1.2316322326660156e+00 3fff 9da6200000000000 293 */ 294 if(x<8.9999961853027344e-01) { 295 r = -logl(x); 296 if(x>=7.3159980773925781e-01) {y = 1-x; i= 0;} 297 else if(x>=2.3163998126983643e-01) {y= x-(tc-1); i=1;} 298 else {y = x; i=2;} 299 } else { 300 r = 0; 301 if(x>=1.7316312789916992e+00) {y=2-x;i=0;} 302 else if(x>=1.2316322326660156e+00) {y=x-tc;i=1;} 303 else {y=x-1;i=2;} 304 } 305 switch(i) { 306 case 0: 307 z = y*y; 308 p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*(a10+z*a12))))); 309 p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*(a11+z*a13)))))); 310 p = y*p1+p2; 311 r += p-y/2; break; 312 case 1: 313 p = t0+y*t1+tt+y*y*(t2+y*(t3+y*(t4+y*(t5+y*(t6+y*(t7+y*(t8+ 314 y*(t9+y*(t10+y*(t11+y*(t12+y*(t13+y*(t14+y*(t15+y*(t16+ 315 y*(t17+y*t18)))))))))))))))); 316 r += tf + p; break; 317 case 2: 318 p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*u6)))))); 319 p2 = 1+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*v6))))); 320 r += p1/p2-y/2; 321 } 322 } 323 /* x < 8.0 */ 324 else if(ix<0x4002) { 325 i = x; 326 y = x-i; 327 p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); 328 q = 1+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*(r6+y*r7)))))); 329 r = y/2+p/q; 330 z = 1; /* lgamma(1+s) = log(s) + lgamma(s) */ 331 switch(i) { 332 case 7: z *= (y+6); /* FALLTHRU */ 333 case 6: z *= (y+5); /* FALLTHRU */ 334 case 5: z *= (y+4); /* FALLTHRU */ 335 case 4: z *= (y+3); /* FALLTHRU */ 336 case 3: z *= (y+2); /* FALLTHRU */ 337 r += logl(z); break; 338 } 339 /* 8.0 <= x < 2**(p+3) */ 340 } else if (ix<0x3fff+67) { 341 t = logl(x); 342 z = one/x; 343 y = z*z; 344 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*(w6+y*(w7+y*w8))))))); 345 r = (x-half)*(t-one)+w; 346 /* 2**(p+3) <= x <= inf */ 347 } else 348 r = x*(logl(x)-1); 349 if(hx&0x8000) r = nadj - r; 350 RETURNI(r); 351} 352