More on numbers and variables

Yutaka Masuda

February 2020

Back to index.html.

More on numbers and variables

Basic usage

Types

Fortran is a language that has data types; the program must declare what kind of variables and data stay in the program. In R (and Python), there are the types, but it is hidden; the user should clarify the types in some cases. In many scripting languages, like Perl and awk, the user may not be aware there is a type. So, why does Fortran have the types?

The merit of data types is high performance. With the data types, the compiler can know which data is conformable in computations and comparisons. If two data types are not conformable, the compilation may fail, or the compiler converts a data type to the other. Because the type conversion is a significant penalty in computing time and memory usage, the program slows down. The programmers can write the code to avoid the data conversion because they know the types.

Literals

There are two types of numeric values in Fortran: integer and real. Additionally, you can use two types of data: logical and characters.

Most intrinsic functions accept a specific type. For example, sin() accepts a real number, not integer; sin(1.0) is okay but sin(1) fails.

Variables

A variable can hold a typed value. Here is a general form of the declaration of variables, which should be placed at the top of a program unit. Also, to turn off the implicit variable typing, you should use implicit none once before the declarations.

implicit none
variable_type,attribute :: variable1, variable2, ...

The following type-statements are available.

Here are a couple of attributes that we have introduced.

Precision

A numerical type has options in precision. You can change the precision both for literals and for variables. An integer code defines the precision; the code is available in 2 ways.

To use inquiry functions, we usually define an integer constant.

integer,parameter :: dp=selected_real_kind(15,307)

It would be best if you used the constants to specify the precision of the literals and variables.

When you mix different precisions in a formula, the result is aligned to the highest precision. For example, 1.0e0 + 1.0d0 results in a double precision value (2.0d0). Also, a formula with integer numbers and real numbers should produce a real number.

Overflow and non-numbers

We have seen that a computer can hold a number in a particular range. If a number is out of range by arithmetics, the value should be overflow. Also, if the program performs mathematically incorrect operations, the result would be unstable. We will see such a corner case here.

Overflow

The standard integer uses 4 bytes that can hold a number between \(-2147483648\) and \(2147483647\). The upper limit of a numeric type is available with huge(). The following program tries to have a number out of range.

program num
   implicit none
   integer :: k
   k = huge(k)
   print *,k
   k = k + 1
   print *,k
end program num
  2147483647
 -2147483648

It does not produce any errors. The largest number becomes the lower bound when 1 is added. It is called overflow, but you can not detect it by default. The overflow happens a negative value like using the integer variable k, k=-2147483648; k=k-1.

The smallest negative integer is just next to the largest positive integer. The reason is due to a representation of integer number in Fortran (and any other languages). Please see the textbooks in computer science for details.

Infinity

With a floating-pint number, e.g., real, the overflow results in being a weird value. The above code is modified to have a real variable (adding another huge instead of adding 1).

program num
   implicit none
   real :: x
   x = huge(x)
   print *,x
   x = x + huge(x)
   print *,x
end program num
   3.40282347E+38
         Infinity

Variable k has a value, Infinity, which is not a standard number but still a valid value in Fortran. There is another value -Infinity from a formula to produce a negative value.

To generate Infinity and -Infinity, you can use an intrinsic module ieee_arithmetic. The function ieee_value in a combination with a constant ieee_positive_inf (or ieee_negative_inf) produces the value. You can use it to check if an arithmetic result is Infinity or not.

program num
   use ieee_arithmetic
   implicit none
   real :: x,y
   x = ieee_value(x, ieee_positive_inf)
   y = ieee_value(y, ieee_negative_inf)
   print *,x,y
end program num
         Infinity        -Infinity

Once you have had Infinity or -Infinity, you can use it in the formula, but the value may propagate to the result. I leave it to an exercise for readers.

There is a intrinsic function to test whether a number has the infinity or not. It is ieee_is_finite() defined in an intrinsic module, ieee_arithmetic. If the value is finite, it returns .true.; otherwise (i.e., infiity), .false..

program num
   use ieee_arithmetic
   implicit none
   real :: x,y
   x = ieee_value(x, ieee_positive_inf)
   y = 1.0
   print *,ieee_is_finite(x),ieee_is_finite(y)
