Precision in Fortran 90

1. Introduction

Recall that digital computers store numbers in fixed and finite size memory locations. For example, the default on a Sun workstation such as k12.colostate.edu for integer and real variables is 4 Bytes, or 32 bits. Each number is stored as a specific number of Bytes (each Byte = 8 bits). The greater the number of Bytes, the more digits can be stored, and the greater the size of the number. Fortran 90 has intrinsic functions which allow us to interrogate the storage of a particular number in a memory location.

Finite memory size results in an imperfect ability to represent numbers, as only a limited amount of information can be stored in a limited memory space. Recall that integers are stored perfectly, while only the most significant digits of the mantissa and the integer exponent are stored for reals. For integers, either the number fits in the memory location and is stored perfectly, or does not fit and is stored erroneously. We rarely have problems with storing integers, as long as we use 4 Bytes (32 bits) for storage. Two complicating factors exist. First, negative numbers are stored internally in two's complement form - we shall not pursue this here, as it is beyond the scope of our presentation. Secondly, and most worrisome, is that sometimes when integers exceed the storage space available for them (say as a result of an addition or multiplication), there is no warning or error message issued. Thus, to protect ourselves against this, we should use integer variables only sparingly, and only when we can be reasonably certain that they will never be large enough to exceed storage. Typically, integers are used only as counters, or to indicate memory locations, i.e. integers are used for bookkeeping. Thus, we shall turn our attention to real numbers, which are used in most if not all calculations.

The number of digits able to be stored in a real mantissa is referred to as precision. High precision refers to the ability to store more information than low precision. In Fortran 90, this precision is referred to using the symbol p, which indicates the minimum number of decimal digits which can be stored. The exponent is also stored for a real number. In Fortran 90, the decimal size of a real number is referred to with the symbol, r, for the range of the number.

A single precision real number on most computers has a precision, p, of between 6 and 7 decimal digits, and a range, r, of 37. Thus, numbers as small as 10**(-37), and as large as 10**37 can be stored, with between 6 and 7 significant decimal digits.

Example. Consider the real number 123456789. To be stored, the number is first shifted to normalized form, making the mantissa less than one, to: 0.123456789 x 10**9. Thus, the exponent stored internally for this number is 9, which fits well within the available range of 37. Since only six decimal digits can be stored, the mantissa is rounded and truncated to six digits - this is termed a "loss of precision." The number stored is, then, 0.123457 x 10**9. Actually, only the mantissa (123457) and exponent (9) are stored internally. Contemplate on your own when the loss of precision might be important (see the earlier notes).

We have already discussed the consequences of insufficient precision when adding or subtracting, and multiplying or dividing. In this section, we shall present the implementation details of precision in Fortran 90. A specific emphasis in Fortran 90 is to permit the user to specify the precision of variables so as to make programs more portable.

2. Implementation

Different kinds of numbers are stored differently internally. Fortran 90 has the "kind" attribute to specify how a number is stored internally. The data storage of a number is specified in Fortran 90 using the kind attribute:
real, kind = 2 :: a, b, c
real, kind = 4 :: e, f, g
integer, kind = 2 :: i, j, k
integer, kind = 3 :: l, m, n

In the above example, the real variables e, f and g have more precision than the real variables a, b and c. Likewise, the integer variables l, m and n, can be larger and have more digits for storage than the integer variables i, j and k. However, the exact storage of any variable of a specific kind may differ from machine to machine, and the default storage (used if the optional kind = attribute is omitted from the type specification statement) of variables may also vary from machine to machine. To deal with this uncertainty, Fortran 90 allows the user to interrogate the precision of variables in a variety of ways. In fact, this is one of the first things we should do when programming on any new machine. To determine the default kind, the kind intrinsic function may be used:

real :: a
integer :: i, kind_a, kind_i
kind_a = kind(a)
kind_i = kind(i)
print *,'default kind for real is', kind_a
print *,'default kind for int is', kind_i

A number of intrinsic functions exist in Fortran 90 for interrogating the sizes of numbers. See Appendix A.4 on pages 703 to 705 of Ellis, Philips and Lahey for a complete listing. For integers, the bit_size(integer_number) intrinsic function specifies the number of bits used for storage. For reals, two intrinsic functions exist: precision(real_number) returns p, the number of decimal digits of precision, while the range(real_number) returns r, the decimal range of the exponent. Some examples, with additional intrinsics, follow.

