/*
								+--------------------------------+
								|                                |
								|  ***  Hexadecimal floats  ***  |
								|                                |
								|  Copyright © -tHE SWINe- 2014  |
								|                                |
								|          HexFloat.inl          |
								|                                |
								+--------------------------------+
*/

#pragma once
#ifndef __HEXADECIMAL_FLOATING_POINT_CONVERSION_INLINE_INCLUDED
#define __HEXADECIMAL_FLOATING_POINT_CONVERSION_INLINE_INCLUDED

/**
 *	@file low/HexFloat.inl
 *	@author -tHE SWINe-
 *	@date 2014
 *	@brief hexadecimal floating point number conversions
 */

/*
 *								=== CFloatRepresentationTraits ===
 */

/**
 *	@brief information about floating point representation (specialization for float)
 */
template <>
class CFloatRepresentationTraits<float> {
public:
	enum {
		n_exponent_bit_num = 8, /**< number of exponent bits */
		n_exponent_bias = 127, /**< exponent bias */
		n_exponent_special_low = -0x7f, /**< special exponent low value (used in IEEE-754 to represent zero) */
		n_exponent_special_high = 0x80, /**< special exponent high value (used in IEEE-754 to represent infinity or NaN) */
		n_fraction_bit_num = 23, /**< number of fraction bits */
		n_fraction_print_offset = (4 - n_fraction_bit_num % 4) % 4
	};

	typedef uint32_t TIntType; /**< @brief integer type, matching size of the float representation */

	/**
	 *	@brief converts the number to a hexadecimal string representaion
	 *	@param[in] f is nuber to be converted
	 *	@return Returns pointer to a null-terminated string, containing
	 *		hexadecimal representaion of the given number.
	 *	@note The buffer for the string is allocated statically, therefore
	 *		it is changed by the next call to this function, it should
	 *		not be deleted, and this is not thread-safe.
	 */
	static const char *p_s_XString(float f);
};

/**
 *	@brief information about floating point representation (specialization for double)
 */
template <>
class CFloatRepresentationTraits<double> {
public:
	enum {
		n_exponent_bit_num = 11, /**< number of exponent bits */
		n_exponent_bias = 1023, /**< exponent bias */
		n_exponent_special_low = -0x3ff, /**< special exponent low value (used in IEEE-754 to represent zero) */
		n_exponent_special_high = 0x400, /**< special exponent high value (used in IEEE-754 to represent infinity or NaN) */
		n_fraction_bit_num = 52, /**< number of fraction bits */
		n_fraction_print_offset = (4 - n_fraction_bit_num % 4) % 4
	};

	typedef uint64_t TIntType; /**< @brief integer type, matching size of the float representation */

	/**
	 *	@brief converts the number to a hexadecimal string representaion
	 *	@param[in] f is nuber to be converted
	 *	@return Returns pointer to a null-terminated string, containing
	 *		hexadecimal representaion of the given number.
	 *	@note The buffer for the string is allocated statically, therefore
	 *		it is changed by the next call to this function, it should
	 *		not be deleted, and this is not thread-safe.
	 */
	static const char *p_s_XString(double f);
};

/*
 *								=== ~CFloatRepresentationTraits ===
 */

/*
 *								=== CFloatUtils ===
 */

#if defined(_MSC_VER) && !defined(__MWERKS__) && !defined(for) && _MSC_VER <= 1200
#define for if(0) {} else for
#endif // _MSC_VER && !__MWERKS__ && !for && _MSC_VER <= 1200
// msvc 'for' scoping hack

template <class CType>
typename CFloatUtils<CType>::TIntType CFloatUtils<CType>::n_Exponent_Mask()
{
	return ((TIntType(1) << n_exponent_bit_num) - 1) << n_fraction_bit_num;
}

template <class CType>
typename CFloatUtils<CType>::TIntType CFloatUtils<CType>::n_Fraction_Mask()
{
	return (TIntType(1) << n_fraction_bit_num) - 1;
}

