0

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.

jpo38
  • 20,821
  • 10
  • 70
  • 151
  • 4
    If you need to make floating point identical, then you are doing something wrong. You need to switch to either fixed-point or rational arithmetic. – stark Jan 18 '19 at 15:00
  • Have you considered doing interval mathematics, and making the test "is the actual result within the interval?" That is both mathematically and logically rigorous. – Yakk - Adam Nevraumont Jan 18 '19 at 15:01
  • 1
    Won't the standard rounding functions work for you? https://en.cppreference.com/w/cpp/numeric/math/round – Jesper Juhl Jan 18 '19 at 15:01
  • @stark: I'm using 3rd party libraries (like nlopt) that uses doubles. – jpo38 Jan 18 '19 at 15:05
  • @Yakk-AdamNevraumont: I'm actually comparing the output to a binary file with expected outputs. I need a 1-1 matching on bytes representing the double result... – jpo38 Jan 18 '19 at 15:06
  • @JesperJuhl: `std::round` rounds to the nearest integer, but I have significant decimals I must preserve. – jpo38 Jan 18 '19 at 15:07
  • round(a*100)/100 – stark Jan 18 '19 at 15:10
  • @stark: OK this works better than by `ceil`+`pow` code. But again, how to decide what precision I should use (`100`, `1000`, `1000000`) and would that get more chance to end up with the same number on all platforms? – jpo38 Jan 18 '19 at 15:12
  • It was response to your "rounds to the nearest integer" just to demonstrate that it is trivial to do any precision required. And no, there is no guarantee of same in my suggestion or yours. – stark Jan 18 '19 at 15:20
  • Question is why do you need such rounding? This looks like classic [XY problem](http://xyproblem.info/), where Y here question if your strange rounding is correct, where X is some unknown function which uses `exp` and rounding issue which is not explained. Please describe details of your function and explain in what is the manifestation of precision problem. There are some ready solutions for that for sure. – Marek R Jan 18 '19 at 15:20
  • 1
    @MarekR: This shows where the original difference is: https://stackoverflow.com/questions/54250126/why-do-i-get-platform-specific-result-for-stdexp. std::exp produces different output. later this small rounding difference is amplified by `nlopt` iterations. – jpo38 Jan 18 '19 at 15:22
  • @jpo38 Without using interval mathematics or similar, your result is "meh, maybe something like this". That is what you are asking for, and what you are getting. You'll have to hand-tune your results for the precision you are getting, because you asked for that, and your test (some exact result) demands you do that. In a general case, `==` is not an operation that can physically be done on actual real numbers. `double`s are mere approximatoins of actual real numbers; but even if you had perfect representations of real numbers you could not in general do what you are asking. – Yakk - Adam Nevraumont Jan 18 '19 at 15:26
  • @jpo38 in other question you do not explained how problem manifests. You ask why you have an issue and someone provide explanation. We need to know how you algorithm gets corrupted what goes wrong when results are a bit different. This is the place where issue should be fixed. One way to fix it is to calculate accuracy of the calculation and compere result using taking that accuracy into account. – Marek R Jan 18 '19 at 15:29
  • @MarekR: He very well explains how the problem manifests: exp gives different presults on different platforms and the small differences, in his calculations, finally result in quite noticeable differences. What is so hard to understand there? His solution is to set the low byte of the significan to 128 on all platforms. I personally think (have a feeling) that is good enough, but can't really tell without testing, and I can't test it, so no answer. Manipulating bytes of a floating point type is probably undefined behaviour, but that should not matter here. – Rudy Velthuis Jan 18 '19 at 18:29
  • @MarekR: he doesn't ask **why** he has an issue: he already found the answer explaining it. He asks if his solution is good enough. IMO, yes, if tests prove that the results are reliable, precise enough and the same on all relevant platforms. If that is so, who cares about undefined behaviour? I assume that soft floating point libraries are probably too slow. – Rudy Velthuis Jan 18 '19 at 18:34
  • Yes but if he could show what he exactly does, why he needs exact matching then it may turn out that his doing something wrong. Strict equality for floating points is quite unusual and it may indicate more general bug in the code. This is common mistake made by programmers (not only newbies). Also his attempt to solve issue is bad since if compared values are very close to rounding border each value can be rounded to different value. As a result he still will have this issue, but accruing very rarely (bug extremely hard to reproduce). – Marek R Jan 18 '19 at 18:35
  • @MarekR: **He** is not doing anything wrong. **Exp()** returns different results on different platofrms. That is what is wrong. He is merely trying to even out the effects of that. Strict equality may be unusual, but that is irrelevant here. One might expect the exact same calculation to give the exact same results with the exact same inputs (assuming all use IEEE754). Reality shows this is not always the case (especially because the implementations might slightly differ). He is trying to even this out and ISTM that if it produces the desired results on all relevent platforms, it is fine. – Rudy Velthuis Jan 18 '19 at 18:42
  • @RudyVelthuis: Thanks for arguing in my name ;-) You are right, even if we all know double results may differ due to floating point precision, I'm just asking if fixing the least significant byte is an acceptable way to fix it or not. Wherever the difference comes from and regardless the reason why you want to end up with the exact same double (meaning exact same bitset value in the end). – jpo38 Jan 18 '19 at 20:29
  • I tried to reformulate my question to make it less looking like a XY problem. – jpo38 Jan 18 '19 at 20:38
  • If you perform operations in a loop so that a small error at the beginning is greatly magnified at the end of the loop, then your algorithm has a serious problem. You need to use a process which averages out the error accumulation, OR, you need to track error accumulation so that at the end you can fairly say that you have only 20 accurate bits. – Zan Lynx Jan 18 '19 at 20:52
  • If a 1e-17 difference in an intermediate calculation changes a result from 2.36 ti 2.47, you need to look for a bug in the code ir examine the numerical stability if your algorithm. You should understand why you are getting the results you are before you try to round away deviations. Rounding may make your results less accurate. – Eric Postpischil Jan 18 '19 at 21:52

