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