|
| 1 | +#= |
| 2 | +/* @(#)e_lgamma_r.c 1.3 95/01/18 */ |
| 3 | +/* |
| 4 | +* ==================================================== |
| 5 | +* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
| 6 | +* |
| 7 | +* Developed at SunSoft, a Sun Microsystems, Inc. business. |
| 8 | +* Permission to use, copy, modify, and distribute this |
| 9 | +* software is freely granted, provided that this notice |
| 10 | +* is preserved. |
| 11 | +* ==================================================== |
| 12 | +* |
| 13 | +*/ |
| 14 | +=# |
| 15 | + |
| 16 | +#= |
| 17 | +
|
| 18 | +/* __ieee754_lgamma_r(x, signgamp) |
| 19 | + * Reentrant version of the logarithm of the Gamma function |
| 20 | + * with user provide pointer for the sign of Gamma(x). |
| 21 | + * |
| 22 | + * Method: |
| 23 | + * 1. Argument Reduction for 0 < x <= 8 |
| 24 | + * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may |
| 25 | + * reduce x to a number in [1.5,2.5] by |
| 26 | + * lgamma(1+s) = log(s) + lgamma(s) |
| 27 | + * for example, |
| 28 | + * lgamma(7.3) = log(6.3) + lgamma(6.3) |
| 29 | + * = log(6.3*5.3) + lgamma(5.3) |
| 30 | + * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) |
| 31 | + * 2. Polynomial approximation of lgamma around its |
| 32 | + * minimun ymin=1.461632144968362245 to maintain monotonicity. |
| 33 | + * On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use |
| 34 | + * Let z = x-ymin; |
| 35 | + * lgamma(x) = -1.214862905358496078218 + z^2*poly(z) |
| 36 | + * where |
| 37 | + * poly(z) is a 14 degree polynomial. |
| 38 | + * 2. Rational approximation in the primary interval [2,3] |
| 39 | + * We use the following approximation: |
| 40 | + * s = x-2.0; |
| 41 | + * lgamma(x) = 0.5*s + s*P(s)/Q(s) |
| 42 | + * with accuracy |
| 43 | + * |P/Q - (lgamma(x)-0.5s)| < 2**-61.71 |
| 44 | + * Our algorithms are based on the following observation |
| 45 | + * |
| 46 | + * zeta(2)-1 2 zeta(3)-1 3 |
| 47 | + * lgamma(2+s) = s*(1-Euler) + --------- * s - --------- * s + ... |
| 48 | + * 2 3 |
| 49 | + * |
| 50 | + * where Euler = 0.5771... is the Euler constant, which is very |
| 51 | + * close to 0.5. |
| 52 | + * |
| 53 | + * 3. For x>=8, we have |
| 54 | + * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... |
| 55 | + * (better formula: |
| 56 | + * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) |
| 57 | + * Let z = 1/x, then we approximation |
| 58 | + * f(z) = lgamma(x) - (x-0.5)(log(x)-1) |
| 59 | + * by |
| 60 | + * 3 5 11 |
| 61 | + * w = w0 + w1*z + w2*z + w3*z + ... + w6*z |
| 62 | + * where |
| 63 | + * |w - f(z)| < 2**-58.74 |
| 64 | + * |
| 65 | + * 4. For negative x, since (G is gamma function) |
| 66 | + * -x*G(-x)*G(x) = pi/sin(pi*x), |
| 67 | + * we have |
| 68 | + * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) |
| 69 | + * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 |
| 70 | + * Hence, for x<0, signgam = sign(sin(pi*x)) and |
| 71 | + * lgamma(x) = log(|Gamma(x)|) |
| 72 | + * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); |
| 73 | + * Note: one should avoid compute pi*(-x) directly in the |
| 74 | + * computation of sin(pi*(-x)). |
| 75 | + * |
| 76 | + * 5. Special Cases |
| 77 | + * lgamma(2+s) ~ s*(1-Euler) for tiny s |
| 78 | + * lgamma(1) = lgamma(2) = 0 |
| 79 | + * lgamma(x) ~ -log(|x|) for tiny x |
| 80 | + * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero |
| 81 | + * lgamma(inf) = inf |
| 82 | + * lgamma(-inf) = inf (bug for bug compatible with C99!?) |
| 83 | + * |
| 84 | + */ |
| 85 | +
|
| 86 | +=# |
| 87 | + |
| 88 | +# Matches OpenLibm behavior (except commented out |x|≥2^52 early exit) |
| 89 | +function _logabsgamma(x::Float64) |
| 90 | + ux = reinterpret(UInt64, x) |
| 91 | + hx = ux >>> 32 % Int32 |
| 92 | + lx = ux % UInt32 |
| 93 | + |
| 94 | + #= purge off +-inf, NaN, +-0, tiny and negative arguments =# |
| 95 | + signgam = 1 |
| 96 | + isneg = hx < Int32(0) |
| 97 | + ix = hx & 0x7fffffff |
| 98 | + ix ≥ 0x7ff00000 && return x * x, signgam |
| 99 | + ix | lx == 0x00000000 && return Inf, signgam |
| 100 | + if ix < 0x3b900000 #= |x|<2**-70, return -log(|x|) =# |
| 101 | + if isneg |
| 102 | + signgam = -1 |
| 103 | + return -log(-x), signgam |
| 104 | + else |
| 105 | + return -log(x), signgam |
| 106 | + end |
| 107 | + end |
| 108 | + if isneg |
| 109 | + # ix ≥ 0x43300000 && return Inf, signgam #= |x|>=2**52, must be -integer =# |
| 110 | + t = sinpi(x) |
| 111 | + iszero(t) && return Inf, signgam #= -integer =# |
| 112 | + nadj = logπ - log(abs(t * x)) |
| 113 | + if t < 0.0; signgam = -1; end |
| 114 | + x = -x |
| 115 | + end |
| 116 | + if ix < 0x40000000 #= x < 2.0 =# |
| 117 | + i = round(x, RoundToZero) |
| 118 | + f = x - i |
| 119 | + if f == 0.0 #= purge off 1; 2 handled by x < 8.0 branch =# |
| 120 | + return 0.0, signgam |
| 121 | + elseif i == 1.0 |
| 122 | + r = 0.0 |
| 123 | + c = 1.0 |
| 124 | + else |
| 125 | + r = -log(x) |
| 126 | + c = 0.0 |
| 127 | + end |
| 128 | + if f ≥ 0.7315998077392578 |
| 129 | + y = 1.0 + c - x |
| 130 | + z = y * y |
| 131 | + p1 = @evalpoly(z, 7.72156649015328655494e-02, 6.73523010531292681824e-02, 7.38555086081402883957e-03, 1.19270763183362067845e-03, 2.20862790713908385557e-04, 2.52144565451257326939e-05) |
| 132 | + p2 = z * @evalpoly(z, 3.22467033424113591611e-01, 2.05808084325167332806e-02, 2.89051383673415629091e-03, 5.10069792153511336608e-04, 1.08011567247583939954e-04, 4.48640949618915160150e-05) |
| 133 | + p = muladd(p1, y, p2) |
| 134 | + r += muladd(y, -0.5, p) |
| 135 | + elseif f ≥ 0.2316399812698364 # or, the lb? 0.2316322326660156 |
| 136 | + y = x - 0.46163214496836225 - c |
| 137 | + z = y * y |
| 138 | + w = z * y |
| 139 | + p1 = @evalpoly(w, 4.83836122723810047042e-01, -3.27885410759859649565e-02, 6.10053870246291332635e-03, -1.40346469989232843813e-03, 3.15632070903625950361e-04) |
| 140 | + p2 = @evalpoly(w, -1.47587722994593911752e-01, 1.79706750811820387126e-02, -3.68452016781138256760e-03, 8.81081882437654011382e-04, -3.12754168375120860518e-04) |
| 141 | + p3 = @evalpoly(w, 6.46249402391333854778e-02, -1.03142241298341437450e-02, 2.25964780900612472250e-03, -5.38595305356740546715e-04, 3.35529192635519073543e-04) |
| 142 | + p = muladd(z, p1, -muladd(w, -muladd(p3, y, p2), -3.63867699703950536541e-18)) |
| 143 | + r += p - 1.21486290535849611461e-1 |
| 144 | + else |
| 145 | + y = x - c |
| 146 | + p1 = y * @evalpoly(y, -7.72156649015328655494e-02, 6.32827064025093366517e-01, 1.45492250137234768737, 9.77717527963372745603e-01, 2.28963728064692451092e-01, 1.33810918536787660377e-02) |
| 147 | + p2 = @evalpoly(y, 1.0, 2.45597793713041134822, 2.12848976379893395361, 7.69285150456672783825e-01, 1.04222645593369134254e-01, 3.21709242282423911810e-03) |
| 148 | + r += muladd(y, -0.5, p1 / p2) |
| 149 | + end |
| 150 | + elseif ix < 0x40200000 #= x < 8.0 =# |
| 151 | + i = round(x, RoundToZero) |
| 152 | + y = x - i |
| 153 | + z = 1.0 |
| 154 | + p = 0.0 |
| 155 | + u = x |
| 156 | + while u ≥ 3.0 |
| 157 | + p -= 1.0 |
| 158 | + u = x + p |
| 159 | + z *= u |
| 160 | + end |
| 161 | + p = y * @evalpoly(y, -7.72156649015328655494e-2, 2.14982415960608852501e-1, 3.25778796408930981787e-1, 1.46350472652464452805e-1, 2.66422703033638609560e-2, 1.84028451407337715652e-3, 3.19475326584100867617e-5) |
| 162 | + q = @evalpoly(y, 1.0, 1.39200533467621045958, 7.21935547567138069525e-1, 1.71933865632803078993e-1, 1.86459191715652901344e-2, 7.77942496381893596434e-4, 7.32668430744625636189e-6) |
| 163 | + r = log(z) + muladd(0.5, y, p / q) |
| 164 | + elseif ix < 0x43900000 #= 8.0 ≤ x < 2^58 =# |
| 165 | + z = 1.0 / x |
| 166 | + y = z * z |
| 167 | + w = muladd(z, @evalpoly(y, 8.33333333333329678849e-2, -2.77777777728775536470e-3, 7.93650558643019558500e-4, -5.95187557450339963135e-4, 8.36339918996282139126e-4, -1.63092934096575273989e-3), 4.18938533204672725052e-1) |
| 168 | + r = muladd(x - 0.5, log(x) - 1.0, w) |
| 169 | + else #= 2^58 ≤ x ≤ Inf =# |
| 170 | + r = muladd(x, log(x), -x) |
| 171 | + end |
| 172 | + if isneg |
| 173 | + r = nadj - r |
| 174 | + end |
| 175 | + return r, signgam |
| 176 | +end |
0 commit comments