FLOATING-POINT PRECISION IN FORTRAN

Fortran is a numeric beast, which is exactly why its still popular amongst scientific programmers. Languages like C are still somewhat constrained when it comes to numerical precision. In C, a long double gets you…

Fortran 90 and beyond provides a kind parameter which can be used to specify precision of both reals and integers. Fortran uses the function selected_real_kind() to create a parameter of real type for kind. Inputs to this are precision and range. Here is an example:

integer, parameter :: dp = selected_real_kind(15, 307)
real (kind=dp) :: x, y

This produces a parameter, dp, used to define a double precision real with 15 significant digits, and an exponent range or 307. Single and quad precision can also be defined in a similar way:

integer, parameter :: sp = selected_real_kind(6, 37)
integer, parameter :: qp = selected_real_kind(33, 4931)

In Fortran 2008, the constants real32real64, and real128 can be used instead of selected_real_kind(). How does this work? Consider the following piece of Fortran code to calculate the repeating decimal 1.0/7.0. The answer should be 0.142857 142857… Note that the floating point numbers used have to be tagged as a specific kind (e.g. 1.0_dp), otherwise the outcome can be erroneous.

integer, parameter :: dp = selected_real_kind(15, 307)
real (kind=dp) :: d
d = 1.0_dp / 7.0_dp

What happens when we run it  in all three precisions?

   SP = 0.14285714924335479736328125000000000
   DP = 0.14285714285714284921269268124888185
   QP = 0.14285714285714285714285714285714285

The BOLD numbers show the precision of the number, with correct digits. The Underline numbers are those correct beyond the precision. Now compare this to C’s long double which is typically 80-bit extended precision. Here’s the C code:

long double d;
d = 1.0L / 7.0L;

Here’s the result:

0.1428571428571428571409210675491330277964

Leave a Reply