template <class CType>
typename CFloatUtils<CType>::TIntType CFloatUtils<CType>::n_Mantissa_One()
{
	return TIntType(1) << n_fraction_bit_num;
}

template <class CType>
int CFloatUtils<CType>::n_Get_Exponent(TType f)
{
	UFloatInt uf;
	uf.f = f;
	return int((uf.n & n_Exponent_Mask()) >> n_fraction_bit_num) - n_exponent_bias;
}

template <class CType>
int CFloatUtils<CType>::n_Get_SignBit(TType f)
{
	UFloatInt uf;
	uf.f = f;
	return int(uf.n >> n_sign_bit_shift);
}

template <class CType>
typename CFloatUtils<CType>::TIntType CFloatUtils<CType>::n_Get_Fraction(TType f)
{
	UFloatInt uf;
	uf.f = f;
	return uf.n & n_Fraction_Mask();
}

template <class CType>
typename CFloatUtils<CType>::TIntType CFloatUtils<CType>::n_Get_Mantissa(TType f)
{
	UFloatInt uf;
	uf.f = f;
	int n_exponent = int((uf.n & n_Exponent_Mask()) >> n_fraction_bit_num) - n_exponent_bias;
	CFloatUtils<CType>::TIntType n_fraction = uf.n & n_Fraction_Mask();
	if(n_exponent > n_exponent_special_low)
		return n_fraction | n_Mantissa_One(); // normal
	else
		return n_fraction; // denormal
}

template <class CType>
typename CFloatUtils<CType>::TType CFloatUtils<CType>::f_MakeFloat(bool b_negative_sign, TIntType n_mantissa, int n_exponent)
{
	_ASSERTE(n_exponent + n_exponent_bias >= 0 &&
		((n_exponent + n_exponent_bias) >> n_exponent_bit_num) == 0);
	// make sure the exponent fits (does not care about the special values)

	_ASSERTE((n_mantissa & ~n_Fraction_Mask()) == n_Mantissa_One()); // todo - this is likely a bug, if one wants to make a denormal number!
	// the high bit of mantissa must be set, the mantissa must not be bigger

	UFloatInt di;
	di.n = ((b_negative_sign)? TIntType(1) << n_sign_bit_shift : 0) |
		(((TIntType)(n_exponent + n_exponent_bias)) << n_fraction_bit_num) |
		(n_mantissa & n_Fraction_Mask());
	// make float

	return di.f;
}

template <class CType>
typename CFloatUtils<CType>::TType CFloatUtils<CType>::f_Infinity()
{
	return f_MakeFloat(false, n_Mantissa_One(), n_exponent_special_high);
}

template <class CType>
typename CFloatUtils<CType>::TType CFloatUtils<CType>::f_Not_a_Number()
{
	return f_MakeFloat(true, n_Mantissa_One() | n_Fraction_Mask(), n_exponent_special_high);
}

template <class CType>
typename CFloatUtils<CType>::TType CFloatUtils<CType>::f_ParseXError()
{
	const int n_shift = 4 - n_fraction_bit_num % 4; // so that when printed as hexadecimal, it would spell correctly
	return f_MakeFloat(false, ((((sizeof(CType) > 4)? ((uint64_t(0x15deadde) << 28) | 0xaddead0) :
		0x15dead0) >> n_shift) & n_Fraction_Mask()) | n_Mantissa_One(), n_exponent_special_low);
}

template <class CType>
bool CFloatUtils<CType>::b_Is_NaN(TType f)
{
	return n_Get_Exponent(f) == n_exponent_special_high && n_Get_Fraction(f) != 0;
}

template <class CType>
bool CFloatUtils<CType>::b_Is_Finite(TType f)
{
	return n_Get_Exponent(f) != n_exponent_special_high;
}

template <class CType>
bool CFloatUtils<CType>::b_Is_Infinity(TType f)
{
	return f == f_Infinity();
}