end program num
 F T

NaN: Not-a-number

There is another status as a result of floating-point arithmetic. If the result is not mathematically defined, like \(0/0\), the result becomes a value, not-a-number or simply NaN. It is not a number but a valid floating-point value. Note that it happens only in floating-point arithmetics; for integer arithmetics, it never happens. See the following example.

program num
   implicit none
   integer :: k
   real :: x
   k = 0
   x = 0.0
   k = k/k
   x = x/x
   print *,k,x
end program num
           1              NaN

The value of NaN is also available with the function, ieee_value in the module, ieee_arithmetic. See the following example for details.

program num
   use ieee_arithmetic
   implicit none
   real :: x
   x = ieee_value(x, ieee_quiet_nan)
   print *,x
end program num
              NaN

The value of NaN often happens in computations involving tiny numeric values close to 0. It indicates an arithmetic error, and you should not proceed with the computations because, once NaN is in a formula, the result is always NaN. Also, it is tricky to check if a value is NaN or not because NaN is not a number, and the regular rule of comparisons does not work.

program num
   use ieee_arithmetic
   implicit none
   real :: y
   ! both NaN
   y = 0.0
   y = y/y
   if(y == ieee_value(x, ieee_quiet_nan)) then
      print *,"It is NaN."
   else
      print *,"It is not NaN."
   end if
end program num
 It is not NaN.

A typical idiom to check NaN is to compare it with itself like y==y. While it should always be true for any numbers, NaN returns .false.. So, the comparison looks like if(y/=y) print *,"It is NaN.".

In Modern Fortran, it is better to use ieee_is_nan() to test NaN. It is also defined in ieee_arithmetic. If the value is NaN, it returns .true.; otherwise .false.. See the following example.

program num
   use ieee_arithmetic
   implicit none
   real :: x,y
   x = ieee_value(x, ieee_positive_inf)
   if(ieee_is_nan(x)) then
      print *,"It is NaN."
   else
      print *,"It is not NaN."
   end if
end program num

Exercises

  1. Using a real variable, compute \(0/0\), \(1/0\), and \(0/1\).
  2. See how Infinity (or -Infinity) works in the various arithmetic formula, including addition, subtraction, division, multiplication, and mathematical functions. See the result from multiplying 0 to Infinity.
  3. See how NaN works in the arithmetic formula as the above exercise.

Floating-point arithmetic

Significant digits

Consider measurements collected in the field: milk yield, body weight, and so on. Theoretically, these measurements potentially have infinite digits, and the same values never come out. However, we usually limit the number of digits of measurement by rounding the last digit. Such a numerical value has inaccurate information, especially in the last digit. When a value has uncertainty caused by the rounding, the number of decimal digits of the value is called the significant digits.

Similarly, in many programming languages, the real number is treated as a floating-point number. The last digit is likely rounded, and therefore, the it is not accurate. In such a case, the significant digits should be the apparent number of digits.

The significant digits can be counted with the significand of the exponential expression. For example, the value 12.34 is expressed as \(1.234\times 10^{3}\) so, the significant digit is 4 implying the 5th digit has been rounded. Similarly, the value 0.01234 has 4 significant digits because it is \(1.234\times 10^{-2}\).

The significant digits may decrease but never increase in numerical computations. Consider the subtraction of two floating-point values: \(1.234-1.233\). The result is \(0.001\), and the significant digit is 1. The significant digits of computation are limited, with the smallest significant digits of the original values. This is true for any computations with numerical values with significant digits.

Recurring decimal in the binary numeral system

In the decimal numeral system, a number is not rational, and it is represented as repeating decimals. For example, \(1/3\) is \(0.333333...=0.\dot{3}\). In this case, you have to add the exact \(1/3\) with all repeating decimals to restore 1 by an addition \(1/3+1/3+1/3\). If you use a finite representation of the value, like \(0.3333\), you lose some precision, and you get an approximation instead of the exact result.

