6

I'm computing sin of pi radians in emacs using this function:

#+begin_src emacs-lisp :eval never-export :exports both
(sin pi)
#+end_src

#+RESULTS:
: 1.2246467991473532e-16

The manual docs mentions that the argument is passed in radians. Just to be sure, I confirmed that sin pi radians is infact zero: https://en.wikipedia.org/wiki/Sine#Special_values

But it seems to be working fine for cosine functions:

#+begin_src emacs-lisp :eval never-export :exports both
(cos pi)
#+end_src

#+RESULTS:
: -1.0

I guess I'm missing something here. Can somebody point what I'm doing wrong ?

12
  • 7
    "Floating point numbers are like piles of sand; every time you move one you lose a little sand and pick up a little dirt." Kernighan & Plauger, The Elements of Programming Style, p.117, attribute it to a "wise old programmer". – NickD Jan 21 at 16:14
  • 3
  • 2
    While these references about floating-point may come in handy, they don't actually answer the specific questions about why (sin pi) and (cos pi) gave the values they did (both are very good approximations to sin(pi) and cos(pi), in fact). If you just throw up your hands and say FP is noisy and stop there (like floating-point-gui.de) you'll fall into traps like continuing to evaluate sin at approximate points near πœ‹, where the function itself will magnify errors irrespective of how you implement it, with or without FP. – Sūrya Siddhānta Jan 22 at 15:21
  • 1
    This isn't just an accident of high-quality math library implementation, though: it's also because, unlike the sine function, the cosine function is well-conditioned near πœ‹: even if your inputs are badly approximated near πœ‹, so you ask for cos(π‘₯β‹…(1 + 𝛿)) when you want cos(π‘₯), what you will get is cos(π‘₯)β‹…(1 + πœ€) where πœ€ is a relative error much smaller than the relative error 𝛿 in the input. (The condition number of cos approaches zero at πœ‹.) – Sūrya Siddhānta Jan 23 at 16:22
  • 1
    @BobJarvis-ReinstateMonica Actually, floating-point rounding errors are only a small part of what's going on here. Any system of numbers and arithmetic that can represent πœ‹ only approximately will yield similar results, even if sine and cosine are evaluated exactly with no rounding errors! – Sūrya Siddhānta Jan 23 at 20:17
19

The Emacs Lisp constant pi is not actually the number πœ‹. Rather, the Emacs Lisp constant pi is an approximation to πœ‹, good to about 16 digits -- pi is exactly 3.141592653589793115997963468544185161590576171875, which is the IEEE 754 binary64 floating-point number closest to πœ‹. This is also the floating-point number closest to 3.141592653589793, which is why, when you evaluate pi, Emacs will show the much shorter 3.141592653589793 as an unambiguous abbreviation for it.

So when you pass pi to the sin function, you are not asking it to compute sin(πœ‹), which is zero. Instead, you are asking it to compute sin(3.141592653589793115997963468544185161590576171875), which is slightly above zero because 3.141592653589793115997963468544185161590576171875 is slightly below πœ‹.

The Emacs Lisp function sin may also return some additional error, but the error is very small: the value of Emacs Lisp (sin pi) differs from the true value of sin(3.141592653589793115997963468544185161590576171875) by less than .000000000000002% (that is, 2 × 10⁻¹⁷) of the magnitude of the true value. In particular, the true value of sin(3.141592653589793115997963468544185161590576171875) is irrational and so doesn't have a finite decimal expansion, but it is roughly

1.2246467991473531772260659322749979970830539012998 × 10⁻¹⁢,

whereas Emacs Lisp returns

1.22464679914735320717376402945839660462569212467758006379625612680683843791484832763671875 × 10⁻¹⁢,

which is wrong only after the sixteenth digit. (Again, this is the nearest floating-point number to 1.2246467991473532 × 10⁻¹⁢, which is why Emacs Lisp prints it unambiguously and more concisely as 1.2246467991473532e-16.)

