#371 Fix round(2^53)

closed-fixed
None
2008-06-03
2008-05-12
No

C99 round() failed with an argument of 2^53, because
adding 0.5 to that number caused overflow. With this
patch, it returns its argument if the argument is
already integral--as the C99 reference implementation
(F.9.6.6 paragraph 2) does. Both roundf() and roundl()
are patched similarly.

Speed is almost unaffected: a few percent slower for
non-integral arguments, and a couple percent faster
for integral arguments. The F.9.6.6 implementation
would run at about one-third this speed because it
saves and restores the whole FP environment.

With this patch, I get the same results as with this
NetBSD code:

http://archive.netbsd.se/?ml=freebsd-ports-bugs&a=2005-04&m=870458

(search for "round(double x)"), but that code runs
about ten percent slower; and a minimal change to
the 2004 mingwex sources seemed preferable anyway.

Discussion

  • Anonymous - 2008-05-12
     
    Attachments
  • Keith Marshall

    Keith Marshall - 2008-05-14

    Alternative reference implementation.

     
    Attachments
  • Keith Marshall

    Keith Marshall - 2008-05-14

    Logged In: YES
    user_id=823908
    Originator: NO

    Danny has proposed alternative implementations, which appear similar to the NetBSD example cited. To me, these seem unnecessarily verbose. I've attached an alternative reference implementation, for `round()' and `llround()' only, together with a test case which seems to yield correct results in each of the corner cases reported in these two related tickets.
    File Added: alt-round.tar.gz

     
  • Anonymous - 2008-05-16

    Logged In: YES
    user_id=231935
    Originator: YES

    I have no objection to anything in your 2008-05-14 12:59 posting.

    It's sound engineering practice to factor out a core "helper" routine.
    The one you've written is clear, fast, and passes all the tests I've
    been developing. You seem to have the matter well in hand, so perhaps
    it's best for me to focus on testing. I've long thought that we should
    have unit tests for libmingwex, and the round() family is a good place
    to start. I'll have little time before September to propose a general
    framework worthy of consideration, but I'll make time to test anything
    you write for round() and its kin.

     
  • Danny Smith

    Danny Smith - 2008-05-16

    Logged In: YES
    user_id=11494
    Originator: NO

    If you look at the implementation of trunc , you will see that it just uses a temporary change in FPU control word to get FE_TOWARDZERO. So rather than fiddle with trunc and FP addition, why not just do the obvious for round, eg:

    double
    round (double _x){
    double res;
    unsigned short saved_cw;
    unsigned short tmp_cw;
    __asm__ ("fnstcw %0;" : "=m" (saved_cw)); /* save FPU control word */
    tmp_cw = (saved_cw & ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO))
    | FE_TONEAREST;
    __asm__ ("fldcw %0;" : : "m" (tmp_cw));
    __asm__ ("frndint;" : "=t" (res) : "0" (_x)); /* round to nearest */
    __asm__ ("fldcw %0;" : : "m" (saved_cw) ); /* restore saved control word */
    return res;
    }

    long
    lround (double _x)
    {
    double res = round (_x);

    if (!isfinite (res)
    || res > (double) LONG_MAX
    || res < (double) LONG_MIN)
    {
    errno = ERANGE;
    /* Undefined behaviour, so we could return anything. */
    /* return res > 0.0 ? LONG_MAX : LONG_MIN; */
    }
    return (long) res;
    }

    long long
    llround (double _x)
    {
    double res = round (_x);

    if (!isfinite (res)
    || res > (double) LONG_LONG_MAX
    || res < (double) LONG_LONG_MIN)
    {
    errno = ERANGE;
    /* Undefined behaviour, so we could return anything. */
    /* return res > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
    }
    return (long long) res;
    }

    This should be more efficient.

    With that, and Keith test case I get:

    Input value: 4.99999999999999940000e-001
    round(x): 0.00000000000000000000e+000
    llround(x): 0

    Input value: -4.99999999999999940000e-001
    round(x): 0.00000000000000000000e+000
    llround(x): 0

    Input value: 5.00000000000000000000e-001
    round(x): 0.00000000000000000000e+000
    llround(x): 0

    Input value: -5.00000000000000000000e-001
    round(x): 0.00000000000000000000e+000
    llround(x): 0

    Input value: 5.00000000000000110000e-001
    round(x): 1.00000000000000000000e+000
    llround(x): 1

    Input value: -5.00000000000000110000e-001
    round(x): -1.00000000000000000000e+000
    llround(x): -1

    Input value: 9.00719925474099100000e+015
    round(x): 9.00719925474099100000e+015
    llround(x): 9007199254740991

    Input value: 4.50359962737049700000e+015
    round(x): 4.50359962737049700000e+015
    llround(x): 4503599627370497

    Input value: -4.50359962737049700000e+015
    round(x): -4.50359962737049700000e+015
    llround(x): -4503599627370497

    NB: round (5.00000000000000000000e-001) == 0.0, using the round-to-even rule for 0.5

    Danny

     
  • Anonymous - 2008-05-16

    Logged In: YES
    user_id=231935
    Originator: YES

    Cloning trunc(), but setting FE_TONEAREST, is probably the fastest
    and most reliable way to round to nearest, because all the work is
    done by FRNDINT. But FE_TONEAREST resolves exact ties in favor of
    the even integer, whereas C99's round() picks the integer farther
    away from zero.

    You could get C99 round() semantics by adding this at the bottom:

    if(0.0 < _x && 0.5 == _x - res) return res + 1.0;
    else if(_x < 0.0 && -0.5 == _x - res) return res - 1.0;
    else return res;

    That seems to be only slightly faster than the trunc() approach.
    Instead of adjusting half the real numbers, it adjusts only a
    much smaller number of points; but it has to test every number,
    and that overhead is about the same as with trunc().

     
  • Danny Smith

    Danny Smith - 2008-05-22

    Logged In: YES
    user_id=11494
    Originator: NO

    Thanks Greg, I had forgotten about C99 vs i386 FE_TONEAREST difference.
    Keith's approach in alt-round.tar.gz looks good then. Is someone willing extend that to float and long double versions?

    Danny

     
  • Keith Marshall

    Keith Marshall - 2008-05-22

    Logged In: YES
    user_id=823908
    Originator: NO

    Danny wrote:
    > Keith's approach in alt-round.tar.gz looks good then. Is someone
    > willing extend that to float and long double versions?

    Thanks Danny. Yes, I'll follow up with that. I already have a generic implementation working for all nine of the affected functions; (it uses your FE_TONEAREST technique, with a modification of Greg's proposed mid-way adjustment). I'll tidy it up, WRT comments, ChangeLog, etc., and post it here for review this weekend; it reproduces the correct results from my preceding test case, (with one caveat; I don't seem to be able to printf a long double value on Win2K -- neither the `l' nor the `L' modifiers documented at http://msdn.microsoft.com/de-de/library/tcxf1dw6.aspx, for use with the `f' format spec, (but conspicuously not apparently with `e' or `g'), seem to work as documented, for me).

     
  • Anonymous - 2008-05-22

    Logged In: YES
    user_id=231935
    Originator: YES

    > I don't seem to be able to printf a long double value on Win2K

    IIRC, msvcrt's printf() assumes 'long double' is a synonym for
    'double'; but snprintf() should work, because we've implemented
    it via __mingw_snprintf().

     
  • Keith Marshall

    Keith Marshall - 2008-05-22

    Logged In: YES
    user_id=823908
    Originator: NO

    Greg wrote:
    > IIRC, msvcrt's printf() assumes 'long double' is a synonym
    > for 'double';

    Yuk! While in mingw32, the two are distinct...

    > but snprintf() should work, because we've implemented it
    > via __mingw_snprintf().

    All the more reason, IMO, why we should extend the __mingw_snprintf() formatting engine for use in replacements for the entire gamut of printf functions, (as I've suggested in https://sourceforge.net/support/tracker.php?aid=1961893\); native MSVCRT printf doesn't appear to adequately support mingw32 data types.

    But, that's a subject for another patch...

     
  • Keith Marshall

    Keith Marshall - 2008-05-25

    Alternative patch file

     
    Attachments
  • Keith Marshall

    Keith Marshall - 2008-05-25

    Logged In: YES
    user_id=823908
    Originator: NO

    I've attached an updated patch, against current CVS, which addresses the rounding issue, as previously discussed, for all of {,l,ll}round{,f,l}() functions. This reproduces the correct results from my prior test case, for each of the {,l,ll}round() functions, and also appears to DTRT for each the other six; (testing on GNU/Linux, rounding is correct for _all_ cases, in the vicinity of 0.5; however, I have not (yet) adjusted the extremes to allow for differences in precision of the input data types).

    File Added: alt-round.patch.gz

     
  • Chris Sutcliffe

    Chris Sutcliffe - 2008-05-25

    Logged In: YES
    user_id=570619
    Originator: NO

    Keith, should I hold off applying the patch and wait for the patch that addesses the adjustment of extremes, or is this good to go?

     
  • Keith Marshall

    Keith Marshall - 2008-05-25

    Logged In: YES
    user_id=823908
    Originator: NO

    Chris,

    I was, and still am, willing to commit this myself, but I wanted to give Danny, Greg, FX and others an opportunity to comment first.

    I'm hoping no further patch will be required. I haven't tested with any other input values than x = +/-0.5 +/- nextafter(0.5, 0.0), (for which all cases yield correct results), Greg's 2^53 value, and the two corner case values for which FX reported a problem. For {,l,ll}roundf() the overflow case equivalent to 2^53 would be 2^24, (I believe), and I think it should be 2^64 for {,l,ll}roundl(), (although I haven't found any definitive confirmation that 80-bit long doubles do use a 64-bit mantissa, but http://tinyurl.com/3jn5hl hints that it may be so).

    For the two corner case values suggested by FX, I don't know his basis for those particular choices, so it isn't clear to me how to choose equivalents for the float and long double cases.

     
  • Anonymous - 2008-05-25

    Logged In: YES
    user_id=231935
    Originator: YES

    > I wanted to give
    > Danny, Greg, FX and others an opportunity to comment first.

    I'll try to make time for testing soon, but it'll probably take until
    Tuesday before I can do that. Monday's a national holiday here, and
    right now I'm wrestling with a production problem in my day job.

    > For {,l,ll}roundf() the overflow
    > case equivalent to 2^53 would be 2^24, (I believe), and I think it
    > should be 2^64 for {,l,ll}roundl()

    I think that's all correct:

    $echo '#include <float.h>' | cpp -dM - |grep MANT |grep '[1-9]'

    #define __FLT_MANT_DIG__ 24
    #define __LDBL_MANT_DIG__ 64
    #define __DBL_MANT_DIG__ 53

    > corner case values suggested by FX, I don't know his basis for
    > those particular choices

    I guess
    4.99999999999999944488848768742172978818416595459e-1L
    is the next-lower neighbor of one-half; I think he got the other one
    4.503599627370497e+15L
    from this message:
    http://gcc.gnu.org/ml/gcc-patches/2006-10/msg00917.html
    which gives no exact rationale for it, but it appears to be 1+2^52.

     
  • Anonymous - 2008-06-03

    Logged In: YES
    user_id=231935
    Originator: YES

    I haven't been able to find anything wrong with the proposed
    implementation, and it fixes all reported problems I'm aware
    of, so I recommend adopting it. Here:

    http://cvs.savannah.gnu.org/viewvc/lmi/lmi/round_test.cpp?view=log

    is my unit test suite. It's long and rambling, and tightly
    coupled to other unit tests I use in my day job, but I think
    it has pretty good coverage for floating-point rounding.
    Twelve of its 1116 tests fail with the old implementation;
    all succeed with the new.

     
  • Keith Marshall

    Keith Marshall - 2008-06-03

    Logged In: YES
    user_id=823908
    Originator: NO

    Thanks, Greg. I've committed it.

     
  • Keith Marshall

    Keith Marshall - 2008-06-03
    • milestone: --> IINR_-_Include_In_Next_Release
    • assigned_to: nobody --> keithmarshall
    • status: open --> closed-fixed
     

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks