+++ /dev/null
-#include <stdio.h>
-#include <stdint.h>
-#include "util.h"
-
-int eulpf(float x)
-{
- union { float f; uint32_t i; } u = { x };
- int e = u.i>>23 & 0xff;
-
- if (!e)
- e++;
- return e - 0x7f - 23;
-}
-
-int eulp(double x)
-{
- union { double f; uint64_t i; } u = { x };
- int e = u.i>>52 & 0x7ff;
-
- if (!e)
- e++;
- return e - 0x3ff - 52;
-}
-
-int eulpl(long double x)
-{
-#if LDBL_MANT_DIG == 53
- return eulp(x);
-#elif LDBL_MANT_DIG == 64
- union { long double f; struct {uint64_t m; uint16_t e; uint16_t pad;} i; } u = { x };
- int e = u.i.e & 0x7fff;
-
- if (!e)
- e++;
- return e - 0x3fff - 63;
-#else
- // TODO
- return 0;
-#endif
-}
-
-float ulperrf(float got, float want, float dwant)
-{
- if (isnan(got) && isnan(want))
- return 0;
- if (got == want) {
- if (signbit(got) == signbit(want))
- return dwant;
- return inf;
- }
- if (isinf(got)) {
- got = copysignf(0x1p127, got);
- want *= 0.5;
- }
- return scalbn(got - want, -eulpf(want)) + dwant;
-}
-
-float ulperr(double got, double want, float dwant)
-{
- if (isnan(got) && isnan(want))
- return 0;
- if (got == want) {
- if (signbit(got) == signbit(want))
- return dwant;
- return inf; // treat 0 sign errors badly
- }
- if (isinf(got)) {
- got = copysign(0x1p1023, got);
- want *= 0.5;
- }
- return scalbn(got - want, -eulp(want)) + dwant;
-}
-
-float ulperrl(long double got, long double want, float dwant)
-{
-#if LDBL_MANT_DIG == 53
- return ulperr(got, want, dwant);
-#elif LDBL_MANT_DIG == 64
- if (isnan(got) && isnan(want))
- return 0;
- if (got == want) {
- if (signbit(got) == signbit(want))
- return dwant;
- return inf;
- }
- if (isinf(got)) {
- got = copysignl(0x1p16383L, got);
- want *= 0.5;
- }
- return scalbnl(got - want, -eulpl(want)) + dwant;
-#else
- // TODO
- return inf;
-#endif
-}
-
-#define length(a) (sizeof(a)/sizeof*(a))
-#define flag(x) {x, #x}
-static struct {
- int flag;
- char *s;
-} eflags[] = {
- flag(INEXACT),
- flag(INVALID),
- flag(DIVBYZERO),
- flag(UNDERFLOW),
- flag(OVERFLOW)
-};
-
-char *estr(int f)
-{
- static char buf[256];
- char *p = buf;
- int i, all = 0;
-
- for (i = 0; i < length(eflags); i++)
- if (f & eflags[i].flag) {
- p += sprintf(p, "%s%s", all ? "|" : "", eflags[i].s);
- all |= eflags[i].flag;
- }
- if (all != f) {
- p += sprintf(p, "%s%d", all ? "|" : "", f & ~all);
- all = f;
- }
- p += sprintf(p, "%s", all ? "" : "0");
- return buf;
-}
-
-char *rstr(int r)
-{
- switch (r) {
- case RN: return "RN";
-#ifdef FE_TOWARDZERO
- case RZ: return "RZ";
-#endif
-#ifdef FE_UPWARD
- case RU: return "RU";
-#endif
-#ifdef FE_DOWNWARD
- case RD: return "RD";
-#endif
- }
- return "R?";
-}
+++ /dev/null
-#include <fenv.h>
-#include <float.h>
-#include <math.h>
-
-#undef RN
-#undef RZ
-#undef RD
-#undef RU
-#ifdef FE_TONEAREST
-#define RN FE_TONEAREST
-#else
-#define RN 0
-#endif
-#ifdef FE_TOWARDZERO
-#define RZ FE_TOWARDZERO
-#else
-#define RZ -1
-#endif
-#ifdef FE_DOWNWARD
-#define RD FE_DOWNWARD
-#else
-#define RD -1
-#endif
-#ifdef FE_UPWARD
-#define RU FE_UPWARD
-#else
-#define RU -1
-#endif
-
-#undef INEXACT
-#undef INVALID
-#undef DIVBYZERO
-#undef UNDERFLOW
-#undef OVERFLOW
-#ifdef FE_INEXACT
-#define INEXACT FE_INEXACT
-#else
-#define INEXACT 0
-#endif
-#ifdef FE_INVALID
-#define INVALID FE_INVALID
-#else
-#define INVALID 0
-#endif
-#ifdef FE_DIVBYZERO
-#define DIVBYZERO FE_DIVBYZERO
-#else
-#define DIVBYZERO 0
-#endif
-#ifdef FE_UNDERFLOW
-#define UNDERFLOW FE_UNDERFLOW
-#else
-#define UNDERFLOW 0
-#endif
-#ifdef FE_OVERFLOW
-#define OVERFLOW FE_OVERFLOW
-#else
-#define OVERFLOW 0
-#endif
-
-#undef inf
-#undef nan
-#define inf INFINITY
-#define nan NAN
-
-#define T(...) {__FILE__, __LINE__, __VA_ARGS__},
-
-#define POS char *file; int line;
-struct d_d {POS int r; double x; double y; float dy; int e; };
-struct f_f {POS int r; float x; float y; float dy; int e; };
-struct l_l {POS int r; long double x; long double y; float dy; int e; };
-struct ff_f {POS int r; float x; float x2; float y; float dy; int e; };
-struct dd_d {POS int r; double x; double x2; double y; float dy; int e; };
-struct ll_l {POS int r; long double x; long double x2; long double y; float dy; int e; };
-struct d_di {POS int r; double x; double y; float dy; long long i; int e; };
-struct f_fi {POS int r; float x; float y; float dy; long long i; int e; };
-struct l_li {POS int r; long double x; long double y; float dy; long long i; int e; };
-struct di_d {POS int r; double x; long long i; double y; float dy; int e; };
-struct fi_f {POS int r; float x; long long i; float y; float dy; int e; };
-struct li_l {POS int r; long double x; long long i; long double y; float dy; int e; };
-struct d_i {POS int r; double x; long long i; int e; };
-struct f_i {POS int r; float x; long long i; int e; };
-struct l_i {POS int r; long double x; long long i; int e; };
-struct d_dd {POS int r; double x; double y; float dy; double y2; float dy2; int e; };
-struct f_ff {POS int r; float x; float y; float dy; float y2; float dy2; int e; };
-struct l_ll {POS int r; long double x; long double y; float dy; long double y2; float dy2; int e; };
-struct ff_fi {POS int r; float x; float x2; float y; float dy; long long i; int e; };
-struct dd_di {POS int r; double x; double x2; double y; float dy; long long i; int e; };
-struct ll_li {POS int r; long double x; long double x2; long double y; float dy; long long i; int e; };
-struct fff_f {POS int r; float x; float x2; float x3; float y; float dy; int e; };
-struct ddd_d {POS int r; double x; double x2; double x3; double y; float dy; int e; };
-struct lll_l {POS int r; long double x; long double x2; long double x3; long double y; float dy; int e; };
-#undef POS
-
-char *estr(int);
-char *rstr(int);
-
-float ulperr(double got, double want, float dwant);
-float ulperrf(float got, float want, float dwant);
-float ulperrl(long double got, long double want, float dwant);
-
-static int checkexcept(int got, int want, int r)
-{
- if (r == RN)
- return got == want || got == (want|INEXACT);
- return (got|INEXACT|UNDERFLOW) == (want|INEXACT|UNDERFLOW);
-}
-
-static int checkulp(float d, int r)
-{
- // TODO: we only care about >=1.5 ulp errors for now, should be 1.0
- if (r == RN)
- return fabsf(d) < 1.5;
- return 1;
-}
-
-static int checkcr(long double y, long double ywant, int r)
-{
- if (isnan(ywant))
- return isnan(y);
- return y == ywant && signbit(y) == signbit(ywant);
-}
-