CODEWITHSUNDEEP
floating-pointfortranprecisiongfortran
user1887919
user1887919

Reputation: 923

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?

Upvotes: 2

Views: 2542

Answers (2)

Eric Postpischil
Eric Postpischil

Reputation: 222846

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.

Upvotes: 7

francescalus
francescalus

Reputation: 32366

Let's look at why the result from sin gives 1. In your code

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

This is a floating point approximation to the value of pi. It's probably pretty close to the real value, but it isn't it. Also, the sin function has errors in the evaluation of sin. We hope those are small errors.

We expect that cos(pi/2) has value zero. Your floating point calculation cos(PI/2) has an error of about 6.1232339957367660E-017 from the mathematical answer. Let's assume that your sin calculation has a similar magnitude of error.

Now look at the value of epsilon(0d). This is the smallest number for which 1d0+epsilon(0d0) is not equal to 1d0. The assumed error is much smaller than this number, in the model (for which you report "~ 2.2e-16").

So, 1 is the closest representable number to the real value of the floating point calculation.

Consider the program

  use, intrinsic :: iso_fortran_env, only : real128, real64
  implicit none

  real(real64), parameter :: PI=4.D0*ATAN(1._real64)
  real(real128), parameter :: PI_approx = PI   ! Not 4*ATAN(1._real128)

  print *, SIN(PI/2), COS(PI/2)
  print *, SIN(PI_approx/2), COS(PI_approx/2)

end

This calculates (possibly) the sin of your PI/2 but at a higher precision (using the same approximation for pi). My compiler reports a value different from 1 in the second case, but the difference being much less than epsilon(0._real64).

As a style point, in general it is best to avoid datan and use the generic atan. double precision may also be replaced by appropriate kind parameters. These are shown in my program above.

Upvotes: 1

Related Questions