Printf Rounding Behavior for Doubles

printf rounding behavior for doubles

It's "round half to even" or "Banker's rounding". The last digit of the rounded representation is chosen to be even if the number is exactly half way between the two.

http://linuxgazette.net/144/misc/lg/a_question_of_rounding_in_issue_143.html:

"For the GNU C library, the rounding rule used by printf() is "bankers rounding" or "round to even". This is more correct than some other C libraries, as the C99 specification says that conversion to decimal should use the currently selected IEEE rounding mode (default bankers rounding)."

Trying to recreate printf's behaviour with doubles and given precisions (rounding) and have a question about handling big numbers

You can't round to a particular decimal precision with binary floating point arithmetic. It's not just possible. At small magnitudes, the errors are small enough that you can still get the right answer, but in general it doesn't work.

The only way to round a floating point number as decimal is to do all the arithmetic in decimal. Basically you start with the mantissa, converting it to decimal like an integer, then scale it by powers of 2 (the exponent) using decimal arithmetic. The amount of (decimal) precision you need to keep at each step is roughly (just a bit over) the final decimal precision you want. If you want an exact result, though, it's on the order of the base-2 exponent range (i.e. very large).

Typically rather than using base 10, implementations will use a base that's some large power of 10, since it's equivalent to work with but much faster. 1000000000 is a nice base because it fits in 32 bits and lets you treat your decimal representation as an array of 32-bit ints (comparable to how BCD lets you treat decimal representations as arrays of 4-bit nibbles).

My implementation in musl is dense but demonstrates this approach near-optimally and may be informative.

Rounding doubles - .5 - sprintf

It seems you have to use math round function for correct rounding.

printf("%.2f %.2f\n", 5.555, round(5.555 * 100.)/100.);

This gives the following output on my machine:

5.55 5.56

C printf float rounding

Ok, so I guess my previous algorithm didn't fully implement rounding so let's take a look how the result could be rounded. I'll use your example number as my example number too. First, we find that 0.000011/(2^(-17)) = 0.000011*(2^17) = 1.441792 so the power is -17. Then, we output "1.", subtract 1 from 1.441792 and multiply it by 16, giving 7.068672. We output 7, subtract 7 from it and multiply by 16, giving 1.09875199999999. We output 1, subtract 1 from it and multiply by 16, giving 1.58003199999985. We output 1, subtract 1 from it and multiply by 16, giving 9.28051199999754. Then output 9, subtract 9, multiply by 16, the result is 4.48819199996069. We output 4, subtract 4, multiply by 16, the result is 7.81107199937105.

Now, we're going to output the last character. Now we do the magic. Because 7.81107199937105 is closer to 8 than to 7, we output "8". This magic is done only for the last character. For the non-last characters, the integral part is always used and the fractional part is not used at all in determining which character to output. Then, after this, we output "p-17" because the power was -17.

Note that the usual rounding rules say that 7.5 which is equally close to both 7 and 8 is rounded to 8, not to 7 and 6.5 would be rounded to 7, not to 6. However, if you want to implement the round half to even and you encounter e.g. 6.5 then it is rounded down to 6 because 6 is even and 7 is not. I'm not sure if what you found about the banker round applies also to %a, the only thing you can do is to test implementing various rounding algorithms and see which gives the same results as the real %a of printf. Shouldn't be that hard because the different rounding algorithms just differ on how half is handled. The rest is rounded to the closest number.

By the way, I was wrong in your previous question (How %a conversion work in printf statement?) in saying that for 3.2 which has the non-rounded representation 1.999999....p+1 and the rounded representation 1.99999ap+1 the last "a" would occur due to limited floating point precision. Of course it occurs due to rounding and not due to limited floating point precision, as you probably have realized by now.

C++ printf Rounding?

It's a little tricky with the field width

#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include <sstream>
#include <limits>

#define INV_SCALE 100000000

static const int      WIDTH   = std::ceil(
                                    std::log10(std::numeric_limits<uint64_t>::max())
                                ) + 1 /* for the decimal dot */;
static const uint64_t INPUT   = 1033468;
static const double   DIVISOR = double(INV_SCALE);
static const int      PREC    = std::ceil(std::log10(DIVISOR));

static const double   DAVIDS_SAMPLE = 1000000.000033;

namespace {
std::string to_string(double d, int prec) {
    std::stringstream s;
    s << std::fixed
      << std::setw(WIDTH)
      << std::setprecision(prec) 
      << d;
    // find where the width padding ends    
    auto start = s.str().find_first_not_of(" ");
    // and trim it left on return
    return start != std::string::npos ? 
                    &(s.str().c_str()[start]) : "" ;
}
}

int main() {
    for (auto& s : 
            {to_string(INPUT/DIVISOR, PREC), to_string(DAVIDS_SAMPLE, 6)} 
        ) std::cout << s << std::endl;

    return /*EXIT_SUCCESS*/ 0;
}

output:

0.01033468
1000000.000033

---

Related Topics