X-Git-Url: http://nsz.repo.hu/git/?p=musl;a=blobdiff_plain;f=src%2Fmath%2Flgammaf_r.c;h=976986aa4aa13e2ec8d1f712787594e02b231457;hp=9955b2f93f1d07a2adb2bcc1aa0bdf0fa0825f32;hb=0cbb65479147ecdaa664e88cc2a5a925f3de502f;hpb=b69f695acedd4ce2798ef9ea28d834ceccc789bd diff --git a/src/math/lgammaf_r.c b/src/math/lgammaf_r.c index 9955b2f9..976986aa 100644 --- a/src/math/lgammaf_r.c +++ b/src/math/lgammaf_r.c @@ -17,8 +17,6 @@ static const float two23= 8.3886080000e+06, /* 0x4b000000 */ -half= 5.0000000000e-01, /* 0x3f000000 */ -one = 1.0000000000e+00, /* 0x3f800000 */ pi = 3.1415927410e+00, /* 0x40490fdb */ a0 = 7.7215664089e-02, /* 0x3d9e233f */ a1 = 3.2246702909e-01, /* 0x3ea51a66 */ @@ -83,8 +81,6 @@ w4 = -5.9518753551e-04, /* 0xba1c065c */ w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ w6 = -1.6309292987e-03; /* 0xbad5c4e8 */ -static const float zero = 0.0000000000e+00; - static float sin_pif(float x) { float y,z; @@ -104,12 +100,12 @@ static float sin_pif(float x) */ z = floorf(y); if (z != y) { /* inexact anyway */ - y *= (float)0.5; - y = (float)2.0*(y - floorf(y)); /* y = |x| mod 2.0 */ - n = (int) (y*(float)4.0); + y *= 0.5f; + y = 2.0f*(y - floorf(y)); /* y = |x| mod 2.0 */ + n = (int)(y*4.0f); } else { if (ix >= 0x4b800000) { - y = zero; /* y must be even */ + y = 0.0f; /* y must be even */ n = 0; } else { if (ix < 0x4b000000) @@ -123,18 +119,18 @@ static float sin_pif(float x) switch (n) { case 0: y = __sindf(pi*y); break; case 1: - case 2: y = __cosdf(pi*((float)0.5-y)); break; + case 2: y = __cosdf(pi*(0.5f - y)); break; case 3: - case 4: y = __sindf(pi*(one-y)); break; + case 4: y = __sindf(pi*(1.0f - y)); break; case 5: - case 6: y = -__cosdf(pi*(y-(float)1.5)); break; - default: y = __sindf(pi*(y-(float)2.0)); break; + case 6: y = -__cosdf(pi*(y - 1.5f)); break; + default: y = __sindf(pi*(y - 2.0f)); break; } return -y; } -float lgammaf_r(float x, int *signgamp) +float __lgammaf_r(float x, int *signgamp) { float t,y,z,nadj,p,p1,p2,p3,q,r,w; int32_t hx; @@ -148,7 +144,7 @@ float lgammaf_r(float x, int *signgamp) if (ix >= 0x7f800000) return x*x; if (ix == 0) - return one/zero; + return 1.0f/0.0f; if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ if (hx < 0) { *signgamp = -1; @@ -158,12 +154,12 @@ float lgammaf_r(float x, int *signgamp) } if (hx < 0) { if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */ - return one/zero; + return 1.0f/0.0f; t = sin_pif(x); - if (t == zero) /* -integer */ - return one/zero; + if (t == 0.0f) /* -integer */ + return 1.0f/0.0f; nadj = logf(pi/fabsf(t*x)); - if (t < zero) + if (t < 0.0f) *signgamp = -1; x = -x; } @@ -176,25 +172,25 @@ float lgammaf_r(float x, int *signgamp) if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ r = -logf(x); if (ix >= 0x3f3b4a20) { - y = one - x; + y = 1.0f - x; i = 0; } else if (ix >= 0x3e6d3308) { - y = x - (tc-one); + y = x - (tc-1.0f); i = 1; } else { y = x; i = 2; } } else { - r = zero; + r = 0.0f; if (ix >= 0x3fdda618) { /* [1.7316,2] */ - y = (float)2.0 - x; + y = 2.0f - x; i = 0; } else if (ix >= 0x3F9da620) { /* [1.23,1.73] */ y = x - tc; i = 1; } else { - y = x - one; + y = x - 1.0f; i = 2; } } @@ -204,7 +200,7 @@ float lgammaf_r(float x, int *signgamp) p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); p = y*p1+p2; - r += (p-(float)0.5*y); + r += p - 0.5f*y; break; case 1: z = y*y; @@ -217,34 +213,36 @@ float lgammaf_r(float x, int *signgamp) break; case 2: p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); - p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); - r += (-(float)0.5*y + p1/p2); + p2 = 1.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); + r += -0.5f*y + p1/p2; } } else if (ix < 0x41000000) { /* x < 8.0 */ i = (int)x; - y = x-(float)i; + y = x - (float)i; p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); - q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); - r = half*y+p/q; - z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ + q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); + r = 0.5f*y+p/q; + z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */ switch (i) { - case 7: z *= y + (float)6.0; /* FALLTHRU */ - case 6: z *= y + (float)5.0; /* FALLTHRU */ - case 5: z *= y + (float)4.0; /* FALLTHRU */ - case 4: z *= y + (float)3.0; /* FALLTHRU */ - case 3: z *= y + (float)2.0; /* FALLTHRU */ + case 7: z *= y + 6.0f; /* FALLTHRU */ + case 6: z *= y + 5.0f; /* FALLTHRU */ + case 5: z *= y + 4.0f; /* FALLTHRU */ + case 4: z *= y + 3.0f; /* FALLTHRU */ + case 3: z *= y + 2.0f; /* FALLTHRU */ r += logf(z); break; } } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */ t = logf(x); - z = one/x; + z = 1.0f/x; y = z*z; w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); - r = (x-half)*(t-one)+w; + r = (x-0.5f)*(t-1.0f)+w; } else /* 2**58 <= x <= inf */ - r = x*(logf(x)-one); + r = x*(logf(x)-1.0f); if (hx < 0) r = nadj - r; return r; } + +weak_alias(__lgammaf_r, lgammaf_r);