This number returned by (sin pi) is not simply some random error! It turns out that it is actually a good approximation to the error in Emacs Lisp's approximation of πœ‹! This happens because when π‘₯ is near πœ‹, sin(π‘₯) = (πœ‹ βˆ’ π‘₯) + 𝑂((πœ‹ βˆ’ π‘₯)Β³); this is the Taylor series of sine about πœ‹. So if pi = πœ‹ βˆ’ 𝑒 with absolute error 𝑒, then (sin pi) computes an approximation to sin(πœ‹ βˆ’ 𝑒) β‰ˆ πœ‹ βˆ’ (πœ‹ βˆ’ 𝑒) = 𝑒, good to about 16 digits.

If you add the value of the Emacs Lisp constant pi (πœ‹ βˆ’ 𝑒 = 3.141592653589793115997963468544185161590576171875) and the value of (sin pi) (1.22464679914735320717376402945839660462569212467758006379625612680683843791484832763671875 × 10⁻¹⁢), using exact arbitrary-precision arithmetic, what you get is

3.141592653589793238462643383279505878966979117714660462569212467758006379625612680683843791484832763671875,

in which, if you are a Standard Nerd, you will start to recognize all the digits up to the 32nd or so, since pi + sin(pi) β‰ˆ pi + 𝑒 = πœ‹ βˆ’ 𝑒 + 𝑒 = πœ‹. The familiarity ends there (or earlier if you did not spend an inordinate time as a Standard Nerd Child memorizing digits of πœ‹) because pi and (sin pi) each only have about 16 digits of precision, and so together they have only 32 digits of precision.


You should be aware, though, that the mathematical sine function is said to be ill-conditioned near πœ‹ (and near any odd multiple of πœ‹). What this means is that it magnifies errors in inputs near πœ‹. Specifically:

  • Suppose you want to find the sine of some number π‘₯ near πœ‹.
  • Suppose you don't actually know π‘₯ itself, but only an approximation π‘₯β‹…(1 + 𝛿), with relative error 𝛿.
    • For example, perhaps π‘₯ is a physical dimension, and you made a small error in measurement because you don't have an infinitely precise ruler or compass to measure with in the real world.
  • If you then compute sin(π‘₯β‹…(1 + 𝛿)), you will get, at best, an approximation sin(π‘₯)β‹…(1 + πœ€) to sin(π‘₯). The relative error πœ€ of the result may be arbitrarily large even if the relative error 𝛿 of the input is very small, and even if the sine function is evaluated exactly.

The last part means that sine's being ill-conditioned near πœ‹ is a property of the function, not of the system of floating-point arithmetic. Even if you used some kind of infinite-precision transcendental arithmetic system with no rounding error whatsoever (arbitrary-precision rational arithmetic wouldn't cut it because sine is transcendental), the sine function itself could magnify your small input error near πœ‹ into a large output error.

Example: Suppose the true value of π‘₯ is the floating-point number nearest to 3.14159000000001, but the approximation you have owing to measurement error is the floating-point number nearest to 3.14159 -- this is in error by only a fraction 𝛿 β‰ˆ 3 × 10⁻¹⁡ of the true value.

Then sin(π‘₯) β‰ˆ 2.653589783138678 × 10⁻⁢ but your approximation to it will be sin(π‘₯β‹…(1 + 𝛿)) β‰ˆ 2.65358979335273 × 10⁻⁢ in which about half the digits are wrong.

So while you have storage for 16 digits of precision, and your input was good to nearly all of them, only 9 digits of the output are good and your output error in sin(π‘₯)β‹…(1 + πœ€) is πœ€ β‰ˆ 4 × 10⁻⁹, a million times worse than the input error, because of the question you asked (and not because you used floating-point arithmetic to answer it).

The lesson is that if you find yourself trying to answer the question of the sine of an input π‘₯ near πœ‹, you should consider trying to rephrase the question another way -- for example, perhaps instead of measuring π‘₯ directly near πœ‹, you can measure the much smaller 𝑦 = πœ‹ βˆ’ π‘₯, and then use the approximation sin(π‘₯) = 𝑦 + 𝑂(𝑦³).

