Double precision with Fortran sin and cos

Question

I'm investigating floating point accuracy in Fortran 95 and I have the following code:

  PROGRAM MAIN
  IMPLICIT NONE
  integer, parameter :: dp = SELECTED_REAL_KIND(15,307)
  REAL(dp), PARAMETER :: DEG2RAD = ATAN(1.d0)/45.d0

! ---- Calculate sine and cosine
  WRITE (*,*) COS(60.d0*DEG2RAD) - DBLE(0.5)
  WRITE (*,*) SIN(60.d0*DEG2RAD) - SQRT(3.d0)/2.d0

  END PROGRAM MAIN

I'm expecting both of the printed lines to be zero, but only the latter one is. Could anyone explain why that is? Of course I understand that the errors are smaller than double precision, but still the first one prints 1e-16, while the second one shows only zeros. Is there a way to make them both "exactly" zero?

1e-16 is within the possibilities of double precision. In single precision the error is bigger. You cannot ask for more. You should be using `0.5d0` instead of `dble(0.5)` but it want make it zero. — Vladimir F Героям слава, Jun 16 '20 at 09:14
Also note that the math functions are coming from a system library, not from gfortran. They will behave differently in different operating systems. In Linux they come from https://www.gnu.org/software/libc/manual/html_node/Trig-Functions.html Other compiler may use different libraries with different series used and with rounding properties. — Vladimir F Героям слава, Jun 16 '20 at 10:11
More in https://stackoverflow.com/questions/2284860/how-does-c-compute-sin-and-other-math-functions — Vladimir F Героям слава, Jun 16 '20 at 10:17

Double precision with Fortran sin and cos

0 Answers0