Blame src/math/util.c

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;
Szabolcs Nagy c58b8e
	if (got == want) {
Szabolcs Nagy c58b8e
		if (signbit(got) == signbit(want))
Szabolcs Nagy c58b8e
			return dwant;
Szabolcs Nagy c58b8e
		return inf;
Szabolcs Nagy c58b8e
	}
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;
Szabolcs Nagy c58b8e
	if (got == want) {
Szabolcs Nagy c58b8e
		if (signbit(got) == signbit(want))
Szabolcs Nagy c58b8e
			return dwant;
Szabolcs Nagy c58b8e
		return inf; // treat 0 sign errors badly
Szabolcs Nagy c58b8e
	}
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;
Szabolcs Nagy c58b8e
	if (got == want) {
Szabolcs Nagy c58b8e
		if (signbit(got) == signbit(want))
Szabolcs Nagy c58b8e
			return dwant;
Szabolcs Nagy c58b8e
		return inf;
Szabolcs Nagy c58b8e
	}
nsz f9d179
	if (isinf(got)) {
nsz f9d179
		got = copysignl(0x1p16383L, got);
nsz f9d179
		want *= 0.5;
nsz f9d179
	}
Szabolcs Nagy b8fbda
	return scalbnl(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 c5feb4
#ifdef FE_TOWARDZERO
nsz f9d179
	case RZ: return "RZ";
nsz c5feb4
#endif
nsz c5feb4
#ifdef FE_UPWARD
nsz f9d179
	case RU: return "RU";
nsz c5feb4
#endif
nsz c5feb4
#ifdef FE_DOWNWARD
nsz f9d179
	case RD: return "RD";
nsz c5feb4
#endif
nsz f9d179
	}
nsz f9d179
	return "R?";
nsz f9d179
}