Search code examples
cmathfloating-pointroundingtrigonometry

sin, cos, tan and rounding error


I'm doing some trigonometry calculations in C/C++ and am running into problems with rounding errors. For example, on my Linux system:

#include <stdio.h>
#include <math.h>

int main(int argc, char *argv[]) {
    printf("%e\n", sin(M_PI));
    return 0;
}

This program gives the following output:

1.224647e-16

when the correct answer is of course 0.

How much rounding error can I expect when using trig functions? How can I best handle that error? I'm familiar with the Units in Last Place technique for comparing floating point numbers, from Bruce Dawson's Comparing Floating Point Numbers, but that doesn't seem to work here, since 0 and 1.22e-16 are quite a few ULPs apart.


Solution

  • An IEEE double stores 52 bits of mantissa, with the "implicit leading one" forming a 53 bit number. An error in the bottom bit of a result therefore makes up about 1/2^53 of the scale of the numbers. Your output is of the same order as 1.0, so that comes out to just about exactly one part in 10^16 (because 53*log(2)/log(10) == 15.9).

    So yes. This is about the limit of the precision you can expect. I'm not sure what the ULP technique you're using is, but I suspect you're applying it wrong.