projects
/
musl
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Merge remote branch 'nsz/master'
[musl]
/
src
/
math
/
jnf.c
diff --git
a/src/math/jnf.c
b/src/math/jnf.c
index
648db32
..
fd29110
100644
(file)
--- a/
src/math/jnf.c
+++ b/
src/math/jnf.c
@@
-13,14
+13,9
@@
* ====================================================
*/
* ====================================================
*/
+#define _GNU_SOURCE
#include "libm.h"
#include "libm.h"
-static const float
-two = 2.0000000000e+00, /* 0x40000000 */
-one = 1.0000000000e+00; /* 0x3F800000 */
-
-static const float zero = 0.0000000000e+00;
-
float jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
float jnf(int n, float x)
{
int32_t i,hx,ix, sgn;
@@
-46,7
+41,7
@@
float jnf(int n, float x)
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
sgn = (n&1)&(hx>>31); /* even n -- 0, odd n -- sign(x) */
x = fabsf(x);
if (ix == 0 || ix >= 0x7f800000) /* if x is 0 or inf */
- b =
zero
;
+ b =
0.0f
;
else if((float)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x);
else if((float)n <= x) {
/* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
a = j0f(x);
@@
-62,11
+57,11
@@
float jnf(int n, float x)
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
* J(n,x) = 1/n!*(x/2)^n - ...
*/
if (n > 33) /* underflow */
- b =
zero
;
+ b =
0.0f
;
else {
temp = 0.5f * x;
b = temp;
else {
temp = 0.5f * x;
b = temp;
- for (a=
one
,i=2; i<=n; i++) {
+ for (a=
1.0f
,i=2; i<=n; i++) {
a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
a *= (float)i; /* a = n! */
b *= temp; /* b = (x/2)^n */
}
@@
-120,10
+115,10
@@
float jnf(int n, float x)
q1 = tmp;
}
m = n+n;
q1 = tmp;
}
m = n+n;
- for (t=
zero
, i = 2*(n+k); i>=m; i -= 2)
- t =
one
/(i/x-t);
+ for (t=
0.0f
, i = 2*(n+k); i>=m; i -= 2)
+ t =
1.0f
/(i/x-t);
a = t;
a = t;
- b =
one
;
+ b =
1.0f
;
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
/* estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
* Hence, if n*(log(2n/x)) > ...
* single 8.8722839355e+01
@@
-133,7
+128,7
@@
float jnf(int n, float x)
* likely underflow to zero
*/
tmp = n;
* likely underflow to zero
*/
tmp = n;
- v =
two
/x;
+ v =
2.0f
/x;
tmp = tmp*logf(fabsf(v*tmp));
if (tmp < 88.721679688f) {
for (i=n-1,di=(float)(i+i); i>0; i--) {
tmp = tmp*logf(fabsf(v*tmp));
if (tmp < 88.721679688f) {
for (i=n-1,di=(float)(i+i); i>0; i--) {
@@
-141,7
+136,7
@@
float jnf(int n, float x)
b *= di;
b = b/x - a;
a = temp;
b *= di;
b = b/x - a;
a = temp;
- di -=
two
;
+ di -=
2.0f
;
}
} else {
for (i=n-1,di=(float)(i+i); i>0; i--){
}
} else {
for (i=n-1,di=(float)(i+i); i>0; i--){
@@
-149,12
+144,12
@@
float jnf(int n, float x)
b *= di;
b = b/x - a;
a = temp;
b *= di;
b = b/x - a;
a = temp;
- di -=
two
;
+ di -=
2.0f
;
/* scale b to avoid spurious overflow */
if (b > 1e10f) {
a /= b;
t /= b;
/* scale b to avoid spurious overflow */
if (b > 1e10f) {
a /= b;
t /= b;
- b =
one
;
+ b =
1.0f
;
}
}
}
}
}
}
@@
-182,9
+177,9
@@
float ynf(int n, float x)
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 (hx < 0)
if (hx < 0)
- return
zero/zero
;
+ return
0.0f/0.0f
;
sign = 1;
if (n < 0) {
n = -n;
sign = 1;
if (n < 0) {
n = -n;
@@
-195,7
+190,7
@@
float ynf(int n, float x)
if (n == 1)
return sign*y1f(x);
if (ix == 0x7f800000)
if (n == 1)
return sign*y1f(x);
if (ix == 0x7f800000)
- return
zero
;
+ return
0.0f
;
a = y0f(x);
b = y1f(x);
a = y0f(x);
b = y1f(x);