sin() returns a slightly different results for the different platforms

Question

I've noticed during CI that this code

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

int main()
{
    const double d = 3.81219767052986080458;
    fprintf(stderr, "%.20f\n", sin(d));
    return 0;
}

compilled with gcc test.c -lm -o test
returns a slightly different results for my 2 test platforms.
It's -0.62146009873891927544
witn gcc version 13.2.0 libc 2.39-0ubuntu8.3 Kubuntu 24.04
and -0.62146009873891938646
with gcc version 12.2.0 libc 2.36-9+deb12u7+ci202405171200+astra5+b1 Astra Linux 1.8

I wonder what affects it? libm? gcc? may be CPU?
This slight difference may be accumulated during integration in a tight loops.
Is there a simple way to get rid of it without performance penalties? A gcc compiler param, for example?

Answer 1

I wonder what affects it? libm?

Differences in sin for identical inputs are generally caused by different implementations of sin in the standard C library. The implementation of sin may be in a “math library” called libm, or it may be in other files, depending on the platform.

(Another possibility is that the input numeral, 3.81219767052986080458, was parsed differently in the different compiler versions, but that is less likely.)

Ideally, a sin implementation returns the number that equals the mathematical sine of its input operand rounded to the nearest value representable in the floating-point format.¹ This would produce the first result you show or, more precisely, −.62146009873891927544065083566238172352313995361328125. However, the C standard does not require this; it allows implementations to produce results that are not correctly rounded.

For some math library functions, humans do not yet know how to compute the entire function with correct rounding in an execution time with a known bound. We can approximate the functions to any desired precision using extended-precision arithmetic, and we can write iterative algorithms that compute more and more precision until the result is known so accurately that we know which way it should be rounded. However, some results of some functions may be so close to the point where rounded changes (halfway between representable values) that a great deal of precision is needed. We do not yet know how much precision is needed for all functions. So we can write an iterative algorithm that works given enough time, but, for some functions, we do not know how much time is required to cover all cases.

Many users would not tolerate math library routines that might, in some cases, run for long periods of time. They want routines that run in short fixed or limited times. Again, it is currently impossible for humans to implement for all math library functions because we do not know how much precision is needed for the worst cases for all of the math library functions.

This is not the case for sine for the IEEE-754 binary32 and binary64 formats (commonly used for float and double). All the worst cases for it have been completely mapped. So a correctly rounded sine can be implemented and it can be done with reasonable speed, although some people prefer faster implementations that are not correctly rounded.

The CORE-MATH project is working on further development of correctly rounded implementations.

Even if somebody has chosen not to implemented a correctly rounded sin, 3.81219767052986080458 is slightly surprising to find a discrepancy in, because it is so small. With periodic functions, the usual design is to reduce arguments outside a base interval to the base interval, then evaluate a polynomial that approximates the function in the base interval. For sine, an implementation usually uses at least a little extra precision for that reduction to the base interval and for the polynomial evaluation, and a little extra precision should suffice for that case. I do note that the sine of the double nearest 3.81219767052986080458 is just about 0.2% away from the midpoint between the two nearest representable values. So eight or nine extra bits ought to suffice. But maybe one of your sine implementations did not do that.

Is there a simple way to get rid of it without performance penalties?

If the problem is in the sin implementation, you would have to use a different math library. Its performance may differ.

Footnote

¹ For functions in general, the rounding would break ties in favor of the nearest representable value with an even low digit (bit for base-two floating-point) in its significand. This is not relevant for sine since its only representable operand with a rational sine is zero, so none of the operands involve ties.