program get_sizes

!
! set storage
!
   real (kind = 1) :: a
   real (kind = 2) :: b
   real (kind = 4) :: c
   integer (kind = 1) :: i
   integer (kind = 2) :: j
   integer (kind = 3) :: k
!
! interrogate real variables
!
   print *,'precision of real(1) =', &
    precision(a)
   print *,'precision of real(2) =', &
    precision(b)
   print *,'precision of real(4) =', &
    precision(c)
!
   print *,'range of real(1) =', range(a)
   print *,'range of real(2) =', range(b)
   print *,'range of real(4) =', range(c)
!
   print *,'maximum exponent of real(1) ='&
    ,maxexponent(a)
   print *,'maximum exponent of real(2) ='&
    ,maxexponent(b)
   print *,'maximum exponent of real(4) ='&
    ,maxexponent(c)

!
   print *,'minimum exponent of real(1) ='&
    , minexponent(a)
   print *,'minimum exponent of real(2) ='&
    , minexponent(b)
   print *,'minimum exponent of real(4) ='&
    , minexponent(c)
!
! now for the integers
!
   print *,'bits in integer(1) ='&
    ,bit_size(i)
   print *,'bits in integer(2) ='&
    ,bit_size(j)
   print *,'bits in integer(3) ='&
    ,bit_size(k)
!
stop
!
end program get_sizes

Fortran 90 has the capability to specify the minimum size of the mantissa and exponent. Formally:

p = the minimum number of decimal (base 10) digits in the mantissa, and
r = the minimum decimal (base 10) size of the exponent,

Then the minimum precision can be specified for real numbers as follows:

real (kind=selected_real_kind(p=6)) :: a
real (kind=selected_real_kind(p=12)) :: b

Note that Fortran guarantees only a minimum precision. The kind returned meets or exceeds the minimum we specify as the argument of the selected_real_kind intrinsic. On most machines, a will be a 4-Byte real, and b will be an 8-Byte real. It is also possible to specify the minimum range, although this is not nearly so critical as specifying the precision:

real (kind=selected_real_kind(r=30)) :: c
real (kind=selected_real_kind(r=100)) :: d

Again, on most machines, c will be a 4_byte real, and d will be an 8_byte real. Alternatively, we can specify both the precision and range, and the kind returned is the minimum which meets both criteria:

real(kind=selected_real_kind(p=6,r=30)) & :: e

Here, e will typically be returned as a 4-Byte real. Or, we can determine kinds directly, using the intrinsic function selected_real_kind() (see appendix A.4)

integer :: real_kind_p6, real_kind_p12
!
   real_kind_p6 = selected_real_kind(p=6)
   real_kind_p12 = selected_real_kind(p=12)

We recommend specifying the precision only, as this is most critical. In the rare cases when an overflow results during computation, this will cause the machine to "crash," and an error message will be issued to the screen. Then, the user can go back and edit the code to obtain additional range where needed.

3. Exercises

Exercise 1. Determine all permissible kinds for: a) integers, and b) reals on your computer. Try all values of kind= from 1 to 10.

Exercise 2. For all permissible kinds of integers, determine the number of binary bits used to represent the numbers.

Exercise 3. For all permissible kinds of reals, determine the precision and ranges used to represent the numbers. (Note: You can do exercises 1 through 3 in a single program.)

Exercise 4. Repeat the Pi exercise using at least 12 digits of precision. Where in n (the number of sides in the polygon) does precision begin to cause inaccuracy?

Exercise 5. Write a program to compute Pi vs. n for both single precision (r = 6), and double precision (r = 8). Compute the error in Pi vs. n for both situations. Construct a Spyglass data file which displays the errors in two columns vs. n for n = 10, 100, 1,000, and 10,000. Plot the error using Spyglass. Remember to write as the first line a dummy heading with two numbers, and a dummy (first) column, along with both answers as columns 2 and 3 for all remaining lines.

Exercise 6. Duplicate the calculations in Figure 10.2 on page 318 of Ellis, Philips and Lahey, except use p = 12 instead of p = 14 as they have done. Discuss how precision affects the computation.