11

I'm in the process of converting a program from Scilab code to C++. One loop in particular is producing a slightly different result than the original Scilab code (it's a long piece of code so I'm not going to include it in the question but I'll try my best to summarise the issue below).

The problem is, each step of the loop uses calculations from the previous step. Additionally, the difference between calculations only becomes apparent around the 100,000th iteration (out of approximately 300,000).

Note: I'm comparing the output of my C++ program with the outputs of Scilab 5.5.2 using the "format(25);" command. Meaning I'm comparing 25 significant digits. I'd also like to point out I understand how precision cannot be guaranteed after a certain number of bits but read the sections below before commenting. So far, all calculations have been identical up to 25 digits between the two languages.

In attempts to get to the bottom of this issue, so far I've tried:

  1. Examining the data type being used:

I've managed to confirm that Scilab is using IEEE 754 doubles (according to the language documentation). Also, according to Wikipedia, C++ isn't required to use IEEE 754 for doubles, but from what I can tell, everywhere I use a double in C++ it has perfectly match Scilab's results.

  1. Examining the use of transcendental functions:

I've also read from What Every Computer Scientist Should Know About Floating-Point Arithmetic that IEEE does not require transcendental functions to be exactly rounded. With that in mind, I've compared the results of these functions (sin(), cos(), exp()) in both languages and again, the results appear to be the same (up to 25 digits).

  1. The use of other functions and predefined values:

