Is hardcoding least significant byte of a double a good rounding strategy?

Question

I have a function doing some mathematical computation and returning a double. It ends up with different results under Windows and Android due to std::exp implementation beging different (Why do I get platform-specific result for std::exp?). The e-17 rounding difference gets propagated and in the end it's not just a rounding difference that I get (results can change 2.36 to 2.47 in the end). As I compare the result to some expected values, I want this function to return the same result on all platform.

So I need to round my result. The simpliest solution to do this is apparently (as far as I could find on the web) to do std::ceil(d*std::pow<double>(10,precision))/std::pow<double>(10,precision). However, I feel like this could still end up with different results depending on the platform (and moreover, it's hard to decide what precision should be).

I was wondering if hard-coding the least significant byte of the double could be a good rounding strategy.

This quick test seems to show that "yes":

#include <iostream>
#include <iomanip>

double roundByCast( double d )
{
    double rounded = d;
    unsigned char* temp = (unsigned char*) &rounded;
    // changing least significant byte to be always the same
    temp[0] = 128;
    return rounded;
}

void showRoundInfo( double d, double rounded )
{
    double diff = std::abs(d-rounded);
    std::cout << "cast: " << d << " rounded to " << rounded << " (diff=" << diff << ")" << std::endl;
}

void roundIt( double d )
{
    showRoundInfo( d, roundByCast(d) );
}

int main( int argc, char* argv[] )
{
    roundIt( 7.87234042553191493141184764681 );
    roundIt( 0.000000000000000000000184764681 );
    roundIt( 78723404.2553191493141184764681 );
}

This outputs:

cast: 7.87234 rounded to 7.87234 (diff=2.66454e-14)
cast: 1.84765e-22 rounded to 1.84765e-22 (diff=9.87415e-37)
cast: 7.87234e+07 rounded to 7.87234e+07 (diff=4.47035e-07)

My question is:

• Is unsigned char* temp = (unsigned char*) &rounded safe or is there an undefined behaviour here, and why?
• If there is no UB (or if there is a better way to do this without UB), is such a round function safe and accurate for all input?

Note: I know floating point numbers are inaccurate. Please don't mark as duplicate of Is floating point math broken? or Why Are Floating Point Numbers Inaccurate?. I understand why results are different, I'm just looking for a way to make them be identical on all targetted platforms.

---

Edit, I may reformulate my question as people are asking why I have different values and why I want them to be the same.

Let's say you get a double from a computation that could end up with a different value due to platform specific implementations (like std::exp). If you want to fix those different double to end up having the exact same memory representation (1) on all platforms, and you want to loose the fewest precision as possible, then, is fixing the least significant byte a good approach? (because I feel that rounding to an arbitrary given precision is likely to loose more information than this trick).

(1) By "same representation", I mean that if you transform it to a std::bitset, you want to see the same bits sequence for all platform.

Answer 1

Is unsigned char* temp = (unsigned char*) &rounded safe or is there an undefined behaviour here, and why?

It is well defined, as aliasing through unsigned char is allowed.

is such a round function safe and accurate for all input?

No. You cannot perfectly fix this problem with truncating/rounding. Consider, that one implementation gives 0x.....0ff, and the other 0x.....100. Setting the lsb to 0x00 will make the original 1 ulp difference to 256 ulps.

No rounding algorithm can fix this.

You have two options:

• don't use floating point, use some other way (for example, fixed point)
• embed a floating point library into your application, which only uses basic floating point arithmetic (+, -, *, /, sqrt), and don't use -ffast-math, or any equivalent option. This way, if you're on a IEEE-754 compatible platform, floating point results should be the same, as IEEE-754 mandates that basic operations should be calculated "perfectly". It means as if the operation calculated at infinite precision, and then rounded to the resulting representation.

Btw, if an input 1e-17 difference means a huge output difference, then your problem/algorithm is ill-conditioned, which generally should be avoided, as it usually doesn't give you meaningful results.

Answer 2

No, rounding is not a strategy for removing small errors, or guaranteeing agreement with calculations performed with errors.

For any slicing of the number line into ranges, you will successfully eliminate most slight deviations (by placing them in the same bucket and clamping to the same value), but you greatly increase the deviation if your original pair of values straddle a boundary.

In your particular case of hardcoding the least significant byte, the very near values

0x1.mmmmmmm100

and

0x1.mmmmmmm0ff

have a deviation of only one ULP... but after your rounding, they differ by 256 ULP. Oops!

Answer 3

Your function won't work because of aliasing.

double roundByCast( double d )
{
    double rounded = d;
    unsigned char* temp = (unsigned char*) &rounded;
    // changing least significant byte to be always the same
    temp[0] = 128;
    return rounded;
}

Casting to unsigned char* for temp is allowed, because char* casts are the exception to the aliasing rules. That's necessary for functions like read, write, memcpy, etc, so that they can copy values to and from byte representations.

However, you aren't allowed to write to temp[0] and then assume that rounded changed. You must create a new double variable (on the stack is fine) and memcpy temp back to it.

Answer 4

What you are doing is totally, totally misguided.

Your problem is not that you are getting different results (2.36 vs. 2.47). Your problem is that at least one of these results, and likely both, have massive errors. Your Windows and Android results are not just different, they are WRONG. (At least one of them, and you have no idea which one).

Find out why you get these massive errors and change your algorithms to not increase tiny rounding errors massively. Or you have a problem that is inherently chaotic, in which case the difference between results is actually very useful information.

What you are trying just makes the rounding errors 256 times bigger, and if two different results end in ....1ff and ....200 hexadecimal, then you change these to ....180 and ....280, so even the difference between slightly different numbers can grow by a factor 256.

And on a bigendian machine your code will just go kaboom!!!