4 Answers4

5

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!

Ben Voigt
  • 277,958
  • 43
  • 419
  • 720
  • But if I round the two values to an arbitrary precision they are likely to differ by much more than 256 ULP, no? – jpo38 Jan 18 '19 at 22:27
  • @jpo38: It doesn't matter. You are not achieving your goal of making the results equal. Using the values as-is is better than rounding. If 256 ULP is acceptable, then so is using the values as-is. If your calculation is magnifying the existing uncertainty, it will magnify 256 ULP even more. – Ben Voigt Jan 18 '19 at 22:41
  • I think all of this is purely academical. We don't know what values the outcome should be, or if the calculations after that particular rounding are actually done badly or if they can be improved. Only OP knows if the values he gets are acceptable or not. If he thinks the values are good enough for his purposes, and the rounding guarantees he gets consistent results on different platforms, then the rounding is valid, even if that creates differences of up to128 ulp (or 7 out of 53 bits). ISTM that only OP can judge if what he gets is good enough/acceptable or not. – Rudy Velthuis Jan 19 '19 at 01:22
  • @Rudy: OP's proposed "fix" will increase the deviation in some cases, therefore anyone can judge that it does not meet the stated goal of making results uniform across platforms. – Ben Voigt Jan 19 '19 at 03:51
  • @BenVoigt: in what way? It is assumed that by his trick, the results are exactly the same across all platforms. That can easily be checked. If that is indeed the case, there is no deviation between platforms at all. – Rudy Velthuis Jan 19 '19 at 16:20
  • @RudyVelthuis: That assumption is wrong. My answer explains why, using an example. Two more recent answers have given almost identical examples. – Ben Voigt Jan 19 '19 at 16:21
3

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.

