# More on numbers and variables

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.

- Integer number
- An exact number in a range.
- Literal: a sign (
`+`

or`-`

) and a sequence of numbers (`0`

to`9`

); e.g.,`123045`

and`-612`

- Default: 4-byte integer between \(-2147483648\) and \(2147483647\).

- Real number
- A floating-point number approximating a value within the defined precision.
- Literal using the decimal separator:
`3.14`

and`-0.002`

for single precision - Literal using exponential i.e., \(x\times 10^{y}\):
`0.314e1`

and`12.234e-3`

for single precision`0.314d1`

and`12.234d-3`

for double precision

- Default:
- 4-byte single precision (approximately between \(\pm 1.18\times 10^{-38}\) and \(\pm 3.40\times 10^{38}\))
- 8-byte double precision (approximately between \(\pm 2.23\times 10^{-308}\) and \(\pm 1.80\times 10^{308}\))

- Logical value
- Alternatative value
- Literal:
`.true.`

or`.false.`

- Characters
- A sequence of letters (string)
- Literal: encircled string by
`'`

or`"`

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.

`integer ::`

for integer variables (default: 4-byte)`real ::`

for real variables (default: 4-byte single precision)`logical ::`

for logical variables`character ::`

for a single character`character(len=n) ::`

for a string with`n`

letters

Here are a couple of attributes that we have introduced.

`parameter`

: the variable becomes a constant, which does not change the value.`allocatable`

: the array variable is allocatable i.e., resizable.

### 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.

- Precision inquiry function
- Integer:
`selected_int_kind(n)`

for the range of \(-10^n<i<10^n\). - Real:
`selected_real_kind(p,r)`

for at least`p`

digits in significand and a decimal-exponent range of at least`r`

.

- Integer:
- Constants defined in intrinsic module
`iso_fortran_env`

- Integer:
`int8`

,`int16`

,`int32`

, and`int64`

for 8-bit (1-byte), 16-bit (2-byte), 32-bit (4-byte), and 64-bit (8-byte) integers. - Real:
`real32`

for 32-bit (4-byte single precision) and`real64`

for 64-bit (8-byte double precision) numbers.

- Integer:

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.

- Literals: literal expression +
`_`

+ precision constant- For integer number:
`128_int64`

- For real number:
`3.14_real64`

- For integer number:
- Real literals: use
`d`

instead`e`

in exponential expression for double precision- Example:
`3.14e0`

for single precision and`3.14d0`

for double precision

- Example:
- Variables: precision code enclosed by
`()`

in the declaration- Example:
`integer(int64) :: i`

- Example:
`real(real64) :: x`

- Example:

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

- Using a real variable, compute \(0/0\), \(1/0\), and \(0/1\).
- 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`

. - 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

- 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?
- Considering the above question, can all the computations be done with the low-precision arithmetic? If not, why?
- 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.

- Use higher precision values. Sometimes it is just a workaround of the issue.
- Modify a formula to reduce errors. It is essential but more technical and complicated to create a program.

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:

- Read \(n\) from the keyboard.
- Initialize the result variable as
`p=1`

. - Repeat to read a value,
`x`

, from keyboard \(n\) times. The input is multiplied by the previous result as`p = p*x`

- Compute the product.
- 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 staring value of
`p`

is 0, instead of 1. - The multiplication
`p*x`

is replaced with`p+log(x)`

.

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

- Modify the above program to support 0.
- Modify the above program to support 0 and negative values.

Back to index.html.