|
nsz |
f9d179 |
#include <stdio.h>
|
|
nsz |
f9d179 |
#include <stdint.h>
|
|
nsz |
f9d179 |
#include "util.h"
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
int eulpf(float x)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
union { float f; uint32_t i; } u = { x };
|
|
nsz |
f9d179 |
int e = u.i>>23 & 0xff;
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
if (!e)
|
|
nsz |
f9d179 |
e++;
|
|
nsz |
f9d179 |
return e - 0x7f - 23;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
int eulp(double x)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
union { double f; uint64_t i; } u = { x };
|
|
nsz |
f9d179 |
int e = u.i>>52 & 0x7ff;
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
if (!e)
|
|
nsz |
f9d179 |
e++;
|
|
nsz |
f9d179 |
return e - 0x3ff - 52;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
int eulpl(long double x)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
#if LDBL_MANT_DIG == 53
|
|
nsz |
f9d179 |
return eulp(x);
|
|
nsz |
f9d179 |
#elif LDBL_MANT_DIG == 64
|
|
nsz |
f9d179 |
union { long double f; struct {uint64_t m; uint16_t e; uint16_t pad;} i; } u = { x };
|
|
nsz |
f9d179 |
int e = u.i.e & 0x7fff;
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
if (!e)
|
|
nsz |
f9d179 |
e++;
|
|
nsz |
f9d179 |
return e - 0x3fff - 63;
|
|
nsz |
f9d179 |
#else
|
|
nsz |
f9d179 |
// TODO
|
|
nsz |
f9d179 |
return 0;
|
|
nsz |
f9d179 |
#endif
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
float ulperrf(float got, float want, float dwant)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
if (isnan(got) && isnan(want))
|
|
nsz |
f9d179 |
return 0;
|
|
nsz |
f9d179 |
if (got == want && signbit(got) == signbit(want))
|
|
nsz |
f9d179 |
return dwant;
|
|
nsz |
f9d179 |
if (isinf(got)) {
|
|
nsz |
f9d179 |
got = copysignf(0x1p127, got);
|
|
nsz |
f9d179 |
want *= 0.5;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
return scalbn(got - want, -eulpf(want)) + dwant;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
float ulperr(double got, double want, float dwant)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
if (isnan(got) && isnan(want))
|
|
nsz |
f9d179 |
return 0;
|
|
nsz |
f9d179 |
if (got == want && signbit(got) == signbit(want))
|
|
nsz |
f9d179 |
return dwant;
|
|
nsz |
f9d179 |
if (isinf(got)) {
|
|
nsz |
f9d179 |
got = copysign(0x1p1023, got);
|
|
nsz |
f9d179 |
want *= 0.5;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
return scalbn(got - want, -eulp(want)) + dwant;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
float ulperrl(long double got, long double want, float dwant)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
#if LDBL_MANT_DIG == 53
|
|
nsz |
f9d179 |
return ulperr(got, want, dwant);
|
|
nsz |
f9d179 |
#elif LDBL_MANT_DIG == 64
|
|
nsz |
f9d179 |
if (isnan(got) && isnan(want))
|
|
nsz |
f9d179 |
return 0;
|
|
nsz |
f9d179 |
if (got == want && signbit(got) == signbit(want))
|
|
nsz |
f9d179 |
return dwant;
|
|
nsz |
f9d179 |
if (isinf(got)) {
|
|
nsz |
f9d179 |
got = copysignl(0x1p16383L, got);
|
|
nsz |
f9d179 |
want *= 0.5;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
return scalbn(got - want, -eulpl(want)) + dwant;
|
|
nsz |
f9d179 |
#else
|
|
nsz |
f9d179 |
// TODO
|
|
nsz |
f9d179 |
return inf;
|
|
nsz |
f9d179 |
#endif
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
#define length(a) (sizeof(a)/sizeof*(a))
|
|
nsz |
f9d179 |
#define flag(x) {x, #x}
|
|
nsz |
f9d179 |
static struct {
|
|
nsz |
f9d179 |
int flag;
|
|
nsz |
f9d179 |
char *s;
|
|
nsz |
f9d179 |
} eflags[] = {
|
|
nsz |
f9d179 |
flag(INEXACT),
|
|
nsz |
f9d179 |
flag(INVALID),
|
|
nsz |
f9d179 |
flag(DIVBYZERO),
|
|
nsz |
f9d179 |
flag(UNDERFLOW),
|
|
nsz |
f9d179 |
flag(OVERFLOW)
|
|
nsz |
f9d179 |
};
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
char *estr(int f)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
static char buf[256];
|
|
nsz |
f9d179 |
char *p = buf;
|
|
nsz |
f9d179 |
int i, all = 0;
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
for (i = 0; i < length(eflags); i++)
|
|
nsz |
f9d179 |
if (f & eflags[i].flag) {
|
|
nsz |
f9d179 |
p += sprintf(p, "%s%s", all ? "|" : "", eflags[i].s);
|
|
nsz |
f9d179 |
all |= eflags[i].flag;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
if (all != f) {
|
|
nsz |
f9d179 |
p += sprintf(p, "%s%d", all ? "|" : "", f & ~all);
|
|
nsz |
f9d179 |
all = f;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
p += sprintf(p, "%s", all ? "" : "0");
|
|
nsz |
f9d179 |
return buf;
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
char *rstr(int r)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
switch (r) {
|
|
nsz |
f9d179 |
case RN: return "RN";
|
|
nsz |
f9d179 |
case RZ: return "RZ";
|
|
nsz |
f9d179 |
case RU: return "RU";
|
|
nsz |
f9d179 |
case RD: return "RD";
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
return "R?";
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
void setupfenv(int r)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
fesetround(r);
|
|
nsz |
f9d179 |
feclearexcept(FE_ALL_EXCEPT);
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|
|
nsz |
f9d179 |
int getexcept(void)
|
|
nsz |
f9d179 |
{
|
|
nsz |
f9d179 |
return fetestexcept(INEXACT|INVALID|DIVBYZERO|UNDERFLOW|OVERFLOW);
|
|
nsz |
f9d179 |
}
|
|
nsz |
f9d179 |
|