My C code works with -O3 but not with -O0
mulle-kybernetik.com
mulle-kybernetik.com
Casting the value to double ends up converting the long value 0x7fffffffffffffff to the nearest double value: 0x8000000000000000. As the -O0 version CORRECTLY reports, this does not round-trip back to the same value in the "long" type. Many other values, though not all, down to about 1/1024 of that value (1 / 2^(63-53)) will also fail to round-trip for similar reasons.
Unless my coffee-deficient brain is missing something at the moment, it should be the case that any integer with 53 bits or fewer between the first and last 1 bit (inclusive) will roundtrip cleanly. Any other integer will not.
Edit: fixed a typo above, and coded up the idea I expected to work and ran it through quickcheck for a few min, and this version seems to be correct ('int' return rather than bool is just because haskell's ffi doesn't natively support C99 bool):
#include <limits.h>
int fits(long x) {
if (x == LONG_MIN) return 0;
unsigned long ux = x < 0 ? -x : x;
while (ux > 0x1fffffffffffffUL && !(ux & 1)) {
ux /= 2;
}
return ux <= 0x1fffffffffffffUL;
}Add enough zeroes, and you’ll run out of exponent range.
From what I understand from the spec it should be the nearest in value, no? Not the nearest in memory representation.
The language spec (or at least a summary of it) is linked in the article, and it is pretty loose: nearest higher or nearest lower, chosen by the implementation (regardless of which is nearer).
IEEE double-precision float cannot store LONG_MAX (which is 2^63-1 with 8-byte longs) precisely, so it gets converted to 2^63; which you cannot cast back to a long, because it doesn't fit (resulting in undefined behaviour).
if( value < (double) LONG_MIN || value > (double) LONG_MAX)
return( 0);
Which then leads to this line:
l_val = (long) value;
When -O0 and -O3 give different results, either at least one of them is incorrect and you've stumbled on a compiler bug, or both of them are correct and you're invoking UB (the far more common situation).
EDIT: no, I think I misunderstood it: it's not that value is larger than (double) LONG_MAX, it IS (double)LONG_MAX, so of course "value > (double)LONG_MAX" is false.
The problem is that the "(long)((double)LONG_MAX))" is undefined behaviour on implementations that round (double)LONG_MAX upwards instead of downwards. Which is allowed by the standard. Ok, cool :)
"When a value of integer type is converted to a real floating type, if the value being converted can be represented exactly in the new type, it is unchanged. If the value being converted is in the range of values that can be represented but cannot be represented exactly, the result is either the nearest higher or nearest lower representable value, chosen in an implementation-defined manner. If the value being converted is outside the range of values that can be represented, the behavior is undefined."
See it says implementation-defined manner, not according to the current rounding mode.
Test case:
#pragma STDC FENV_ACCESS ON
#include <stdio.h>
#include <stdint.h>
#include <fenv.h>
int main()
{
fesetround(FE_DOWNWARD);
printf("%f\n", (double)UINT64_MAX);
fesetround(FE_UPWARD);
printf("%f\n", (double)UINT64_MAX);
return 0;
}
$ gcc -std=c99 -frounding-math x.c -lm
$ ./a.out
18446744073709551616.000000
18446744073709551616.000000
bool double_to_uint64 (double x, uint64_t *out)
{
double y = trunc(x);
if (y >= 0.0 && y < ldexp(1.0, 64)) {
*out = y;
return true;
} else {
return false;
}
}
Explanation:
- we first round the value to an integer (but still a floating point value),
- we then check that this integer is in the valid range of the target integer type by comparing it with exact integer values (0 and 2^N),
- if the check passes, then converting this integer to the target integer type is safe, and if the check fails, then conversion is not possible.
Of course if you literally need to convert to "long" you have a problem because the width of "long" is not known, but that is a rather different concern. I argue types with unknown width like "long" should almost never be used anyway.
(based on my answer here: https://stackoverflow.com/questions/8905246/how-to-check-if-...)
The answer by ambrop7 below solves a different problem but appears to be correct. The reason escaped me at first and is super subtle. The trick is that LONG_MAX is 2^63 - 1, not 2^63. And the subtlety is that 2^63 is guaranteed to be exactly representable in IEEE double because it is an even power of 2, which 2^63-1 is not.
I don't care much for runtime ldexp() anyhow. So I'd be tempted to just pre-compute the exact limits -2^63 and nextafter(+2^63, 0) and encode them as doubles manually (omitting some #if method of portably determining the width of LONG):
#define FLONG_MIN 9223372036854775808.0 // exact -2^63
#define FLONG_MAX 9223372036854774784.0 // exact nextbefore(2^63)
if (value < FLONG_MIN || value > FLONG_MAX)
return(0);
When you work with floating point, you need to remember you work with a tolerance to epsilon for comparisons because you are rounding to 1/n^2 precision and different floating point units perform the conversion in different ways.
You must abandon the idea of '==' for floats.
This is why his code is unpredictable, because you cannot guarantee the conversion of any integer to and from float is the same number. Period. The LSBs of the mantissa can and do change, which is why we mask to a precision or use signal-to-noise comparisons when evaluation bit drift between FP computations.
He has the first part correct, < and > are your friend with FP. But to get past the '==' hurdle, he needs to define his tolerance, the code should be something like:
if (fabs(f1 - f2) > TOLERANCE) ... fits = true.
I was irked by his arrogance when he asks, "Intel CPUs have a history of bugs. Did I hit one of those?" First, learn about floating point, then, work on an FPUnit team for 10 years, and even then, don't assume you're smarter than a team of floating point architects, you're not.
Be aware that +0.0 and -0.0 are different floating point values but represent the same real number, so +0.0 == -0.0 follows.
People who say == means nothing for floating point and you always need epsilon checks are wrong, plain and simple. == is very well defined. Don't confuse the definition of floating point operations with common practices for using them effectively.
You can iterate through all non-NaN values and check that successive ones are indeed not equal:
#include <math.h>
#include <stdint.h>
#include <assert.h>
#include <stdio.h>
#include <inttypes.h>
int main()
{
float x = (float)-INFINITY;
uint64_t count = 1;
while (x != (float)INFINITY) {
float y = nextafterf(x, (float)INFINITY);
assert(y != x);
x = y;
++count;
}
printf("Found %" PRIu64 " floats.\n", count);
return 0;
}
$ gcc -std=c99 -O3 a.c -lm -o a
$ ./a
Found 4278190081 floats.
Interestingly, this only finds one zero (-0.0), hence the assert doesn't actually fail around zero.
(a) get x casted to a float
(b) get `div` instead (like in Python 2) and obtain a false result
(c) get cursed out with a type error.
These three options are available because it's possible to determine whether a given real number is representable as an int or not. This is not possible with floats.
What you can do is:
- determine if a float is an integer (trunc(x) == x),
- convert a float to a certain integer type with some kind of rounding, or get an error if it's out of range (see my comment with double_to_uint64),
- convert a float to a certain integer type exactly, or get an error if it's not representable (e.g. by doing both of the above).
The basic reason that so many people fail to use floats correctly is that they act like operations on floats are equivalent to operations on the real numbers they represent, when in fact they are usually defined as the operation on real numbers rounded to a representable value.
No system of finite representations is able to express arbitrary real numbers.
But paradoxically you can store “any” real number, eg the number made up of all future lottery draws.
It's kind of a problem since most programmers aren't writing programs to do scientific/modeling. They're writing programs that are doing decimal math and currency where IEEE 754 is inappropriate.
And I'll repeat what the parent said, again: decimal math.
>Most programs aren't doing "decimal math"
And yet they do. Most programs do decimal math with floats, either because the programmers don't know that's not a thing (and consider floats the same as decimals), or because the language doesn't give them anything else (e.g. Javascript up until recently), or because workarounds are cumbersome (e.g. representing as fractions or working with integers for money with a scale factor and downscaling in the end, etc).
But in almost all programs, the kind of math the program needs to be done, from money handling, to estimating distances, to layout, etc are better done as decimal.
Floating point don't represent a way to solve any kind of actual problem, except the problem of speed and memory representation. So their use is because of computer constraints, not because the problem domain requires floating calculations.
Decimal math isn't built into most of the languages or processors we use, so that's why I keep asking what you'll use. It's a kind of hand-wavy response. You'll need to seek out an appropriate library and use it consistently, which is a lot more trouble than just using the built-in IEEE. Not impossible, just usually more trouble than it's worth.
Whether it changes everything is neither here not there.
It might still result in the same output (e.g. the lack of precision might not matter for layout purposes), but there's absolutely no domain need to do the calculations in FP except for performance or lack of a better number type in the language.
In all of these cases, the exact match would be better than lack of precision, aside from the performance and similar (memory, etc) costs.
That is, if decimal was computing-wise as fast and as available as FP, there would be absolutely no reason to use FP. It's not a choice of mathematical need (e.g. not imposed by the calculation) but a choice of computing constraints.
>Decimal math isn't built into most of the languages or processors we use, so that's why I keep asking what you'll use.
In my experience many languages have it (C#, Java, Python, Javascript is getting it, etc) and people don't use it.
And in most cases, being in the processors wouldn't matter - we use FP for all kinds of non-critical performance calculations, not because it's needed for the speedup, but because it's there.
When there is indeed a need for the CPU-support, I'd use FP myself, sure. But not in most other cases, and absolutely not for money and other cases where precision matters (but people still use FP).
bfloat.
incompatible with and far superior to IEEE 754 for one very obvious application. surprised you missed it.
These days I write my floating point code assuming float is a 32-bit IEEE 754 floating point number, and double is 64-bit. You can get those guarantees on any desktop hardware with the right flags, and the same semantics are commonly available on other platforms too.
Picking a well-defined implementation makes it much easier to reason about conversions between integer and floating point types. In fact, it allows you to reason about a lot of operations, e.g. 1.0 + 2.0 == 3.0 is true; (float)183 == 183.0 is true; 0.1 == (double)0.1f is false; etc.
Totally false. Any integer between -(2^53) and 2^53 can be converted to IEEE 64-bit double without any loss of information.
Those are a specific integer range, not "any integer".
clang test.c -O0 -fsanitize=undefined
./a.out
[...]
test.c:17:12: runtime error: 9.22337e+18 is outside the range of representable values of type 'long'
Interestingly gcc doesn't throw that warning.
There also exists low-cost random-sampling-based ASAN implementation that can be enabled in production: Google uses GWP-ASAN for all server-side applications as well as Chrome on Windows/Mac. See https://www.youtube.com/watch?v=RQGWMLkwrKc for details.
Here 14%, https://www.jetbrains.com/lp/devecosystem-2019/cpp/
Here 40 - 55%, https://www.bfilipek.com/2019/12/cpp-status-2019.html
At CppCon 2015 or something, at Herb's question during his keynote, about 1% of the audience as per his comment on the video.
The reason they aren't enabled by default is that that's not what they're designed for. They have a significant performance impact, you can't enable them all at once, they conflict with other security features and they may introduce security issues.
These are developer features. They aren't there to run your production code, they are there to test during development and bug finding.
$ g++ -x c++ -fsanitize=undefined -
#include <iostream>
int main() {
int *a = nullptr;
std::cout << std::addressof(*a) << std::endl;
}
^D
$ ./a.out
<stdin>:5:29: runtime error: reference binding to null pointer of type 'int'
0> neither that operator nor the & operator is evaluated and the result is as if both were omitted, except that the constraints on the operators still apply and the result is not an lvalue
Can you cite me something for that? The very first mention of dereferencing in the C++03 standard - ISO/IEC 14882:2003 1.9 Program exection ¶ 4 (page 5) - would seem to disagree:
> Certain other operations are described in this International Standard as undefined (for example, the effect of dereferencing the null pointer). [Note: this International Standard imposes no requirements on the behavior of programs that contain undefined behavior. ]
EDIT: As does 8.3.2 References ¶ 4 (page 136):
> [...] [Note: in particular, a null reference cannot exist in a well-defined program, because the only way to create such a reference would be to bind it to the “object” obtained by dereferencing a null pointer, which causes undefined behavior. As described in 9.6, a reference cannot be bound directly to a bit-field. ]
> Notes from the October 2003 meeting:
> ...
> We agreed that the approach in the standard seems okay: p = 0; *p; is not inherently an error. An lvalue-to-rvalue conversion would give it undefined behavior.
http://open-std.org/JTC1/SC22/WG21/docs/cwg_active.html#232
WRT your edit: no, that says dereferencing a null pointer and binding a reference to the result produces undefined behaviour. That agrees with what I was saying. "which" refers to the whole of "to bind it to the "object" obtained by dereferencing a null pointer", not just to "dereferencing a null pointer".
sizeof on invalid pointer dereferences is also pretty common.
This sounds like a good UB learning experience. The parent comments sounds like the author is running into UB more than they really want to.
Note that ruining sanitizers in prod might be insecure. They're for development.
On the other hand, the clang static analyzer may have a shocking amount of messages when run first on a large existing code base, and some of those warnings can be considered more "opinions" than warnings. It's still makes sense and is very rewarding to make a code base "static analyzer clean".
The runtime sanitizers in comparison are very precise and always pointed to actual "sleeper bugs", it's almost definitely a good idea to use them and take their warnings serious.
But anyway, clang ASAN, UBSAN, TSAN and the static analyzer are all really excellent and important tools for everybody writing C or C++ code.
PS: the reason why those checks are optional is that they increase compilation time (sometimes dramatically, like 10x slower compilation or more), and they add runtime instrumentation code which both increases the executables size and decreases performance dramatically (also 2..10x times or more, although the clang sanitizers are really quite fast compared to other solutions).
The best projects are the ones that start out from day 1 with all the pedantic warnings turned on, warnings-as-errors, static analysis and sanitizers run as part of the automated build and any peep out of them gets treated just like an error would get treated. When you start the project out that way, there isn't that de-motivating initial hump to get over.
(TIL - thanks!)
Current codebase generates around 3000 casting errors. I doubt I'm the only one with this kind of "history" to deal with AND the application is crashing with memory access errors.
What's my cleanup plan for this? Multiple compilers with every warning enabled on multiple platforms. Even a platform we're not targeting. We've mapped out which files and which lines have the biggest code smell. Now we have a giant map on a 65" tv to guide us.
Why all this attention? We're moving this nightmare from 32 to 64 bits. Parts of it were originally 16bit which have already been updated. Cast errors alone are now signposts to other bad code.
It helps that you then almost certainly have buy-in and are allowed to treat this as a priority because it's tied to work that the business is willing to prioritize. You're already over the hurdle of "customers don't pay to fix compiler warnings".
Getting rid of the warnings and more importantly the related other debris means being able to use enhancements we make each week. Well... after qa/test have done their thing.
Been good practice for my own indie dev. Treat every warning as an error and make a habit of running a detailed reporting build each day or week. The trick of using other platforms and compilers I brought from years back.
* Probably because the Ada practitioners are all trapped in Scifs somewhere in the inner mantle of the Earth.
It's pretty much established that even expert C and C++ programmers, especially for larger code bases, will end up making some sort of mistake that will cause a security vulnerability or undefined behavior.
ripgrep cannot be a drop-in replacement. Despite that, it can certainly replace grep in a wide variety of use cases. See: https://github.com/BurntSushi/ripgrep/blob/master/FAQ.md#pos...
I used Delphi on Windows. Most of the code I wrote in the last 20 years is in Pascal.
And that works really well on Windows. You have a stable API, and it does not matter what language the API is written in. Any language can use the API in the same way.
I tried to run some old projects this month. My 20 year old Delphi Windows programs run better in WINE on Linux than most programs I wrote 5 years using Linux tools, because the libraries have changed, but the API has not
20 years ago it was advertised as the safe C alternative.
This should be fairly obvious with knowledge about how floating point numbers are represented internally IMO.
Edit: Be more precise about what can be represented.
A good way to get around this is reading https://floating-point-gui.de/ to weed out any preconceptions, but yeah it's difficult to steer novices there without them stepping on one of the pitfalls first.
It's also worth noting that every finite double with magnitude larger than 2^52 has a precise integer value; it's just that once you get beyond 2^53, not every integer is representable.
[edit]
Since the extra value is precisely a power of 2 (-2^52), then it will round correctly, however the value is arguably not precisely -2^52 since it has an epsilon of greater than 1.
To be precise: every integer, positive or negative, with magnitude less than 2^53+1, is exactly represented in double-precision. “Extra values in two’s complement” don’t (and couldn’t possibly) effect this at all, since it is a statement about abstract integers and floating-point numbers, neither of which depends on two’s complement representations.
In particular, -2^52 has a sign field of 1, an exponent field of 1023+52=1075, and an all-zero significand field. This number is exactly -2^52.
int fits_long(double value)
{
double max = (double) (1L << 52)
double min = (double) -max;
return min <= value && value <= max;
}
I think this is a good solution without using the math library and having access to the FP status registers. Otherwise the lrint solution seems better.
int fits_long(double value)
{
double max = (double) (1L << 52)
double min = (double) -max;
if( min <= value && value <= max)
return( (double) (long) value == value);
return( 0);
}The value 2^80 can be exactly represented in double-precision as well, but it is also the most precise representation of 2^28 other integers, so it is ambiguous which integer is being represented.
[edit]
To clarify, the representation you suggest would also be the best possible representation of -2^52-1.
Analogously, `float x = 2.25; int a = x;` assigns the value two to a, but this does not imply that the integer two also represents 2.25. Two is just two.
Integers are inherently different because calculations with integers are naturally discrete, while floating-point calculations are a discrete approximation of the reals, which are not discrete.
There are two general purposes for floating-point numbers: scientific computing, where you start with an imprecise value and the precision multiplicatively accumulates (the original use). And as a hardware optimization for calculation of non-integer values (a common use-case today). In neither case does it make sense to treat a floating-point value as a precise number, which matches common advice to not compare floating-point values by equality.
But I also see that StehpenCannon has spoken :-); there is going to be very little chance he is wrong when it comes to arguing about floating point! He also notes the following, but not quite so explicitly...
The problem is poorly posed because the question is, which mantissas, exponents and sign bits fit into a long.
The algorithm is backwards because simple comparisons in integer space cannot compute this; but I think the algorithm should be,
int fits_long(double value) {
unsigned int trailingZeros = countOfTrailingZeros(mantissaOf(value));
bool fits = (bitsInAMantissa + 2 - trailingZeros + exponentOf(value) < bitsInALong);
return fits;
}
if double is always > 0 and going to an unsigned long, the 2 above would be 1. The 2 represents the sign bit in the double. Given that sizeof(double) == sizeof(long) all bits in the two representations are accounted for, so there is no information loss.
checking:
mantissa = 0 (with a hidden msb of 1) means, trailingZeros = bitsInAMantissa -> fits will be true when the exponent value can be 0..62. So this represents each of the +/- 2^exponent values
mantissa = 1 x bitsInAMantissa means, trailingZeros = 0 -> fits will be true when the exponent value can be 0..(bitsInALong-BitsInAMantissa-2). So, +/- (all ones) * 2^exponent value.
The representation for going from double to long is not "smooth".
]
countOfTrailingZeros() is a favorite of bit twiddlers. "Hackers Delight" or "bithacks" https://graphics.stanford.edu/~seander/bithacks.html
[edited to improve readability]
My experience has been different -- forcing SSE instructions gives me a different result on some math calculations. Core2 cpu, boost odeint calculations. Clang or gcc.
Do you have a reference for why it's rare?
Edit: check for me whether just calling lrint(x) works. The manpage doesn't specify that lrint() will set FE_INEXACT, but it seems weird to me that it wouldn't.
> The floating-point environment has thread storage duration. The initial state for a thread's floating-point environment is the current state of the floating-point environment of the thread that creates it at the time of creation.
Great, thanks, now I have to go back and restart some of those code reviews I've been doing of certain third party matrix math libraries...
Annex F:
The lrint and llrint functions provide floating-to-integer conversion as prescribed by IEC 60559. They round according to the current rounding direction. If the rounded value is outside the range of the return type, the numeric result is unspecified and the ''invalid'' floating-point exception is raised. When they raise no other floating-point exception and the result differs from the argument, they raise the ''inexact'' floating-point exception.
I don't know if this is supposed to be a joke or part of the setup for an explanatory post about undefined behaviour, but that list is in exactly the wrong order.
What you're seeing is not excess precision due to wide registers but excess precision due to optimization and constant propagation, which means GCC calculates a fast path for (argc == 1) that doesn't round correctly and ends up with "it fits".
Interestingly it does optimize to the correct "doesn't fit" with -mfpmath=387 -fexcess-precision=standard, so I guess this is a bug in how GCC treats SSE math. The sanitizer (-fsanitize=float-cast-overflow) also notices the problem.
if( ! fits)
#include <math.h>
#include <fenv.h>
int fits_long( double d)
{
long l_val;
double d_val;
// may be needed ?
// #pragma STDC FENV_ACCESS ON
feclearexcept( FE_INVALID);
l_val = lrint( d);
d_val = (double) l_val;
if( fetestexcept( FE_INVALID))
return( 0);
return( d_val == d);
}
Comparing floats is more subtle than most programmers realize, and there really isn't a one-size-fits-all solution.
Things to consider due to the nature of fp representation - comparing results close to zero is different (i.e. "is a small" needs a differen test than "are a & b close"
- the distance between fp numbers depends on their magnitude, so comparing two large numbers to each other shouldn't have the same bounds as comparing two numbers near 1, say[1].
- if you aren't quite careful you can easily create tests where a == b but b != a , which can cause sorting issues, etc.
Hand-wavily speaking if you want to do this "right", you should probably look at doing the analysis in ULP (units in last place) rather than directly on the floats. Don't do it for values near zero though. And have a fast path for differently signed values.
The above doesn't even get into denormalized values.
[1] note that what people usually mean for epsilon is the version of machine epsilon that is the difference between 1 and the next representable float above 1 [2]. So by definition this is smaller than the representable difference between any two numbers in larger decades
[2] MS .NET somewhat confusingly defines Epsilon as the smallest representable normalized number.
I'd say the cleanest would be to decode exponent and mantissa, check if the exponent is within the 64-bit limit of long, then check if there's any bits set below the decimal point. (+plus some extra care for two's complement negative numbers)
The problem with this is of course that this would be platform dependent.
Is there a way to get the largest double smaller or equal than some positive integer?
I really dislike the arrogant programmer trope. Can we all stop?