Blame src/math/gen/util.c

nsz f9d179
#include <stdlib.h>
nsz f9d179
#include <stdio.h>
nsz f9d179
#include <string.h>
nsz f9d179
#include <stdint.h>
nsz f9d179
#include "gen.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
double ulperr(double y, double ycr, double dy)
nsz f9d179
{
nsz f9d179
	return dy + scalbn(ycr - y, -eulp(ycr));
nsz f9d179
}
nsz f9d179
nsz f9d179
char *skipstr(char *buf, char *sep)
nsz f9d179
{
nsz f9d179
	while (*buf && strchr(sep, *buf))
nsz f9d179
		buf++;
nsz f9d179
	return buf;
nsz f9d179
}
nsz f9d179
nsz f9d179
int splitstr(char **a, int n, char *buf, char *sep)
nsz f9d179
{
nsz f9d179
	int i, j;
nsz f9d179
nsz f9d179
	for (i = j = 0; j < n; j++) {
nsz f9d179
		for (; buf[i] && strchr(sep, buf[i]); i++)
nsz f9d179
				buf[i] = 0;
nsz f9d179
		a[j] = buf + i;
nsz f9d179
		if (buf[i] == 0)
nsz f9d179
			break;
nsz f9d179
		for (i++; buf[i] && !strchr(sep, buf[i]); i++);
nsz f9d179
	}
nsz f9d179
	return j;
nsz f9d179
}
nsz f9d179
nsz f9d179
char *dropcomm(char *buf)
nsz f9d179
{
nsz f9d179
	char *p = buf;
nsz f9d179
nsz f9d179
	for (; *p; p++)
nsz f9d179
		if (*p == '/' && p[1] == '/') {
nsz f9d179
			*p = 0;
nsz f9d179
			break;
nsz f9d179
		}
nsz f9d179
	return buf;
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
int econv(int *f, char *s)
nsz f9d179
{
nsz f9d179
	char *a[16];
nsz f9d179
	char *e;
nsz f9d179
	int i,j,k,n;
nsz f9d179
nsz f9d179
	*f = 0;
nsz f9d179
	n = splitstr(a, length(a), s, "|");
nsz f9d179
	for (i = 0; i < n; i++) {
nsz f9d179
		for (j = 0; j < length(eflags); j++)
nsz f9d179
			if (strcmp(a[i], eflags[j].s) == 0) {
nsz f9d179
				*f |= eflags[j].flag;
nsz f9d179
				break;
nsz f9d179
			}
nsz f9d179
		if (j == length(eflags)) {
nsz f9d179
			k = strtol(a[i], &e, 0);
nsz f9d179
			if (*e)
nsz f9d179
				return -1;
nsz f9d179
			*f |= k;
nsz f9d179
		}
nsz f9d179
	}
nsz f9d179
	return 0;
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
int rconv(int *r, char *s)
nsz f9d179
{
nsz f9d179
	if (strcmp(s, "RN") == 0)
nsz f9d179
		*r = RN;
nsz f9d179
	else if (strcmp(s, "RZ") == 0)
nsz f9d179
		*r = RZ;
nsz f9d179
	else if (strcmp(s, "RD") == 0)
nsz f9d179
		*r = RD;
nsz f9d179
	else if (strcmp(s, "RU") == 0)
nsz f9d179
		*r = RU;
nsz f9d179
	else
nsz f9d179
		return -1;
nsz f9d179
	return 0;
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