/* 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.
		 *
		 * For too big numbers the float may have no decimals.
		 * That case is detected with x+1.0 == x+0.5; when that
		 * happens, the number is returned unchanged.
		 *
		 * r is the result of the rounding.
		 */
		r = x;
		x += 0.5;
		if (x == r + 1.0)
			return r;
		r = floor(x);
		/*
		 * 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 with signs inverted and 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 = x;
	x -= 0.5;
	if (x == r - 1.0)
		return r;
	r = ceil(x);
	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));

	/* bigger numbers */
	/* 2^52, +1, +2, +3 */
	assert(rint( 4503599627370496.0) ==  4503599627370496.0);
	assert(rint(-4503599627370496.0) == -4503599627370496.0);
	assert(rint( 4503599627370497.0) ==  4503599627370497.0);
	assert(rint(-4503599627370497.0) == -4503599627370497.0);
	assert(rint( 4503599627370498.0) ==  4503599627370498.0);
	assert(rint(-4503599627370498.0) == -4503599627370498.0);
	assert(rint( 4503599627370499.0) ==  4503599627370499.0);
	assert(rint(-4503599627370499.0) == -4503599627370499.0);

	/* 2^53-2, -1, +0, +2 */
	assert(rint( 9007199254740990.0) ==  9007199254740990.0);
	assert(rint(-9007199254740990.0) == -9007199254740990.0);
	assert(rint( 9007199254740991.0) ==  9007199254740991.0);
	assert(rint(-9007199254740991.0) == -9007199254740991.0);
	assert(rint( 9007199254740992.0) ==  9007199254740992.0);
	assert(rint(-9007199254740992.0) == -9007199254740992.0);
	assert(rint( 9007199254740994.0) ==  9007199254740994.0);
	assert(rint(-9007199254740994.0) == -9007199254740994.0);
	/* 2^54-2, +0, +4 */
	assert(rint( 18014398509481982.0) ==  18014398509481982.0);
	assert(rint(-18014398509481982.0) == -18014398509481982.0);
	assert(rint( 18014398509481984.0) ==  18014398509481984.0);
	assert(rint(-18014398509481984.0) == -18014398509481984.0);
	assert(rint( 18014398509481988.0) ==  18014398509481988.0);
	assert(rint(-18014398509481988.0) == -18014398509481988.0);
	/* 2^55-4 */
	assert(rint( 36028797018963964.0) ==  36028797018963964.0);
	assert(rint(-36028797018963964.0) == -36028797018963964.0);

	return 0;
}