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?
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.
ieee_selected_real_kind() to select a kind parameter for given desired characteristics.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.
episilon(1d0)for your compiler?sinandcosare accurate to within7E-017. That gives you the answer forcos.1-7E-017is equal to1.[You say you understand the "usual floating point issues". Do you need this aspect explaining, or do you call this a slight oversight?]