Search code examples
floating-pointfortranprecisiongfortran

Precision of Sin vs Cos functions in Fortran

Consider the code,

PROGRAM TRIG_TEST
IMPLICIT NONE

DOUBLE PRECISION, PARAMETER :: PI=4.D0*DATAN(1.0D)

print *, sin(PI/2.0), cos(PI/2.0)

END PROGRAM TRIG_TEST

Compiling with gfortran outputs,

1.0000000000000000 6.1232339957367660E-017

I am aware of the usual floating point issues but is there a reason why the sin function is identically 1, but the cos function is not identically zero?

Solution

  • The following assumes double is the IEEE 754 basic 64-bit binary format. Common implementations of the trigonometric routines are less accurate than the format supports. However, for this answer, let’s assume they return the most accurate results possible.

    π cannot be exactly represented in double. The closest possible value is 884279719003555 / 281474976710656 or 3.141592653589793115997963468544185161590576171875. Let’s call this p.

    The sine of p/2 is about 1 − 1.8747•10−33. The two values representable in double on either side of that are 1 and 0.99999999999999988897769753748434595763683319091796875, which is about 1 − 1.11•10−16. The closer of those is 1, so the closest representable value to sine of p/2 is exactly 1.

    The cosine of p/2 is about 6.123233995736765886•10−17. The closest value representable in double to that is 6.12323399573676603586882014729198302312846062338790031898128063403419218957424163818359375•10−17.

    So the results you observed are the closest possible results to the true mathematical values.