In Fortran, a numeric value is represented as binary, and some numbers are rational in decimal but not in binary. A typical example is \(0.1\), which becomes a repeating decimal number in the binary numeral system. It may result in unexpected results.

Exercises

  1. Some people show 4 or more digits in their regression coefficients or heritability from statistical analyses. Is it advantageous to show all the small digits?
  2. Considering the above question, can all the computations be done with the low-precision arithmetic? If not, why?
  3. Sum up \(0.1\) ten times and see what the result is.

Error

Computations using floating-point numbers may produce non-intuitive results i.e., erroneous values. It is mainly caused by a finite expression of real numbers. The error is categorized into several groups.

Rounding error

Rounding error (or round-off error) occurs when a number with high precision (more decimal digits) is rounded to a number with low precision (fewer digits). Typically, it happens when you mix single and double precision numbers in a formula and round the result into single-precision or integer number. It is also caused when a formula has a high-precision literal, but it is defined as lower precision.

program error
   print *,3.14159265359         ! treated as single precision
   print *,3.14159265359*2       ! same as above
   print *,dble(3.14159265359*2) ! lower precision arithmetic
end program error

Cancellation

When two values with a similar magnitude are added or subtracted together, the significant digits will be decreased. See the following example.

program error
   print *,sqrt(111.00)
   print *,sqrt(111.01)
   print *,sqrt(111.01) - sqrt(111.00)
end program error
   10.5356541
   10.5361280
   4.73976135E-04

For the first two values, the last digit is rounded, and therefore it has an error. The effective digit of each of them is 8. The result of subtraction is shown as 0.00473976135 but, the numbers below the 8th digit under the decimal point are erroneous. In the end, the result has only 6 effective digits, which has less information than the original numbers. If you repeat such computations many times, an amount of information in the result should be greatly missing.

Loss of trailing digit

It is a combination of rounding error and cancellation, but it happens when the absolute difference between two values is huge.

program error
   print *,12345678.0 - 0.1
end program error
   12345678.0

The subtraction is ignored because the latter is too small, and the former has no rooms to handle such a small number.

Summary

There is some possible source to produce numerical errors. To avoid errors, we have two options.

In cases where a scientist without a computer-science-background writes numerical programs, the reduction of errors is not necessarily a priority in their software. Sometimes the error is serious but sometimes not.

Exercises

Exercises in this section are not provided.

Some algorithms with numerics

Product of many numbers

Algorithm

There is a case where you have to multiply many numbers many times. In quantitative genetics, this method is used in the computation of likelihood in derivative-free REML (DF REML). As you can imagine, when many numbers are involved, the product may become huge (or tiny), and it may accumulate errors during the computation. A typical method to handle the issue is to convert multiplication to addition by log for higher-precision variables. For positive values, the mathematical formula of this algorithm is as follows.

\[ p = \prod_{i=1}^{n}x_{i} = \exp\left[\sum_{i=1}^{n}\log x_i\right]\quad\text{for}\quad x_i>0 \]

Naive program

Here is a simple program to compute the product of elements in an array. The program works as:

  1. Read \(n\) from the keyboard.
  2. Initialize the result variable as p=1.
  3. Repeat to read a value, x, from keyboard \(n\) times. The input is multiplied by the previous result as p = p*x
  4. Compute the product.
  5. print the result on the screen.
program num
   implicit none
   integer :: n,i
   real :: p,x

   read *,n
   if(n<=0) stop "error: negative n"

   ! read values from keyboard
   p = 1
   do i=1,n
      read(*,*) x
      if(x<=0) stop "error: negative value"
      p = p*x
   end do

   print *,"product = ",p
end program num

Modified program

To implement the modified algorithm, you have to change some code.

The modified code is as follows.

program num
   implicit none
   integer :: n,i
   real :: p,x

   read *,n
   if(n<=0) stop "error: negative n"

   ! read values from keyboard
   p = 0
   do i=1,n
      read(*,*) x
      if(x<=0) stop "error: negative value"
      p = p + log(x)
   end do

   p = exp(p)
   print *,"product = ",p
end program num

Exercises

  1. Modify the above program to support 0.
  2. Modify the above program to support 0 and negative values.

Back to index.html.