Implementing C99's roundf(), round(), and roundl()

[David Schultz]

On Tue, Dec 02, 2003, Bruce Evans wrote:
> On Mon, 1 Dec 2003, Steve Kargl wrote:
> 
> > On Mon, Dec 01, 2003 at 07:05:18PM +1100, Bruce Evans wrote:
> > > All the other corner cases need to be checked.  It's possibly to check
> > > all 2^32 cases for floats (once you know the correct results).
> >
> > Do you have code to do this check?
> >
> > > Other things to check: setting of exception flags.  I'm not sure if the
> > > settings by ceil() are the right ones and the only ones.
> 
> I thought of a good way after righting the above: roundf() is supposed to
> be equivalent to rintf() with certain rounding, so set the rounding mode
> and call rintf() to determine the correct value.  Unfortunately there
> might not be a correct rounding mode (what does round to nearest do
> for integers?  I think it rounds to even, but roundf() requires rounding
> up half-way cases).

You're right on the money.  The C99 standard says:

	The double version of round behaves as though implemented by
	
	#include <math.h>
	#include <fenv.h>
	#pragma STDC FENV_ACCESS ON
	double round(double x)
	{
		double result;
		fenv_t save_env;
		feholdexcept(&save_env);
		result = rint(x);
		if (fetestexcept(FE_INEXACT)) {
			fesetround(FE_TOWARDZERO);
			result = rint(copysign(0.5 + fabs(x), x));
		}
		feupdateenv(&save_env);
		return result;
	}

	The round functions may, but are not required to, raise the
	inexact exception for noninteger numeric arguments, as this
	implementation does.

This should be a faster implementation in that it only has one
branch.  The one with ceil() has a branch for the inf/NaN test,
another that I added to avoid spurious inexact exceptions, and a
third to test the sign of the input.

However, I think we may be able to use floor() instead of
fesetround()/rint()/feupdateenv() and do something similar without
changing the rounding mode:

	double
	round(double x)
	{
		double result;
		fp_except_t fe;

		fe = fpresetsticky(0);
		result = rint(x);
		if (fpgetsticky() & FP_X_IMP) {
			result = copysign(floor(fabs(x) + 0.5), x);
			fe |= FP_X_IMP;
		}
		fpresetsticky(fe);
		return (result);
	}

Does this seem reasonable?

By the way, the man pages refer to fpresetsticky(3), but not
fpsetsticky(3).  The MI ieee.h header declares only fpsetsticky(),
but some architectures override these with equivalent macros for
both.  If you have a good idea of how this is supposed to be, it
would be great if you could fix it.

Another option is to just implement the relevant parts of fenv.h.
They are fairly straightforward, and there's no sense in hiding
behind the brokenness of gcc's i386 backend forever.  (By the way,
if we did implement fenv.h, does anyone know if there are any
legal issues or attribution requirements involved in using the
exact round() code proposed by the C standard?)