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[] = "@(#)gamma.c 5.1 (Berkeley) 12/02/92"; 10 #endif /* not lint */ 11 12 #include <math.h> 13 #include <errno.h> 14 15 /* METHOD: 16 * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x)) 17 * At negative integers, return +Inf, and set errno. 18 * 19 * x < 6.5: use one of two rational approximation, 20 * to log(G(x)), expanded around 2 for x near integral; 21 * around the minimum for x near half-integral. The two 22 * regions overlap. 23 * In the range [2.0, 2.5], it is necessary to expand around 2. 24 * In the range ~[1.462-.22, 1.462+.25], the expansion around 25 * the minimum is necessary to get <1ulp accuracy. 26 * 27 * x >= 6.5: Use the asymptotic approximation (Stirling's formula.) 28 * 29 * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x)) 30 * 31 * Keep extra precision in multiplying (x-.5)(log(x)-1), to 32 * avoid premature round-off. 33 * 34 * Special values: NaN, +/-Inf, 0, Negative Integers. 35 * Neg integer: Set overflow trap; return +Inf; errno = EDOM 36 * +Inf: Return +Inf; errno = ERANGE; 37 * -Inf: Return +Inf; errno = EDOM; 38 * NaN: Return NaN; errno = EDOM; 39 */ 40 #define x0 .461632144968362356785 41 #define LEFT -.3955078125 42 #define a0_hi 0.88560319441088874992 43 #define a0_lo -.00000000000000004996427036469019695 44 #define lns2pi_hi 0.418945312500000 45 #define lns2pi_lo -.000006779295327258219670263595 46 47 #define UNDERFL (1e-1020 * 1e-1020) 48 double small_gam(double); 49 double smaller_gam(double); 50 struct Double large_gam(double); 51 double neg_gam(double); 52 struct Double ratfun_gam(double, double); 53 /** 54 #define NP 5 55 static double P[] = { 56 0.57410538218150719558252603747, 57 0.24541366696467897812183878159, 58 0.00513973619299223308948265654, 59 0.00129391688253677823901288679, 60 0.00222188711638167000692045683 61 }; 62 #define NQ 9 63 static double Q[] = { 64 1.33984375, 65 0.981446340605471312379393111769, 66 -0.19172028764767945485658628968, 67 -0.13543838178180836462338731962, 68 0.028432780865671299780350622655, 69 0.004720852857293747484312973484, 70 -0.00162320758342873413572482466, 71 8.63879091186865255905247274e-05, 72 5.67776543645974456238616906e-06, 73 -1.1130244665113561369974706e-08 74 }; 75 /**/ 76 #define P0 .621389571821820863029017800727g 77 #define P1 .265757198651533466104979197553, 78 #define P2 .00553859446429917461063308081748, 79 #define P3 .00138456698304096573887145282811, 80 #define P4 .00240659950032711365819348969808 81 82 #define Q0 1.4501953125 83 #define Q1 1.06258521948016171343454061571 84 #define Q2 -.207474561943859936441469926649 85 #define Q3 -.146734131782005422506287573015 86 #define Q4 .0307878176156175520361557573779 87 #define Q5 .00512449347980666221336054633184 88 #define Q6 -.00176012741431666995019222898833 89 #define Q7 .0000935021023573788935372153030556 90 #define Q8 .00000613275507472443958924745652239 91 92 #define Pa0 .0833333333333333148296162562474 93 #define Pa1 -.00277777777774548123579378966497 94 #define Pa2 .000793650778754435631476282786423 95 #define Pa3 -.000595235082566672847950717262222 96 #define Pa4 .000841428560346653702135821806252 97 #define Pa5 -.00189773526463879200348872089421 98 #define Pa6 .00569394463439411649408050664078 99 #define Pa7 -.0144705562421428915453880392761 100 101 static struct Double large_gam __P((double)); 102 static double neg_gam __P((double)); 103 static struct Double ratfun_gam __P((double, double)); 104 static double small_gam __P((double)); 105 static double smaller_gam __P((double)); 106 107 double 108 gamma(x) 109 double x; 110 { 111 double zero = 0; 112 struct Double u; 113 if (x > 6 + x0 + LEFT) { 114 if(x > 171.63) 115 return(1.0/zero); 116 u = large_gam(x); 117 return(exp__D(u.a, u.b)); 118 } else if (x >= 1.0 + LEFT + x0) 119 return (small_gam(x)); 120 else if (x > 1e-18) 121 return (smaller_gam(x)); 122 else if(x > 0) { 123 1 + 1e-20; /* raise inexact flag */ 124 return(1/x); 125 } 126 else if (x <= 0) 127 return (neg_gam(x)); 128 else { /* x = NaN */ 129 errno = EDOM; 130 return (x); 131 } 132 } 133 134 /* TRUNC sets trailing bits in a floating-point number to zero. 135 * is a temporary variable. 136 */ 137 138 #if defined(vax) || defined(tahoe) 139 #define _IEEE 0 140 #define TRUNC(x) (double) (float) (x) 141 #else 142 #define _IEEE 1 143 #define TRUNC(x) *(((int *) &x) + 1) &= 0xf8000000 144 #define infnan(x) 0.0 145 #endif 146 147 /* Accurate to max(ulp(1/128) absolute, 2^-75 relative) error. */ 148 static struct Double 149 large_gam(x) 150 double x; 151 { 152 double z, p; 153 int i; 154 struct Double t, u, v; 155 pua = (long *) &u.a, pua++; 156 pva = (long *) &v.a, pva++; 157 if (x == infinity()) { 158 u.b = 0, u.a = x; 159 return u; 160 } 161 z = 1.0/(x*x); 162 p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*Pa5+z*(Pa6+z*Pa7))))); 163 p = p/x; /* |e| < 2.8e-18 */ 164 /* 0 < p < 1/64 (at x = 5.5) */ 165 u = log__D(x); 166 u.a -= 1.0; 167 v.a = (x -= .5); 168 TRUNC(v.a); 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 /* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */ 173 t.b += lns2pi_lo; t.b += p; /* small pieces ( < 1/64, assuming t < 1e14) */ 174 u.a = lns2pi_hi + t.b; u.a += t.a; 175 u.b = t.a - u.a; 176 u.b += lns2pi_hi; u.b += t.b; 177 return (u); 178 } 179 /* Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.) 180 It also has correct monotonicity. 181 */ 182 static double 183 small_gam(x) 184 double x; 185 { 186 double y, t; 187 struct Double yy, r; 188 pt = ((long *) &t)+1; 189 pra = ((long *) &r.a)+1; 190 y = x - 1; 191 if (y <= 1.0 + (LEFT + x0)) { 192 yy = ratfun_gam(y - x0, 0); 193 return (yy.a + yy.b); 194 } 195 r.a = y--; 196 TRUNC(r.a); 197 yy.a = r.a - 1.0; 198 yy.b = r.b = y - yy.a; 199 for (; --y > 1.0 + (LEFT + x0); yy.a--) { 200 t = r.a*yy.a; 201 r.b = r.a*yy.b + y*r.b; 202 r.a = t + r.b; 203 TRUNC(r.a); 204 t -= r.a; 205 r.b += t; 206 } /* now want r * gamma(y); */ 207 yy = ratfun_gam(y - x0, 0); 208 y = r.b*(yy.a+yy.b) + r.a*yy.b; 209 y += yy.a*r.a; 210 return (y); 211 } 212 /* Good on (0, 1+x0+LEFT]. Accurate to 1ulp on [.25+x0+LEFT, 1+x0+LEFT]. 213 * Below this, x+LEFT-x0 introduces additional rounding errror of .5ulp. 214 */ 215 static double 216 smaller_gam(x) 217 double x; 218 { 219 double t, d; 220 struct Double r, xx; 221 if (x < x0 + LEFT) { 222 t = x, TRUNC(t); 223 d = (t+x)*(x-t); 224 t *= t; 225 xx.a = ((d+t) + x), TRUNC(xx.a); 226 xx.b = x - xx.a; xx.b += t; xx.b += d; 227 t = (1.0-x0); t += x; 228 d = (1.0-x0); d -= t; d += x; 229 x = xx.a + xx.b; 230 } else { 231 xx.a = x, TRUNC(xx.a); 232 xx.b = x - xx.a; 233 t = x - x0; 234 d = (-x0 -t); d += x; 235 } 236 r = ratfun_gam(t, d); 237 d = r.a/x, TRUNC(d); 238 r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b; 239 return (d + r.a/x); 240 } 241 /* returns (z+c)^2 * P(z)/Q(z) + a0 */ 242 static struct Double 243 ratfun_gam(z, c) 244 double z, c; 245 { 246 int i; 247 double p, q, hi, lo; 248 struct Double r; 249 250 q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8))))))); 251 p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4))); 252 253 /* return r.a + r.b = a0 + (z+c)^2 * p/q, with r.a truncated to 26 bits. */ 254 p = p/q; 255 hi = z, TRUNC(hi); /* hi+lo ~= z + c */ 256 lo = z - hi; lo += c; 257 lo *= (hi+z); /* q+lo = (z+c)*(z+c) */ 258 q = hi*hi; 259 hi = (q + lo), TRUNC(hi); /* hi+lo = q+lo */ 260 q -= hi; 261 lo += q; 262 z = p, TRUNC(z); /* z+q = p */ 263 q = p - z; 264 lo = lo*p + hi*q + a0_lo; 265 hi *= z; 266 r.a = hi + a0_hi, TRUNC(r.a); 267 r.b = ((a0_hi-r.a) + hi) + lo; 268 return (r); 269 } 270 #define lpi_hi 1.1447298858494001638 271 #define lpi_lo .0000000000000000102659511627078262 272 /* Error: within 3.5 ulp for x < 171. For large x, see lgamma. */ 273 static double 274 neg_gam(x) 275 double x; 276 { 277 int sgn = 1; 278 struct Double lg, lsine; 279 double y, z, one = 1.0, zero = 0.0; 280 y = floor(x + .5); 281 if (y == x) { 282 if (-x == infinity()) { /* G(-inf) = NaN */ 283 errno = EDOM; 284 return (signaling_nan(0)); 285 } 286 errno = ERANGE; 287 return (one/zero); /* G(-(integer) -> oo */ 288 } 289 z = fabs(x - y); 290 y = ceil(x); 291 if (y*.5 == ceil(.5*y)) { 292 sgn = -1; 293 printf("neg\n"); 294 } 295 if (x < 1 - (6 + x0 + LEFT)) { 296 if(x < -190) { 297 UNDERFL; 298 z = .5*ceil(x); 299 if (z==ceil(z)) return (-0); 300 else return (0); 301 } 302 y = 1 - x; 303 if (1 - y == x) { 304 lg = large_gam(y); 305 lsine = log__D(sin(M_PI*z)); 306 } else { 307 x = -x; 308 lg = large_gam(x); 309 lsine = log__D(x*sin(M_PI*z)); 310 } 311 lg.b += lsine.b - lpi_lo; 312 y = (-(lg.b + lsine.a) + lpi_hi) - lg.a; 313 z = -lg.a - y; z+= lpi_hi; z -= lsine.a; z -= lg.b; 314 y = exp__D(y, z); 315 if(sgn < 0) y = -y; 316 return (y); 317 } 318 y = 1-x; 319 if (1-y == x) 320 y = ngamma(y); 321 else /* 1-x is inexact */ 322 y = -x*ngamma(-x); 323 if (sgn < 0) y = -y; 324 return (M_PI / (y*sin(M_PI*z))); 325 } 326