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.7 (Berkeley) 12/14/92"; 10 #endif /* not lint */ 11 12 #include <math.h> 13 #include <errno.h> 14 15 #include "mathimpl.h" 16 17 /* Log gamma function. 18 * Error: x > 0 error < 1.3ulp. 19 * x > 4, error < 1ulp. 20 * x > 9, error < .6ulp. 21 * x < 0, all bets are off. 22 * Method: 23 * x > 6: 24 * Use the asymptotic expansion (Stirling's Formula) 25 * 0 < x < 6: 26 * Use gamma(x+1) = x*gamma(x) 27 * Use rational approximation in 28 * the range 1.2, 2.5 29 * x < 0: 30 * Use the reflection formula, 31 * G(1-x)G(x) = PI/sin(PI*x) 32 * Special values: 33 * non-positive integer returns +Inf. 34 * NaN returns NaN 35 */ 36 #if defined(vax) || defined(tahoe) 37 /* double and float have same size exponent field */ 38 #define TRUNC(x) (double) (float) (x) 39 #define _IEEE 0 40 #else 41 #define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 42 #define _IEEE 1 43 #define infnan(x) (zero/zero) 44 #endif 45 46 extern double log1p(double); 47 static double small_lgam(double); 48 static double large_lgam(double); 49 static double neg_lgam(double); 50 static double zero = 0.0, one = 1.0; 51 int signgam; 52 53 #define lns2pi .418938533204672741780329736405 54 #define UNDERFL (1e-1020 * 1e-1020) 55 56 #define LEFT (1.0 - (x0 + .25)) 57 #define RIGHT (x0 - .218) 58 /* 59 /* Constants for approximation in [1.244,1.712] 60 */ 61 #define x0 0.461632144968362356785 62 #define x0_lo -.000000000000000015522348162858676890521 63 #define a0_hi -0.12148629128932952880859 64 #define a0_lo .0000000007534799204229502 65 #define r0 -2.771227512955130520e-002 66 #define r1 -2.980729795228150847e-001 67 #define r2 -3.257411333183093394e-001 68 #define r3 -1.126814387531706041e-001 69 #define r4 -1.129130057170225562e-002 70 #define r5 -2.259650588213369095e-005 71 #define s0 1.714457160001714442e+000 72 #define s1 2.786469504618194648e+000 73 #define s2 1.564546365519179805e+000 74 #define s3 3.485846389981109850e-001 75 #define s4 2.467759345363656348e-002 76 /* 77 * Constants for approximation in [1.71, 2.5] 78 */ 79 #define a1_hi 4.227843350984671344505727574870e-01 80 #define a1_lo 4.670126436531227189e-18 81 #define p0 3.224670334241133695662995251041e-01 82 #define p1 3.569659696950364669021382724168e-01 83 #define p2 1.342918716072560025853732668111e-01 84 #define p3 1.950702176409779831089963408886e-02 85 #define p4 8.546740251667538090796227834289e-04 86 #define q0 1.000000000000000444089209850062e+00 87 #define q1 1.315850076960161985084596381057e+00 88 #define q2 6.274644311862156431658377186977e-01 89 #define q3 1.304706631926259297049597307705e-01 90 #define q4 1.102815279606722369265536798366e-02 91 #define q5 2.512690594856678929537585620579e-04 92 #define q6 -1.003597548112371003358107325598e-06 93 /* 94 * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf]. 95 */ 96 #define pb0 .0833333333333333148296162562474 97 #define pb1 -.00277777777774548123579378966497 98 #define pb2 .000793650778754435631476282786423 99 #define pb3 -.000595235082566672847950717262222 100 #define pb4 .000841428560346653702135821806252 101 #define pb5 -.00189773526463879200348872089421 102 #define pb6 .00569394463439411649408050664078 103 #define pb7 -.0144705562421428915453880392761 104 105 double 106 lgamma(double x) 107 { 108 double r; 109 signgam = 1; 110 if (!finite(x)) 111 if (_IEEE) 112 return (x+x); 113 else return (infnan(EDOM)); 114 115 if (x > 6 + RIGHT) { 116 r = large_lgam(x); 117 return (r); 118 } else if (x > 1e-16) 119 return (small_lgam(x)); 120 else if (x > -1e-16) { 121 if (x < 0) 122 signgam = -1, x = -x; 123 return (-log(x)); 124 } else 125 return (neg_lgam(x)); 126 } 127 128 static double 129 large_lgam(double x) 130 { 131 double z, p, x1; 132 int i; 133 struct Double t, u, v; 134 u = log__D(x); 135 u.a -= 1.0; 136 if (x > 1e15) { 137 v.a = x - 0.5; 138 TRUNC(v.a); 139 v.b = (x - v.a) - 0.5; 140 t.a = u.a*v.a; 141 t.b = x*u.b + v.b*u.a; 142 if (_IEEE == 0 && !finite(t.a)) 143 return(infnan(ERANGE)); 144 return(t.a + t.b); 145 } 146 x1 = 1./x; 147 z = x1*x1; 148 p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7)))))); 149 /* error in approximation = 2.8e-19 */ 150 151 p = p*x1; /* error < 2.3e-18 absolute */ 152 /* 0 < p < 1/64 (at x = 5.5) */ 153 x = x - 0.5; 154 TRUNC(v.a); /* truncate v.a to 26 bits. */ 155 v.b = x - v.a; 156 t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ 157 t.b = v.b*u.a + z*u.b; 158 t.b += p; t.b += lns2pi; /* return t + lns2pi + p */ 159 return (t.a + t.b); 160 } 161 162 static double 163 small_lgam(double x) 164 { 165 int x_int; 166 double y, z, t, r = 0, p, q, hi, lo; 167 struct Double rr; 168 x_int = (x + .5); 169 y = x - x_int; 170 if (x_int <= 2 && y > RIGHT) { 171 t = y - x0; 172 y--; x_int++; 173 goto CONTINUE; 174 } else if (y < -LEFT) { 175 t = y +(1.0-x0); 176 CONTINUE: 177 z = t - x0_lo; 178 p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5)))); 179 q = s0+z*(s1+z*(s2+z*(s3+z*s4))); 180 r = t*(z*(p/q) - x0_lo); 181 t = .5*t*t; 182 z = 1.0; 183 switch (x_int) { 184 case 6: z = (y + 5); 185 case 5: z *= (y + 4); 186 case 4: z *= (y + 3); 187 case 3: z *= (y + 2); 188 rr = log__D(z); 189 rr.b += a0_lo; rr.a += a0_hi; 190 return(((r+rr.b)+t+rr.a)); 191 case 2: return(((r+a0_lo)+t)+a0_hi); 192 case 0: r -= log1p(x); 193 default: rr = log__D(x); 194 rr.a -= a0_hi; rr.b -= a0_lo; 195 return(((r - rr.b) + t) - rr.a); 196 } 197 } else { 198 p = p0+y*(p1+y*(p2+y*(p3+y*p4))); 199 q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6))))); 200 p = p*(y/q); 201 t = (double)(float) y; 202 z = y-t; 203 hi = (double)(float) (p+a1_hi); 204 lo = a1_hi - hi; lo += p; lo += a1_lo; 205 r = lo*y + z*hi; /* q + r = y*(a0+p/q) */ 206 q = hi*t; 207 z = 1.0; 208 switch (x_int) { 209 case 6: z = (y + 5); 210 case 5: z *= (y + 4); 211 case 4: z *= (y + 3); 212 case 3: z *= (y + 2); 213 rr = log__D(z); 214 r += rr.b; r += q; 215 return(rr.a + r); 216 case 2: return (q+ r); 217 case 0: rr = log__D(x); 218 r -= rr.b; r -= log1p(x); 219 r += q; r-= rr.a; 220 return(r); 221 default: rr = log__D(x); 222 r -= rr.b; 223 q -= rr.a; 224 return (r+q); 225 } 226 } 227 } 228 229 #define lpi_hi 1.1447298858494001638 230 #define lpi_lo .0000000000000000102659511627078262 231 /* Error: within 3.5 ulp for x < 171. For large x, see lgamma. */ 232 static double 233 neg_lgam(double x) 234 { 235 double y, z, one = 1.0, zero = 0.0; 236 237 z = floor(x + .5); 238 if (z == x) { /* convention: G(-(integer)) -> +Inf */ 239 if (_IEEE) 240 return (one/zero); 241 else 242 return (infnan(ERANGE)); 243 } 244 y = .5*ceil(x); 245 if (y == ceil(y)) 246 signgam = -1; 247 x = -x; 248 z = fabs(x + z); /* 0 < z <= .5 */ 249 if (z < .25) 250 z = sin(M_PI*z); 251 else 252 z = cos(M_PI*(0.5-z)); 253 z = -log(z*x/M_PI); 254 255 if (x > 6. + RIGHT) 256 y -= large_lgam(x); 257 else 258 y = -small_lgam (x); 259 return (y + z); 260 } 261