template <class CType>
bool CFloatUtils<CType>::b_Is_NegativeInfinity(TType f)
{
	return f == -f_Infinity();
}

template <class CType>
typename CFloatUtils<CType>::TType CFloatUtils<CType>::f_ParseXFloat(const char *p_s_float)
{
	TType f_resutlt;
	if(!ParseXFloat(f_resutlt, p_s_float))
		return f_ParseXError();
	return f_resutlt;
}

template <class CType>
bool CFloatUtils<CType>::ParseXFloat(TType &r_f_result, const char *p_s_float)
{
	CFloatCommon::TParseInfo t_info = CFloatCommon::t_ParseXLiteral(p_s_float);
	if(!t_info.b_parsed)
		return false;
	// parse the literal (common part for doubles and floats)

	TIntType n_mantissa;
	{
		int n_highest_mantissa_bit = CFloatCommon::n_HighestBit_Set(t_info.n_mantissa);
		t_info.n_exponent += n_highest_mantissa_bit - 1 - t_info.n_mantissa_dp_position;
		// count number of bits above the decimal point (may be negative)
		// note that the bounds of exponent will be checked in f_MakeDouble()

		if(n_highest_mantissa_bit > n_mantissa_bit_num) {
			n_mantissa = TIntType(t_info.n_mantissa >> (n_highest_mantissa_bit - n_mantissa_bit_num));
			// data lost, least significant digits shifted out to the right
		} else if(n_highest_mantissa_bit < n_mantissa_bit_num) {
			n_mantissa = TIntType(t_info.n_mantissa << (n_mantissa_bit_num - n_highest_mantissa_bit));
			// no data lost, just zeroes shifted in from right
		} else
			n_mantissa = TIntType(t_info.n_mantissa); // !!
		// the position of the highest (1.0) bit of mantissa will be checked again in f_MakeDouble()
	}
	// adjust mantissa shift, compensate in exponent

	{
		if(!(t_info.n_exponent + n_exponent_bias >= 0 &&
		   ((t_info.n_exponent + n_exponent_bias) >> n_exponent_bit_num) == 0))
			return false;
		// make sure the exponent fits (does not care about the special values)

		if((n_mantissa & ~n_Fraction_Mask()) != n_Mantissa_One())
			return false;
		// the high bit of mantissa must be set, the mantissa must not be bigger
	}
	// make the checks in f_MakeFloat() runtime checks

	r_f_result = f_MakeFloat(t_info.b_sign_bit, n_mantissa, t_info.n_exponent);
	// put it all together

	return true;
}

template <class CType>
typename CFloatUtils<CType>::TType CFloatUtils<CType>::f_ULP(TType f)
{
	bool b_sign = (n_Get_SignBit(f) != 0);
	int n_exponent = n_Get_Exponent(f) - n_fraction_bit_num;
	return (n_exponent + n_exponent_bias >= 0)? f_MakeFloat(b_sign, n_Mantissa_One(), n_exponent) : 0; // handle underflow
}

template <class CType>
bool CFloatUtils<CType>::b_TieBreak_ToEven()
{
	/*
	//										 <- 16B double ->
	//											<- fraction->
	const double _special = f_ParseXDouble("0x0.00000000000008p+0"); // one, at the position one past the LSB
	const double _oddone =  f_ParseXDouble("0x1.0000000000001p+0"); // one, ending with a single one at LSB
	const double _evenone = f_ParseXDouble("0x1.0000000000002p+0"); // one, ending with a single one to the left of LSB
	*/
	const TType special = f_ULP(1) / 2; // one, at the position one past the LSB
	const TType oddone = 1 + f_ULP(1); // one, ending with a single one at LSB
	const TType evenone = 1 + 2 * f_ULP(1); // one, ending with a single one to the left of LSB
	const TType one = 1;

	volatile TType v;
	v = one; v += special;
	if(v != one)
		return false;
	v = oddone; v += special;
	if(v != evenone) // odd + half rounds to even
		return false;
	v = evenone; v += special;
	if(v != evenone) // even + half rounds to the same even
		return false;
	v = -one; v -= special;
	if(v != -one)
		return false;
	v = -oddone; v -= special;
	if(v != -evenone) // -odd - half rounds to -even
		return false;
	v = -evenone; v -= special;
	if(v != -evenone) // -even - half rounds to the same -even
		return false;
	v -= special; // odd number of ops to detect alternate mode
	if(v != -evenone) // -even - half rounds to the same -even
		return false;

	return true;
}