This phenomenon is, incidentally, related to why the historic trigonometric functions with crazy-sounding names like versine and exsecant exist -- it wasn't to torture students in school, but rather to allow navigational and astronomical questions to be phrased so they admit computational approximations without magnifying errors.


Why does it seem to work with cosine? The true value of cos(3.141592653589793115997963468544185161590576171875) is an irrational number: βˆ’0.9999999999999999999999999999999925012010866907120… Emacs Lisp returns a floating-point number (always a rational number), and the two nearest floating-point numbers are βˆ’1 and βˆ’0.99999999999999988897769753748434595763683319091796875.

Although Emacs Lisp doesn't always guarantee to return the nearest floating-point number to the true answer for transcendental functions like cosine, in this case it happens to do so, giving βˆ’1. (Usually when the nearest floating-point number is that much closer (relative error below 10⁻³²) than the next-nearest (relative error above 10⁻¹⁢), math libraries will return the nearest one.) So Emacs Lisp happens to return βˆ’1 as an approximation to cos(3.141592653589793115997963468544185161590576171875), even though the question wasn't to compute cos(πœ‹).

This isn't just an accident of high-quality math library implementation, though! It's also because, unlike the sine function, the cosine function is well-conditioned near πœ‹: even if your inputs are badly approximated near πœ‹, so you ask for cos(π‘₯β‹…(1 + 𝛿)) when you want cos(π‘₯), what you will get is cos(π‘₯)β‹…(1 + πœ€) where πœ€ is a relative error much smaller than the relative error 𝛿 in the input.

Example: Suppose you have a much worse approximation to πœ‹ given by 3.14159264305308116860260270186699926853179931640625 (abbreviated 3.141592643053081), which is πœ‹β‹…(1 + 𝛿) whose relative error 𝛿 from πœ‹ is more than 10⁻⁹.

If you ask for the cosine of this approximation when you really want cos(πœ‹), it turns out that cos(πœ‹β‹…(1 + 𝛿)) = cos(3.14159264305308116860260270186699926853179931640625) = βˆ’0.99999999999999994448884937843286910504136142018913… = cos(πœ‹)β‹…(1 + πœ€) whose relative error from the result you wanted (βˆ’1) is πœ€ < 10⁻¹⁢.

That is, the relative error of the result is a million times better than the relative error of the input.

So cosine is much more tolerant to bad approximations of inputs near πœ‹ than sine is.

In more formal terms, the condition number of the function cos(π‘₯) is π‘₯β‹…cosβ€²(π‘₯)/cos(π‘₯) = βˆ’π‘₯β‹…tan(π‘₯) which approaches zero at πœ‹ (well-conditioned), whereas the condition number of the function sin(π‘₯) is π‘₯β‹…sinβ€²(π‘₯)/sin(π‘₯) = π‘₯/tan(π‘₯) which approaches infinity at πœ‹ (ill-conditioned).

In general, you should avoid trying to ask about functions at approximations to points where they are ill-conditioned, and phrase your questions in terms of functions at approximations to points where they are well-conditioned -- the condition number is roughly the factor by which the function itself magnifies errors (before any additional rounding errors in the computation take effect).

2
  • I think the sign is wrong in your Taylor approximation, but a very nice answer otherwise – doetoe Jan 22 at 7:02
  • 1
    @doetoe I have corrected an odd number of sign errors. I hope I have now reduced it to zero mod 2. – Sūrya Siddhānta Jan 22 at 15:10
6

I think this is simply the expected answer. As a float, 1.2246467991473532e-16 is not different from zero, given the floating point error Emacs Lisp can handle. One get the same answer in R for instance:

sin(pi)
#> [1] 1.224647e-16

and also in Python, as shown in this thread.

Actually, this is said in the Emacs Lisp manual:

Floating-point computations often involve rounding errors, as the numbers have a fixed amount of precision.

Most programming languages, other than computer algebra systems (e.g., Maxima or Maple) will give such approximate answers.

1
5

You are doing nothing wrong. This function is working with floating point numbers which have good, but limited, accuracy. 1.2246467991473532e-16 is pretty close to 0 (0.0... with 15 zeros before the 1224... kicks in).

Your Answer

By clicking β€œPost Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.