Harder than it looks: rounding a float to the nearest integer

Floats are outrageously complicated. Tread carefully! If you ever deal with the nitty gritty of floating point numbers you should read every post on the topic on Bruce Dawson's blog. It's up to 15 posts long.

TLDR - if you think you know all there is to know about floats then you're wrong. Read on!

http://randomascii.wordpress.com/2013/02/07/float-precision-... (has links to all posts in series)

I wrote a very similar article about how testing for integer overflow is also harder than it looks: http://blog.reverberate.org/2012/12/testing-for-integer-over...

Wouldn't this be easier if both numbers were just represented bitwise, and the float were converted to integer that way? Seems more straightforward to me.

Floating point representation works by having "exponent" bits, and "mantissa" or fraction bits. Also a sign bit. For explanation's sake let's say we have a single precision 32-bit FPU. IEEE 754 single precision is 1 sign bit, 8 exponent bits, and 23 mantissa bits. The exponent bits are subtracted by 127 to yield the exponent.

So, the basic formula is:

  (-1)*sign * 2^exponent-127 * 1.mantissa

So right away, if exponent < 0x7F, round to 0.

If exponent > 9E (for 32-bit int), round to infinity.

Multiply

  1.mantissa * 2^exponent-127

This will be a power of 2 so there's probably a fast bit-shifting method for this. I'm just too lazy to look up how this should be done.

If exponent > 0x96, the exponent's too big to represent fractions of 1, so you're done.

Otherwise, you're within 1 number. If the remaining decimal < 0.5, truncate. Otherwise, add 1 and truncate.

The title piqued my interest in how to implement rounding of floats, on the bit-fiddling level - certainly not an easy task! Disappointingly though, the article starts out by using type casting. If it's a solution on that level we're after, why not just use roundf()?

Doesn't the modulus function avoid most of the issues, here?

    float RoundToNearestWholeNumber(float f) {
        if (double.IsNanOrInfinity(f)) return f;
        if (f < 0) return -RoundToNearestWholeNumber(-f);
        
        // make sure we're working with a consistent value
        float v = DiscardExtraPrecision(f);
        // the fractional part is guaranteed to be representable
        float r = v % 1;
        // largest whole value up to v (the floor)
        float w = v - r;
        // switch to the ceiling if in the upper half
        if (r > 0.5) return w + 1;
        return w;
    }

Feel free to disillusion me.

Hi,

I don't understand the problem with a) rounding floats to integer and I don't understand b) why you would ever have to overflow an integers.

a) Wouldn't it solve the problem to write an algorithm that can use mathematical axioms like +,-,/,* virtually using strings? I mean you want it to calculate 99999999999999999999999999999999999999999999999999999 + 9999999999999999999999999999, why not teach it how we do this instead of having to use bit-hacking tricks? Maybe my example isn't good, but I mean writing an algorithm that reads 999 for example and 1000 knows that it now has to put a 1 infront of 999. You see, that's what I mean by manually.

b) And can you please explain "why" it isn't enough to use ie. a linked list of 8bit virtual address spaces and allocate using virtual address spaces everything that fits into RAM and everything that goes beyond the addressable space onto swap?

Why would you have to worry about flipping bits in an overflowed integer/float, when you could use that mechanism. Don't you just recursively run through the virtual address map until you're able modify the allocated data in that memory segments again? I don't see a problem with this implementation except for performance maybe.

When the algirhtms knows about the 32bit limit for example (constant), you can tell it to not allocate real memory with that size, but allocate virtual memory with that size dependant on the input data size. Yes I know that the CPU has built-in mechanism for the mathematical axioms, but doing it manually would solve the problem, from my limited perspective.

Sorry if I'm blind for seeing the problem, I'm not ignorant, just really curious in how this can be solved once and forever.

Just like there is a relatively clear separation between electrical engineering and software engineering, I'm starting to think that we need separate words for this kind of stuff and for hacking up a CRUD-some-tables MVP using Ruby on Rails. They're nearly entirely unrelated skills.

I, for one, am very happy that I do not need to think about this stuff. I think that rounding a float to an integer is super easy:

    5.3.round

Now that we know floats have coarse granularity on large scale, why not reduce the scale first?

