summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorRich Felker <dalias@aerifal.cx>2012-11-15 08:37:58 -0500
committerRich Felker <dalias@aerifal.cx>2012-11-15 08:37:58 -0500
commit8bb181622222f2ee3462c8b021bcae4fcdbbd37a (patch)
tree166a1be834c8cd5a4c2521d41d7ee506bd8c8e51 /src
parent22781b4d8ecaf97ac52d2a501a2e7d3bc1920ce0 (diff)
parent68847ecd3a644d3ffd9be58603be1a3fa437dff0 (diff)
downloadmusl-8bb181622222f2ee3462c8b021bcae4fcdbbd37a.tar.gz
Merge remote-tracking branch 'nsz/math'
Diffstat (limited to 'src')
-rw-r--r--src/complex/__cexp.c2
-rw-r--r--src/complex/__cexpf.c2
-rw-r--r--src/complex/cacos.c2
-rw-r--r--src/complex/cacosf.c2
-rw-r--r--src/complex/cacosh.c2
-rw-r--r--src/complex/cacoshf.c2
-rw-r--r--src/complex/cacoshl.c2
-rw-r--r--src/complex/cacosl.c2
-rw-r--r--src/complex/casin.c4
-rw-r--r--src/complex/casinf.c4
-rw-r--r--src/complex/casinh.c4
-rw-r--r--src/complex/casinhf.c4
-rw-r--r--src/complex/casinhl.c4
-rw-r--r--src/complex/casinl.c4
-rw-r--r--src/complex/catanh.c4
-rw-r--r--src/complex/catanhf.c4
-rw-r--r--src/complex/catanhl.c4
-rw-r--r--src/complex/ccos.c2
-rw-r--r--src/complex/ccosf.c2
-rw-r--r--src/complex/ccosh.c26
-rw-r--r--src/complex/ccoshf.c26
-rw-r--r--src/complex/ccosl.c2
-rw-r--r--src/complex/cexp.c12
-rw-r--r--src/complex/cexpf.c12
-rw-r--r--src/complex/clog.c2
-rw-r--r--src/complex/clogf.c2
-rw-r--r--src/complex/clogl.c2
-rw-r--r--src/complex/conj.c2
-rw-r--r--src/complex/conjf.c2
-rw-r--r--src/complex/conjl.c2
-rw-r--r--src/complex/cproj.c2
-rw-r--r--src/complex/cprojf.c2
-rw-r--r--src/complex/cprojl.c2
-rw-r--r--src/complex/csin.c4
-rw-r--r--src/complex/csinf.c4
-rw-r--r--src/complex/csinh.c26
-rw-r--r--src/complex/csinhf.c26
-rw-r--r--src/complex/csinl.c4
-rw-r--r--src/complex/csqrt.c14
-rw-r--r--src/complex/csqrtf.c14
-rw-r--r--src/complex/ctan.c4
-rw-r--r--src/complex/ctanf.c4
-rw-r--r--src/complex/ctanh.c10
-rw-r--r--src/complex/ctanhf.c10
-rw-r--r--src/complex/ctanl.c4
-rw-r--r--src/fenv/fenv.c2
-rw-r--r--src/fenv/feupdateenv.c1
-rw-r--r--src/internal/libm.h31
-rw-r--r--src/math/__invtrigl.c42
-rw-r--r--src/math/__invtrigl.h63
-rw-r--r--src/math/acos.c1
-rw-r--r--src/math/acosl.c9
-rw-r--r--src/math/asinl.c4
-rw-r--r--src/math/atan2.c1
-rw-r--r--src/math/atan2l.c32
-rw-r--r--src/math/atanl.c45
-rw-r--r--src/math/exp10l.c4
-rw-r--r--src/math/expm1l.c20
-rw-r--r--src/math/fma.c2
-rw-r--r--src/math/fmaf.c1
-rw-r--r--src/math/fmal.c1
-rw-r--r--src/math/hypot.c9
-rw-r--r--src/math/hypotf.c6
-rw-r--r--src/math/ilogb.c8
-rw-r--r--src/math/ilogbf.c8
-rw-r--r--src/math/ilogbl.c8
-rw-r--r--src/math/llrintl.c1
-rw-r--r--src/math/log1pl.c9
-rw-r--r--src/math/log2l.c12
-rw-r--r--src/math/logb.c17
-rw-r--r--src/math/logbf.c14
-rw-r--r--src/math/logbl.c11
-rw-r--r--src/math/logl.c9
-rw-r--r--src/math/lrint.c1
-rw-r--r--src/math/lrintl.c1
-rw-r--r--src/math/modf.c6
-rw-r--r--src/math/modff.c6
-rw-r--r--src/math/nearbyint.c1
-rw-r--r--src/math/nearbyintf.c1
-rw-r--r--src/math/nearbyintl.c1
-rw-r--r--src/math/nextafter.c2
-rw-r--r--src/math/nextafterf.c2
-rw-r--r--src/math/nexttoward.c2
-rw-r--r--src/math/nexttowardf.c2
-rw-r--r--src/math/scalbn.c15
-rw-r--r--src/math/scalbnf.c15
86 files changed, 330 insertions, 372 deletions
diff --git a/src/complex/__cexp.c b/src/complex/__cexp.c
index f603e2be..05ac28c7 100644
--- a/src/complex/__cexp.c
+++ b/src/complex/__cexp.c
@@ -83,5 +83,5 @@ double complex __ldexp_cexp(double complex z, int expt)
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
- return cpack(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2);
+ return CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2);
}
diff --git a/src/complex/__cexpf.c b/src/complex/__cexpf.c
index 47168e8f..69b54045 100644
--- a/src/complex/__cexpf.c
+++ b/src/complex/__cexpf.c
@@ -63,6 +63,6 @@ float complex __ldexp_cexpf(float complex z, int expt)
half_expt = expt - half_expt;
SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
- return cpackf(cosf(y) * exp_x * scale1 * scale2,
+ return CMPLXF(cosf(y) * exp_x * scale1 * scale2,
sinf(y) * exp_x * scale1 * scale2);
}
diff --git a/src/complex/cacos.c b/src/complex/cacos.c
index 3aca0519..27c35636 100644
--- a/src/complex/cacos.c
+++ b/src/complex/cacos.c
@@ -7,5 +7,5 @@
double complex cacos(double complex z)
{
z = casin(z);
- return cpack(M_PI_2 - creal(z), -cimag(z));
+ return CMPLX(M_PI_2 - creal(z), -cimag(z));
}
diff --git a/src/complex/cacosf.c b/src/complex/cacosf.c
index 563766e7..11852659 100644
--- a/src/complex/cacosf.c
+++ b/src/complex/cacosf.c
@@ -5,5 +5,5 @@
float complex cacosf(float complex z)
{
z = casinf(z);
- return cpackf((float)M_PI_2 - crealf(z), -cimagf(z));
+ return CMPLXF((float)M_PI_2 - crealf(z), -cimagf(z));
}
diff --git a/src/complex/cacosh.c b/src/complex/cacosh.c
index c2dfc1ba..8c68cb01 100644
--- a/src/complex/cacosh.c
+++ b/src/complex/cacosh.c
@@ -5,5 +5,5 @@
double complex cacosh(double complex z)
{
z = cacos(z);
- return cpack(-cimag(z), creal(z));
+ return CMPLX(-cimag(z), creal(z));
}
diff --git a/src/complex/cacoshf.c b/src/complex/cacoshf.c
index 37ff8800..ade01c09 100644
--- a/src/complex/cacoshf.c
+++ b/src/complex/cacoshf.c
@@ -3,5 +3,5 @@
float complex cacoshf(float complex z)
{
z = cacosf(z);
- return cpackf(-cimagf(z), crealf(z));
+ return CMPLXF(-cimagf(z), crealf(z));
}
diff --git a/src/complex/cacoshl.c b/src/complex/cacoshl.c
index 2a04e27b..65342557 100644
--- a/src/complex/cacoshl.c
+++ b/src/complex/cacoshl.c
@@ -9,6 +9,6 @@ long double complex cacoshl(long double complex z)
long double complex cacoshl(long double complex z)
{
z = cacosl(z);
- return cpackl(-cimagl(z), creall(z));
+ return CMPLXL(-cimagl(z), creall(z));
}
#endif
diff --git a/src/complex/cacosl.c b/src/complex/cacosl.c
index 5992e056..7fd4a2f6 100644
--- a/src/complex/cacosl.c
+++ b/src/complex/cacosl.c
@@ -11,6 +11,6 @@ long double complex cacosl(long double complex z)
long double complex cacosl(long double complex z)
{
z = casinl(z);
- return cpackl(PI_2 - creall(z), -cimagl(z));
+ return CMPLXL(PI_2 - creall(z), -cimagl(z));
}
#endif
diff --git a/src/complex/casin.c b/src/complex/casin.c
index 79aff278..dfdda988 100644
--- a/src/complex/casin.c
+++ b/src/complex/casin.c
@@ -11,6 +11,6 @@ double complex casin(double complex z)
x = creal(z);
y = cimag(z);
- w = cpack(1.0 - (x - y)*(x + y), -2.0*x*y);
- return clog(cpack(-y, x) + csqrt(w));
+ w = CMPLX(1.0 - (x - y)*(x + y), -2.0*x*y);
+ return clog(CMPLX(-y, x) + csqrt(w));
}
diff --git a/src/complex/casinf.c b/src/complex/casinf.c
index cb9863f6..93f0e335 100644
--- a/src/complex/casinf.c
+++ b/src/complex/casinf.c
@@ -9,6 +9,6 @@ float complex casinf(float complex z)
x = crealf(z);
y = cimagf(z);
- w = cpackf(1.0 - (x - y)*(x + y), -2.0*x*y);
- return clogf(cpackf(-y, x) + csqrtf(w));
+ w = CMPLXF(1.0 - (x - y)*(x + y), -2.0*x*y);
+ return clogf(CMPLXF(-y, x) + csqrtf(w));
}
diff --git a/src/complex/casinh.c b/src/complex/casinh.c
index f2b3fef8..b57fe8c4 100644
--- a/src/complex/casinh.c
+++ b/src/complex/casinh.c
@@ -4,6 +4,6 @@
double complex casinh(double complex z)
{
- z = casin(cpack(-cimag(z), creal(z)));
- return cpack(cimag(z), -creal(z));
+ z = casin(CMPLX(-cimag(z), creal(z)));
+ return CMPLX(cimag(z), -creal(z));
}
diff --git a/src/complex/casinhf.c b/src/complex/casinhf.c
index ed4af643..a11bf902 100644
--- a/src/complex/casinhf.c
+++ b/src/complex/casinhf.c
@@ -2,6 +2,6 @@
float complex casinhf(float complex z)
{
- z = casinf(cpackf(-cimagf(z), crealf(z)));
- return cpackf(cimagf(z), -crealf(z));
+ z = casinf(CMPLXF(-cimagf(z), crealf(z)));
+ return CMPLXF(cimagf(z), -crealf(z));
}
diff --git a/src/complex/casinhl.c b/src/complex/casinhl.c
index e5d80cef..714f1893 100644
--- a/src/complex/casinhl.c
+++ b/src/complex/casinhl.c
@@ -8,7 +8,7 @@ long double complex casinhl(long double complex z)
#else
long double complex casinhl(long double complex z)
{
- z = casinl(cpackl(-cimagl(z), creall(z)));
- return cpackl(cimagl(z), -creall(z));
+ z = casinl(CMPLXL(-cimagl(z), creall(z)));
+ return CMPLXL(cimagl(z), -creall(z));
}
#endif
diff --git a/src/complex/casinl.c b/src/complex/casinl.c
index f9aa8ded..0916c60f 100644
--- a/src/complex/casinl.c
+++ b/src/complex/casinl.c
@@ -14,7 +14,7 @@ long double complex casinl(long double complex z)
x = creall(z);
y = cimagl(z);
- w = cpackl(1.0 - (x - y)*(x + y), -2.0*x*y);
- return clogl(cpackl(-y, x) + csqrtl(w));
+ w = CMPLXL(1.0 - (x - y)*(x + y), -2.0*x*y);
+ return clogl(CMPLXL(-y, x) + csqrtl(w));
}
#endif
diff --git a/src/complex/catanh.c b/src/complex/catanh.c
index b1628022..e248d9b9 100644
--- a/src/complex/catanh.c
+++ b/src/complex/catanh.c
@@ -4,6 +4,6 @@
double complex catanh(double complex z)
{
- z = catan(cpack(-cimag(z), creal(z)));
- return cpack(cimag(z), -creal(z));
+ z = catan(CMPLX(-cimag(z), creal(z)));
+ return CMPLX(cimag(z), -creal(z));
}
diff --git a/src/complex/catanhf.c b/src/complex/catanhf.c
index e1d1e648..4a5eb040 100644
--- a/src/complex/catanhf.c
+++ b/src/complex/catanhf.c
@@ -2,6 +2,6 @@
float complex catanhf(float complex z)
{
- z = catanf(cpackf(-cimagf(z), crealf(z)));
- return cpackf(cimagf(z), -crealf(z));
+ z = catanf(CMPLXF(-cimagf(z), crealf(z)));
+ return CMPLXF(cimagf(z), -crealf(z));
}
diff --git a/src/complex/catanhl.c b/src/complex/catanhl.c
index 0a9374a3..a5dd538e 100644
--- a/src/complex/catanhl.c
+++ b/src/complex/catanhl.c
@@ -8,7 +8,7 @@ long double complex catanhl(long double complex z)
#else
long double complex catanhl(long double complex z)
{
- z = catanl(cpackl(-cimagl(z), creall(z)));
- return cpackl(cimagl(z), -creall(z));
+ z = catanl(CMPLXL(-cimagl(z), creall(z)));
+ return CMPLXL(cimagl(z), -creall(z));
}
#endif
diff --git a/src/complex/ccos.c b/src/complex/ccos.c
index 5754c238..645aec29 100644
--- a/src/complex/ccos.c
+++ b/src/complex/ccos.c
@@ -4,5 +4,5 @@
double complex ccos(double complex z)
{
- return ccosh(cpack(-cimag(z), creal(z)));
+ return ccosh(CMPLX(-cimag(z), creal(z)));
}
diff --git a/src/complex/ccosf.c b/src/complex/ccosf.c
index 9b72c4f4..9a67241f 100644
--- a/src/complex/ccosf.c
+++ b/src/complex/ccosf.c
@@ -2,5 +2,5 @@
float complex ccosf(float complex z)
{
- return ccoshf(cpackf(-cimagf(z), crealf(z)));
+ return ccoshf(CMPLXF(-cimagf(z), crealf(z)));
}
diff --git a/src/complex/ccosh.c b/src/complex/ccosh.c
index 81f29432..401f3c60 100644
--- a/src/complex/ccosh.c
+++ b/src/complex/ccosh.c
@@ -55,23 +55,23 @@ double complex ccosh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
- return cpack(cosh(x), x * y);
+ return CMPLX(cosh(x), x * y);
if (ix < 0x40360000) /* small x: normal case */
- return cpack(cosh(x) * cos(y), sinh(x) * sin(y));
+ return CMPLX(cosh(x) * cos(y), sinh(x) * sin(y));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
- return cpack(h * cos(y), copysign(h, x) * sin(y));
+ return CMPLX(h * cos(y), copysign(h, x) * sin(y));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
- z = __ldexp_cexp(cpack(fabs(x), y), -1);
- return cpack(creal(z), cimag(z) * copysign(1, x));
+ z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
+ return CMPLX(creal(z), cimag(z) * copysign(1, x));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
- return cpack(h * h * cos(y), h * sin(y));
+ return CMPLX(h * h * cos(y), h * sin(y));
}
}
@@ -85,7 +85,7 @@ double complex ccosh(double complex z)
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
- return cpack(y - y, copysign(0, x * (y - y)));
+ return CMPLX(y - y, copysign(0, x * (y - y)));
/*
* cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
@@ -95,8 +95,8 @@ double complex ccosh(double complex z)
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
- return cpack(x * x, copysign(0, x) * y);
- return cpack(x * x, copysign(0, (x + x) * y));
+ return CMPLX(x * x, copysign(0, x) * y);
+ return CMPLX(x * x, copysign(0, (x + x) * y));
}
/*
@@ -108,7 +108,7 @@ double complex ccosh(double complex z)
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
- return cpack(y - y, x * (y - y));
+ return CMPLX(y - y, x * (y - y));
/*
* cosh(+-Inf + I NaN) = +Inf + I d(NaN).
@@ -121,8 +121,8 @@ double complex ccosh(double complex z)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000)
- return cpack(x * x, x * (y - y));
- return cpack((x * x) * cos(y), x * sin(y));
+ return CMPLX(x * x, x * (y - y));
+ return CMPLX((x * x) * cos(y), x * sin(y));
}
/*
@@ -136,5 +136,5 @@ double complex ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
- return cpack((x * x) * (y - y), (x + x) * (y - y));
+ return CMPLX((x * x) * (y - y), (x + x) * (y - y));
}
diff --git a/src/complex/ccoshf.c b/src/complex/ccoshf.c
index 683e77fa..90acfe05 100644
--- a/src/complex/ccoshf.c
+++ b/src/complex/ccoshf.c
@@ -48,43 +48,43 @@ float complex ccoshf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
- return cpackf(coshf(x), x * y);
+ return CMPLXF(coshf(x), x * y);
if (ix < 0x41100000) /* small x: normal case */
- return cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y));
+ return CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
- return cpackf(h * cosf(y), copysignf(h, x) * sinf(y));
+ return CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
- z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
- return cpackf(crealf(z), cimagf(z) * copysignf(1, x));
+ z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
+ return CMPLXF(crealf(z), cimagf(z) * copysignf(1, x));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
- return cpackf(h * h * cosf(y), h * sinf(y));
+ return CMPLXF(h * h * cosf(y), h * sinf(y));
}
}
if (ix == 0 && iy >= 0x7f800000)
- return cpackf(y - y, copysignf(0, x * (y - y)));
+ return CMPLXF(y - y, copysignf(0, x * (y - y)));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
- return cpackf(x * x, copysignf(0, x) * y);
- return cpackf(x * x, copysignf(0, (x + x) * y));
+ return CMPLXF(x * x, copysignf(0, x) * y);
+ return CMPLXF(x * x, copysignf(0, (x + x) * y));
}
if (ix < 0x7f800000 && iy >= 0x7f800000)
- return cpackf(y - y, x * (y - y));
+ return CMPLXF(y - y, x * (y - y));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000)
- return cpackf(x * x, x * (y - y));
- return cpackf((x * x) * cosf(y), x * sinf(y));
+ return CMPLXF(x * x, x * (y - y));
+ return CMPLXF((x * x) * cosf(y), x * sinf(y));
}
- return cpackf((x * x) * (y - y), (x + x) * (y - y));
+ return CMPLXF((x * x) * (y - y), (x + x) * (y - y));
}
diff --git a/src/complex/ccosl.c b/src/complex/ccosl.c
index e37825a9..d787047f 100644
--- a/src/complex/ccosl.c
+++ b/src/complex/ccosl.c
@@ -8,6 +8,6 @@ long double complex ccosl(long double complex z)
#else
long double complex ccosl(long double complex z)
{
- return ccoshl(cpackl(-cimagl(z), creall(z)));
+ return ccoshl(CMPLXL(-cimagl(z), creall(z)));
}
#endif
diff --git a/src/complex/cexp.c b/src/complex/cexp.c
index 3b8bb752..5118e00e 100644
--- a/src/complex/cexp.c
+++ b/src/complex/cexp.c
@@ -44,22 +44,22 @@ double complex cexp(double complex z)
/* cexp(x + I 0) = exp(x) + I 0 */
if ((hy | ly) == 0)
- return cpack(exp(x), y);
+ return CMPLX(exp(x), y);
EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if (((hx & 0x7fffffff) | lx) == 0)
- return cpack(cos(y), sin(y));
+ return CMPLX(cos(y), sin(y));
if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
- return cpack(y - y, y - y);
+ return CMPLX(y - y, y - y);
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
- return cpack(0.0, 0.0);
+ return CMPLX(0.0, 0.0);
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
- return cpack(x, y - y);
+ return CMPLX(x, y - y);
}
}
@@ -78,6 +78,6 @@ double complex cexp(double complex z)
* - x = NaN (spurious inexact exception from y)
*/
exp_x = exp(x);
- return cpack(exp_x * cos(y), exp_x * sin(y));
+ return CMPLX(exp_x * cos(y), exp_x * sin(y));
}
}
diff --git a/src/complex/cexpf.c b/src/complex/cexpf.c
index 0cf13a3d..1a09964c 100644
--- a/src/complex/cexpf.c
+++ b/src/complex/cexpf.c
@@ -44,22 +44,22 @@ float complex cexpf(float complex z)
/* cexp(x + I 0) = exp(x) + I 0 */
if (hy == 0)
- return cpackf(expf(x), y);
+ return CMPLXF(expf(x), y);
GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if ((hx & 0x7fffffff) == 0)
- return cpackf(cosf(y), sinf(y));
+ return CMPLXF(cosf(y), sinf(y));
if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
- return cpackf(y - y, y - y);
+ return CMPLXF(y - y, y - y);
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
- return cpackf(0.0, 0.0);
+ return CMPLXF(0.0, 0.0);
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
- return cpackf(x, y - y);
+ return CMPLXF(x, y - y);
}
}
@@ -78,6 +78,6 @@ float complex cexpf(float complex z)
* - x = NaN (spurious inexact exception from y)
*/
exp_x = expf(x);
- return cpackf(exp_x * cosf(y), exp_x * sinf(y));
+ return CMPLXF(exp_x * cosf(y), exp_x * sinf(y));
}
}
diff --git a/src/complex/clog.c b/src/complex/clog.c
index 6f10a115..12aae9c7 100644
--- a/src/complex/clog.c
+++ b/src/complex/clog.c
@@ -10,5 +10,5 @@ double complex clog(double complex z)
r = cabs(z);
phi = carg(z);
- return cpack(log(r), phi);
+ return CMPLX(log(r), phi);
}
diff --git a/src/complex/clogf.c b/src/complex/clogf.c
index f3aec54d..e9b32e60 100644
--- a/src/complex/clogf.c
+++ b/src/complex/clogf.c
@@ -8,5 +8,5 @@ float complex clogf(float complex z)
r = cabsf(z);
phi = cargf(z);
- return cpackf(logf(r), phi);
+ return CMPLXF(logf(r), phi);
}
diff --git a/src/complex/clogl.c b/src/complex/clogl.c
index 5b84ba59..18f16088 100644
--- a/src/complex/clogl.c
+++ b/src/complex/clogl.c
@@ -13,6 +13,6 @@ long double complex clogl(long double complex z)
r = cabsl(z);
phi = cargl(z);
- return cpackl(logl(r), phi);
+ return CMPLXL(logl(r), phi);
}
#endif
diff --git a/src/complex/conj.c b/src/complex/conj.c
index 4aceea7b..0b3f5f46 100644
--- a/src/complex/conj.c
+++ b/src/complex/conj.c
@@ -2,5 +2,5 @@
double complex conj(double complex z)
{
- return cpack(creal(z), -cimag(z));
+ return CMPLX(creal(z), -cimag(z));
}
diff --git a/src/complex/conjf.c b/src/complex/conjf.c
index 31556800..9af6b2c3 100644
--- a/src/complex/conjf.c
+++ b/src/complex/conjf.c
@@ -2,5 +2,5 @@
float complex conjf(float complex z)
{
- return cpackf(crealf(z), -cimagf(z));
+ return CMPLXF(crealf(z), -cimagf(z));
}
diff --git a/src/complex/conjl.c b/src/complex/conjl.c
index 01332262..67f11b9d 100644
--- a/src/complex/conjl.c
+++ b/src/complex/conjl.c
@@ -2,5 +2,5 @@
long double complex conjl(long double complex z)
{
- return cpackl(creall(z), -cimagl(z));
+ return CMPLXL(creall(z), -cimagl(z));
}
diff --git a/src/complex/cproj.c b/src/complex/cproj.c
index 1cf9bb94..15f358a1 100644
--- a/src/complex/cproj.c
+++ b/src/complex/cproj.c
@@ -3,6 +3,6 @@
double complex cproj(double complex z)
{
if (isinf(creal(z)) || isinf(cimag(z)))
- return cpack(INFINITY, copysign(0.0, creal(z)));
+ return CMPLX(INFINITY, copysign(0.0, creal(z)));
return z;
}
diff --git a/src/complex/cprojf.c b/src/complex/cprojf.c
index 71129743..653be5e8 100644
--- a/src/complex/cprojf.c
+++ b/src/complex/cprojf.c
@@ -3,6 +3,6 @@
float complex cprojf(float complex z)
{
if (isinf(crealf(z)) || isinf(cimagf(z)))
- return cpackf(INFINITY, copysignf(0.0, crealf(z)));
+ return CMPLXF(INFINITY, copysignf(0.0, crealf(z)));
return z;
}
diff --git a/src/complex/cprojl.c b/src/complex/cprojl.c
index 72e50cf5..6731aaa2 100644
--- a/src/complex/cprojl.c
+++ b/src/complex/cprojl.c
@@ -9,7 +9,7 @@ long double complex cprojl(long double complex z)
long double complex cprojl(long double complex z)
{
if (isinf(creall(z)) || isinf(cimagl(z)))
- return cpackl(INFINITY, copysignl(0.0, creall(z)));
+ return CMPLXL(INFINITY, copysignl(0.0, creall(z)));
return z;
}
#endif
diff --git a/src/complex/csin.c b/src/complex/csin.c
index 246a4595..ad8ae67a 100644
--- a/src/complex/csin.c
+++ b/src/complex/csin.c
@@ -4,6 +4,6 @@
double complex csin(double complex z)
{
- z = csinh(cpack(-cimag(z), creal(z)));
- return cpack(cimag(z), -creal(z));
+ z = csinh(CMPLX(-cimag(z), creal(z)));
+ return CMPLX(cimag(z), -creal(z));
}
diff --git a/src/complex/csinf.c b/src/complex/csinf.c
index 3aabe8f8..60b3cbaa 100644
--- a/src/complex/csinf.c
+++ b/src/complex/csinf.c
@@ -2,6 +2,6 @@
float complex csinf(float complex z)
{
- z = csinhf(cpackf(-cimagf(z), crealf(z)));
- return cpackf(cimagf(z), -crealf(z));
+ z = csinhf(CMPLXF(-cimagf(z), crealf(z)));
+ return CMPLXF(cimagf(z), -crealf(z));
}
diff --git a/src/complex/csinh.c b/src/complex/csinh.c
index fe16f06b..0f8035d1 100644
--- a/src/complex/csinh.c
+++ b/src/complex/csinh.c
@@ -55,23 +55,23 @@ double complex csinh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
- return cpack(sinh(x), y);
+ return CMPLX(sinh(x), y);
if (ix < 0x40360000) /* small x: normal case */
- return cpack(sinh(x) * cos(y), cosh(x) * sin(y));
+ return CMPLX(sinh(x) * cos(y), cosh(x) * sin(y));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
- return cpack(copysign(h, x) * cos(y), h * sin(y));
+ return CMPLX(copysign(h, x) * cos(y), h * sin(y));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
- z = __ldexp_cexp(cpack(fabs(x), y), -1);
- return cpack(creal(z) * copysign(1, x), cimag(z));
+ z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
+ return CMPLX(creal(z) * copysign(1, x), cimag(z));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
- return cpack(h * cos(y), h * h * sin(y));
+ return CMPLX(h * cos(y), h * h * sin(y));
}
}
@@ -85,7 +85,7 @@ double complex csinh(double complex z)
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
- return cpack(copysign(0, x * (y - y)), y - y);
+ return CMPLX(copysign(0, x * (y - y)), y - y);
/*
* sinh(+-Inf +- I 0) = +-Inf + I +-0.
@@ -94,8 +94,8 @@ double complex csinh(double complex z)
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
- return cpack(x, y);
- return cpack(x, copysign(0, y));
+ return CMPLX(x, y);
+ return CMPLX(x, copysign(0, y));
}
/*
@@ -107,7 +107,7 @@ double complex csinh(double complex z)
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
- return cpack(y - y, x * (y - y));
+ return CMPLX(y - y, x * (y - y));
/*
* sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
@@ -122,8 +122,8 @@ double complex csinh(double complex z)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000)
- return cpack(x * x, x * (y - y));
- return cpack(x * cos(y), INFINITY * sin(y));
+ return CMPLX(x * x, x * (y - y));
+ return CMPLX(x * cos(y), INFINITY * sin(y));
}
/*
@@ -137,5 +137,5 @@ double complex csinh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
- return cpack((x * x) * (y - y), (x + x) * (y - y));
+ return CMPLX((x * x) * (y - y), (x + x) * (y - y));
}
diff --git a/src/complex/csinhf.c b/src/complex/csinhf.c
index bbb116c2..49697f02 100644
--- a/src/complex/csinhf.c
+++ b/src/complex/csinhf.c
@@ -48,43 +48,43 @@ float complex csinhf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
- return cpackf(sinhf(x), y);
+ return CMPLXF(sinhf(x), y);
if (ix < 0x41100000) /* small x: normal case */
- return cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y));
+ return CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
- return cpackf(copysignf(h, x) * cosf(y), h * sinf(y));
+ return CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
- z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
- return cpackf(crealf(z) * copysignf(1, x), cimagf(z));
+ z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
+ return CMPLXF(crealf(z) * copysignf(1, x), cimagf(z));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
- return cpackf(h * cosf(y), h * h * sinf(y));
+ return CMPLXF(h * cosf(y), h * h * sinf(y));
}
}
if (ix == 0 && iy >= 0x7f800000)
- return cpackf(copysignf(0, x * (y - y)), y - y);
+ return CMPLXF(copysignf(0, x * (y - y)), y - y);
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
- return cpackf(x, y);
- return cpackf(x, copysignf(0, y));
+ return CMPLXF(x, y);
+ return CMPLXF(x, copysignf(0, y));
}
if (ix < 0x7f800000 && iy >= 0x7f800000)
- return cpackf(y - y, x * (y - y));
+ return CMPLXF(y - y, x * (y - y));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000)
- return cpackf(x * x, x * (y - y));
- return cpackf(x * cosf(y), INFINITY * sinf(y));
+ return CMPLXF(x * x, x * (y - y));
+ return CMPLXF(x * cosf(y), INFINITY * sinf(y));
}
- return cpackf((x * x) * (y - y), (x + x) * (y - y));
+ return CMPLXF((x * x) * (y - y), (x + x) * (y - y));
}
diff --git a/src/complex/csinl.c b/src/complex/csinl.c
index 4ad86745..4e9f86c3 100644
--- a/src/complex/csinl.c
+++ b/src/complex/csinl.c
@@ -8,7 +8,7 @@ long double complex csinl(long double complex z)
#else
long double complex csinl(long double complex z)
{
- z = csinhl(cpackl(-cimagl(z), creall(z)));
- return cpackl(cimagl(z), -creall(z));
+ z = csinhl(CMPLXL(-cimagl(z), creall(z)));
+ return CMPLXL(cimagl(z), -creall(z));
}
#endif
diff --git a/src/complex/csqrt.c b/src/complex/csqrt.c
index 21fb879d..8a2ba608 100644
--- a/src/complex/csqrt.c
+++ b/src/complex/csqrt.c
@@ -51,12 +51,12 @@ double complex csqrt(double complex z)
/* Handle special cases. */
if (z == 0)
- return cpack(0, b);
+ return CMPLX(0, b);
if (isinf(b))
- return cpack(INFINITY, b);
+ return CMPLX(INFINITY, b);
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
- return cpack(a, t); /* return NaN + NaN i */
+ return CMPLX(a, t); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@@ -66,9 +66,9 @@ double complex csqrt(double complex z)
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
- return cpack(fabs(b - b), copysign(a, b));
+ return CMPLX(fabs(b - b), copysign(a, b));
else
- return cpack(a, copysign(b - b, b));
+ return CMPLX(a, copysign(b - b, b));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@@ -87,10 +87,10 @@ double complex csqrt(double complex z)
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
- result = cpack(t, b / (2 * t));
+ result = CMPLX(t, b / (2 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
- result = cpack(fabs(b) / (2 * t), copysign(t, b));
+ result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
}
/* Rescale. */
diff --git a/src/complex/csqrtf.c b/src/complex/csqrtf.c
index 16487c23..ab5102f0 100644
--- a/src/complex/csqrtf.c
+++ b/src/complex/csqrtf.c
@@ -43,12 +43,12 @@ float complex csqrtf(float complex z)
/* Handle special cases. */
if (z == 0)
- return cpackf(0, b);
+ return CMPLXF(0, b);
if (isinf(b))
- return cpackf(INFINITY, b);
+ return CMPLXF(INFINITY, b);
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
- return cpackf(a, t); /* return NaN + NaN i */
+ return CMPLXF(a, t); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@@ -58,9 +58,9 @@ float complex csqrtf(float complex z)
* csqrtf(-inf + y i) = 0 + inf i
*/
if (signbit(a))
- return cpackf(fabsf(b - b), copysignf(a, b));
+ return CMPLXF(fabsf(b - b), copysignf(a, b));
else
- return cpackf(a, copysignf(b - b, b));
+ return CMPLXF(a, copysignf(b - b, b));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@@ -74,9 +74,9 @@ float complex csqrtf(float complex z)
*/
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
- return cpackf(t, b / (2.0 * t));
+ return CMPLXF(t, b / (2.0 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
- return cpackf(fabsf(b) / (2.0 * t), copysignf(t, b));
+ return CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b));
}
}
diff --git a/src/complex/ctan.c b/src/complex/ctan.c
index 4741a4df..c0926374 100644
--- a/src/complex/ctan.c
+++ b/src/complex/ctan.c
@@ -4,6 +4,6 @@
double complex ctan(double complex z)
{
- z = ctanh(cpack(-cimag(z), creal(z)));
- return cpack(cimag(z), -creal(z));
+ z = ctanh(CMPLX(-cimag(z), creal(z)));
+ return CMPLX(cimag(z), -creal(z));
}
diff --git a/src/complex/ctanf.c b/src/complex/ctanf.c
index 9bbeb051..009b1921 100644
--- a/src/complex/ctanf.c
+++ b/src/complex/ctanf.c
@@ -2,6 +2,6 @@
float complex ctanf(float complex z)
{
- z = ctanhf(cpackf(-cimagf(z), crealf(z)));
- return cpackf(cimagf(z), -crealf(z));
+ z = ctanhf(CMPLXF(-cimagf(z), crealf(z)));
+ return CMPLXF(cimagf(z), -crealf(z));
}
diff --git a/src/complex/ctanh.c b/src/complex/ctanh.c
index dd569fc3..0461050d 100644
--- a/src/complex/ctanh.c
+++ b/src/complex/ctanh.c
@@ -95,9 +95,9 @@ double complex ctanh(double complex z)
*/
if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */
- return cpack(x, (y == 0 ? y : x * y));
+ return CMPLX(x, (y == 0 ? y : x * y));
SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
- return cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)));
+ return CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)));
}
/*
@@ -105,7 +105,7 @@ double complex ctanh(double complex z)
* ctanh(x +- i Inf) = NaN + i NaN
*/
if (!isfinite(y))
- return cpack(y - y, y - y);
+ return CMPLX(y - y, y - y);
/*
* ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
@@ -114,7 +114,7 @@ double complex ctanh(double complex z)
*/
if (ix >= 0x40360000) { /* x >= 22 */
double exp_mx = exp(-fabs(x));
- return cpack(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
+ return CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
}
/* Kahan's algorithm */
@@ -123,5 +123,5 @@ double complex ctanh(double complex z)
s = sinh(x);
rho = sqrt(1 + s * s); /* = cosh(x) */
denom = 1 + beta * s * s;
- return cpack((beta * rho * s) / denom, t / denom);
+ return CMPLX((beta * rho * s) / denom, t / denom);
}
diff --git a/src/complex/ctanhf.c b/src/complex/ctanhf.c
index 7d746134..a7e1a5fc 100644
--- a/src/complex/ctanhf.c
+++ b/src/complex/ctanhf.c
@@ -44,17 +44,17 @@ float complex ctanhf(float complex z)
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
- return cpackf(x, (y == 0 ? y : x * y));
+ return CMPLXF(x, (y == 0 ? y : x * y));
SET_FLOAT_WORD(x, hx - 0x40000000);
- return cpackf(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
+ return CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
}
if (!isfinite(y))
- return cpackf(y - y, y - y);
+ return CMPLXF(y - y, y - y);
if (ix >= 0x41300000) { /* x >= 11 */
float exp_mx = expf(-fabsf(x));
- return cpackf(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
+ return CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
}
t = tanf(y);
@@ -62,5 +62,5 @@ float complex ctanhf(float complex z)
s = sinhf(x);
rho = sqrtf(1 + s * s);
denom = 1 + beta * s * s;
- return cpackf((beta * rho * s) / denom, t / denom);
+ return CMPLXF((beta * rho * s) / denom, t / denom);
}
diff --git a/src/complex/ctanl.c b/src/complex/ctanl.c
index 4b4c99b6..ac1c3e0a 100644
--- a/src/complex/ctanl.c
+++ b/src/complex/ctanl.c
@@ -8,7 +8,7 @@ long double complex ctanl(long double complex z)
#else
long double complex ctanl(long double complex z)
{
- z = ctanhl(cpackl(-cimagl(z), creall(z)));
- return cpackl(cimagl(z), -creall(z));
+ z = ctanhl(CMPLXL(-cimagl(z), creall(z)));
+ return CMPLXL(cimagl(z), -creall(z));
}
#endif
diff --git a/src/fenv/fenv.c b/src/fenv/fenv.c
index f77599bc..b2dfd254 100644
--- a/src/fenv/fenv.c
+++ b/src/fenv/fenv.c
@@ -19,7 +19,7 @@ int fetestexcept(int mask)
int fegetround(void)
{
- return 0;
+ return FE_TONEAREST;
}
int fesetround(int r)
diff --git a/src/fenv/feupdateenv.c b/src/fenv/feupdateenv.c
index 50cef8e5..f45ed7c0 100644
--- a/src/fenv/feupdateenv.c
+++ b/src/fenv/feupdateenv.c
@@ -2,6 +2,7 @@
int feupdateenv(const fenv_t *envp)
{
+ #pragma STDC FENV_ACCESS ON
int ex = fetestexcept(FE_ALL_EXCEPT);
fesetenv(envp);
feraiseexcept(ex);
diff --git a/src/internal/libm.h b/src/internal/libm.h
index a71c4c05..46c4b564 100644
--- a/src/internal/libm.h
+++ b/src/internal/libm.h
@@ -157,38 +157,23 @@ long double __tanl(long double, long double, int);
long double __polevll(long double, const long double *, int);
long double __p1evll(long double, const long double *, int);
-// FIXME: not needed when -fexcess-precision=standard is supported (>=gcc4.5)
-/*
- * Attempt to get strict C99 semantics for assignment with non-C99 compilers.
- */
-#if 1
+#if 0
+/* Attempt to get strict C99 semantics for assignment with non-C99 compilers. */
#define STRICT_ASSIGN(type, lval, rval) do { \
volatile type __v = (rval); \
(lval) = __v; \
} while (0)
#else
+/* Should work with -fexcess-precision=standard (>=gcc-4.5) or -ffloat-store */
#define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
#endif
-
/* complex */
-union dcomplex {
- double complex z;
- double a[2];
-};
-union fcomplex {
- float complex z;
- float a[2];
-};
-union lcomplex {
- long double complex z;
- long double a[2];
-};
-
-/* x + y*I is not supported properly by gcc */
-#define cpack(x,y) ((union dcomplex){.a={(x),(y)}}.z)
-#define cpackf(x,y) ((union fcomplex){.a={(x),(y)}}.z)
-#define cpackl(x,y) ((union lcomplex){.a={(x),(y)}}.z)
+#ifndef CMPLX
+#define CMPLX(x, y) __CMPLX(x, y, double)
+#define CMPLXF(x, y) __CMPLX(x, y, float)
+#define CMPLXL(x, y) __CMPLX(x, y, long double)
+#endif
#endif
diff --git a/src/math/__invtrigl.c b/src/math/__invtrigl.c
index a821842c..a51330e2 100644
--- a/src/math/__invtrigl.c
+++ b/src/math/__invtrigl.c
@@ -28,9 +28,8 @@
#include "__invtrigl.h"
#if LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-/*
- * asinl() and acosl()
- */
+
+/* coefficients used in asinl() and acosl() */
const long double
pS0 = 1.66666666666666666631e-01L,
pS1 = -4.16313987993683104320e-01L,
@@ -45,38 +44,9 @@ qS3 = -1.68285799854822427013e+00L,
qS4 = 3.90699412641738801874e-01L,
qS5 = -3.14365703596053263322e-02L;
-/*
- * atanl()
- */
-const long double atanhi[] = {
- 4.63647609000806116202e-01L,
- 7.85398163397448309628e-01L,
- 9.82793723247329067960e-01L,
- 1.57079632679489661926e+00L,
-};
-
-const long double atanlo[] = {
- 1.18469937025062860669e-20L,
- -1.25413940316708300586e-20L,
- 2.55232234165405176172e-20L,
- -2.50827880633416601173e-20L,
-};
-
-const long double aT[] = {
- 3.33333333333333333017e-01L,
- -1.99999999999999632011e-01L,
- 1.42857142857046531280e-01L,
- -1.11111111100562372733e-01L,
- 9.09090902935647302252e-02L,
- -7.69230552476207730353e-02L,
- 6.66661718042406260546e-02L,
- -5.88158892835030888692e-02L,
- 5.25499891539726639379e-02L,
- -4.70119845393155721494e-02L,
- 4.03539201366454414072e-02L,
- -2.91303858419364158725e-02L,
- 1.24822046299269234080e-02L,
-};
-
+const long double pi_hi = 3.1415926535897932384626433832795L;
const long double pi_lo = -5.01655761266833202345e-20L;
+const long double pio2_hi = 1.57079632679489661926L;
+const long double pio2_lo = -2.50827880633416601173e-20L;
+
#endif
diff --git a/src/math/__invtrigl.h b/src/math/__invtrigl.h
index c3ad3c49..22748b68 100644
--- a/src/math/__invtrigl.h
+++ b/src/math/__invtrigl.h
@@ -32,15 +32,6 @@
#define BIAS (LDBL_MAX_EXP - 1)
#define MANH_SIZE LDBL_MANH_SIZE
-/* Approximation thresholds. */
-#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */
-#define ACOS_CONST (BIAS - 65) /* 2**-65 */
-#define ATAN_CONST (BIAS + 65) /* 2**65 */
-#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */
-
-/* 0.95 */
-#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
-
/* Constants shared by the long double inverse trig functions. */
#define pS0 __pS0
#define pS1 __pS1
@@ -54,56 +45,24 @@
#define qS3 __qS3
#define qS4 __qS4
#define qS5 __qS5
-#define atanhi __atanhi
-#define atanlo __atanlo
-#define aT __aT
+#define pi_hi __pi_hi
#define pi_lo __pi_lo
+#define pio2_hi __pio2_hi
+#define pio2_lo __pio2_lo
-#define pio2_hi atanhi[3]
-#define pio2_lo atanlo[3]
-#define pio4_hi atanhi[1]
-
-#ifdef STRUCT_DECLS
-typedef struct longdouble {
- uint64_t mant;
- uint16_t expsign;
-} LONGDOUBLE;
-#else
-typedef long double LONGDOUBLE;
-#endif
-
-extern const LONGDOUBLE pS0, pS1, pS2, pS3, pS4, pS5, pS6;
-extern const LONGDOUBLE qS1, qS2, qS3, qS4, qS5;
-extern const LONGDOUBLE atanhi[], atanlo[], aT[];
-extern const LONGDOUBLE pi_lo;
-
-#ifndef STRUCT_DECLS
-static inline long double
-P(long double x)
-{
- return (x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 + x * \
- (pS4 + x * (pS5 + x * pS6)))))));
-}
-
-static inline long double
-Q(long double x)
-{
- return (1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5)))));
-}
+extern const long double pS0, pS1, pS2, pS3, pS4, pS5, pS6;
+extern const long double qS1, qS2, qS3, qS4, qS5;
+extern const long double pi_hi, pi_lo, pio2_hi, pio2_lo;
-static inline long double
-T_even(long double x)
+static long double P(long double x)
{
- return (aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] + x * \
- (aT[8] + x * (aT[10] + x * aT[12]))))));
+ return x * (pS0 + x * (pS1 + x * (pS2 + x * (pS3 +
+ x * (pS4 + x * (pS5 + x * pS6))))));
}
-static inline long double
-T_odd(long double x)
+static long double Q(long double x)
{
- return (aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] + x * \
- (aT[9] + x * aT[11])))));
+ return 1.0 + x * (qS1 + x * (qS2 + x * (qS3 + x * (qS4 + x * qS5))));
}
-#endif
#endif
diff --git a/src/math/acos.c b/src/math/acos.c
index 54d266ee..0eb15bed 100644
--- a/src/math/acos.c
+++ b/src/math/acos.c
@@ -38,6 +38,7 @@
static const double
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
+// FIXME
static const volatile double
pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
static const double
diff --git a/src/math/acosl.c b/src/math/acosl.c
index cc565336..83857d49 100644
--- a/src/math/acosl.c
+++ b/src/math/acosl.c
@@ -23,9 +23,7 @@ long double acosl(long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
-
-static const long double
-pi = 3.14159265358979323846264338327950280e+00L;
+#define ACOS_CONST (BIAS - 65) /* 2**-65 */
long double acosl(long double x)
{
@@ -41,7 +39,8 @@ long double acosl(long double x)
if (expsign > 0)
return 0.0; /* acos(1) = 0 */
else
- return pi + 2.0 * pio2_lo; /* acos(-1)= pi */
+ // FIXME
+ return pi_hi + 2.0 * pio2_lo; /* acos(-1)= pi */
}
return (x - x) / (x - x); /* acos(|x|>1) is NaN */
}
@@ -60,7 +59,7 @@ long double acosl(long double x)
s = sqrtl(z);
r = p / q;
w = r * s - pio2_lo;
- return pi - 2.0 * (s + w);
+ return pi_hi - 2.0 * (s + w);
} else { /* x > 0.5 */
z = (1.0 - x) * 0.5;
s = sqrtl(z);
diff --git a/src/math/asinl.c b/src/math/asinl.c
index ddd807e2..7572767c 100644
--- a/src/math/asinl.c
+++ b/src/math/asinl.c
@@ -24,6 +24,10 @@ long double asinl(long double x)
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
static const long double huge = 1.000e+300;
+static const long double pio4_hi = 7.85398163397448309628e-01L;
+#define ASIN_LINEAR (BIAS - 32) /* 2**-32 */
+/* 0.95 */
+#define THRESH ((0xe666666666666666ULL>>(64-(MANH_SIZE-1)))|LDBL_NBIT)
long double asinl(long double x)
{
diff --git a/src/math/atan2.c b/src/math/atan2.c
index 143c3834..825a38fe 100644
--- a/src/math/atan2.c
+++ b/src/math/atan2.c
@@ -39,6 +39,7 @@
#include "libm.h"
+// FIXME
static const volatile double
tiny = 1.0e-300;
static const double
diff --git a/src/math/atan2l.c b/src/math/atan2l.c
index 0fc901c8..45cbfcc1 100644
--- a/src/math/atan2l.c
+++ b/src/math/atan2l.c
@@ -24,10 +24,8 @@ long double atan2l(long double y, long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
-static const volatile long double
-tiny = 1.0e-300;
-static const long double
-pi = 3.14159265358979323846264338327950280e+00L;
+// FIXME:
+static const volatile long double tiny = 1.0e-300;
long double atan2l(long double y, long double x)
{
@@ -55,9 +53,9 @@ long double atan2l(long double y, long double x)
if (expty==0 && ((uy.bits.manh&~LDBL_NBIT)|uy.bits.manl)==0) {
switch(m) {
case 0:
- case 1: return y; /* atan(+-0,+anything)=+-0 */
- case 2: return pi+tiny; /* atan(+0,-anything) = pi */
- case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */
+ case 1: return y; /* atan(+-0,+anything)=+-0 */
+ case 2: return pi_hi+tiny; /* atan(+0,-anything) = pi */
+ case 3: return -pi_hi-tiny; /* atan(-0,-anything) =-pi */
}
}
/* when x = 0 */
@@ -69,15 +67,15 @@ long double atan2l(long double y, long double x)
switch(m) {
case 0: return pio2_hi*0.5+tiny; /* atan(+INF,+INF) */
case 1: return -pio2_hi*0.5-tiny; /* atan(-INF,+INF) */
- case 2: return 1.5*pio2_hi+tiny; /*atan(+INF,-INF)*/
- case 3: return -1.5*pio2_hi-tiny; /*atan(-INF,-INF)*/
+ case 2: return 1.5*pio2_hi+tiny; /* atan(+INF,-INF) */
+ case 3: return -1.5*pio2_hi-tiny; /* atan(-INF,-INF) */
}
} else {
switch(m) {
- case 0: return 0.0; /* atan(+...,+INF) */
- case 1: return -0.0; /* atan(-...,+INF) */
- case 2: return pi+tiny; /* atan(+...,-INF) */
- case 3: return -pi-tiny; /* atan(-...,-INF) */
+ case 0: return 0.0; /* atan(+...,+INF) */
+ case 1: return -0.0; /* atan(-...,+INF) */
+ case 2: return pi_hi+tiny; /* atan(+...,-INF) */
+ case 3: return -pi_hi-tiny; /* atan(-...,-INF) */
}
}
}
@@ -95,11 +93,11 @@ long double atan2l(long double y, long double x)
else /* safe to do y/x */
z = atanl(fabsl(y/x));
switch (m) {
- case 0: return z; /* atan(+,+) */
- case 1: return -z; /* atan(-,+) */
- case 2: return pi - (z-pi_lo); /* atan(+,-) */
+ case 0: return z; /* atan(+,+) */
+ case 1: return -z; /* atan(-,+) */
+ case 2: return pi_hi-(z-pi_lo); /* atan(+,-) */
default: /* case 3 */
- return (z-pi_lo) - pi; /* atan(-,-) */
+ return (z-pi_lo)-pi_hi; /* atan(-,-) */
}
}
#endif
diff --git a/src/math/atanl.c b/src/math/atanl.c
index 36072c17..33ecff0d 100644
--- a/src/math/atanl.c
+++ b/src/math/atanl.c
@@ -23,8 +23,53 @@ long double atanl(long double x)
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include "__invtrigl.h"
+
+#define ATAN_CONST (BIAS + 65) /* 2**65 */
+#define ATAN_LINEAR (BIAS - 32) /* 2**-32 */
static const long double huge = 1.0e300;
+static const long double atanhi[] = {
+ 4.63647609000806116202e-01L,
+ 7.85398163397448309628e-01L,
+ 9.82793723247329067960e-01L,
+ 1.57079632679489661926e+00L,
+};
+
+static const long double atanlo[] = {
+ 1.18469937025062860669e-20L,
+ -1.25413940316708300586e-20L,
+ 2.55232234165405176172e-20L,
+ -2.50827880633416601173e-20L,
+};
+
+static const long double aT[] = {
+ 3.33333333333333333017e-01L,
+ -1.99999999999999632011e-01L,
+ 1.42857142857046531280e-01L,
+ -1.11111111100562372733e-01L,
+ 9.09090902935647302252e-02L,
+ -7.69230552476207730353e-02L,
+ 6.66661718042406260546e-02L,
+ -5.88158892835030888692e-02L,
+ 5.25499891539726639379e-02L,
+ -4.70119845393155721494e-02L,
+ 4.03539201366454414072e-02L,
+ -2.91303858419364158725e-02L,
+ 1.24822046299269234080e-02L,
+};
+
+static long double T_even(long double x)
+{
+ return aT[0] + x * (aT[2] + x * (aT[4] + x * (aT[6] +
+ x * (aT[8] + x * (aT[10] + x * aT[12])))));
+}
+
+static long double T_odd(long double x)
+{
+ return aT[1] + x * (aT[3] + x * (aT[5] + x * (aT[7] +
+ x * (aT[9] + x * aT[11]))));
+}
+
long double atanl(long double x)
{
union IEEEl2bits u;
diff --git a/src/math/exp10l.c b/src/math/exp10l.c
index 4b2b7ef6..22a4636a 100644
--- a/src/math/exp10l.c
+++ b/src/math/exp10l.c
@@ -5,8 +5,8 @@
long double exp10l(long double x)
{
static const long double p10[] = {
- 1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
- 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+ 1e-15L, 1e-14L, 1e-13L, 1e-12L, 1e-11L, 1e-10L,
+ 1e-9L, 1e-8L, 1e-7L, 1e-6L, 1e-5L, 1e-4L, 1e-3L, 1e-2L, 1e-1L,
1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15
};
diff --git a/src/math/expm1l.c b/src/math/expm1l.c
index ca0d141f..21a86c00 100644
--- a/src/math/expm1l.c
+++ b/src/math/expm1l.c
@@ -44,13 +44,7 @@
*
* Relative error:
* arithmetic domain # trials peak rms
- * IEEE -45,+MAXLOG 200,000 1.2e-19 2.5e-20
- *
- * ERROR MESSAGES:
- *
- * message condition value returned
- * expm1l overflow x > MAXLOG MAXNUM
- *
+ * IEEE -45,+maxarg 200,000 1.2e-19 2.5e-20
*/
#include "libm.h"
@@ -61,7 +55,6 @@ long double expm1l(long double x)
return expm1(x);
}
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
-static const long double MAXLOGL = 1.1356523406294143949492E4L;
/* exp(x) - 1 = x + 0.5 x^2 + x^3 P(x)/Q(x)
-.5 ln 2 < x < .5 ln 2
@@ -83,19 +76,20 @@ C1 = 6.93145751953125E-1L,
C2 = 1.428606820309417232121458176568075500134E-6L,
/* ln 2^-65 */
minarg = -4.5054566736396445112120088E1L,
-huge = 0x1p10000L;
+/* ln 2^16384 */
+maxarg = 1.1356523406294143949492E4L;
long double expm1l(long double x)
{
long double px, qx, xx;
int k;
- /* Overflow. */
- if (x > MAXLOGL)
- return huge*huge; /* overflow */
+ if (isnan(x))
+ return x;
+ if (x > maxarg)
+ return x*0x1p16383L; /* overflow, unless x==inf */
if (x == 0.0)
return x;
- /* Minimum value.*/
if (x < minarg)
return -1.0;
diff --git a/src/math/fma.c b/src/math/fma.c
index 7076d4b9..17f1cdcc 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -89,6 +89,7 @@ static int getexp(long double x)
double fma(double x, double y, double z)
{
+ #pragma STDC FENV_ACCESS ON
long double hi, lo1, lo2, xy;
int round, ez, exy;
@@ -306,6 +307,7 @@ static inline struct dd dd_mul(double a, double b)
*/
double fma(double x, double y, double z)
{
+ #pragma STDC FENV_ACCESS ON
double xs, ys, zs, adj;
struct dd xy, r;
int oround;
diff --git a/src/math/fmaf.c b/src/math/fmaf.c
index cc427fb5..a1c7f4f8 100644
--- a/src/math/fmaf.c
+++ b/src/math/fmaf.c
@@ -37,6 +37,7 @@
*/
float fmaf(float x, float y, float z)
{
+ #pragma STDC FENV_ACCESS ON
double xy, result;
uint32_t hr, lr;
diff --git a/src/math/fmal.c b/src/math/fmal.c
index be64f145..ccbe434d 100644
--- a/src/math/fmal.c
+++ b/src/math/fmal.c
@@ -162,6 +162,7 @@ static inline struct dd dd_mul(long double a, long double b)
*/
long double fmal(long double x, long double y, long double z)
{
+ #pragma STDC FENV_ACCESS ON
long double xs, ys, zs, adj;
struct dd xy, r;
int oround;
diff --git a/src/math/hypot.c b/src/math/hypot.c
index ba4c7575..9a4cbdb3 100644
--- a/src/math/hypot.c
+++ b/src/math/hypot.c
@@ -117,12 +117,7 @@ double hypot(double x, double y)
t2 = a - t1;
w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
}
- if (k != 0) {
- uint32_t high;
- t1 = 1.0;
- GET_HIGH_WORD(high, t1);
- SET_HIGH_WORD(t1, high+(k<<20));
- return t1*w;
- }
+ if (k)
+ w = scalbn(w, k);
return w;
}
diff --git a/src/math/hypotf.c b/src/math/hypotf.c
index 9fd77e6a..4d80178d 100644
--- a/src/math/hypotf.c
+++ b/src/math/hypotf.c
@@ -80,9 +80,7 @@ float hypotf(float x, float y)
t2 = a - t1;
w = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
}
- if (k != 0) {
- SET_FLOAT_WORD(t1, 0x3f800000+(k<<23));
- return t1*w;
- }
+ if (k)
+ w = scalbnf(w, k);
return w;
}
diff --git a/src/math/ilogb.c b/src/math/ilogb.c
index 0a3a6a46..5a1819d8 100644
--- a/src/math/ilogb.c
+++ b/src/math/ilogb.c
@@ -8,13 +8,17 @@ int ilogb(double x)
if (!e) {
u.bits <<= 12;
- if (u.bits == 0)
+ if (u.bits == 0) {
+ FORCE_EVAL(0/0.0f);
return FP_ILOGB0;
+ }
/* subnormal x */
for (e = -0x3ff; u.bits < (uint64_t)1<<63; e--, u.bits<<=1);
return e;
}
- if (e == 0x7ff)
+ if (e == 0x7ff) {
+ FORCE_EVAL(0/0.0f);
return u.bits<<12 ? FP_ILOGBNAN : INT_MAX;
+ }
return e - 0x3ff;
}
diff --git a/src/math/ilogbf.c b/src/math/ilogbf.c
index 272cbdac..42cd62e2 100644
--- a/src/math/ilogbf.c
+++ b/src/math/ilogbf.c
@@ -8,13 +8,17 @@ int ilogbf(float x)
if (!e) {
u.bits <<= 9;
- if (u.bits == 0)
+ if (u.bits == 0) {
+ FORCE_EVAL(0/0.0f);
return FP_ILOGB0;
+ }
/* subnormal x */
for (e = -0x7f; u.bits < (uint32_t)1<<31; e--, u.bits<<=1);
return e;
}
- if (e == 0xff)
+ if (e == 0xff) {
+ FORCE_EVAL(0/0.0f);
return u.bits<<9 ? FP_ILOGBNAN : INT_MAX;
+ }
return e - 0x7f;
}
diff --git a/src/math/ilogbl.c b/src/math/ilogbl.c
index ed9ddcbc..1512934f 100644
--- a/src/math/ilogbl.c
+++ b/src/math/ilogbl.c
@@ -14,15 +14,19 @@ int ilogbl(long double x)
int e = u.bits.exp;
if (!e) {
- if (m == 0)
+ if (m == 0) {
+ FORCE_EVAL(0/0.0f);
return FP_ILOGB0;
+ }
/* subnormal x */
for (e = -0x3fff+1; m < (uint64_t)1<<63; e--, m<<=1);
return e;
}
- if (e == 0x7fff)
+ if (e == 0x7fff) {
+ FORCE_EVAL(0/0.0f);
/* in ld80 msb is set in inf */
return m & (uint64_t)-1>>1 ? FP_ILOGBNAN : INT_MAX;
+ }
return e - 0x3fff;
}
#endif
diff --git a/src/math/llrintl.c b/src/math/llrintl.c
index 32bb8b03..3449f6f2 100644
--- a/src/math/llrintl.c
+++ b/src/math/llrintl.c
@@ -18,6 +18,7 @@ raises inexact (with tonearest or upward rounding mode)
*/
long long llrintl(long double x)
{
+ #pragma STDC FENV_ACCESS ON
int e;
e = fetestexcept(FE_INEXACT);
diff --git a/src/math/log1pl.c b/src/math/log1pl.c
index 1400d365..edb48df1 100644
--- a/src/math/log1pl.c
+++ b/src/math/log1pl.c
@@ -46,11 +46,6 @@
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1.0, 9.0 100000 8.2e-20 2.5e-20
- *
- * ERROR MESSAGES:
- *
- * log singularity: x-1 = 0; returns -INFINITY
- * log domain: x-1 < 0; returns NAN
*/
#include "libm.h"
@@ -123,8 +118,8 @@ long double log1pl(long double xm1)
/* Test for domain errors. */
if (x <= 0.0) {
if (x == 0.0)
- return -INFINITY;
- return NAN;
+ return -1/x; /* -inf with divbyzero */
+ return 0/0.0f; /* nan with invalid */
}
/* Separate mantissa from exponent.
diff --git a/src/math/log2l.c b/src/math/log2l.c
index 8ebce9c4..345b395d 100644
--- a/src/math/log2l.c
+++ b/src/math/log2l.c
@@ -50,11 +50,6 @@
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
- *
- * ERROR MESSAGES:
- *
- * log singularity: x = 0; returns -INFINITY
- * log domain: x < 0; returns NAN
*/
#include "libm.h"
@@ -113,8 +108,7 @@ static const long double S[4] = {
long double log2l(long double x)
{
- volatile long double z;
- long double y;
+ long double y, z;
int e;
if (isnan(x))
@@ -123,8 +117,8 @@ long double log2l(long double x)
return x;
if (x <= 0.0) {
if (x == 0.0)
- return -INFINITY;
- return NAN;
+ return -1/(x+0); /* -inf with divbyzero */
+ return 0/0.0f; /* nan with invalid */
}
/* separate mantissa from exponent */
diff --git a/src/math/logb.c b/src/math/logb.c
index f7cd7613..624425a8 100644
--- a/src/math/logb.c
+++ b/src/math/logb.c
@@ -1,20 +1,17 @@
-#include <limits.h>
#include "libm.h"
/*
special cases:
- logb(+-0) = -inf
+ logb(+-0) = -inf, and raise divbyzero
logb(+-inf) = +inf
logb(nan) = nan
-these are calculated at runtime to raise fp exceptions
*/
-double logb(double x) {
- int i = ilogb(x);
-
- if (i == FP_ILOGB0)
- return -1.0/fabs(x);
- if (i == FP_ILOGBNAN || i == INT_MAX)
+double logb(double x)
+{
+ if (!isfinite(x))
return x * x;
- return i;
+ if (x == 0)
+ return -1/(x+0);
+ return ilogb(x);
}
diff --git a/src/math/logbf.c b/src/math/logbf.c
index 934827f8..950d3569 100644
--- a/src/math/logbf.c
+++ b/src/math/logbf.c
@@ -1,12 +1,10 @@
-#include <limits.h>
#include "libm.h"
-float logbf(float x) {
- int i = ilogbf(x);
-
- if (i == FP_ILOGB0)
- return -1.0f/fabsf(x);
- if (i == FP_ILOGBNAN || i == INT_MAX)
+float logbf(float x)
+{
+ if (!isfinite(x))
return x * x;
- return i;
+ if (x == 0)
+ return -1/(x+0);
+ return ilogbf(x);
}
diff --git a/src/math/logbl.c b/src/math/logbl.c
index 5d04abd7..f3850745 100644
--- a/src/math/logbl.c
+++ b/src/math/logbl.c
@@ -1,4 +1,3 @@
-#include <limits.h>
#include "libm.h"
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double logbl(long double x)
@@ -8,12 +7,10 @@ long double logbl(long double x)
#else
long double logbl(long double x)
{
- int i = ilogbl(x);
-
- if (i == FP_ILOGB0)
- return -1.0/fabsl(x);
- if (i == FP_ILOGBNAN || i == INT_MAX)
+ if (!isfinite(x))
return x * x;
- return i;
+ if (x == 0)
+ return -1/(x+0);
+ return ilogbl(x);
}
#endif
diff --git a/src/math/logl.c b/src/math/logl.c
index ffd81038..ef2b5515 100644
--- a/src/math/logl.c
+++ b/src/math/logl.c
@@ -50,11 +50,6 @@
* In the tests over the interval exp(+-10000), the logarithms
* of the random arguments were uniformly distributed over
* [-10000, +10000].
- *
- * ERROR MESSAGES:
- *
- * log singularity: x = 0; returns -INFINITY
- * log domain: x < 0; returns NAN
*/
#include "libm.h"
@@ -121,8 +116,8 @@ long double logl(long double x)
return x;
if (x <= 0.0) {
if (x == 0.0)
- return -INFINITY;
- return NAN;
+ return -1/(x+0); /* -inf with divbyzero */
+ return 0/0.0f; /* nan with invalid */
}
/* separate mantissa from exponent */
diff --git a/src/math/lrint.c b/src/math/lrint.c
index 67091835..bdca8b7c 100644
--- a/src/math/lrint.c
+++ b/src/math/lrint.c
@@ -28,6 +28,7 @@ as a double.
#if LONG_MAX < 1U<<53 && defined(FE_INEXACT)
long lrint(double x)
{
+ #pragma STDC FENV_ACCESS ON
int e;
e = fetestexcept(FE_INEXACT);
diff --git a/src/math/lrintl.c b/src/math/lrintl.c
index 5eb1ba7e..b2a8106d 100644
--- a/src/math/lrintl.c
+++ b/src/math/lrintl.c
@@ -18,6 +18,7 @@ raises inexact (with tonearest or upward rounding mode)
*/
long lrintl(long double x)
{
+ #pragma STDC FENV_ACCESS ON
int e;
e = fetestexcept(FE_INEXACT);
diff --git a/src/math/modf.c b/src/math/modf.c
index de45069f..8f337ef0 100644
--- a/src/math/modf.c
+++ b/src/math/modf.c
@@ -1,5 +1,4 @@
-#include <math.h>
-#include <stdint.h>
+#include "libm.h"
double modf(double x, double *iptr)
{
@@ -33,5 +32,6 @@ double modf(double x, double *iptr)
}
u.n &= ~mask;
*iptr = u.x;
- return x - *iptr;
+ STRICT_ASSIGN(double, x, x - *iptr);
+ return x;
}
diff --git a/src/math/modff.c b/src/math/modff.c
index 84d0b82a..90caf316 100644
--- a/src/math/modff.c
+++ b/src/math/modff.c
@@ -1,5 +1,4 @@
-#include <math.h>
-#include <stdint.h>
+#include "libm.h"
float modff(float x, float *iptr)
{
@@ -33,5 +32,6 @@ float modff(float x, float *iptr)
}
u.n &= ~mask;
*iptr = u.x;
- return x - *iptr;
+ STRICT_ASSIGN(float, x, x - *iptr);
+ return x;
}
diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c
index 7a4c58cf..f4e8aac4 100644
--- a/src/math/nearbyint.c
+++ b/src/math/nearbyint.c
@@ -6,6 +6,7 @@
double nearbyint(double x)
{
#ifdef FE_INEXACT
+ #pragma STDC FENV_ACCESS ON
int e;
e = fetestexcept(FE_INEXACT);
diff --git a/src/math/nearbyintf.c b/src/math/nearbyintf.c
index 39c3d73b..092e9ffa 100644
--- a/src/math/nearbyintf.c
+++ b/src/math/nearbyintf.c
@@ -4,6 +4,7 @@
float nearbyintf(float x)
{
#ifdef FE_INEXACT
+ #pragma STDC FENV_ACCESS ON
int e;
e = fetestexcept(FE_INEXACT);
diff --git a/src/math/nearbyintl.c b/src/math/nearbyintl.c
index 0ff4b1f9..82852492 100644
--- a/src/math/nearbyintl.c
+++ b/src/math/nearbyintl.c
@@ -11,6 +11,7 @@ long double nearbyintl(long double x)
long double nearbyintl(long double x)
{
#ifdef FE_INEXACT
+ #pragma STDC FENV_ACCESS ON
int e;
e = fetestexcept(FE_INEXACT);
diff --git a/src/math/nextafter.c b/src/math/nextafter.c
index a3b42c99..9ee82518 100644
--- a/src/math/nextafter.c
+++ b/src/math/nextafter.c
@@ -27,7 +27,7 @@ double nextafter(double x, double y)
e = ux.bits >> 52 & 0x7ff;
/* raise overflow if ux.value is infinite and x is finite */
if (e == 0x7ff)
- return x + x;
+ FORCE_EVAL(x+x);
/* raise underflow if ux.value is subnormal or zero */
if (e == 0)
FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/nextafterf.c b/src/math/nextafterf.c
index b703487b..22b61dce 100644
--- a/src/math/nextafterf.c
+++ b/src/math/nextafterf.c
@@ -26,7 +26,7 @@ float nextafterf(float x, float y)
e = ux.bits & 0x7f800000;
/* raise overflow if ux.value is infinite and x is finite */
if (e == 0x7f800000)
- return x + x;
+ FORCE_EVAL(x+x);
/* raise underflow if ux.value is subnormal or zero */
if (e == 0)
FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/nexttoward.c b/src/math/nexttoward.c
index 7355f2f1..6f32eca4 100644
--- a/src/math/nexttoward.c
+++ b/src/math/nexttoward.c
@@ -36,7 +36,7 @@ double nexttoward(double x, long double y)
e = ux.bits>>52 & 0x7ff;
/* raise overflow if ux.value is infinite and x is finite */
if (e == 0x7ff)
- return x + x;
+ FORCE_EVAL(x+x);
/* raise underflow if ux.value is subnormal or zero */
if (e == 0)
FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/nexttowardf.c b/src/math/nexttowardf.c
index 8648be6a..9a693b1a 100644
--- a/src/math/nexttowardf.c
+++ b/src/math/nexttowardf.c
@@ -28,7 +28,7 @@ float nexttowardf(float x, long double y)
e = ux.bits & 0x7f800000;
/* raise overflow if ux.value is infinite and x is finite */
if (e == 0x7f800000)
- return x + x;
+ FORCE_EVAL(x+x);
/* raise underflow if ux.value is subnormal or zero */
if (e == 0)
FORCE_EVAL(x*x + ux.value*ux.value);
diff --git a/src/math/scalbn.c b/src/math/scalbn.c
index c9c7af80..003141e3 100644
--- a/src/math/scalbn.c
+++ b/src/math/scalbn.c
@@ -10,8 +10,10 @@ double scalbn(double x, int n)
if (n > 1023) {
x *= 0x1p1023;
n -= 1023;
- if (n > 1023)
- return x * 0x1p1023;
+ if (n > 1023) {
+ STRICT_ASSIGN(double, x, x * 0x1p1023);
+ return x;
+ }
}
} else if (n < -1022) {
x *= 0x1p-1022;
@@ -19,10 +21,13 @@ double scalbn(double x, int n)
if (n < -1022) {
x *= 0x1p-1022;
n += 1022;
- if (n < -1022)
- return x * 0x1p-1022;
+ if (n < -1022) {
+ STRICT_ASSIGN(double, x, x * 0x1p-1022);
+ return x;
+ }
}
}
INSERT_WORDS(scale, (uint32_t)(0x3ff+n)<<20, 0);
- return x * scale;
+ STRICT_ASSIGN(double, x, x * scale);
+ return x;
}
diff --git a/src/math/scalbnf.c b/src/math/scalbnf.c
index 243dafdf..f94b5d59 100644
--- a/src/math/scalbnf.c
+++ b/src/math/scalbnf.c
@@ -10,8 +10,10 @@ float scalbnf(float x, int n)
if (n > 127) {
x *= 0x1p127f;
n -= 127;
- if (n > 127)
- return x * 0x1p127f;
+ if (n > 127) {
+ STRICT_ASSIGN(float, x, x * 0x1p127f);
+ return x;
+ }
}
} else if (n < -126) {
x *= 0x1p-126f;
@@ -19,10 +21,13 @@ float scalbnf(float x, int n)
if (n < -126) {
x *= 0x1p-126f;
n += 126;
- if (n < -126)
- return x * 0x1p-126f;
+ if (n < -126) {
+ STRICT_ASSIGN(float, x, x * 0x1p-126f);
+ return x;
+ }
}
}
SET_FLOAT_WORD(scale, (uint32_t)(0x7f+n)<<23);
- return x * scale;
+ STRICT_ASSIGN(float, x, x * scale);
+ return x;
}