4

If we use algorithms with double and float arithmetic, how can we guarantee that the results are the same running it in Python and C, in x86 and x64 Linux and Windows computers and ARM microcontrollers?

We re using an algorithm that uses:

  • double + double
  • double + float
  • double exp(double)
  • float * float

On the same computer, compiling it for x86 and x64 MinGW gives different results. The algorithm makes a lot of math so any small error will make a difference in the end.

Right now the ARM mcu implementation gives the same results as the x86, but after seeing this I'm not sure whether is right or not.

EDIT

Precision loss is not a problem in this case, as long as it's the same in all implementations

EDIT 2

I'm finding these links very helpful, some hints are already in the comments:

rnunes
  • 2,785
  • 7
  • 28
  • 56
  • I think you can only rely on a very specific configuration, since C doesn't mandate a specific implementation of floating point (dunno about python) – Eugene Sh. Oct 17 '17 at 17:51
  • Are you doing anything with integers as well or do you keep everything as double and floats? – Foon Oct 17 '17 at 17:53
  • 3
    Generally the issue on x86 is that the x87 registers are 80 bit, so, as long as the compiler doesn't spill them, you get higher precision than what would be "correct" for 64-bit double precision arithmetic. Try to use `-ffloat-store` and see if x86 and x86_64 give the same results. – Matteo Italia Oct 17 '17 at 17:58
  • Why is determinism so important? – klutt Oct 17 '17 at 17:59
  • @Foon just floats and doubles – rnunes Oct 17 '17 at 18:08
  • @klutt we implemennt most of our algorithms in Python but some of them need to run in an ARM mcu as well. the output needs to be the same in every implementation, for the same input – rnunes Oct 17 '17 at 18:09
  • 1
    If you need the consistency that badly, you will have to give up on native FP implementations and use some homebrew/thirdparty portable libraries. And I would recommend fixed point... – Eugene Sh. Oct 17 '17 at 18:11
  • Different platforms define different sizes for different data types. Use stdint library and things like `uint8_t` to ensure you have a desired size. – Burstful Oct 17 '17 at 18:17
  • 2
    Make all FP variables `volatile` and never use a compound equation: `sum = a+b*c` --> `product = b*c; y = a + product;` Insure rounding modes match and disable FP optimizations. For math libraries - roll your own. This will significantly reduce variation and reduce performance. – chux - Reinstate Monica Oct 17 '17 at 18:17
  • @CoreyLakey we already guaranteed that – rnunes Oct 17 '17 at 18:19
  • Hmm.. although this cannot be done in some embedded environments, turning off optimizations may help on your Windows and Linux platforms. This sounds like it could be compiler specific. – Burstful Oct 17 '17 at 18:24
  • Consider using Java and `strictfp`/`StrictMath`. – user2357112 Oct 17 '17 at 18:26
  • "consider using Java".. Horrible answer to an issue involving cross-platform implementation to EMBEDDED. – Burstful Oct 17 '17 at 18:27
  • @CoreyLakey I didn't want do to that since it doubles processing time. Also, I believe it gives the same results (in the MCU) with or without optimization, tomorrow I'll test it out – rnunes Oct 17 '17 at 18:30
  • I would see if the results are different despite processing time, and if you can "fix" the issue with optimization turned off then start your journey into assembly differences :) – Burstful Oct 17 '17 at 18:33
  • 1
    Could it be down to differences in rounding defaults, or compiler “optimizations”? A good way to test your algorithm for instabilities is to rerun it with different rounding settings. In C99 you have `fenv` for controlling this, for Python try my [pyfenv](https://github.com/ldo/pyfenv/) module. – Lawrence D'Oliveiro Oct 18 '17 at 03:11
  • @Broman Without determinism, can you even define the semantic of a program? – curiousguy Nov 07 '18 at 05:01

2 Answers2

3

If we use algorithms with double and float arithmetic, how can we guarantee that the results are the same running it in Python and C, in x86 and x64 Linux and Windows computers and ARM microcontrollers?

Generally speaking, you cannot do so except possibly by carefully implementing your own FP operations. If you are using the various languages' standard operators and libraries and the underlying floating-point hardware, then you cannot be ensured of exact reproducability of results across different implementations.

In the first place, there is an issue with the internal representation of floating-point numbers. C does not specify the representation to be used, and even if all else were equal, that means you cannot rely on the same C program running on different implementations (e.g. x86_64 and ARM) to compute identical results.

In practice, most everyone uses IEEE 754 floating-point formats these days, and CPython uses the underlying C implementation's double type to back its floats. Even then, however, IEEE allows for a certain small amount of variation between conforming implementations. Even directives and compilation options that request strict conformance to IEEE specifications cannot fully work around that.

Additionally, you specify that you want to handle both double and float, in both C and Python, but Python doesn't have a native analog of float. Its native floating-point format (probably) corresponds to a C double. Operations performed on different floating-point data types necessarily produce different results, even when the operands are numerically equivalent, and the difference can persist across type conversions, such as converting a double result to float.

There are additional details to consider as well, at the (machine) code-generation level, such as whether or when intermediate results are copied out of FPU registers into main memory (which may involve rounding) and the order in which operations are performed.

We re using an algorithm that uses:

double + double
double + float
double exp(double)
float * float

If you want to minimize the differences in computed values, then start by choosing one floating-point data type, and using it everywhere. For consistency between Python and C implementations, that should be double.

You should also consider disabling any and all optimizations that might change the order in which FP operations are evaluated. That might be all optimizations whatever. If you have options available in your C compiler to enforce strict IEEE conformance, then turn those on.

You should furthermore test the equivalency of the exp() functions on all relevant platforms. You may need to provide your own implementation.


Whatever you do, you should recognize that if your various implementations produce different results despite all being correct in some algorithmic sense, then that's a result in itself. It tells you something about the true precision of the computation, as implemented.

You must never forget that most computer FP operations produce approximate results, so even if you did manage to get all the implementations to produce identical results, that doesn't mean that those results are necessarily more right in an absolute sense than other nearby FP values. If numeric consistency is a requirement, then you ought to quantify that in terms of a specific precision of the results, implement your algorithm in a way that will deliver that precision, and ignore differences at precision higher than the one chosen.

John Bollinger
  • 160,171
  • 8
  • 81
  • 157
  • 2
    To add "In practice, most everyone uses IEEE 754 floating-point formats these days" except for various corners of the spec that were hard to implement or mis-implemented. Many strive for adherence yet do not make it. Example: I was disappointed on the weak/broken implementations of [`fma()`](https://stackoverflow.com/q/42166563/2410359) – chux - Reinstate Monica Oct 18 '17 at 03:02
0

It is hard. Doubles and floats are not formalised in C or C++ standard and its accuracy is up to complier/cpu implementation. For example both float and double are allowed to be same.

From C++17 draft (similar in other papers) basic.fundamental

There are three floating-point types: float, double, and long double. The type double provides at least as much precision as float, and the type long double provides at least as much precision as double. The set of values of the type float is a subset of the set of values of the type double; the set of values of the type double is a subset of the set of values of the type long double. The value representation of floating-point types is implementation-defined. [ Note: This document imposes no requirements on the accuracy of floating-point operations; see also [support.limits]. — end note ]

I don't think there is IEEE 754 mentioned in C or C++ standard.

Python derives this, where floating types are referenced on C implementation, where formalisation is also up to implementaion

There are three distinct numeric types: integers, floating point numbers, and complex numbers. In addition, Booleans are a subtype of integers. Integers have unlimited precision. Floating point numbers are usually implemented using double in C;

Luka Rahne
  • 10,336
  • 3
  • 34
  • 56