-3

I just checked numpy's sine function. Apparently, it produce highly inaccurate results around pi. 

In [26]: import numpy as np

In [27]: np.sin(np.pi)
Out[27]: 1.2246467991473532e-16

The expected result is 0. Why is numpy so inaccurate there?

To some extend, I feel uncertain whether it is acceptable to regard the calculated result as inaccurate: Its absolute error comes within one machine epsilon (for binary64), whereas the relative error is +inf -- reason why I feel somewhat confused. Any idea?

[Edit] I fully understand that floating-point calculation can be inaccurate. But most of the floating-point libraries can manage to deliver results within a small range of error. Here, the relative error is +inf, which seems unacceptable. Just imagine that we want to calculate 

1/(1e-16 + sin(pi))

The results would be disastrously wrong if we use numpy's implementation.

zell
  • 9,830
  • 10
  • 62
  • 115
  • 3
    You did notice the `e-16` there, right? – Fred Larson Sep 07 '16 at 18:35
  • 1
    "So inaccurate" is relative I guess...the difference between 0 and 0.00000000000000012246468 is relatively small considering the computers ability to accurately represent the infinite sequence pi – Lost Sep 07 '16 at 18:35
  • @FredLarson Yes. I did notice the e-16 part. As mentioned above, its absolute error comes within one machine epsilon (for binary64), whereas the relative error is +inf – zell Sep 07 '16 at 18:36
  • 4
    [This](http://stackoverflow.com/questions/20903384/numpy-sinpi-returns-negative-value) is a very similar SO thread - floats can only be so accurate and the [allclose](http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html) method is useful in arrays to check if values are 'close enough'. – Aurora0001 Sep 07 '16 at 18:37
  • Possible duplicate of [Why Are Floating Point Numbers Inaccurate?](http://stackoverflow.com/questions/21895756/why-are-floating-point-numbers-inaccurate) – Fred Larson Sep 07 '16 at 18:37
  • 3
    Also, it's inaccurate to say that relative error is "+Inf" - "[relative error is undefined when the true value is zero](https://en.wikipedia.org/wiki/Approximation_error#Uses_of_relative_error)". – Aurora0001 Sep 07 '16 at 18:51

2 Answers2

7

The main problem here is that np.pi is not exactly π, it's a finite binary floating point number that is close to the true irrational real number π but still off by ~1e-16. np.sin(np.pi) is actually returning a value closer to the true infinite-precision result for sin(np.pi) (i.e. the ideal mathematical sin() function being given the approximated np.pi value) than 0 would be.

Robert Kern
  • 13,118
  • 3
  • 35
  • 32
  • a = np.sin(np.pi) ... np.allclose(a, np.PZERO) # True ... np.allclose(a, np.NZERO) # True can be used to produce a 'zero' if you don't like current representation –  Sep 07 '16 at 19:24
-1

The value is dependent upon the algorithm used to compute it. A typical implementation will use some quickly-converging infinite series, carried out until it converges within one machine epsilon. Many modern chips (starting with the Intel 960, I think) had such functions in the instruction set.

To get 0 returned for this, we would need either a notably more accurate algorithm, one that ran extra-precision arithmetic to guarantee the closest-match result, or something that recognizes special cases: detect a multiple of PI and return the exact value.

Prune
  • 76,765
  • 14
  • 60
  • 81
  • 3
    I'm pretty sure it's just using a look-up table, not an iterative algorithm. Also `np.pi` is just a float so it's only approximately equal to pi. – wim Sep 07 '16 at 18:44
  • I tried looking up what NumPy does; can't find it. Do you have a reference, perhaps? I'd like to know whether it doesn't bother with the native chip functions; the software I work with "trusts" the chip, so a good reference could teach me a few things and utterly torpedo this answer. – Prune Sep 07 '16 at 19:04
  • 2
    It might use the [C library implementation](http://stackoverflow.com/a/2285277/10077). – Fred Larson Sep 07 '16 at 19:06
  • “detect a multiple of PI and return the exact value” is easy. The only multiple of π representable as a floating-point number is 0. There is no reason to invoke algorithm inaccuracy here: the result shown in the question seems about right for sine applied to the nearest binary64 approximation of π. – Pascal Cuoq Sep 07 '16 at 21:59
  • 1
    To get 0 returned for this, we'd need a much *less* accurate algorithm. In this case, `np.sin` is returning the best possible answer that it could possibly give for the given input of `np.pi` (which is exactly equal to `3.141592653589793115997963468544185161590576171875`). It's returning the closest representable float to the true mathematical answer, with an absolute error of around `3e-33`, or about `0.12` ulps. In contrast, an answer of `0` would have an absolute error of around `1e-16`, and would be billions of ulps out. – Mark Dickinson Sep 08 '16 at 12:07
  • @PascalCuoq This is a user-surprise issue: the user put in the package-defined representation of π, and has reasonable expectation of getting an answer of 0. As a QA engineer, I'd file a bug; as an implementation engineer, I'd close it as "will not fix". :-) – Prune Sep 08 '16 at 22:14