>From c4cf94d37e2c556a2f77b1debe2fe9ec269fdb3e Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Sat, 21 Nov 2015 21:23:30 +0000 Subject: [PATCH] math: make -ffloat-store work on i386 old i386 toolchains don't support -fexcess-precision=standard and then musl configure defaults to -ffloat-store to avoid problems because of excess precision. however there are cases when excess precision must not be rounded away so -ffloat-store breaks the code. this patch adds (double_t) casts to double precision arithmetics when the excess precision is required for correct results. (without -ffloat-store or on non-i386 targets the generated code is not changed.) --- src/math/__rem_pio2.c | 2 +- src/math/__rem_pio2f.c | 2 +- src/math/hypot.c | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/math/__rem_pio2.c b/src/math/__rem_pio2.c index a40db9f..d403f81 100644 --- a/src/math/__rem_pio2.c +++ b/src/math/__rem_pio2.c @@ -118,7 +118,7 @@ int __rem_pio2(double x, double *y) if (ix < 0x413921fb) { /* |x| ~< 2^20*(pi/2), medium size */ medium: /* rint(x/(pi/2)), Assume round-to-nearest. */ - fn = x*invpio2 + toint - toint; + fn = (double_t)x*invpio2 + toint - toint; n = (int32_t)fn; r = x - fn*pio2_1; w = fn*pio2_1t; /* 1st round, good to 85 bits */ diff --git a/src/math/__rem_pio2f.c b/src/math/__rem_pio2f.c index f516385..4473c1c 100644 --- a/src/math/__rem_pio2f.c +++ b/src/math/__rem_pio2f.c @@ -51,7 +51,7 @@ int __rem_pio2f(float x, double *y) /* 25+53 bit pi is good enough for medium size */ if (ix < 0x4dc90fdb) { /* |x| ~< 2^28*(pi/2), medium size */ /* Use a specialized rint() to get fn. Assume round-to-nearest. */ - fn = x*invpio2 + toint - toint; + fn = (double_t)x*invpio2 + toint - toint; n = (int32_t)fn; *y = x - fn*pio2_1 - fn*pio2_1t; return n; diff --git a/src/math/hypot.c b/src/math/hypot.c index 29ec6a4..6071bf1 100644 --- a/src/math/hypot.c +++ b/src/math/hypot.c @@ -12,10 +12,10 @@ static void sq(double_t *hi, double_t *lo, double x) { double_t xh, xl, xc; - xc = x*SPLIT; + xc = (double_t)x*SPLIT; xh = x - xc + xc; xl = x - xh; - *hi = x*x; + *hi = (double_t)x*x; *lo = xh*xh - *hi + 2*xh*xl + xl*xl; } -- 2.4.1