projects
/
musl
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
make dynamic linker accept : or \n as path separator
[musl]
/
src
/
math
/
lgammaf_r.c
diff --git
a/src/math/lgammaf_r.c
b/src/math/lgammaf_r.c
index
9955b2f
..
976986a
100644
(file)
--- a/
src/math/lgammaf_r.c
+++ b/
src/math/lgammaf_r.c
@@
-17,8
+17,6
@@
static const float
two23= 8.3886080000e+06, /* 0x4b000000 */
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 */
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 */
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;
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 */
*/
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) {
} else {
if (ix >= 0x4b800000) {
- y =
zero
; /* y must be even */
+ y =
0.0f
; /* y must be even */
n = 0;
} else {
if (ix < 0x4b000000)
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:
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 3:
- case 4: y = __sindf(pi*(
one-
y)); break;
+ case 4: y = __sindf(pi*(
1.0f -
y)); break;
case 5:
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;
}
}
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;
{
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)
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;
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 */
}
if (hx < 0) {
if (ix >= 0x4b000000) /* |x| >= 2**23, must be -integer */
- return
one/zero
;
+ return
1.0f/0.0f
;
t = sin_pif(x);
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));
nadj = logf(pi/fabsf(t*x));
- if (t <
zero
)
+ if (t <
0.0f
)
*signgamp = -1;
x = -x;
}
*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) {
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) {
i = 0;
} else if (ix >= 0x3e6d3308) {
- y = x - (tc-
one
);
+ y = x - (tc-
1.0f
);
i = 1;
} else {
y = x;
i = 2;
}
} else {
i = 1;
} else {
y = x;
i = 2;
}
} else {
- r =
zero
;
+ r =
0.0f
;
if (ix >= 0x3fdda618) { /* [1.7316,2] */
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 {
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;
}
}
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;
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;
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)))));
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;
}
} 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))))));
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 =
hal
f*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.5
f*y+p/q;
+ z =
1.0f
; /* lgamma(1+s) = log(s) + lgamma(s) */
switch (i) {
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);
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)))));
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 */
} else /* 2**58 <= x <= inf */
- r = x*(logf(x)-
one
);
+ r = x*(logf(x)-
1.0f
);
if (hx < 0)
r = nadj - r;
return r;
}
if (hx < 0)
r = nadj - r;
return r;
}
+
+weak_alias(__lgammaf_r, lgammaf_r);