#include <math.h>
#include <fenv.h>
#include <errno.h>

_MSL_IMP_EXP_C double _MSL_MATH_CDECL nexttoward(double x, long double y) _MSL_CANT_THROW
{
	int increase;
	int positive;
	unsigned long* lx = (unsigned long*)&x;
	unsigned long* ly = (unsigned long*)&y;
#if _MSL_LITTLE_ENDIAN
	unsigned long x_exp = lx[1] & 0x7FF00000;
	if ((ly[1] & 0x7FF00000) == 0x7FF00000 && (ly[1] & 0x000FFFFF) | ly[0] ||
	    x == y)
		return (double)y;
	if (x_exp == 0x7FF00000 && lx[1] & 0x000FFFFF | lx[0])
		return x;
#else
	unsigned long x_exp = lx[0] & 0x7FF00000;
	if ((ly[0] & 0x7FF00000) == 0x7FF00000 && (ly[0] & 0x000FFFFF) | ly[1] ||
	    x == y)
		return (double)y;
	if (x_exp == 0x7FF00000 && lx[0] & 0x000FFFFF | lx[1])
		return x;
#endif
	increase = y > x;
	positive = x > 0;
	if (x_exp != 0x7FF00000 || increase != positive)
	{
		unsigned long long* llx = (unsigned long long*)&x;
		if (x == 0)
		{
			if (increase)
				*llx = 0x0000000000000001ULL;
			else
				*llx = 0x8000000000000001ULL;
		}
		else
		{
			if (increase == positive)
				++(*llx);
			else
				--(*llx);
		}
	#if _MSL_LITTLE_ENDIAN
		x_exp = lx[1] & 0x7FF00000;
	#else
		x_exp = lx[0] & 0x7FF00000;
	#endif
		if (x_exp == 0 && (math_errhandling & MATH_ERREXCEPT))
			feraiseexcept(FE_UNDERFLOW | FE_INEXACT);
		else if (x_exp == 0x7FF00000)
		{
			if (math_errhandling & MATH_ERRNO)
				errno = ERANGE;
			if (math_errhandling & MATH_ERREXCEPT)
				feraiseexcept(FE_OVERFLOW | FE_INEXACT);
		}
	}
	return x;
}

_MSL_IMP_EXP_C double _MSL_MATH_CDECL nextafter(double x, double y) _MSL_CANT_THROW
{
	return nexttoward(x, (long double)y);
}

_MSL_IMP_EXP_C float _MSL_MATH_CDECL nexttowardf(float x, long double y) _MSL_CANT_THROW
{
	int increase;
	int positive;
	unsigned long* lx = (unsigned long*)&x;
	unsigned long* ly = (unsigned long*)&y;
	unsigned long x_exp = *lx & 0x7F800000;
#if _MSL_LITTLE_ENDIAN
	if ((ly[1] & 0x7FF00000) == 0x7FF00000 && (ly[1] & 0x000FFFFF) | ly[0] ||
	    x == y)
		return (float)y;
#else
	if ((ly[0] & 0x7FF00000) == 0x7FF00000 && (ly[0] & 0x000FFFFF) | ly[1] ||
	    x == y)
		return (float)y;
#endif
	if (x_exp == 0x7F800000 && *lx & 0x007FFFFF)
		return x;
	increase = y > x;
	positive = x > 0;
	if (x_exp != 0x7F800000 || increase != positive)
	{
		if (x == 0)
		{
			if (increase)
				*lx = 0x00000001UL;
			else
				*lx = 0x80000001UL;
		}
		else
		{
			if (increase == positive)
				++(*lx);
			else
				--(*lx);
		}
		x_exp = *lx & 0x7F800000;
		if (x_exp == 0 && (math_errhandling & MATH_ERREXCEPT))
			feraiseexcept(FE_UNDERFLOW | FE_INEXACT);
		else if (x_exp == 0x7F800000)
		{
			if (math_errhandling & MATH_ERRNO)
				errno = ERANGE;
			if (math_errhandling & MATH_ERREXCEPT)
				feraiseexcept(FE_OVERFLOW | FE_INEXACT);
		}
	}
	return x;
}

_MSL_IMP_EXP_C float _MSL_MATH_CDECL nextafterf(float x, float y) _MSL_CANT_THROW
{
	return nexttowardf(x, (long double)y);
}
