[cairo] Banker's rounding problems and Nautilus

[Daniel Amelang]

On 12/1/06, Bill Spitzak <spitzak at d2.com> wrote:
> The following code appears to do floor(x+.5) quickly. However two
> additions are necessary, as it is impossible to get the extra correction
> of .25 and the large factor into the same constant becuase it requires
> higher precision than the double can hold. Somebody needs to test if
> this is still faster, and also that the optimizer does not merge the two
> constants. Some technique to produce an 80-bit constant would also work.
>
> #include <stdio.h>
> #include <stdlib.h>
> #include <math.h>
>
> #if __BIG_ENDIAN__
> #define iman_ 1
> #else
> #define iman_ 0
> #endif
>
> typedef union {
>    double v;
>    long word[2];
>    long long big;
> } Double_Long;
>
> // Does floor(x) but quickly. Works for -1073741824 to 1073748823.9977785
> inline long fast_floor(double val) {
>    Double_Long v; v.v = (val-.25) + (68719476736.0*32768.0*1.5);
>    return v.word[iman_] >> 1;
> }
>
> // Does floor(x+.5) but quickly. Works for -1073741824.5 to
> 1073741823.49987785
> inline long fast_rfloor(double val) {
>    Double_Long v; v.v = (val+.25) + (68719476736.0*32768.0*1.5);
>    return v.word[iman_] >> 1;
> }
>
> int main(int argc, char** argv) {
>    int i, j;
>    double f;
>    if (argc < 2) {
>      fprintf(stderr, "Return fast_rfloor(x) for each argument\n");
>      return 1;
>    }
>    for (i = 1; i < argc; i++) {
>      f = strtod(argv[i],0);
>      j = fast_rfloor(f);
>      printf("fast_rfloor(%g) = %d\n", f,j);
>    }
>    return 0;
> }

I ran some tests on this approach, and found that it actually doesn't
match the floor(+.5) behavior very closely. For example, in my random
hammering test, I got several thousand like this:

Rounding 0.499882361618747161902120978993:
Floor+0.5 says: 0
Bill says: 1
Rounding 316076.499983850168064236640930175781:
Floor+0.5 says: 316076
Bill says: 316077
Rounding 693704.499969894299283623695373535156:
Floor+0.5 says: 693704
Bill says: 693705
Rounding 1046961.499895715387538075447082519531:
Floor+0.5 says: 1046961
Bill says: 1046962

After messing with it for quite some time, I think I'm ready to say
that you can't turn banker's rounding into arithmetic rounding by just
add a number and shifting. BUT, I discovered that you can get pretty
close by forcing the rounding to be performed as far down away from
the decimal point as possible. So, I think Bill and Behdad were on the
right track with the 31.1 fixed point number approach.

Here's one that converts to a 33.19 fixed-point number (using the
entire mantissa). My random hammer hasn't yet found a case where it
doesn't perform away-from-zero rounding properly, and it's pretty
fast. Tomorrow, I'll clean it up, turn it into a full-fledged patch
with documentation and perf profiles. In the meantime, here's a sneak
peek:

#define CAIRO_MAGIC_NUMBER_FIXED_33_19 (12884901888.0)
int
_cairo_lround (double d)
{
    union {
        double d;
        uint32_t ui32[2];
    } u;

    u.d = (d + 0.499999) + CAIRO_MAGIC_NUMBER_FIXED_33_19;
#ifdef FLOAT_WORDS_BIGENDIAN
    return (u.ui32[0] << 13) | (u.ui32[1] >> 19);
#else
    return (u.ui32[1] << 13) | (u.ui32[0] >> 19);
#endif
}

Dan