Re: Re: questionnable nint() behavior

[Steve Kargl <sgk at troutmask dot apl dot washington dot edu>]

On Tue, Mar 21, 2006 at 12:52:47AM +0000, Timothy C Prince wrote:
> 
> On Mon, Mar 20, 2006 at 10:24:12PM +0100, Dominique Dhumieres wrote:
> > I have found in my file the following program
> > 
> >       program main
> >       real x, y
> >       x = 8388609.0
> >       y = 0.4999999701976776123046875
> >       print '(A6,F9.1,A3,I8)', 'nint (', x, ') =', nint (x)
> >       print '(A6,F10.8,A3,I2,A18,L1)', 'nint (', y, ') =', nint (y),    &
> >      &      ', where y < 0.5 = ', y < 0.5
> >       end
> > 
> > xlf, ifc, and g95 return
> > 
> > nint (8388609.0) = 8388609
> > nint (0.49999997) = 0, where y < 0.5 = T
> > 
> > IMHO the correct result, g77, pgf, and gfc return
> > 
> > nint (8388609.0) = 8388610
> > nint (0.49999997) = 1, where y < 0.5 = T
> > 
> > as if nint(a) was implemented as int(a+0.5) without
> > taking care of rounding mode (rounding to nearest
> > rounds to even on tie). I don't remember where I borrowed
> > the code, but I thought the problem was fixed in gfortran.
> >
> 
> In trans-intrinsics.c, you'll find
> 
> /* This is needed because the gcc backend only implements FIX_TRUNC_EXPR
>    NINT(x) = INT(x + ((x > 0) ? 0.5 : -0.5)).  */
> 
> static tree
> build_round_expr (stmtblock_t * pblock, tree arg, tree type)
> 
> A much better algorithm may be found in FreeBSD's round[fl].c. :-)
> 
> ______________________________________
> We could get away with this, by promoting float to double, and double
> to long double, as happens already in gfortran i386 for some -march=i
> values. g77 (libf2c) had a more careful and slow version, using floor().
> Tim Prince

FreeBSD's roundf:

        if (!isfinite(x))
                return (x);

        if (x >= 0.0) {
                t = floorf(x);
                if (t - x <= -0.5)
                        t += 1.0;
                return (t);
        } else {
                t = floorf(-x);
                if (t + x <= -0.5)
                        t += 1.0;
                return (-t);
        }

I'll being on an airplane for 10 hours this week.  Maybe I'll look at this.

-- 
Steve