I repeated the above steps for the use of sqrt() and pow(). As well as the value of Pi (I'm using M_PI in C++ and %pi in Scilab). Again, the results were the same.

  1. Lastly, I've rewritten the loop (very carefully) in order to ensure that the code is identical between the two languages.

Note: Interestingly, I noticed that for all the above calculations the results between the two languages match farther than the actual result of the calculations (outside of floating point arithmetic). For example:

Value of sin(x) using Wolfram Alpha = 0.123456789.....

Value of sin(x) using Scilab & C++ = 0.12345yyyyy.....

Where even once the value computed using Scilab or C++ started to differ from the actual result (from Wolfram). Each language's result still matched each other. This leads me to believe that most of the values are being calculated (between the two languages) in the same way. Even though they're not required to by IEEE 754.

My original thinking was one of the first three points above are implemented differently between the two languages. But from what I can tell everything seems to produce identical results.

Is it possible that even though all the inputs to these loops are identical, the results can be different? Possibly because a very small error (past what I can see with 25 digits) is occurring that accumulates over time? If so, how can I go about fixing this issue?

curiousguy
  • 8,038
  • 2
  • 40
  • 58
Paul Warnick
  • 903
  • 2
  • 13
  • 26
  • Can you simply persist the floating-point values without conversion or formatting? It looks like the HDF5 management library in Scilab should be able to write attributes in "H5T_IEEE_F64LE" format (assuming your architecture is little-endian). Then, you just need to persist the same HDF5 format from C++ and you can compare them exactly. – Useless Jun 20 '16 at 17:15
  • 3
    "So far, all calculations have been identical up to 25 digits". Considering that IEEE double has 15-17 decimal digits of precision, that's a better than average result. – n. m. could be an AI Jun 20 '16 at 17:20
  • Neither C++ nor C# even have the same results as *themselves* if you compile differently (eg x87 vs SSE, for C# that comes down to x86 vs x64). – harold Jun 20 '16 at 17:25
  • @Useless I'm not exactly sure what you mean by persisting the floats. I'll look into HDF5 but unfortunately, I'm not able to edit the Scilab code in any way (if it affects the output). If that's what you're suggesting. – Paul Warnick Jun 20 '16 at 17:46
  • @harold To what extent? Meaning if you were to execute 1 + 1 in C++ it would result in the same output on most compilers. Is this something that would effect most calculations only slightly? And since I'm finding that 25 digits between languages are matching. Could this error occur sometime after those digits and still affect the program? – Paul Warnick Jun 20 '16 at 17:50
  • @harold: You can control the behavior through compiler flags. If you compile with `/fp:precise`, you will get identical results, irrespective of the floating point unit in use. The effect is that a value stored in an 80bits FPU register is copied to memory (thereby rounding it to 64bits), and then read back before the next operation is performed. Performance plummets, but it ensures identical results across floating point units (FPU vs. SSE/AVX). – IInspectable Jun 20 '16 at 18:01
  • Why do you even need that many digits of precision? – Jesper Juhl Jun 20 '16 at 18:32
  • @JesperJuhl I don't exactly need such a high level of precision. I just need the results from the two languages to be identical. Thus, I tested the inputs to see how similar they were in order to determine if that was the problem. They just happened to be similar for 25 digits (the most I could display in Scilab). – Paul Warnick Jun 20 '16 at 18:35
  • If they are identical up to 25 digits and you don't even need that much, why don't you just use - say - the first 10 digits? And if they are equal than all is good and you just ignore the rest and move on...? – Jesper Juhl Jun 20 '16 at 18:37
  • 1
    I echo the suggestion to experiment with compiler flags. Some compilers default to aggressive optimization settings that allow the re-arrangement floating-point expressions (e.g. in ways that are equivalent mathematically but not in limited-precision floating-point arithmetic), may contract multiplies and adds into FMAs, or use SIMD reductions. Check the compiler documentation for the switch(es) that promise the strictest compliance with IEEE-754 and other relevant standards. For example on my Intel compiler that is `/fp:strict`. You may lose a lot of performance when using such a flag. – njuffa Jun 20 '16 at 22:18
  • Persist = save to disk. And, writing intermediate values to disk shouldn't change the result unless it affects optimization. – Useless Jun 21 '16 at 09:04
  • "I've compared the results of `sin(), cos(), exp()` in both languages and the results appear to be the same" - does that apply to _all_ inputs or just _some_ inputs? For instance, a common implementation technique is to use `sin(x)==x` for small x, but exactly how small should `x` be? – MSalters Jun 21 '16 at 10:28

4 Answers4

7

No, the format of the numbering system does not guarantee equivalent answers from functions in different languages.

Functions, such as sin(x), can be implemented in different ways, using the same language (as well as different languages). The sin(x) function is an excellent example. Many implementations will use a look-up table or look-up table with interpolation. This has speed advantages. However, some implementations may use a Taylor Series to evaluate the function. Some implementations may use polynomials to come up with a close approximation.

Having the same numeric format is one hurdle to solve between languages. Function implementation is another.

Remember, you need to consider the platform as well. A program that uses an 80-bit floating point processor will have different results than a program that uses a 64-bit floating point software implementation.

Thomas Matthews
  • 56,849
  • 17
  • 98
  • 154
  • As I suggested in one of your earlier posts, you should write your own functions to provide more consistency in results especially between two different languages. For example, write a `sin` function in FORTRAN and implement it the same way in C++. This gives you better accuracy when comparing results. Otherwise you are comparing apples to oranges. – Thomas Matthews Jun 20 '16 at 17:22
  • Actually, that was what I planned to do next. The only problem is the results need to be identical. Therefore, I need to use whatever functions Scilab is currently using to calculate the result and I'm having trouble locating them in the source code. – Paul Warnick Jun 20 '16 at 17:56
  • Also in regards to your answer, when you mention having the same numeric format. Are you aware if I already do or not? I know Scilab is following IEEE for double precision, and I'm assuming C++ is as well. Is it possible that C++ is just similar enough to appear the same, but really it's a different implementation? – Paul Warnick Jun 20 '16 at 17:59
  • Can you write your own transcendental functions in Scilab? This is what I'm suggesting. Otherwise, try to get the source code for the Scilab compiler or math functions. – Thomas Matthews Jun 20 '16 at 17:59
  • Lastly, I'm running both of these programs on the same platform (i.e. my personal computer). Would that avoid the problem of the last point in your answer? – Paul Warnick Jun 20 '16 at 18:00
  • 1
    No, the code is what matters. I can have a floating point processor on my computer and have a program that doesn't use it. You have to verify that all dependencies are the same between the two functions. Again, this is simplified when you write your own functions in both languages. – Thomas Matthews Jun 20 '16 at 18:04
  • 3
    @PaulWarnick: Whether or not your C++ compiler implements IEEE 754 floating point representation is documented in the compiler's manual. There is no requirement for an implementation of the C++ language to use any given floating point representation. – IInspectable Jun 20 '16 at 18:50
6

Some architectures provide the capability of using extended precision floating point registers (e.g. 80 bits internally, versus 64-bit values in RAM). So, it's possible to get slightly different results for the same calculation, depending on how the computations are structured, and the optimization level used to compile the code.

Jim Lewis
  • 43,505
  • 7
  • 82
  • 96
  • So even if I was to use the exact same functions (meaning, copy whatever Scilab is using for sin(), cos(), etc. into C++) there is a chance that the results could still be different? – Paul Warnick Jun 20 '16 at 18:09
  • If one language/implementation is using x87 and the other is using, say, SSE or running on SPARC or Power or something else, then yes. You are not guaranteed identical results to the last bit. – Jesper Juhl Jun 20 '16 at 18:25
  • @JesperJuhl How would I go about determine if the two languages are using the same architecture? – Paul Warnick Jun 20 '16 at 19:54
  • 1
    If they're running on the same CPU, they're using the same architecture. As to whether they're using the same _part_ of that architecture - you'd need to know what scilab uses, and then figure out how to make your C++ compiler use the same instructions (eg, disabling AVX, disabling fast-math, etc.) – Useless Jun 21 '16 at 09:06
3

Yes, it's possible to have a different results. It's possible even if you are using exactly the same source code in the same programming language for the same platform. Sometimes it's enough to have a different compiler switch; for example -ffastmath would lead the compiler to optimize your code for speed rather than accuracy, and, if your computational problem is not well-conditioned to begin with, the result may be significantly different.

For example, suppose you have this code:

 x_8th = x*x*x*x*x*x*x*x;

One way to compute this is to perform 7 multiplications. This would be the default behavior for most compilers. However, you may want to speed this up by specifying compiler option -ffastmath and the resulting code would have only 3 multiplications:

temp1 = x*x; temp2 = temp1*temp1; x_8th = temp2*temp2;

The result would be slightly different because finite precision arithmetic is not associative, but sufficiently close for most applications and much faster. However, if your computation is not well-conditioned that small error can quickly get amplified into a large one.

Michael
  • 5,775
  • 2
  • 34
  • 53
1

Note that it is possible that the Scilab and C++ are not using the exact same instruction sequence, or that one uses FPU and the other uses SSE, so there may not be a way to get them to be exactly the same.

As commented by IInspectable, if your compiler has _control87() or something similar, you can use it to change the precision and/or rounding settings. You could try combinations of this to see if it has any effect, but again, even you manage to get the settings identical for Scilab and C++, differences in the actual instruction sequences may be the issue.

http://msdn.microsoft.com/en-us/library/e9b52ceh.aspx

If SSE is used, I"m not sure what can be adjusted as I don't think SSE has an 80 bit precision mode.

In the case of using FPU in 32 bit mode, and if your compiler doesn't have something like _control87, you could use assembly code. If inline assembly is not allowed, you would need to call an assembly function. This example is from an old test program:

static short fcw; /* 16 bit floating point control word */
/* ... */
/* set precision control to extended precision */
__asm{
    fnstcw fcw
    or fcw,0300h
    fldcw fcw
}
rcgldr
  • 27,407
  • 3
  • 36
  • 61
  • What exactly would this achieve? Would it increase the precision of my C++ calculations? If so I'm not entirely sure that would solve the issue as I'm trying to exactly match the results of the Scilab code (which apparently already has the same precision as my C++ code). – Paul Warnick Jun 20 '16 at 18:03
  • 64bit compilers on x64 usually do not use the FPU at all, but implement all floating point math in the SSE unit. That means 64-bit floating point numbers, without the option to switch to extended precision. In case you do want to set the control word, you don't need assembly either, when compiling with Visual Studio. It provides the [_control87](https://msdn.microsoft.com/en-us/library/e9b52ceh.aspx) function for that. – IInspectable Jun 20 '16 at 18:55
  • Alright, this makes sense (although it's all new to me so I'll have to do some googling). If it makes any difference I'm compiling my C++ via gcc and I'm just running Scilab by hitting execute in the Scilab console (5.5.2). Any useful suggestions that I might want to know for these methods of execution? – Paul Warnick Jun 20 '16 at 21:08
  • @PaulWarnick - I updated my answer to clarify the issue and to reflect the comment by IInspectable. – rcgldr Jun 20 '16 at 21:29