/* By Pedro Gimeno Fortea, donated to the public domain */

#include <assert.h>
#include <math.h>
#include <stdio.h>


double
rint(double x)
{
	/* Perform nearest-or-even rounding. */
	double r;
	if (x <= 0.0)
	{
		/* Positive/Negative zero should be returned unchanged. */
		if (x == 0.0)
			return x;
		/*
		 * Subtracting 0.5 from a number too close to -0.5 can
		 * round up to exactly -1.0, producing incorrect results,
		 * so this takes the opposite approach: to add 0.5 to the
		 * negative number, so that it goes closer to zero (or at
		 * most to +0.5, which is dealt with later), avoiding the
		 * precision issue.
		 *
		 * r is the result of the rounding.
		 */
		r = floor(x += 0.5);
		/*
		 * If the rounding did not produce exactly input+0.5 then
		 * we're mostly done. Just be careful to return minus zero
		 * when input+0.5 >= 0, as that's what rint should return
		 * with negative input.
		 */
		if (r != x || x == 0.0)
			return x >= 0.0 ? -0.0 : r;
		/*
		 * We're in the case where the original fractional part was
		 * exactly 0.5 (detected with floor(input+0.5) == input+0.5).
		 * We need to round to even. Dividing input+0.5 by 2, taking
		 * the floor and multiplying by 2 yields the closest even
		 * number. This part assumes that division by 2 is exact
		 * (no underflow possible here).
		 */
		return floor(x * 0.5) * 2.0;
	}

	/*
	 * The positive case is similar but using ceil instead of floor,
	 * except that the case input-0.5 == 0 does not need to be checked,
	 * as the correct sign is already produced.
	 */
	r = ceil(x -= 0.5);
	if (r != x)
		return x <= 0.0 ? 0.0 : r;
	return ceil(x * 0.5) * 2.0;
}



int test_equal_z(double a, double b)
{
	char a1[20];
	char b1[20];
	sprintf(a1, "%+g", a);
	sprintf(b1, "%+g", b);
	return a == b && a1[0] == b1[0];
}

int main()
{
	assert(rint( 2.5) ==  2.0);
	assert(rint( 2.4) ==  2.0);
	assert(rint( 2.6) ==  3.0);
	assert(rint( 1.5) ==  2.0);
	assert(rint( 1.4) ==  1.0);
	assert(rint( 1.6) ==  2.0);
	assert(rint(-2.5) == -2.0);
	assert(rint(-2.4) == -2.0);
	assert(rint(-2.6) == -3.0);
	assert(rint(-1.5) == -2.0);
	assert(rint(-1.4) == -1.0);
	assert(rint(-1.6) == -2.0);

	/* Test minus zero */
	assert(test_equal_z(rint( 0.0),  0.0));
	assert(test_equal_z(rint(-0.0), -0.0));
	assert(test_equal_z(rint( 0.5),  0.0));
	assert(test_equal_z(rint(-0.5), -0.0));

	/* Corner cases for IEEE double */
	/* Last float having fractional part */
	assert(rint(-4503599627370495.5) == -4503599627370496.0);
	assert(rint( 4503599627370495.5) ==  4503599627370496.0);
	assert(rint(-4503599627370494.5) == -4503599627370494.0);
	assert(rint( 4503599627370494.5) ==  4503599627370494.0);
	/* abs value strictly greater than 0.5 should be rounded away from 0 */
	/* (these numbers are 0.5 + 1 ulp) */
	assert(rint( 0.5000000000000001) ==  1.0);
	assert(rint(-0.5000000000000001) == -1.0);
	/* abs value strictly less than 0.5 should be rounded towards 0 */
	/* (these values plus 1 ulp equal 0.5) */
	assert(test_equal_z(rint( 0.49999999999999994),  0.0));
	assert(test_equal_z(rint(-0.49999999999999994), -0.0));
	/* denormals */
	assert(test_equal_z(rint( 5e-324),  0.0));
	assert(test_equal_z(rint(-5e-324), -0.0));
	assert(test_equal_z(rint( 1e-322),  0.0));
	assert(test_equal_z(rint(-1e-322), -0.0));

	return 0;
}