Sin, cos, tan inaccurate

Question

Why sinl give incorrect results when the argument is close to the non-zero edge of pi? Why sinl give incorrect results with a large argument? The following code illustrates this.

Note that the digits used to initialize the variable pi do not match the exact 64-bit double value. The compiler selects the closest value, which is 3.14159265358979323851280895940618620443274267017841339111328125 . The expected sine value can be found using libquadmath, gnu MPFR lib, or an online calculator such as http://www.ttmath.org/online_calculator .

 #include <stdio.h> #include <math.h> int main (int argc, char *argv []) { volatile long double pi = 3.14159265358979323846L; volatile long double big = 9223372035086174241L; volatile long double expected1 = -5.0165576126683320235E-20L; volatile long double expected2 = -4.2053336735954077951E-10L; double result; double ex1 = expected1, ex2 = expected2; result = sinl (pi); printf("expected: %g, \nreturned: %g\n\n", ex1, result); result = sinl (big); printf("expected: %g, \nreturned: %g\n\n", ex2, result); return 0; } 

I am using gcc 4.7.3. Using volatile allows the compiler to replace the sinl() call with the result of hard coding. My computer has an Intel Core i7 processor and is running Windows. I print the results as double instead of long double, because the gcc mingw port I use does not support printing a long double. Here is the output of the program:

 expected: -5.01656e-020, returned: -5.42101e-020 expected: -4.20533e-010, returned: -0.011874 

Answer 1

Inaccuracy can be attributed to the fsin processor command used by the sinl library code. The fsin, fcos, and fptan instructions do not match 1.0 ulp, as Intel claims: http://notabs.org/fpuaccuracy/

Answer 2

To achieve 1 ULP accuracy for multiples of pi, the internal M_PI constant should have approximately 106-bit accuracy (or 128 for long doubles).

At the recovery stage, the ideal implementation would somehow generate the missing 53 or 64 bits of accuracy after subtraction (x - M_PI) , since a naive implementation would calculate this intermediate value as zero. The problem, of course, becomes more and more when the argument is a large integer multiplied by a nonzero one.

66 bit internal accuracy for M_PI is not enough for 1 ULP accuracy. Then again, you can re-read the statements and check whether the accuracy of 1 ULP has been indicated with respect to the result or argument.

Answer 3

The GNU libc documentation (available when you run info libc math errors ) contains 1 known "ul_" error for cosl on x86 and "x86_64 / fpu". It does not document anything for sinl . I can reproduce similar huge errors for cosl around pi / 2 on my x86_64 machine.

Perhaps you should report this as a documentation error for users of man glibc and Linux; I can’t imagine that it is worth implementing the “correct” fix.

If you really want fast and accurate sinl , I'm not too sure where to look. CRlibm does sin (option for double s). MPFR will process it, but it will be many times slower than fsin .