148402Sbostic /*- 248402Sbostic * Copyright (c) 1985 The Regents of the University of California. 348402Sbostic * All rights reserved. 448402Sbostic * 548402Sbostic * %sccs.include.proprietary.c% 634120Sbostic */ 734120Sbostic 824600Szliu #ifndef lint 9*56952Sbostic static char sccsid[] = "@(#)lgamma.c 5.5 (Berkeley) 12/02/92"; 1034120Sbostic #endif /* not lint */ 1124600Szliu 12*56952Sbostic #include <math.h> 13*56952Sbostic #include <errno.h> 1424600Szliu 15*56952Sbostic #include "mathimpl.h" 1624600Szliu 17*56952Sbostic /* TRUNC sets trailing bits in a floating-point number to zero. 18*56952Sbostic * ptrx points to the second half of the floating-point number. 19*56952Sbostic * x1 is a temporary variable. 2024600Szliu */ 2124600Szliu 22*56952Sbostic #if defined(vax) || defined(tahoe) 23*56952Sbostic #define LSHFT_PLUS_1 (134217728.+1.) 24*56952Sbostic #define TRUNC(x, dummy, x1) x1 = x*LSHFT_PLUS_1, x -= x1, x += x1 25*56952Sbostic #else 26*56952Sbostic #define MASK 0xf8000000 27*56952Sbostic #define TRUNC(x, ptrx, dummy) *ptrx &= MASK 28*56952Sbostic #endif 2935679Sbostic 30*56952Sbostic #define x0 0.461632144968362356785 31*56952Sbostic #define x0_lo -.000000000000000015522348162858676890521 32*56952Sbostic #define a0_hi -0.12148629053584961146 33*56952Sbostic #define a0_lo 3.07435645275902737e-18 34*56952Sbostic #define LEFT (1.0 - (x0 + .25)) 35*56952Sbostic #define RIGHT (x0 - .218) 36*56952Sbostic #define lns2pi_hi 0.418945312500000 37*56952Sbostic #define lns2pi_lo -.000006779295327258219670263595 3824600Szliu 39*56952Sbostic #define UNDERFL (1e-1020 * 1e-1020) 40*56952Sbostic int signgam; 41*56952Sbostic 42*56952Sbostic #define r0 -0.02771227512955130520 43*56952Sbostic #define r1 -0.2980729795228150847 44*56952Sbostic #define r2 -0.3257411333183093394 45*56952Sbostic #define r3 -0.1126814387531706041 46*56952Sbostic #define r4 -0.01129130057170225562 47*56952Sbostic #define r5 -2.259650588213369095e-05 48*56952Sbostic #define s0 1.7144571600017144419 49*56952Sbostic #define s1 2.7864695046181946481 50*56952Sbostic #define s2 1.5645463655191798047 51*56952Sbostic #define s3 0.34858463899811098496 52*56952Sbostic #define s4 0.024677593453636563481 53*56952Sbostic 54*56952Sbostic #define p0 -7.721566490153286087127140e-02 55*56952Sbostic #define p1 2.077324848654884653970785e-01 56*56952Sbostic #define p2 3.474331160945523535787194e-01 57*56952Sbostic #define p3 1.724375677840324429295524e-01 58*56952Sbostic #define p4 3.546181984297784658205969e-02 59*56952Sbostic #define p5 2.866163630124151557532506e-03 60*56952Sbostic #define p6 6.143168512963655570532770e-05 61*56952Sbostic #define q0 1.000000000000000000000000e+00 62*56952Sbostic #define q1 1.485897307300750567469226e+00 63*56952Sbostic #define q2 8.336064915387758261556045e-01 64*56952Sbostic #define q3 2.185469782512070977585556e-01 65*56952Sbostic #define q4 2.671060746990048983840354e-02 66*56952Sbostic #define q5 1.296631961305647947057711e-03 67*56952Sbostic #define q6 1.484566091079246905938151e-05 68*56952Sbostic 69*56952Sbostic #define NP2 8 70*56952Sbostic double P2[] = { 71*56952Sbostic .0833333333333333148296162562474, 72*56952Sbostic -.00277777777774548123579378966497, 73*56952Sbostic .000793650778754435631476282786423, 74*56952Sbostic -.000595235082566672847950717262222, 75*56952Sbostic .000841428560346653702135821806252, 76*56952Sbostic -.00189773526463879200348872089421, 77*56952Sbostic .00569394463439411649408050664078, 78*56952Sbostic -.0144705562421428915453880392761 7924600Szliu }; 8024600Szliu 81*56952Sbostic static double neg_lgam __P((double)); 82*56952Sbostic static double small_lgam __P((double)); 83*56952Sbostic static double large_lgam __P((double)); 8435679Sbostic 8524600Szliu double 86*56952Sbostic lgamma(x) 87*56952Sbostic double x; 8824600Szliu { 89*56952Sbostic double zero = 0.0, one = 1.0; 90*56952Sbostic double r; 91*56952Sbostic signgam = 1; 92*56952Sbostic if (!finite(x)) { 93*56952Sbostic errno = EDOM; 94*56952Sbostic if (x < 0) 95*56952Sbostic x = -x; 96*56952Sbostic else if (x > 0) 97*56952Sbostic errno = ERANGE; 98*56952Sbostic return (x); 99*56952Sbostic } 100*56952Sbostic if (x > 6 + RIGHT) { 101*56952Sbostic if (x > 1.0e20) 102*56952Sbostic return (x*(log(x)-one)); 103*56952Sbostic r = large_lgam(x); 104*56952Sbostic return (r); 105*56952Sbostic } else if (x > 1e-17) 106*56952Sbostic return (small_lgam(x)); 107*56952Sbostic else if(x > -1e-17) { 108*56952Sbostic if (x < 0) 109*56952Sbostic signgam = -1, x = -x; 110*56952Sbostic return (-log(x)); 111*56952Sbostic } else 112*56952Sbostic return (neg_lgam(x)); 11324600Szliu } 11424600Szliu 115*56952Sbostic /* Accurate to max(ulp(1/128) absolute, 2^-75 relative) error. */ 11624600Szliu static double 117*56952Sbostic large_lgam(x) 118*56952Sbostic double x; 11924600Szliu { 120*56952Sbostic double z, p, x1; 12124600Szliu int i; 122*56952Sbostic long *pva, *pua; 123*56952Sbostic struct Double t, u, v; 124*56952Sbostic pua = (long *) &u.a, pua++; 125*56952Sbostic pva = (long *) &v.a, pva++; 126*56952Sbostic z = 1.0/(x*x); 127*56952Sbostic for (p = P2[i = NP2-1]; --i >= 0;) 128*56952Sbostic p = P2[i] + p*z; /* error in approximation = 2.8e-18 */ 12924600Szliu 130*56952Sbostic p = p/x; /* ulp = 1.7e-18; error < 1.6ulp */ 131*56952Sbostic /* 0 < frac < 1/64 (at x = 5.5) */ 132*56952Sbostic t = log__D(x); 133*56952Sbostic t.a -= 1.0; 134*56952Sbostic u.a = t.a + t.b; 135*56952Sbostic TRUNC (u.a, pua, x1); /* truncate u.a */ 136*56952Sbostic u.b = (t.a - u.a); 137*56952Sbostic u.b += t.b; 138*56952Sbostic x -= .5; 139*56952Sbostic v.a = x; 140*56952Sbostic TRUNC(v.a, pva, x1); /* truncate v.a */ 141*56952Sbostic v.b = x - v.a; 142*56952Sbostic t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ 143*56952Sbostic t.b = v.b*u.a + x*u.b; 144*56952Sbostic z = t.b + lns2pi_lo; /* return t + lns2pi + p */ 145*56952Sbostic z += p; z += lns2pi_hi; 146*56952Sbostic z += t.a; 147*56952Sbostic return (z); 14824600Szliu } 149*56952Sbostic /* Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.) 150*56952Sbostic It also has correct monotonicity. 151*56952Sbostic */ 15224600Szliu static double 153*56952Sbostic small_lgam(x) 154*56952Sbostic double x; 15524600Szliu { 156*56952Sbostic int xi; 157*56952Sbostic double y, z, t, r = 0, p, q; 158*56952Sbostic struct Double rr; 15924600Szliu 160*56952Sbostic /* Do nasty area near the minimum. No good for 1.25 < x < 2.5 */ 161*56952Sbostic if (x < 2.0 - LEFT && x > 1.0 + RIGHT) { 162*56952Sbostic t = x - 1.0; t -= x0; 163*56952Sbostic z = t - x0_lo; 164*56952Sbostic p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5)))); 165*56952Sbostic q = s0+z*(s1+z*(s2+z*(s3+z*s4))); 166*56952Sbostic r = t*(z*(p/q) - x0_lo) + a0_lo; 167*56952Sbostic r += .5*t*t; r += a0_hi; 168*56952Sbostic return (r); 16924600Szliu } 170*56952Sbostic xi = (x + .5); 171*56952Sbostic y = x - xi; 172*56952Sbostic rr.a = rr.b = 0; 173*56952Sbostic if (y < -LEFT) { /* necessary for 2.5 < x < 2.72.. */ 174*56952Sbostic t = y + (1.0 - x0); 175*56952Sbostic z = t - x0_lo; 176*56952Sbostic p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5)))); 177*56952Sbostic q = s0+z*(s1+z*(s2+z*(s3+z*s4))); 178*56952Sbostic r = t*(z*(p/q) - x0_lo) + a0_lo; 179*56952Sbostic r += .5*t*t; 180*56952Sbostic q = a0_hi; 181*56952Sbostic printf("(0)q = %.18lg r = %.18lg\n", q, r); 182*56952Sbostic } else { 183*56952Sbostic p = y*(p0+y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*p6)))))); 184*56952Sbostic q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6))))); 185*56952Sbostic r = p/q; 186*56952Sbostic q = .5*y; 187*56952Sbostic printf("(1)q = %.18lg r = %.18lg\n", q, r); 18824600Szliu } 189*56952Sbostic printf("y = %lg, r = %.18lg\n", y, r); 190*56952Sbostic z = 1.0; 191*56952Sbostic switch (xi) { 192*56952Sbostic case 6: z = (y + 5); 193*56952Sbostic case 5: z *= (y + 4); 194*56952Sbostic case 4: z *= (y + 3); 195*56952Sbostic case 3: z *= (y + 2); 196*56952Sbostic rr = log__D(z); 197*56952Sbostic printf("%.21lg, %lg\n", rr.a, rr.b); 198*56952Sbostic r += rr.b; q += r; 199*56952Sbostic return(q + rr.a); 200*56952Sbostic case 2: return (q + r); 201*56952Sbostic 202*56952Sbostic case 0: printf("r = %lg\n", r); 203*56952Sbostic r -= log1p(x); printf("\t%lg\n", r); 204*56952Sbostic default: rr = log__D(x); 205*56952Sbostic r -= rr.b; printf("\t%lg\n", r); 206*56952Sbostic } 207*56952Sbostic if (q > .5 *rr.a) { 208*56952Sbostic printf("q = %lg, rr = %lg, r = %lg\n", q, rr.a, r); 209*56952Sbostic q -= rr.a; 210*56952Sbostic return(r + q); 211*56952Sbostic } else 212*56952Sbostic printf("q = %lg, rr = %lg, r = %lg\n", q, rr.a, r); 213*56952Sbostic return((r + q) - rr.a); 21424600Szliu } 21524600Szliu 216*56952Sbostic #define lpi_hi 1.1447298858494001638 217*56952Sbostic #define lpi_lo .0000000000000000102659511627078262 218*56952Sbostic /* Error: within 3.5 ulp for x < 171. For large x, see lgamma. */ 21924600Szliu static double 220*56952Sbostic neg_lgam(x) 221*56952Sbostic double x; 22224600Szliu { 223*56952Sbostic struct Double lg, lsine; 224*56952Sbostic double y, z, one = 1.0, zero = 0.0; 22524600Szliu 226*56952Sbostic z = floor(x + .5); 227*56952Sbostic if (z == x) { 228*56952Sbostic errno = EDOM; 229*56952Sbostic return (one/zero); /* convention: G(-(integer)) -> +oo */ 230*56952Sbostic } 231*56952Sbostic y = ceil(x); 232*56952Sbostic if (y*.5 == ceil(.5*y)) 233*56952Sbostic signgam = -1; 23424600Szliu 235*56952Sbostic x = -x; 236*56952Sbostic z = fabs(x + z); /* 0 < z <= .5 */ 237*56952Sbostic if (z < .25) 238*56952Sbostic z = sin(M_PI*z); 239*56952Sbostic else 240*56952Sbostic z = cos(M_PI*(.5-z)); 241*56952Sbostic z = -log(z*x/M_PI); 242*56952Sbostic 243*56952Sbostic if (x > 6.5) 244*56952Sbostic y -= large_lgam(x); 245*56952Sbostic else 246*56952Sbostic y = -small_lgam (x); 247*56952Sbostic return (y + z); 24824600Szliu } 249