X-Git-Url: http://nsz.repo.hu/git/?p=musl;a=blobdiff_plain;f=src%2Fmath%2Ftgammal.c;h=5c1a02a6399692ca60037f346e4e1be1f8163f19;hp=86782a967226f0f6fe50a61ae86cf92576eadf8e;hb=d8a7619e371ff0f226200f6316abb46dd1192f3d;hpb=a9014ac1b9bc813a21752cd40e0a73441860e6e6 diff --git a/src/math/tgammal.c b/src/math/tgammal.c index 86782a96..5c1a02a6 100644 --- a/src/math/tgammal.c +++ b/src/math/tgammal.c @@ -205,42 +205,36 @@ static long double stirf(long double x) long double tgammal(long double x) { long double p, q, z; - int i, sign; - if (isnan(x)) - return NAN; - if (x == INFINITY) - return INFINITY; - if (x == -INFINITY) - return x - x; + if (!isfinite(x)) + return x + INFINITY; + q = fabsl(x); if (q > 13.0) { - sign = 1; - if (q > MAXGAML) - goto goverf; if (x < 0.0) { p = floorl(q); - if (p == q) - return (x - x) / (x - x); - i = p; - if ((i & 1) == 0) - sign = -1; z = q - p; - if (z > 0.5L) { - p += 1.0; - z = q - p; - } - z = q * sinl(PIL * z); - z = fabsl(z) * stirf(q); - if (z <= PIL/LDBL_MAX) { -goverf: - return sign * INFINITY; + if (z == 0) + return 0 / z; + if (q > MAXGAML) { + z = 0; + } else { + if (z > 0.5) { + p += 1.0; + z = q - p; + } + z = q * sinl(PIL * z); + z = fabsl(z) * stirf(q); + z = PIL/z; } - z = PIL/z; + if (0.5 * p == floorl(q * 0.5)) + z = -z; + } else if (x > MAXGAML) { + z = x * 0x1p16383L; } else { z = stirf(x); } - return sign * z; + return z; } z = 1.0; @@ -268,8 +262,9 @@ goverf: return z; small: - if (x == 0.0) - return (x - x) / (x - x); + /* z==1 if x was originally +-0 */ + if (x == 0 && z != 1) + return x / x; if (x < 0.0) { x = -x; q = z / (x * __polevll(x, SN, 8));