}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#include <fenv.h>
+#if LDBL_MANT_DIG == 64
+#define LASTBIT(u) (u.i.m & 1)
+#define SPLIT (0x1p32L + 1)
+#elif LDBL_MANT_DIG == 113
+#define LASTBIT(u) (u.i.lo & 1)
+#define SPLIT (0x1p57L + 1)
+#endif
/*
* A struct dd represents a floating-point number with twice the precision
static inline long double add_adjusted(long double a, long double b)
{
struct dd sum;
- union IEEEl2bits u;
+ union ldshape u;
sum = dd_add(a, b);
if (sum.lo != 0) {
- u.e = sum.hi;
- if ((u.bits.manl & 1) == 0)
+ u.f = sum.hi;
+ if (!LASTBIT(u))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (sum.hi);
{
struct dd sum;
int bits_lost;
- union IEEEl2bits u;
+ union ldshape u;
sum = dd_add(a, b);
* break the ties manually.
*/
if (sum.lo != 0) {
- u.e = sum.hi;
- bits_lost = -u.bits.exp - scale + 1;
- if (bits_lost != 1 ^ (int)(u.bits.manl & 1))
+ u.f = sum.hi;
+ bits_lost = -u.i.se - scale + 1;
+ if ((bits_lost != 1) ^ LASTBIT(u))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return scalbnl(sum.hi, scale);
*/
static inline struct dd dd_mul(long double a, long double b)
{
-#if LDBL_MANT_DIG == 64
- static const long double split = 0x1p32L + 1.0;
-#elif LDBL_MANT_DIG == 113
- static const long double split = 0x1p57L + 1.0;
-#endif
struct dd ret;
long double ha, hb, la, lb, p, q;
- p = a * split;
+ p = a * SPLIT;
ha = a - p;
ha += p;
la = a - ha;
- p = b * split;
+ p = b * SPLIT;
hb = b - p;
hb += p;
lb = b - hb;
*/
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;
/*
* There is no need to worry about double rounding in directed
* rounding modes.
+ * But underflow may not be raised correctly, example in downward rounding:
+ * fmal(0x1.0000000001p-16000L, 0x1.0000000001p-400L, -0x1p-16440L)
*/
+ long double ret;
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+ int e = fetestexcept(FE_INEXACT);
+ feclearexcept(FE_INEXACT);
+#endif
fesetround(oround);
adj = r.lo + xy.lo;
- return scalbnl(r.hi + adj, spread);
+ ret = scalbnl(r.hi + adj, spread);
+#if defined(FE_INEXACT) && defined(FE_UNDERFLOW)
+ if (ilogbl(ret) < -16382 && fetestexcept(FE_INEXACT))
+ feraiseexcept(FE_UNDERFLOW);
+ else if (e)
+ feraiseexcept(FE_INEXACT);
+#endif
+ return ret;
}
adj = add_adjusted(r.lo, xy.lo);