projects
/
musl
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
some assemblers don't like fistpq; use the alt. mnemonic fistpll
[musl]
/
src
/
math
/
jnf.c
diff --git
a/src/math/jnf.c
b/src/math/jnf.c
index
7db93ae
..
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 {
else {
- temp =
x*(float)0.5
;
+ temp =
0.5f * x
;
b = temp;
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 */
}
@@
-106,13
+101,13
@@
float jnf(int n, float x)
float q0,q1,h,tmp;
int32_t k,m;
float q0,q1,h,tmp;
int32_t k,m;
- w = (n+n)/
(float)
x;
- h =
(float)2.0/(float)
x;
+ w = (n+n)/x;
+ h =
2.0f/
x;
z = w+h;
q0 = w;
z = w+h;
q0 = w;
- q1 = w*z -
(float)1.0
;
+ q1 = w*z -
1.0f
;
k = 1;
k = 1;
- while (q1 <
(float)1.0e9
) {
+ while (q1 <
1.0e9f
) {
k += 1;
z += h;
tmp = z*q1 - q0;
k += 1;
z += h;
tmp = z*q1 - q0;
@@
-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,15
+128,15
@@
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));
tmp = tmp*logf(fabsf(v*tmp));
- if (tmp <
(float)8.8721679688e+01
) {
+ if (tmp <
88.721679688f
) {
for (i=n-1,di=(float)(i+i); i>0; i--) {
temp = b;
b *= di;
b = b/x - a;
a = temp;
for (i=n-1,di=(float)(i+i); i>0; i--) {
temp = b;
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 */
/* scale b to avoid spurious overflow */
- if (b >
(float)1e10
) {
+ if (b >
1e10f
) {
a /= b;
t /= b;
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);