A huge roadblock I recently ran into: It's possible, that floating operations may yield different results depending on the compiler or CPU architecture.
As far as I know it's not possible to run a simulation depending on floats on different machines, expecting the same output. Some people recommend switching to integer operations and to keep floating point arithmetic to a minimum. Is this the only way to get deterministic math across machines?
I ran into this while trying to procedurally generate a map based on a seed value. Heightmap, zones/biomes, rivers, towns, etc are all based on a seed. The result must be completely reproducible on different machines, otherwise players would see different maps. The terrain is large and therefore chunked. There's a seed per chunk and only seeds of chunks visible to the client are sent to him.
I now believe that the only way to implement this correctly is passing the full, complete state of a chunk to the player after the generation. There's no sane way to implement the whole generation with integer math only. But there's a positive side-effect: More traffic but less performance needed to regenerate the chunk on the client.
Actually it's simpler to work around this problem in a physics simulation, because the authoritative server can simply sync the properties across the players frequently.
Fixed-point integer is generally easier to handle but you have to do it "manually" as language support is scarce. And then watch out for overflow. It also may not be as fast as FP on modern systems (especially GPUs).
Arguably fixed-point's limitations are easier to understand. Floating will give you within epsilon of the "right" answer, but as we have seen this may vary in weird ways across platforms. Fixed will run into precision loss problems a lot earlier, forcing you to confront and deal with them. You also have to confront truncate vs. round earlier.
Mostly nonsense. You are making the same mistake as Matlab, throwing away all the internal Intel precision, just because you don't know the boundaries and limitations or you don't care.
In short Matlab did as you said in v6.0 and it was wrong.
Extended-precision mode is a disaster because it's confined to the FPU. If it were an 80-bit datatype it wouldn't be a problem. Because of this property, it's only predictable if you write all your math in assembler. If you write C, the generated code is very vulnerable to perturbation when the compiler changes when variables are written back to memory.
It may well give you more accurate results, but if you're doing something that's inherently numerically unstable (in our case, ASIC timing feeding into k-means clustering) it can give you different results in response to tiny code changes. Force-truncating trades accuracy for stability across changes to the software.
And obviously if you want C code that gives the same results on Intel and non-Intel hardware, you can't use Intel 80-bit FPU.
(Edit: I'll concede somewhat on "within epsilon of the right answer", because there are some examples in the Matlab case where 80 bit gets much closer to the right answer to 64-bit.)
As it stands, pjc50 has some on-point advice that directly addresses OPs concerns with consistent results from seed values, and testing. Your reply is: "that is nonsense."
If you want compatibility between x86 CPUs and ARM/PowerPC/MIPS then you have to not use the extended 80-bit results. The 80-bit results can be more accurate (although not always, it depends on exactly how they are used) but the OP wasn't concerned with maximum accuracy but with exact repeatability.
Similarly, if you want cross-architecture consistency you need to not use fmadd instructions because fmadd gives different results from separate multiply and add. The fmadd results are more accurate, but that is not the point in this case.
There's an additional caveat that may be of interest.
If you execute code in parallel, you may pick up different round-off errors unless you take care to combine the results from parallel tasks in the same order every time.
This can be a problem when simulating non-linear dynamical systems (such as the Solar System).
The important question to ask is: do you really need bit-accurate results? If so, then floating point arithmetic isn't probably what you want. However, running the same binary on different (compatible, ie. not ARM vs. x86) CPU architectures should produce identical results almost always (can anyone give a counterexample?).
But in reality, it is often acceptable to have one or two least significant bits differing. Especially when computing with doubles, a relative error in the order of 2^-53 = ~10^-16 isn't really significant. Even 32 bit floats give a precision of 2^-24 = 5*10^-8.
This is usually good enough for e.g. games, although multiplayer games need some kind of authoritative server model to avoid cheating.
Finally, the x87 floating point unit shouldn't be used any more because SSE2 is ubiquitous at this point so issues with 80 bit internal precision should be a thing of the past.
I'm not sure 80-bit floats are a thing of the past. Many calculations still benefit from the extended precision compared to 64-bit. The issue also crops up agian on some platforms that can do 128-bit float, but are then truncated down to 64-bits for storage.
I agree there are ways to avoid that problem, but the general issue of truncation across storage/calculation I think remains an important concept.
> Many calculations still benefit from the extended precision compared to 64-bit.
Out of curiosity, can you name some examples?
64 bit doubles are really precise, to give a practical figure: a 64 bit float is good enough to calculate the orbit of Neptune to a precision of half a millimiter (0.497 mm). That's a distance of about 30 AU = 4.5 * 10^12 m = four and a half billion kilometers.
> The issue also crops up agian on some platforms that can do 128-bit float, but are then truncated down to 64-bits for storage.
Which CPU architectures use 128 bit floats internally? Is it common to use such platforms?
What comes to x87, I'd imagine that the difference in speed in comparison to SSE is already a good enough reason to not use 80 bit floats in most practical applications.
Perhaps some nuclear simulations or particle accelerator physics deal with numbers that vary in magnitude so greatly that increased precision may be useful but these cases are really rare (albeit important).
In short Java went ahead and declared imprecise and architecture-independent float results better than more precise and architecture dependent results.
Matlab followed suite with v6.0 and found big errors in mxmuleps compared to Matlab v5.0.
Such a little rounding problem in the 6th digit e.g. caused a big problem in one architectural problem of mine, where it did cost several millions. Check for the costs of rounding problem in the real world, comp.risks e.g.
https://www.ima.umn.edu/~arnold/disasters/
Basically everything that can be done with 80-bit floats can be done with techniques such as compensated summation and double-double arithmetic. This requires several 64-bit operations per operation, but those operations can take advantage of SIMD capabilities when processing several inputs in parallel.
This is the case for example with the x86 .NET jit. The problem is usually when values are moved to 64bit locations rather than keeping them in 80bit fpu registers. This was a huge headache our regression tests for a .NET engineering/FEM application. .NET has no option to force "strict" FPU mode so there is no way of getting reproducible results. What saved us was that the x64 JIT uses only 64 bit, which makes it deterministic.
As long as you avoid unspecified instructions (reciprocal square-root estimate, fsin) and avoid library functions or code-gen that takes different paths when SSE2/AVX/whatever are available and as long as you run the same binary on every machine (compatibility between PowerPC and x86 is much harder) you can definitely do it. Many games have.
So, if you can restrict yourself to x86 CPUs it's manageable, at least for modest sized chunks of code.
I once ran into this with python. 2 interpreters giving different results. Turns out they compiled 1 interpreter with some --I-don-t-care-about-correct-flops-as-long-as-they-are-fast flags. Still traumatised.
I (the author) agree. I'm mainly trying to point out all the issues with floating point rather than give a full solution. I believe the problem stems from people just not realizing a lot of these problems. Certainly I'm not offering the complete solution here, and for any one of those points could spawn off several articles just to give ideas.
So yes, links are appreciated so readers can get more into the details.
True, but in this case it was to test for near-zero values so absolute error comparison is good enough.
Here's a few macros I have been using to compare floating point values for approximate equality. It should be true if the values are equal to about 6-7 decimal places (in base-10). I use this in testing some code that deals with algorithms that have known numerical precision issues. The errors are inherent to the algorithms, not specifically floating point errors (e.g. solving Kepler's equation for near-parabolic trajectories).
Note that it isn't specified whether "a" or "b" is the "gold" value to compare against, so it's measuring the relative error against a rough average of the two magnitudes.
but in this case it was to test for near-zero values so absolute error comparison is good enough
Unless you are working with small numbers? Suppose you are dividing 1e-30 by 1e-30 than the result should still be 1 so if you return 0 when the second argument is < 1e-6 then you're screwed.
I think boost's close_at_tolerance does what your macro tries to do; basically you say 'here's two numbers, check if their values differ no more than x% of their value' With the advantage that it is not a macro and guards against overflow etc: http://www.boost.org/doc/libs/1_37_0/boost/test/floating_poi....
> Unless you are working with small numbers? Suppose you are dividing 1e-30 by 1e-30 than the result should still be 1 so if you return 0 when the second argument is < 1e-6 then you're screwed.
A very good point which I didn't think of! Obviously the 1e-6 figure is specific to the application.
How should one avoid division by zero then? You can't figure out the relative error without a division, can you? For normal comparisons, using a relative error is the way to go but the zero case is different.
On the other hand, using quantities of magnitudes in 1e-30 is already using half the exponent, a change of units may be due here.
Good question and in this case I really can't tell, it depends on the problem domain (i.e. is your normal data set numbers greater than 1 or rather smaller then 1e-100) and I'm not sure if is possible (in a feasible way) to write one generic function which works for all cases.
I've seen those "abs(fp_number)<some_tiny_value" in front of divisions, because people fear that the division will "blow up" in some way if a zero or sufficiently small number is passed to it.
But it's quite pointless, as the semantics of floating point calculations give you the indication about overflow right in the division result, it will be marked as positive or negative infinity.
I am not sure I understand what you are saying exactly, but I believe in some circumstances you may be working with numerical values which you know to have errors > ulps.
For example, suppose you working with data that you've obtained from some physical measurement, and your knowledge of the measurement process means you know it will have around 1% relative error.
Or another example: suppose you are solving a linear system of equations Ax=b with highly precise right-hand-side data b (< \epsilon floating-point error), but the condition number of the matrix is large. Then perhaps it is only meaningful to expect the solution to be accurate to within k(A) \epsilon, the condition number of A.
(I haven't thought about this stuff for a while, please correct me if this is wrong).
edit: here's another reason. Suppose there's a heavy tradeoff between runtime and accuracy for some algorithm, so for some application it makes sense to use some cheap but somewhat inaccurate approximation. Then you're probably going to have and expect far larger errors. But, you might still be able to get some kind of (worst case?) bounds on the result, and want to assert something about that in test cases or sanity checks.
Good article! Just a couple corrections / comments:
> n * x != x + … + x
Not true if n == 2! For n = 2 they are exactly the same.
Why is this? Because "n * 2" and "n + n" are both exactly one operation, which means they only get rounded once, and since both are mathematically identical they both get rounded to the same value.
> I’m not sure why, but the default string formatting of floating point numbers in most languages, even domain-specific math languages, does not have high enough precision to store the complete number.
I think it's probably related to the fact that printf() doesn't have any format specifiers for float that do What You Really Want, which is: print the shortest number that will parse back to the right value.
So unless people want to implement this algorithm themselves (and, fascinatingly, the best known algorithm to do this was published only five years ago: http://dl.acm.org/citation.cfm?doid=1809028.1806623), you have to choose between having all your numbers be annoyingly long or losing precision for some numbers.
> Not true if n == 2! For n = 2 they are exactly the same.
Not true in many other cases either; a trivial example is n=1, or x=1 and n can be represented exactly. In all of the cases the "!=" in the headings is meant to be read "is not necessarily equal" or "there exist expressions e, e' for which e != e'"
> I think it's probably related to the fact that printf() doesn't have any format specifiers for float that do What You Really Want, which is: print the shortest number that will parse back to the right value.
Lots of language support the proper behaviour these days. Haskell/ghc does, for example.
> Historically, the Python prompt and built-in repr() function would choose the one with 17 significant digits, 0.10000000000000001. Starting with Python 3.1, Python (on most systems) is now able to choose the shortest of these and simply display 0.1.
I can't count the amount of times I see this text being posted. Yet it takes something like a degree in math/CS/... to read and understand it completely; it is really not suitable for everyone.
As such I bookmarked Essential facts about floating point calculations so next time a psychologist or biologists comes asking me if I already knew matlab has a bug when comparing 0 to 0, I can point him/her to this document. And some proper code to compare floating point numbers using a relative error or so.
I do mention if you simply increase the length of the formatted string you can do round-trip. It should be noted however (as I do) that the decimal number is still not equivalent to the actual floating point number, only that when converted back it is the same. This point could, in theory, be an issue if your serialized form crosses systems that use different floating point systems.
If you are trying to get exact compatibility between systems with different floating-point systems then serializing to text is the least of your problems.
It should be noted however (as I do) that the decimal number is still not equivalent to the actual floating point number, only that when converted back it is the same.
Since 2 divides evenly into 10, every binary number has an exact decimal representation, and you can produce a decimal string which is equivalent to any floating point number. You may need quite a few digits, though.
> It is possible to get the same number back though, even when using decimal. If the truncation length is simply extended long enough the parsing will yield the same number back.
But I agree that actual code would have been clearer.
In cases where the iteration can’t be avoided, it’s vital to understand you won’t get the expected result. A certain amount of rounding error will always accumulate.
Maybe I don't fully understand the issues at play here but I remember reading about so-called 'unums'[0], an alternative binary number representation that seems quite amazing. I don't understand why it doesn't get more attention, apart from the fact it would need to be supported in CPU hardware to be of any real use.
As far as I know it's not possible to run a simulation depending on floats on different machines, expecting the same output. Some people recommend switching to integer operations and to keep floating point arithmetic to a minimum. Is this the only way to get deterministic math across machines?
http://gafferongames.com/networking-for-game-programmers/flo...
Edit:
I ran into this while trying to procedurally generate a map based on a seed value. Heightmap, zones/biomes, rivers, towns, etc are all based on a seed. The result must be completely reproducible on different machines, otherwise players would see different maps. The terrain is large and therefore chunked. There's a seed per chunk and only seeds of chunks visible to the client are sent to him.
I now believe that the only way to implement this correctly is passing the full, complete state of a chunk to the player after the generation. There's no sane way to implement the whole generation with integer math only. But there's a positive side-effect: More traffic but less performance needed to regenerate the chunk on the client.
Actually it's simpler to work around this problem in a physics simulation, because the authoritative server can simply sync the properties across the players frequently.