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