sin, cos, tan and rounding error

Question

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.

Answer 1

The answer is only 0 for sin(pi) - did you include all the digits of Pi ?

-Has anyone else noticed a distinct lack of, irony/sense of humour around here?

Answer 2

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.

Answer 3

Sine of π is 0.0.
Sine of M_PI is about 1.224647e-16.

M_PI is not π.

program gives ... 1.224647e-16 when the correct answer is of course 0.

Code gave a correct answer to 7 significant places.

---

The following does not print the sine of π. It prints the sine of a number close to π. See below pic.

π                            // 3.141592653589793 2384626433832795...
printf("%.21\n", M_PI);      // 3.141592653589793 115998
printf("%.21f\n", sin(M_PI));// 0.000000000000000 122465

Note: With the math function sine(x), the slope of the curve is -1.0 at x = π. The difference of π and M_PI is about the sin(M_PI) - as expected.

---

am running into problems with rounding errors

The rounding problem occurs when using M_PI to represent π. M_PI is the double constant closest to π, yet since π is irrational and all finite double are rational, they must differ - even by a small amount. So not a direct rounding issue with sin(), cos(), tan(). sin(M_PI) simple exposed the issue started with using M_PI - an inexact π.

---

This problem, with different non-zero results of sin(M_PI), occurs if code used a different FP type like float, long double or double with something other than 53 binary bits of precision. This is not a precision issue so much as a irrational/rational one.

Answer 4

@Josh Kelley - ok serious answer.
In general you should never compare the results of any operation involving floats or doubles with each other.

The only exceptions is assignment.
float a=10.0;
float b=10.0;
then a==b

Otherwise you always have to write some function like bool IsClose(float a,float b, float error) to allow you to check if two numbers are within 'error' of each other.
Remember to also check signs/use fabs - you could have -1.224647e-16

Answer 5

There are two sources of error. The sin() function and the approximated value of M_PI. Even if the sin() function were 'perfect', it would not return zero unless the value of M_PI were also perfect - which it is not.

Answer 6

I rather think that will be system-dependent. I don't think the Standard has anything to say on how accurate the transcendental functions will be. Unfortunately, I don't remember seeing any discussion of function precision, so you'll probably have to figure it out yourself.

Answer 7

Unless your program requires significant digits out to the 16th decimal place or more, you probably can do the rounding manually. From my experience programming games we always rounded our decimals to a tolerable significant digit. For example:

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

#define HALF 0.5
#define GREATER_EQUAL_HALF(X) (X) >= HALF

double const M_PI = 2 * acos(0.0);

double round(double val, unsigned  places = 1) 
{
    val = val * pow(10.0f, (float)places);
    long longval = (long)val;
    if ( GREATER_EQUAL_HALF(val - longval) ) {
       return ceil(val) / pow(10.0f, (float)places);
    } else {
      return floor(val) / pow(10.0f, (float)places);
    }
}

int main() 
{
    printf("\nValue %lf", round(sin(M_PI), 10));
    return 0;
}

Answer 8

I get the exact same result on my system - I'd say it is close enough

I would solve the problem by changing the format string to "%f\n" :)

However, this gives you a "better" result, or at least on my system it does give -3.661369e-245

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

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

Answer 9

Maybe too low accuracy of implementation

M_PI = 3.14159265358979323846 (M_PI is not π)

http://fresh2refresh.com/c/c-function/c-math-h-library-functions/

It is an inaccuracy in implementation, see Stephen C. Steel's comment under Andy Ross` answer above and chux's answer.