template <class CType>
bool CFloatUtils<CType>::b_TieBreak_AwayFromZero()
{
	/*
	//										 <- 16B double ->
	//											<- fraction->
	const double _special = f_ParseXDouble("0x0.00000000000008p+0"); // one, at the position one past the LSB
	const double _oddone =  f_ParseXDouble("0x1.0000000000001p+0"); // one, ending with a single one at LSB
	const double _evenone = f_ParseXDouble("0x1.0000000000002p+0"); // one, ending with a single one to the left of LSB
	*/
	const TType special = f_ULP(1) / 2; // one, at the position one past the LSB
	const TType oddone = 1 + f_ULP(1); // one, ending with a single one at LSB
	const TType evenone = 1 + 2 * f_ULP(1); // one, ending with a single one to the left of LSB
	const TType one = 1;

	volatile TType v;
	v = one; v += special;
	if(v != oddone) // 1.0 + half rounds to odd
		return false;
	v = oddone; v += special;
	if(v != evenone) // odd + half rounds to even
		return false;
	v = evenone; v += special;
	if(v <= evenone) // even + half rounds to greater than even
		return false;
	v = -one; v -= special;
	if(v != -oddone) // -1.0 - half rounds to -odd
		return false;
	v = -oddone; v -= special;
	if(v != -evenone) // -odd - half rounds to -even
		return false;
	v = -evenone; v -= special;
	if(v >= -evenone) // -even - half rounds to smaller than -even
		return false;
	v -= special; // odd number of ops to detect alternate mode
	if(v >= -evenone) // -even - half rounds to smaller than -even
		return false;

	return true;
}

template <class CType>
int CFloatUtils<CType>::n_FPU_RoundMode()
{
	const TType one = 1;
	const TType special = f_ULP(1) / 256; // much smaller than the least representable digit for 1.0

	int n_comparison = 0, n_prev;
	for(int i = 0; i < 10; ++ i) {
		volatile TType v;
		v = one; v += special;
		if(v == one)
			n_comparison |= 1;
		v = one; v -= special;
		if(v == one)
			n_comparison |= 2;
		v = -one; v += special;
		if(v == -one)
			n_comparison |= 4;
		v = -one; v -= special;
		if(v == -one)
			n_comparison |= 8;
		v = -one; v -= special; // odd number of operations to detect alternate mode
		if(v == -one)
			n_comparison |= 8;
		if(!i)
			n_prev = n_comparison;
		else if(n_comparison != n_prev)
			return CFloatCommon::fmode_Stochastic; // if not the same results, it must be stochastic / alternate
	}

	switch(n_comparison) {
	case 15:
		{
			int n_ties = (b_TieBreak_ToEven())? ((b_TieBreak_AwayFromZero())?
				CFloatCommon::fmode_Nearest_TiesUnknown : CFloatCommon::fmode_Nearest_TiesToEven) :
				((b_TieBreak_AwayFromZero())? CFloatCommon::fmode_Nearest_TiesAway :
				CFloatCommon::fmode_Nearest_TiesUnknown);
			for(int i = 0; i < 10; ++ i) {
				int n_ties1 = (b_TieBreak_ToEven())? ((b_TieBreak_AwayFromZero())?
					CFloatCommon::fmode_Nearest_TiesUnknown : CFloatCommon::fmode_Nearest_TiesToEven) :
					(b_TieBreak_AwayFromZero())? CFloatCommon::fmode_Nearest_TiesAway :
					CFloatCommon::fmode_Nearest_TiesUnknown;
				if(n_ties1 != n_ties)
					return CFloatCommon::fmode_Nearest_TiesStochastic;
			}
			return n_ties;
		}
		break; // to get rid of warnings on some compilers
	case 10:
		return CFloatCommon::fmode_TowardsPosInf;
	case 5:
		return CFloatCommon::fmode_TowardsNegInf;
	case 9:
		return CFloatCommon::fmode_TowardsZero;
	}
	return CFloatCommon::fmode_Unknown;
}