geza
  • 28,403
  • 6
  • 61
  • 135
  • Thanx for your answer. I investigated why I get a huge difference in the end, and it's not the rounding difference that is magnified. The algorithm does thousands and thousands of iterations and constatly takes decisions that modify the inputs of next iteration. At some point, it's on the edge of a threshold for a given decision, and then, the very small difference makes it cross the border....and than changes decision from "do something" to "do something else", and this will end up with huge difference on input for next iterations and then results (that remain both acceptable, but different) – jpo38 Jan 19 '19 at 08:11
  • @jpo38: What your comment describes is exactly "magnification of the difference". You have a non-linearity in your system (when you switch between "do something" to "do something else") and that non-linearity magnifies small differences. The thing you proposed to do in this question was to add another non-linearity, which will magnify small differences, and most importantly, it will not remove your existing non-linearity. – Ben Voigt Jan 19 '19 at 16:24
  • @jpo38: btw, what are you implementing? Fractal renderer? – geza Jan 19 '19 at 16:44
  • @geza: No, part of the program is nonlinear optimization which goes over a lot of iterations. Can't tell what the application is on a public forum ;-) – jpo38 Jan 19 '19 at 17:12
  • @jpo38: ok, I was just curious :) – geza Jan 19 '19 at 17:44
  • @geza: Sure. By the way, your proposal to use a special floating point library or not use floating point is smart but not suitable as I use a third marty for nonlinear optimization and this works with regular `double`. – jpo38 Jan 19 '19 at 17:51
  • @jpo38: Hmm, I'm not sure I understand you. I meant, that you keep using regular `double`, but for all non-basic operations (`exp`, `sin`, etc.), use a custom library. So, don't use `exp` which is in the standard library, but use a custom `libExp`. `libExp` can use `double`+basic operations, of course. – geza Jan 19 '19 at 18:30
  • @geza: I see. Actually I hoped that `boost` had an `exp` function that would have the same implementation on all platform but they did not. – jpo38 Jan 19 '19 at 18:56
2

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!!!

gnasher729
  • 51,477
  • 5
  • 75
  • 98
0

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.

Zan Lynx
  • 53,022
  • 10
  • 79
  • 131
  • There's no aliasing violation here. – Ben Voigt Jan 18 '19 at 22:00
  • @BenVoigt: I'm pretty sure the compiler is allowed to assume that, because a double was in use as a double, it isn't obligated to stop using it because you created a char* to it. I'd have to look it up in the standard but ONCE YOU START USING AN OBJECT AS ITS TYPE, it IS that type, and it doesn't have to be removed from registers and such, as long as an in memory byte copy is generated for char* casts. – Zan Lynx Jan 18 '19 at 22:06
  • Taking the address of an object does require that that object is flushed to memory. The compiler can avoid re-reading memory only if, under the as-if rule, the change cannot be detected by any legal program. This is a legal program that could detect the change, ergo the compiler has to reload `rounded` from memory so long as an alias to it exists, and write operations intervene that potentially alias. Since `char` can alias `double`, the operation `temp[0]` can potentially alias `rounded` (and here it certainly does), so the proposed behavior is not allowed under "as-if". – Ben Voigt Jan 18 '19 at 22:11
  • @BenVoigt I am *almost* entirely certain that I've seen clang generate code that ignored a write through a char* alias. – Zan Lynx Jan 18 '19 at 22:15
  • Well, the code in question is not actually `char*`, it's `unsigned char*`. ISRT the latest standard drafts allow aliasing only through unsigned narrow character and `std::byte`. – Ben Voigt Jan 18 '19 at 22:17
  • The only reading of the Standard I could see in which this code could go wrong, would be if `temp[0] = 128;` was considered to reuse the storage of `rounded`, thus ending the lifetime of the `double`. The programmer would then be obligated to recreate a `double` in the same location before it went out of scope, and unable to access `rounded` as a double until doing so. But that would make memcpy useless unless the dynamic type of the source was already the desired type of the destination, so that interpretation should be rejected. – Ben Voigt Jan 18 '19 at 22:23
  • @ZanLynx Then what would the right way to write this function? Do I simply need to `memcp` the `double` to an `unsignes char` array, modify it and later `memcp` it back to the `double` object? – jpo38 Jan 18 '19 at 22:32