0

I'm copying a script from matlab into a c++ function. However, I constantly get different result for the exp function. For example, following snippet:

std::complex<double> final_b = std::exp(std::complex<double>(0, 1 * pi));

should be equivalent to the MATLAB code

final_b = exp(1i * pi);

But it isn't. For MATLAB, I receive -1 + 0i (which is correct) and for c++, I get -1 + -2.068231e-013*i.

Now I thought at the beginning this is just a rounding error of sorts, but for the actual script I'm using, which has bigger complex exponentials, I get completely different numbers. What is the cause of this? How do I fix this?

Edit: I've manually tried calculating the exponential with eulers formula

exp(x+iy) = exp(x) * (cos(y) + i*sin(y))

and get the same wonky results in c++

spurra
  • 1,007
  • 2
  • 13
  • 38
  • 2
    Are you familiar with the term floating point precision? – 101010 Apr 30 '14 at 15:54
  • 1
    @40two judging by the content of the question, I think the answer is yes. – Tom Fenech Apr 30 '14 at 15:55
  • 2
    just an idea, did you check the precision of pi? – Alessandro Teruzzi Apr 30 '14 at 16:03
  • 3
    It's not guaranteed that MatLab and C++ are computing the `exp` in the same manner. Thus, a small round-off error is tolerable. Have in mind that in the background calculating trigonometric functions like `cos` and `sin` involves iterative procedures like `Newton-Raphson` etc. Please state the amound of error you are getting in-order for us to have a more clear picture of your problem. P.S state also how you are printing those values. – 101010 Apr 30 '14 at 16:04
  • Turns out you are right 40two. I've been reading the results wrong...whenever MATLAB outputs 0, C++ outputs something very close to 0. Sorry for the confusion. – spurra Apr 30 '14 at 16:23

1 Answers1

1

That is called floating point approximation (or imprecision):

If you include the header cfloat there are some definitions. In particular, DBL_EPSILON, which is the smallest number that 1.0 + DBL_EPSILON != 1.0, which is usually 1e-9 (and -2.068231e-013 is much smaller than that. If you do the following piece of code, you can check if it is zero or not:

// The complete formula is std::abs(a - b), but since b is zero, I am ommiting it
if (std::abs(number.imag()) < DBL_EPSILON) {
    // The number is either zero or very close to zero
}

For example, you can see the working code here: http://ideone.com/2OzNZm

Bruno Ferreira
  • 1,621
  • 12
  • 19
  • 2
    Just a note on a detail: DBL_EPSILON is 2.2204460492503131e-016 in my float.h. Idem in Matlab ... maybe you are confusing with FLT_EPSILON, which is about 1.192092896e-07 usually. Beside that tiny detail, it is a good explanation. – Bentoy13 Apr 30 '14 at 16:21
  • Would just simply checking if the number is smaller than DBL_EPSILON and then rounding to zero suffice? – spurra Apr 30 '14 at 16:27
  • @Bentoy13 correct, that value depends on libc implementation and might (I am not sure) depend on compiler flags. – Bruno Ferreira Apr 30 '14 at 17:07
  • @BananaCode In fact, that is the correct way to do. Let me fix the code. – Bruno Ferreira Apr 30 '14 at 17:08
  • Uhm, DBL_EPSILON is the minimum value that summed to `1.0` produces a different value (i.e. it's a measure of the granularity of `double`), how is it related to the precision of the result of trigonometric/exponential functions? – Matteo Italia Apr 30 '14 at 17:45
  • If you check [WolframAlpha](https://www.wolframalpha.com/input/?i=e+%5E+%28i+*+pi%29&lk=4&num=1), the result should be `-1 + 0i` or simply `-1`. But the result is something that is not zero, but really close to it, when doing that calculation using the standard C++ api, that is why it is related to the imprecision of floating point numbers. – Bruno Ferreira Apr 30 '14 at 17:59