/*
 * Floating point functions that are difficult or impossible to implement
 * in pure Haskell.
 *
 * Copyright (C) 2009-2010 Nick Bowler.
 *
 * License BSD2:  2-clause BSD license.  See LICENSE for full terms.
 * This is free software: you are free to change and redistribute it.
 * There is NO WARRANTY, to the extent permitted by law.
 */
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <fenv.h>

#include "cfloat.h"

#pragma STDC FENV_ACCESS ON

int double_format(char *buf, char spec, int precision, double val)
{
	char fmt[] = "%.*f";
	fmt[3] = spec;

	if (buf == NULL)
		return snprintf(NULL, 0, fmt, precision, val);
	return sprintf(buf, fmt, precision, val);
}

double double_signum(double val)
{
	if (signbit(val))
		return -1;
	return 1;
}

float float_signum(float val)
{
	if (signbit(val))
		return -1;
	return 1;
}

int double_classify(double val)
{
	switch (fpclassify(val)) {
	case FP_INFINITE:
		return 0;
	case FP_NAN:
		return 1;
	case FP_NORMAL:
		return 2;
	case FP_SUBNORMAL:
		return 3;
	case FP_ZERO:
		return 4;
	}

	return -1;
}

int float_classify(float val)
{
	switch (fpclassify(val)) {
	case FP_INFINITE:
		return 0;
	case FP_NAN:
		return 1;
	case FP_NORMAL:
		return 2;
	case FP_SUBNORMAL:
		return 3;
	case FP_ZERO:
		return 4;
	}

	return -1;
}

int double_compare(double a, double b)
{
	if (isless(a, b))
		return 0;
	if (a == b)
		return 1;
	if (isgreater(a, b))
		return 2;
	if (isunordered(a, b))
		return 3;
	return -1;
}

int float_compare(float a, float b)
{
	if (isless(a, b))
		return 0;
	if (a == b)
		return 1;
	if (isgreater(a, b))
		return 2;
	if (isunordered(a, b))
		return 3;
	return -1;
}

int set_roundmode(int mode)
{
	int cmode;

	switch (mode) {
	case 0:
		cmode = FE_TONEAREST;
		break;
	case 1:
		cmode = FE_UPWARD;
		break;
	case 2:
		cmode = FE_DOWNWARD;
		break;
	case 3:
		cmode = FE_TOWARDZERO;
		break;
	default:
		return -1;
	}

	return fesetround(cmode);
}

int get_roundmode(void)
{
	int cmode = fegetround();

	switch (cmode) {
	case FE_TONEAREST:
		return 0;
	case FE_UPWARD:
		return 1;
	case FE_DOWNWARD:
		return 2;
	case FE_TOWARDZERO:
		return 3;
	default:
		return -1;
	}
}

int fenv_restore(fenv_t *env, unsigned *excepts)
{
	int raw_excepts = fetestexcept(FE_ALL_EXCEPT);

	if (excepts) {
		*excepts = 0;

#ifdef FE_DIVBYZERO
		if (raw_excepts & FE_DIVBYZERO) *excepts |= 0x01;
#endif
#ifdef FE_INEXACT
		if (raw_excepts & FE_INEXACT)   *excepts |= 0x02;
#endif
#ifdef FE_INVALID
		if (raw_excepts & FE_INVALID)   *excepts |= 0x04;
#endif
#ifdef FE_OVERFLOW
		if (raw_excepts & FE_OVERFLOW)  *excepts |= 0x08;
#endif
#ifdef FE_UNDERFLOW
		if (raw_excepts & FE_UNDERFLOW) *excepts |= 0x10;
#endif
	}

	return fesetenv(env);
}

int fenv_raise_excepts(unsigned excepts)
{
	int raw_excepts = 0;

#ifdef FE_DIVBYZERO
	if (excepts & 0x01) raw_excepts |= FE_DIVBYZERO;
#endif
#ifdef FE_INEXACT
	if (excepts & 0x02) raw_excepts |= FE_INEXACT;
#endif
#ifdef FE_INVALID
	if (excepts & 0x04) raw_excepts |= FE_INVALID;
#endif
#ifdef FE_OVERFLOW
	if (excepts & 0x08) raw_excepts |= FE_OVERFLOW;
#endif
#ifdef FE_UNDERFLOW
	if (excepts & 0x10) raw_excepts |= FE_UNDERFLOW;
#endif
	return feraiseexcept(raw_excepts);
}