float fr = (int)f; if(f-fr>=0.5f)fr+=1.0f; return fr;

It works for all the examples in the article. I'm not saying that I've created a good solution. I'm saying that the author is trying to create a problem out of nowhere.

Thought it would be an article to teach people how to manually do float rounding according to a specific fp encoding scheme... :-(

I don't think that better function that works actually works flawlessly; when passed a NaN or -Infinity, it triggers undefined behavior.

Still, this is a nice article.

There are lots of cute tricks for rounding to integer and to integral values. A few of my favorites:

1. Round to an integral value according to the current rounding mode, allowing inexact to be set (rintf( ) in C, pretty close to what Pascal is going for in his example):

    // All suitably large floats are integral already.
    if (fabsf(x) >= 0x1.0p23f) return x;
    // Otherwise we can force rounding to occur in the
    // units position by simply adding and subtracting
    // 2**23 with the correct sign.
    const float forceRound = copysignf(0x1.0p23f, x);
    return x + forceRound - forceRound;

If you want the result as an integer and it’s in the [-2^23, 2^23] range, you can just rip out bits 0:22 of the floating-point representation of x + forceRound via a mask, which is really simple (faster than hardware convert-to-integer instructions on some architectures!).

If the calls to fabsf and copysignf look like hazards to you, keep in mind that good compilers lower them to one or two bitwise operations on most architectures (though I tend to write these routines in assembly to avoid any risk of missed optimization).

2. Round to nearest, ties away from zero (roundf( ) in C, which can be annoying because it doesn’t correspond to a hardware rounding mode on most architectures):

    // All suitably large floats are integral already.
    if (fabsf(x) >= 0x1.0p23f) return x;
    // Now that we know (int32_t)x cannot raise invalid,
    // we can use that to get truncf(x).
    const float truncated = (int32_t)x;
    // Compare the residual to 0.5 to determine if the
    // result needs to be adjusted away from zero.
    if (fabsf(x - truncated) >= 0.5f)
        return truncated + copysignf(1.0f, x);
    return truncated;

This is also a pretty close analogue to what Pascal was considering. I like my implementation more because it doesn’t spuriously raise the inexact flag (whereas Pascal’s implementation does). This still isn’t a complete implementation of roundf( ), as it gets the sign of zero wrong when the result should be negative zero. Left as an exercise to the reader.

3. Since people asked about doing the rounding using integer operations elsewhere in the discussion: assume the input is known to be in the interesting range already, and that we have some utility functions for converting between a floating-point number and its encoding (which is a NOP architecturally, but requires a union or memcpy or other type-punning in C):

    const uint32_t xRep = toRep(x);
    const uint32_t xAbs = xRep & 0x7fffffff;
    const uint32_t xSgn = xRep & 0x80000000;

    // If |x| > 2**23, or x is Inf or NaN the result is x.
    if (xAbs > 0x4b000000) return x;

    // If |x| < 1.0f, the result is +/-0 if |x| < 0.5f and
    // +/-1 otherwise.
    if (xAbs < 0x3f800000) {
        if (xAbs < 0x3f000000) return fromRep(xSgn);
        return fromRep(xSgn | 0x3f800000); 
    }

    // Otherwise, use the exponent of x to identify the
    // units bit of x.
    const int exponent = xAbs - 0x3f800000 >> 23;
    const uint32_t unitsBit = 0x00800000 >> exponent;
    // "Add 0.5 to x” by adding half the units bit to the
    // encoding.  If this doesn’t change the exponent it
    // really is equivalent to adding 0.5.  If it carries
    // into the exponent, then it’s not actually the same,
    // but it still does what we need it to do.  After
    // “adding 0.5”, we simply mask off the fractional
    // bits to get the result.
    return fromRep(xAbs + (unitsBit >> 1) & -unitsBit);

This doesn’t set inexact, obviously, which makes it unsuitable for some uses (but actually nearbyint( ) is required not to set inexact, so this is pretty useful there). It’s actually even less code than it looks like; C tends to make this sort of thing relatively verbose. On ARM a complete implementation is about 12 instructions, for example.

N.B. although I write this sort of stuff professionally, I’m just sketching code here without even checking if it compiles. Test before use.

Aside from what other commenters have said, this article seems to totally ignore negative numbers.