Date: Mon, 1 Dec 2003 19:05:18 +1100 (EST) From: Bruce Evans <bde@zeta.org.au> To: David Schultz <das@freebsd.org> Cc: Steve Kargl <sgk@troutmask.apl.washington.edu> Subject: Re: Implementing C99's roundf(), round(), and roundl() Message-ID: <20031201182219.O4431@gamplex.bde.org> In-Reply-To: <20031130213951.GA37082@VARK.homeunix.com> References: <20031129000133.GA30662@troutmask.apl.washington.edu> <20031129163105.GA32651@troutmask.apl.washington.edu> <20031130213951.GA37082@VARK.homeunix.com>
index | next in thread | previous in thread | raw e-mail
On Sun, 30 Nov 2003, David Schultz wrote:
> On Sat, Nov 29, 2003, Steve Kargl wrote:
> > On Sat, Nov 29, 2003 at 12:09:11AM -0800, David Schultz wrote:
> > > On Fri, Nov 28, 2003, Steve Kargl wrote:
> > > > Can the math functions round[fl]() be implemented in
> > > > terms of other math(3) functions and still conform to the
> > > > C99 and POSIX standard? For example,
> > > > [code moved later]
> > >
> > > This looks correct to me at first glance, modulo possible problems
> > > with overflow. ...
> >
> > I don't undrestand your overflow comment. ceil[f]() can return Inf
> > and nan, but in those cases round[f]() should also return Inf and nan.
> > The two operations, (t-x) and (t+x), should yield a value in the
> > range [0,1). I'll submit a PR with a man page.
>
> The concern was that ceil() could round a number up to infinity
> when round() is supposed to round the number down. But now that I
> think about it at a reasonable hour, this concern is clearly
> bogus. In base two floating point representations, there isn't
> enough precision to get numbers that large with nonzero fractional
> parts.
It's not completely obvious. I thought of it soon but wondered if I
thought of all the cases. Steve's remark about Infs and NaNs points
to possible problems:
> > > > #include <math.h>
> > > >
> > > > float roundf(float x) {
> > > > float t;
> > > > if (x >= 0.0) {
Suppose x is a NaN. Then it will compare strangely with everything and
we won't get here.
> > > > t = ceilf(x);
> > > > if ((t - x) > 0.5) t -= 1.0;
> > > > return t;
> > > > } else {
We get here for NaNs.
> > > > t = ceilf(-x);
And we really shouldn't do arithmetic on the NaNs. ceilf() should return
its arg for a NaN, but it's not clear what happens for -x. Well, I checked
what happens starting with the QNaN x= 0.0 / 0.0. Almost everything is
broken:
- gcc miscomputes 0.0 / 0.0 at compile time. It gives a positive NaN, but
the npx would give a negative NaN ("Real Indefinite" = the same one except
with the opposite sign).
- gcc invalidly optimizes -x by ORing in the sign bit (even without -O).
It should do the same as the npx (which is to not change anything).
- printf() doesn't print the sign bit for NaNs.
- C99 requires printf() to misspell "NaN" as "nan"
> > > > if ((t + x) > 0.5) t -= 1.0;
> > > > return -t;
Now if x is NaN, t is NaN and we do lots of arithmetic on it. I think
gcc's invalild optimizations cancel out so NaNs will be returned unchanged.
This is not clear :-).
> > > > }
> > > > }
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).
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.
Bruce
help
Want to link to this message? Use this
URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20031201182219.O4431>