#ifdef __HEXFLOAT_UNIT_TESTS

#define FloatUtils__RuntimeAssert(b_exp) \
	do { \
		if(!(b_exp)) { \
			RT_AssertHandler(#b_exp, __LINE__, __FILE__); \
			++ n_error_num; \
		} \
	} while(false)
// want it in release as well

template <class CType>
bool CFloatUtils<CType>::UnitTests()
{
	int n_error_num = 0;

	const int p_special[] = {
		0, 1, 2, 3, 5, 9, 17, 33, 65, 129, 257, 513
	};
	const char *p_s_special[] = {
		0, "0x1.0p+0", "0x1.0p+1", "0x1.8p+1", "0x1.4p+2", "0x1.2p+3", "0x1.1p+4",
		"0x1.08p+5", "0x1.04p+6", "0x1.02p+7", "0x1.01p+8", "0x1.008p+9"
	};
	for(size_t i = 0, n = sizeof(p_special) / sizeof(p_special[0]); i < n; ++ i) {
		TType f = (TType)p_special[i]; // int to float triggers warnings
		const char *p_s_string = TTraits::p_s_XString(f);
		if(p_s_special[i])
			FloatUtils__RuntimeAssert(!strcmp(p_s_special[i], p_s_string)); // make sure the string is printed correctly
		TType p = f_ParseXFloat(p_s_string);
		FloatUtils__RuntimeAssert(p == f); // make sure we get back the same value
	}
	for(int i = 0; i < 1000; ++ i) {
		TType f = ((TType)rand()) / RAND_MAX * 100 - 50;
		const char *p_s_string = TTraits::p_s_XString(f);
		TType p = f_ParseXFloat(p_s_string);
		FloatUtils__RuntimeAssert(p == f); // make sure we get back the same value
	}
	// test printing and parsing

	FloatUtils__RuntimeAssert(f_ParseXFloat("0x1") == 1);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x1.0") == 1);
	FloatUtils__RuntimeAssert(f_ParseXFloat("+0x1") == 1);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x1.0p+0") == 1);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x0.1p+4") == 1);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x10.0p-4") == 1);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x1.0p+1") == 2);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x10.0p-3") == 2);
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x1.0p-1") == TType(.5));
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x0.1p+3") == TType(.5));
	//FloatUtils__RuntimeAssert(f_ParseXFloat("0x1.0p-7f") == 0);
	//float fdenormal = f_ParseXFloat("0x1.123456p-7f"); // well, hard to check
	TType fnot_a_number = f_Not_a_Number();
	//FloatUtils__RuntimeAssert(b_Is_NaN(fnot_a_number = f_ParseXFloat("0x1.123456p+80")));
	FloatUtils__RuntimeAssert(fnot_a_number != fnot_a_number);
	TType fpositive_infinity = f_Infinity();//f_ParseXFloat("0x1.0p+80");
	//FloatUtils__RuntimeAssert(fpositive_infinity > FLT_MAX);
	FloatUtils__RuntimeAssert(b_Is_Infinity(fpositive_infinity));
	TType fnegative_infinity = -f_Infinity();//f_ParseXFloat("-0x1.0p+80");
	//FloatUtils__RuntimeAssert(fnegative_infinity < -FLT_MAX);
	FloatUtils__RuntimeAssert(b_Is_NegativeInfinity(fnegative_infinity));
	FloatUtils__RuntimeAssert(f_ParseXFloat("thrash") == f_ParseXError()); // no "0x"
	FloatUtils__RuntimeAssert(f_ParseXFloat("0xtrash") == f_ParseXError()); // no digits
	FloatUtils__RuntimeAssert(f_ParseXFloat("--0x1.0") == f_ParseXError()); // double sign
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x0.0") == f_ParseXError()); // zero mantissa
	FloatUtils__RuntimeAssert(f_ParseXFloat("0x1.000000000000000fp+0") == f_ParseXError()); // more than 64-bit mantissa
	// simple tests

	FloatUtils__RuntimeAssert(b_Is_NaN(f_Not_a_Number()));
	FloatUtils__RuntimeAssert(b_Is_Infinity(f_Infinity()));
	FloatUtils__RuntimeAssert(b_Is_NegativeInfinity(-f_Infinity()));
	// specials tests

	const TType p_specials[] = {
		0, 1, f_Infinity(), -f_Infinity()
	};
	for(int i = 0; i < 1000; ++ i) {
		volatile TType f;
		if(i < sizeof(p_specials) / sizeof(p_specials[0]))
			f = p_specials[i];
		else
			f = ((TType)rand()) / RAND_MAX * 100 - 50;
		FloatUtils__RuntimeAssert(f_MakeFloat(n_Get_SignBit(f) != 0, n_Get_Mantissa(f), n_Get_Exponent(f)) == f);
	}
	// test de-assembling floats

	for(int i = 0; i < 1000; ++ i) {
		volatile TType f = ((TType)rand()) / RAND_MAX * 100 - 50;
		if(i > 100)
			f = TType(exp(double(f / 100) * 30)); // make the range bigger (but only up to 10^30 to not cause infinities)
		volatile TType f_ulp = f_ULP(f);
		FloatUtils__RuntimeAssert(f_ulp < 0 == f < 0); // sign of ULP is the same as that of f
		FloatUtils__RuntimeAssert(!f || ((f_ulp < 0)? -f_ulp : f_ulp) < ((f < 0)? -f : f)); // ULP is smaller than f
		volatile TType fp = f, fm = f; fp += f_ulp; fm -= f_ulp;
		FloatUtils__RuntimeAssert(fp != f); // adding one ULP changes the value
		FloatUtils__RuntimeAssert(fm != f); // subtracting one ULP changes the value
		volatile TType f_step = TType(f + TType(f_ulp / 2)); f_step -= f;
		FloatUtils__RuntimeAssert(f_step == 0 || f_step == f_ulp); // adding less than half ULP is too small to change the value, if rounding toward nearest
		volatile TType f_step2 = TType(f + TType(f_ulp / 3)); f_step2 -= f;
		FloatUtils__RuntimeAssert(f_step2 == 0 || f_step2 == f_ulp); // adding half ULP is too small to change the value, if rounding toward zero or infinities (half ULP triggers tiebreaking and may not round to the same number)
	}
	// test the ULP function

	return !n_error_num;
}

template <class CType>
void CFloatUtils<CType>::RT_AssertHandler(const char *p_s_expr, int n_line, const char *p_s_file)
{
	fprintf(stderr, "=== Assertion failed ===\n"
		"expression: '%s'\nfile: '%s'\nline: %d\n",
		p_s_expr, p_s_file, n_line);
	// prints message header

#if defined(_DEBUG) && defined(_MSC_VER) && !defined(__MWERKS__)
	__asm {
		int 3
	};
#endif // _DEBUG && _MSC_VER && !__MWERKS__
	// do not abort, though ...
}

#undef FloatUtils__RuntimeAssert

#endif // __HEXFLOAT_UNIT_TESTS

/*
 *								=== ~CFloatUtils ===
 */

#endif // !__HEXADECIMAL_FLOATING_POINT_CONVERSION_INLINE_INCLUDED
