| Index: fusl/src/math/tgammal.c
|
| diff --git a/fusl/src/math/tgammal.c b/fusl/src/math/tgammal.c
|
| index 5336c5b194e86c57797562191ad39cd4af7fafe4..970892920810a384a67a27762b5375dac6ddd45d 100644
|
| --- a/fusl/src/math/tgammal.c
|
| +++ b/fusl/src/math/tgammal.c
|
| @@ -51,9 +51,8 @@
|
| #include "libm.h"
|
|
|
| #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
| -long double tgammal(long double x)
|
| -{
|
| - return tgamma(x);
|
| +long double tgammal(long double x) {
|
| + return tgamma(x);
|
| }
|
| #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
|
| /*
|
| @@ -65,25 +64,17 @@ Peak error = 1.83e-20
|
| Relative error spread = 8.4e-23
|
| */
|
| static const long double P[8] = {
|
| - 4.212760487471622013093E-5L,
|
| - 4.542931960608009155600E-4L,
|
| - 4.092666828394035500949E-3L,
|
| - 2.385363243461108252554E-2L,
|
| - 1.113062816019361559013E-1L,
|
| - 3.629515436640239168939E-1L,
|
| - 8.378004301573126728826E-1L,
|
| - 1.000000000000000000009E0L,
|
| + 4.212760487471622013093E-5L, 4.542931960608009155600E-4L,
|
| + 4.092666828394035500949E-3L, 2.385363243461108252554E-2L,
|
| + 1.113062816019361559013E-1L, 3.629515436640239168939E-1L,
|
| + 8.378004301573126728826E-1L, 1.000000000000000000009E0L,
|
| };
|
| static const long double Q[9] = {
|
| --1.397148517476170440917E-5L,
|
| - 2.346584059160635244282E-4L,
|
| --1.237799246653152231188E-3L,
|
| --7.955933682494738320586E-4L,
|
| - 2.773706565840072979165E-2L,
|
| --4.633887671244534213831E-2L,
|
| --2.243510905670329164562E-1L,
|
| - 4.150160950588455434583E-1L,
|
| - 9.999999999999999999908E-1L,
|
| + -1.397148517476170440917E-5L, 2.346584059160635244282E-4L,
|
| + -1.237799246653152231188E-3L, -7.955933682494738320586E-4L,
|
| + 2.773706565840072979165E-2L, -4.633887671244534213831E-2L,
|
| + -2.243510905670329164562E-1L, 4.150160950588455434583E-1L,
|
| + 9.999999999999999999908E-1L,
|
| };
|
|
|
| /*
|
| @@ -122,15 +113,11 @@ Peak error = 9.44e-21
|
| Relative error spread = 8.8e-4
|
| */
|
| static const long double STIR[9] = {
|
| - 7.147391378143610789273E-4L,
|
| --2.363848809501759061727E-5L,
|
| --5.950237554056330156018E-4L,
|
| - 6.989332260623193171870E-5L,
|
| - 7.840334842744753003862E-4L,
|
| --2.294719747873185405699E-4L,
|
| --2.681327161876304418288E-3L,
|
| - 3.472222222230075327854E-3L,
|
| - 8.333333333333331800504E-2L,
|
| + 7.147391378143610789273E-4L, -2.363848809501759061727E-5L,
|
| + -5.950237554056330156018E-4L, 6.989332260623193171870E-5L,
|
| + 7.840334842744753003862E-4L, -2.294719747873185405699E-4L,
|
| + -2.681327161876304418288E-3L, 3.472222222230075327854E-3L,
|
| + 8.333333333333331800504E-2L,
|
| };
|
|
|
| #define MAXSTIR 1024.0L
|
| @@ -142,15 +129,11 @@ static const long double SQTPI = 2.50662827463100050242E0L;
|
| * Peak relative error 4.2e-23
|
| */
|
| static const long double S[9] = {
|
| --1.193945051381510095614E-3L,
|
| - 7.220599478036909672331E-3L,
|
| --9.622023360406271645744E-3L,
|
| --4.219773360705915470089E-2L,
|
| - 1.665386113720805206758E-1L,
|
| --4.200263503403344054473E-2L,
|
| --6.558780715202540684668E-1L,
|
| - 5.772156649015328608253E-1L,
|
| - 1.000000000000000000000E0L,
|
| + -1.193945051381510095614E-3L, 7.220599478036909672331E-3L,
|
| + -9.622023360406271645744E-3L, -4.219773360705915470089E-2L,
|
| + 1.665386113720805206758E-1L, -4.200263503403344054473E-2L,
|
| + -6.558780715202540684668E-1L, 5.772156649015328608253E-1L,
|
| + 1.000000000000000000000E0L,
|
| };
|
|
|
| /* 1/tgamma(-x) = z P(z)
|
| @@ -160,122 +143,118 @@ static const long double S[9] = {
|
| * Relative error spread = 2.5e-24
|
| */
|
| static const long double SN[9] = {
|
| - 1.133374167243894382010E-3L,
|
| - 7.220837261893170325704E-3L,
|
| - 9.621911155035976733706E-3L,
|
| --4.219773343731191721664E-2L,
|
| --1.665386113944413519335E-1L,
|
| --4.200263503402112910504E-2L,
|
| - 6.558780715202536547116E-1L,
|
| - 5.772156649015328608727E-1L,
|
| --1.000000000000000000000E0L,
|
| + 1.133374167243894382010E-3L, 7.220837261893170325704E-3L,
|
| + 9.621911155035976733706E-3L, -4.219773343731191721664E-2L,
|
| + -1.665386113944413519335E-1L, -4.200263503402112910504E-2L,
|
| + 6.558780715202536547116E-1L, 5.772156649015328608727E-1L,
|
| + -1.000000000000000000000E0L,
|
| };
|
|
|
| static const long double PIL = 3.1415926535897932384626L;
|
|
|
| /* Gamma function computed by Stirling's formula.
|
| */
|
| -static long double stirf(long double x)
|
| -{
|
| - long double y, w, v;
|
| +static long double stirf(long double x) {
|
| + long double y, w, v;
|
|
|
| - w = 1.0/x;
|
| - /* For large x, use rational coefficients from the analytical expansion. */
|
| - if (x > 1024.0)
|
| - w = (((((6.97281375836585777429E-5L * w
|
| - + 7.84039221720066627474E-4L) * w
|
| - - 2.29472093621399176955E-4L) * w
|
| - - 2.68132716049382716049E-3L) * w
|
| - + 3.47222222222222222222E-3L) * w
|
| - + 8.33333333333333333333E-2L) * w
|
| - + 1.0;
|
| - else
|
| - w = 1.0 + w * __polevll(w, STIR, 8);
|
| - y = expl(x);
|
| - if (x > MAXSTIR) { /* Avoid overflow in pow() */
|
| - v = powl(x, 0.5L * x - 0.25L);
|
| - y = v * (v / y);
|
| - } else {
|
| - y = powl(x, x - 0.5L) / y;
|
| - }
|
| - y = SQTPI * y * w;
|
| - return y;
|
| + w = 1.0 / x;
|
| + /* For large x, use rational coefficients from the analytical expansion. */
|
| + if (x > 1024.0)
|
| + w = (((((6.97281375836585777429E-5L * w + 7.84039221720066627474E-4L) * w -
|
| + 2.29472093621399176955E-4L) *
|
| + w -
|
| + 2.68132716049382716049E-3L) *
|
| + w +
|
| + 3.47222222222222222222E-3L) *
|
| + w +
|
| + 8.33333333333333333333E-2L) *
|
| + w +
|
| + 1.0;
|
| + else
|
| + w = 1.0 + w * __polevll(w, STIR, 8);
|
| + y = expl(x);
|
| + if (x > MAXSTIR) { /* Avoid overflow in pow() */
|
| + v = powl(x, 0.5L * x - 0.25L);
|
| + y = v * (v / y);
|
| + } else {
|
| + y = powl(x, x - 0.5L) / y;
|
| + }
|
| + y = SQTPI * y * w;
|
| + return y;
|
| }
|
|
|
| -long double tgammal(long double x)
|
| -{
|
| - long double p, q, z;
|
| +long double tgammal(long double x) {
|
| + long double p, q, z;
|
|
|
| - if (!isfinite(x))
|
| - return x + INFINITY;
|
| + if (!isfinite(x))
|
| + return x + INFINITY;
|
|
|
| - q = fabsl(x);
|
| - if (q > 13.0) {
|
| - if (x < 0.0) {
|
| - p = floorl(q);
|
| - z = q - p;
|
| - 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;
|
| - }
|
| - if (0.5 * p == floorl(q * 0.5))
|
| - z = -z;
|
| - } else if (x > MAXGAML) {
|
| - z = x * 0x1p16383L;
|
| - } else {
|
| - z = stirf(x);
|
| - }
|
| - return z;
|
| - }
|
| + q = fabsl(x);
|
| + if (q > 13.0) {
|
| + if (x < 0.0) {
|
| + p = floorl(q);
|
| + z = q - p;
|
| + 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;
|
| + }
|
| + if (0.5 * p == floorl(q * 0.5))
|
| + z = -z;
|
| + } else if (x > MAXGAML) {
|
| + z = x * 0x1p16383L;
|
| + } else {
|
| + z = stirf(x);
|
| + }
|
| + return z;
|
| + }
|
|
|
| - z = 1.0;
|
| - while (x >= 3.0) {
|
| - x -= 1.0;
|
| - z *= x;
|
| - }
|
| - while (x < -0.03125L) {
|
| - z /= x;
|
| - x += 1.0;
|
| - }
|
| - if (x <= 0.03125L)
|
| - goto small;
|
| - while (x < 2.0) {
|
| - z /= x;
|
| - x += 1.0;
|
| - }
|
| - if (x == 2.0)
|
| - return z;
|
| + z = 1.0;
|
| + while (x >= 3.0) {
|
| + x -= 1.0;
|
| + z *= x;
|
| + }
|
| + while (x < -0.03125L) {
|
| + z /= x;
|
| + x += 1.0;
|
| + }
|
| + if (x <= 0.03125L)
|
| + goto small;
|
| + while (x < 2.0) {
|
| + z /= x;
|
| + x += 1.0;
|
| + }
|
| + if (x == 2.0)
|
| + return z;
|
|
|
| - x -= 2.0;
|
| - p = __polevll(x, P, 7);
|
| - q = __polevll(x, Q, 8);
|
| - z = z * p / q;
|
| - return z;
|
| + x -= 2.0;
|
| + p = __polevll(x, P, 7);
|
| + q = __polevll(x, Q, 8);
|
| + z = z * p / q;
|
| + return z;
|
|
|
| small:
|
| - /* 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));
|
| - } else
|
| - q = z / (x * __polevll(x, S, 8));
|
| - return q;
|
| + /* 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));
|
| + } else
|
| + q = z / (x * __polevll(x, S, 8));
|
| + return q;
|
| }
|
| #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
|
| // TODO: broken implementation to make things compile
|
| -long double tgammal(long double x)
|
| -{
|
| - return tgamma(x);
|
| +long double tgammal(long double x) {
|
| + return tgamma(x);
|
| }
|
| #endif
|
|
|