Fortran floating-point arithmetic
last modified November 14, 2026
A floating-point number stores a value as a sign, an exponent and a fraction, in the same way that scientific notation writes 6.02 times 10 to the 23rd. The fraction is called the significand or mantissa, and its length decides how many digits the number can carry; the exponent decides how large or small it can be. The spacing between representable values is not constant: it doubles with every power of two, so the same type that resolves steps of 1e-7 near 1.0 resolves steps of 8 near 1e8.
This article inspects those limits with epsilon,
spacing and the kind inquiry functions, decodes the bits of a
single-precision value, shows why 0.1 + 0.2 is not
0.3, compares summation orders and the Kahan compensation, uses
the ieee_arithmetic module for NaN, infinity, exceptions and
rounding modes, and finishes with a quadratic formula whose textbook version
loses half of its answer.
Precision, range and spacing
Every real kind answers six questions about itself. digits is the
number of binary digits in the significand, precision the number of
decimal digits that survive a conversion back and forth, range
the decimal exponent range, and tiny and huge the
smallest and largest positive values. epsilon is the gap between
1.0 and the next representable number, and spacing(x) is the gap
at x, which is the more useful quantity once values grow.
program f1
use iso_fortran_env, only: real32, real64
implicit none
real(real32) :: s
s = 1.0_real32
print '(a, i0)', 'digits(real32) = ', digits(s)
print '(a, i0)', 'digits(real64) = ', digits(1.0_real64)
print '(a, i0)', 'precision(real32) = ', precision(s)
print '(a, i0)', 'precision(real64) = ', precision(1.0_real64)
print '(a, i0)', 'range(real32) = ', range(s)
print '(a, i0)', 'range(real64) = ', range(1.0_real64)
print '(a, i0)', 'maxexponent(real32) = ', maxexponent(s)
print '(a, i0)', 'minexponent(real32) = ', minexponent(s)
print '(a, es12.4)', 'epsilon(real32) = ', epsilon(s)
print '(a, es12.4)', 'tiny(real32) = ', tiny(s)
print '(a, es12.4)', 'huge(real32) = ', huge(s)
print '(a, es12.4)', 'spacing(1.0) = ', spacing(s)
print '(a, es12.4)', 'spacing(1.0e8) = ', spacing(1.0e8_real32)
print '(a, es12.4)', 'nearest(1.0) - 1.0 = ', nearest(s, 1.0_real32) - s
print '(a, es12.4)', 'epsilon(real64) = ', epsilon(1.0_real64)
print '(a, es12.4)', 'spacing(1.0_real64) = ', spacing(1.0_real64)
end program f1
A default real has 24 binary digits, which is about 6 decimal digits, and a
real64 has 53, about 15. The exponent of a default real reaches
128 and drops to -125, so huge is 3.4028e38 and tiny
1.1755e-38. epsilon equals spacing(1.0) at 1.19e-7,
and nearest(1.0, 1.0) is exactly one epsilon above 1. The spacing
at 1e8 is 8, which is why a single-precision program cannot count by ones once
its values reach a hundred million.
$ gfortran f1.f90 -o f1 $ ./f1 digits(real32) = 24 digits(real64) = 53 precision(real32) = 6 precision(real64) = 15 range(real32) = 37 range(real64) = 307 maxexponent(real32) = 128 minexponent(real32) = -125 epsilon(real32) = 1.1921E-07 tiny(real32) = 1.1755E-38 huge(real32) = 3.4028E+38 spacing(1.0) = 1.1921E-07 spacing(1.0e8) = 8.0000E+00 nearest(1.0) - 1.0 = 1.1921E-07 epsilon(real64) = 2.2204E-16 spacing(1.0_real64) = 2.2204E-16
The bits of a real
A single-precision value is 32 bits: one sign bit, eight exponent bits and 23
significand bits. The exponent is stored with a bias of 127, and the
significand has an implicit leading 1 that is not stored, so the value is
(1 + m / 2**23) * 2**(e - 127). The intrinsic transfer
reinterprets the bits of a real as an integer without converting anything, and
the bit procedures of the previous article take it apart.
program f6
use iso_fortran_env, only: real32, int32
implicit none
real(real32) :: values(6), decoded
integer(int32) :: bits, sign, exponent, mantissa
integer :: i
values = [1.0_real32, 0.5_real32, 3.0_real32, -2.0_real32, &
0.1_real32, 0.0_real32]
print '(a)', 'value sign exponent mantissa'
do i = 1, size(values)
bits = transfer(values(i), bits)
sign = iand(ishft(bits, -31), 1)
exponent = iand(ishft(bits, -23), 255)
mantissa = iand(bits, z'7FFFFF')
print '(f6.2, 3x, b1, 7x, i4, 5x, z6.6)', values(i), sign, &
exponent, mantissa
end do
bits = transfer(0.1_real32, bits)
print '(a, b32.32)', '0.1 in binary: ', bits
exponent = iand(ishft(bits, -23), 255)
mantissa = iand(bits, z'7FFFFF')
decoded = (1.0_real32 + real(mantissa, real32) / 2.0_real32**23) * &
2.0_real32**(exponent - 127)
print '(a, f0.12)', 'decoded value: ', decoded
end program f6
The table reads like a page from a hardware manual. The value 1.0 has exponent
127, which is 2**0, and a zero significand; 0.5 halves the
exponent to 126; 3.0 has exponent 128 and a significand of
400000, which is 1.5, and 1.5 times 2 is 3. The negative value
differs only in the sign bit, and zero has every field zero. The value 0.1
cannot be written in binary, so its significand is the repeating pattern
4CCCCD rounded up, and decoding it gives 0.100000001490.
$ gfortran f6.f90 -o f6 $ ./f6 value sign exponent mantissa 1.00 0 127 000000 0.50 0 126 000000 3.00 0 128 400000 -2.00 1 128 000000 0.10 0 123 4CCCCD 0.00 0 0 000000 0.1 in binary: 00111101110011001100110011001101 decoded value: .100000001490
Why 0.1 plus 0.2 is not 0.3
One tenth has no finite binary expansion, so the constant 0.1 is already the nearest representable value, and every operation on it carries that small error. Adding 0.2 doubles the error and rounding it to the nearest value gives 0.30000000000000004 in double precision, while the constant 0.3 is stored as 0.29999999999999999. The difference is 5.55e-17, a single bit at that magnitude, and it is enough to make an equality test fail.
program f2
use iso_fortran_env, only: real32, real64
implicit none
real(real64) :: x, y, z
real(real32) :: one, small, absorbed
x = 0.1_real64
y = 0.2_real64
z = 0.3_real64
print '(a, f0.17)', 'stored 0.1: ', x
print '(a, f0.17)', 'x + y: ', x + y
print '(a, f0.17)', 'stored 0.3: ', z
print '(a, es12.4)', 'x + y - z: ', (x + y) - z
print '(a, es12.4)', '0.5+0.25-0.75: ', (0.5_real64 + 0.25_real64) - &
0.75_real64
print '(a, es12.4)', 'spacing(1.0_real64):', spacing(1.0_real64)
one = 1.0_real32
small = 1.0e-8_real32
absorbed = (one + small) - one
print '(a, es12.4)', 'single 1.0+1e-8-1.0:', absorbed
print '(a, es12.4)', 'spacing(1.0_real32):', spacing(1.0_real32)
print '(a, l1)', 'tolerance 1e-6: ', nearly_equal(x + y, z, 1.0e-6_real64)
print '(a, l1)', 'tolerance 1e-16: ', nearly_equal(x + y, z, 1.0e-16_real64)
print '(a, l1)', 'half and quarter: ', nearly_equal(0.5_real64 + &
0.25_real64, 0.75_real64, 1.0e-15_real64)
print '(a, l1)', 'close, 1e-6: ', nearly_equal(1.0_real64, &
1.000001_real64, 1.0e-6_real64)
print '(a, l1)', 'close, 1e-9: ', nearly_equal(1.0_real64, &
1.000001_real64, 1.0e-9_real64)
contains
logical function nearly_equal(a, b, tolerance) result(same)
real(real64), intent(in) :: a, b, tolerance
real(real64) :: scale
scale = max(abs(a), abs(b))
if (scale <= 0.0_real64) then
same = .true.
else
same = abs(a - b) <= tolerance * scale
end if
end function nearly_equal
end program f2
The first three lines show the stored values with seventeen decimals; the fourth shows the residue after subtracting 0.3. Values that are sums of powers of two, such as 0.5, 0.25 and 0.75, are exact, and their combination leaves a clean zero. Single precision is coarser still: adding 1e-8 to 1.0 changes nothing at all, because the spacing at 1.0 is 1.19e-7. The last four lines replace equality with a relative tolerance, which accepts a difference of 5.55e-17 at 1e-6 and rejects it at 1e-16.
$ gfortran f2.f90 -o f2 $ ./f2 stored 0.1: .10000000000000001 x + y: .30000000000000004 stored 0.3: .29999999999999999 x + y - z: 5.5511E-17 0.5+0.25-0.75: 0.0000E+00 spacing(1.0_real64): 2.2204E-16 single 1.0+1e-8-1.0: 0.0000E+00 spacing(1.0_real32): 1.1921E-07 tolerance 1e-6: T tolerance 1e-16: F half and quarter: T close, 1e-6: T close, 1e-9: F
Summation order and the Kahan compensation
Every addition rounds, and the rounding errors of a long sum accumulate in the direction of the running total. Adding 0.1 a million times in single precision is the standard demonstration: the true total is 100000, and the loop drifts by almost a thousand. The Kahan summation keeps a second variable that holds the part of each addition that was rounded away, and subtracts it from the next term, which recovers most of the lost precision at the cost of three extra operations per term.
program f3
use iso_fortran_env, only: real32
implicit none
integer, parameter :: n = 1000000
real(real32) :: naive, compensated, c, y, t
real(real32) :: parts(3)
integer :: i
naive = 0.0
do i = 1, n
naive = naive + 0.1_real32
end do
print '(a, f14.3)', 'one million times 0.1:', naive
print '(a, f14.3)', 'the exact value: ', 100000.0_real32
compensated = 0.0
c = 0.0
do i = 1, n
y = 0.1_real32 - c
t = compensated + y
c = (t - compensated) - y
compensated = t
end do
print '(a, f14.3)', 'Kahan summation: ', compensated
print '(a, es12.4)', 'plain error: ', naive - 100000.0_real32
print '(a, es12.4)', 'Kahan error: ', compensated - 100000.0_real32
parts = [1.0e8_real32, 1.0_real32, -1.0e8_real32]
print '(a, f10.1)', 'left to right: ', (parts(1) + parts(2)) + parts(3)
print '(a, f10.1)', 'cancel first: ', (parts(1) + parts(3)) + parts(2)
print '(a, f10.1)', 'intrinsic sum: ', sum(parts)
end program f3
The naive loop ends at 100958.344, an error of 958, while the compensated sum
lands exactly on 100000. The three values 1e8, 1 and -1e8 show the same effect
in miniature: added from left to right the small value is absorbed and the
result is 0, but adding the cancelling pair first keeps the 1. The intrinsic
sum follows the array order and also returns 0, so a program that
needs the exact result must either sort the terms or compensate.
$ gfortran f3.f90 -o f3 $ ./f3 one million times 0.1: 100958.344 the exact value: 100000.000 Kahan summation: 100000.000 plain error: 9.5834E+02 Kahan error: 0.0000E+00 left to right: 0.0 cancel first: 1.0 intrinsic sum: 0.0
NaN and infinity
IEEE arithmetic defines two values that are not numbers. Infinity is
the result of an overflow or of a division by zero, and it behaves like a
number larger than any other. NaN, short for not a number, is the
result of an operation with no sensible answer, such as zero divided by zero or
infinity minus infinity. NaN propagates: every operation that touches one
produces another, and every comparison with one is false. The module
ieee_arithmetic supplies the tests and the constructors, and the
g0 descriptor prints both values by name.
program f4
use ieee_arithmetic
implicit none
real :: zero, one, infinity, nan, x
zero = 0.0
one = 1.0
infinity = one / zero
nan = zero / zero
print '(a, g0)', '1.0 / 0.0: ', infinity
print '(a, g0)', '0.0 / 0.0: ', nan
print '(a, g0)', 'nan + 1.0: ', nan + 1.0
print '(a, l1)', 'is infinity finite: ', ieee_is_finite(infinity)
print '(a, l1)', 'is nan a nan: ', ieee_is_nan(nan)
print '(a, l1)', 'is infinity a nan: ', ieee_is_nan(infinity)
print '(a, l1)', 'is 1.0 finite: ', ieee_is_finite(1.0)
print '(a, l1)', 'is -0.0 negative: ', ieee_is_negative(-0.0)
print '(a, l1)', 'is 0.0 negative: ', ieee_is_negative(0.0)
x = ieee_value(x, ieee_quiet_nan)
print '(a, g0)', 'ieee_value quiet nan: ', x
x = ieee_value(x, ieee_negative_inf)
print '(a, g0)', 'ieee_value -inf: ', x
print '(a, es12.4)', 'distance to 1.0: ', abs(x + 1.0)
end program f4
The first two lines produce Inf and NaN, and the
third shows the propagation: NaN plus one is still NaN.
ieee_is_finite rejects the infinity and accepts an ordinary value,
and ieee_is_nan identifies the NaN without ever comparing it.
ieee_is_negative reports the sign of zero, so -0.0 is negative
while 0.0 is not, which matters when a sign is carried through a computation.
The last lines build the two special values with ieee_value and
show that the distance from minus infinity to one is infinite.
$ gfortran f4.f90 -o f4 $ ./f4 1.0 / 0.0: Inf 0.0 / 0.0: NaN nan + 1.0: NaN is infinity finite: F is nan a nan: T is infinity a nan: F is 1.0 finite: T is -0.0 negative: T is 0.0 negative: F ieee_value quiet nan: NaN ieee_value -inf: -Inf distance to 1.0: Infinity
Exceptions and rounding modes
The processor records flags for the exceptional events of a computation: invalid operation, division by zero, overflow, underflow and inexact result. A program can clear a flag, run an operation and read the flag back, which is how a numerical library decides whether its result is trustworthy. The same module sets the rounding mode, which decides whether the last bit of a result is rounded to the nearest value, toward zero, toward plus infinity or toward minus infinity.
program f5
use ieee_arithmetic
use iso_fortran_env, only: real64
implicit none
real(real64) :: zero, one, two, result
logical :: flag
zero = 0.0_real64
one = 1.0_real64
two = 2.0_real64
print '(a, l1)', 'invalid flag supported: ', ieee_support_flag(ieee_invalid)
print '(a, l1)', 'rounding modes supported: ', &
ieee_support_rounding(ieee_down)
call ieee_set_flag(ieee_invalid, .false.)
call ieee_set_flag(ieee_divide_by_zero, .false.)
result = zero / zero
call ieee_get_flag(ieee_invalid, flag)
print '(a, l1)', '0/0 raises invalid: ', flag
call ieee_set_flag(ieee_invalid, .false.)
result = one / zero
call ieee_get_flag(ieee_invalid, flag)
print '(a, l1)', '1/0 raises invalid: ', flag
call ieee_get_flag(ieee_divide_by_zero, flag)
print '(a, l1)', '1/0 raises divide by zero: ', flag
result = two / 3.0_real64
print '(a, f0.17)', 'nearest: ', result
call ieee_set_rounding_mode(ieee_down)
result = two / 3.0_real64
print '(a, f0.17)', 'down: ', result
call ieee_set_rounding_mode(ieee_up)
result = two / 3.0_real64
print '(a, f0.17)', 'up: ', result
call ieee_set_rounding_mode(ieee_nearest)
result = two / 3.0_real64
print '(a, f0.17)', 'nearest: ', result
end program f5
Zero divided by zero raises the invalid flag, while one divided by zero raises divide-by-zero and leaves the invalid flag clear, which is exactly the distinction between a meaningless expression and an infinite limit. The four printings of two thirds differ in the last bit: the nearest mode gives 0.66666666666666663, rounding down gives 0.66666666666666662, and rounding up gives 0.66666666666666675, one spacing above. The mode is restored at the end, because it affects every later operation in the program.
$ gfortran f5.f90 -o f5 $ ./f5 invalid flag supported: T rounding modes supported: T 0/0 raises invalid: T 1/0 raises invalid: F 1/0 raises divide by zero: T nearest: .66666666666666663 down: .66666666666666662 up: .66666666666666675 nearest: .66666666666666663
A formula that cancels itself
The quadratic formula is the classic victim of catastrophic
cancellation. When b is much larger than the other
coefficients, -b + sqrt(b**2 - 4ac) subtracts two nearly equal
numbers, and the significant digits of the small root disappear. The cure is
to compute the root that does not cancel, q = -0.5 * (b + sign(sqrt(d),
b)), and obtain the other root from the product of the roots,
x1 * x2 = c / a.
program f7
use iso_fortran_env, only: real64
implicit none
real(real64) :: a, b, c, disc, root, q, x1, x2
a = 1.0_real64
b = 1.0e8_real64
c = 1.0_real64
disc = b * b - 4.0_real64 * a * c
root = sqrt(disc)
print '(a, es14.7)', 'naive root 1: ', (-b + root) / (2.0_real64 * a)
print '(a, es14.7)', 'naive root 2: ', (-b - root) / (2.0_real64 * a)
print '(a, es14.7)', 'residual of root 1:', residual((-b + root) / &
(2.0_real64 * a))
q = -0.5_real64 * (b + sign(root, b))
x1 = q / a
x2 = c / q
print '(a, es14.7)', 'stable root 1: ', x1
print '(a, es14.7)', 'stable root 2: ', x2
print '(a, es14.7)', 'residual of root 2:', residual(x2)
print '(a, es14.7)', 'product of roots: ', x1 * x2
print '(a, es14.7)', 'c / a: ', c / a
contains
real(real64) function residual(x) result(r)
real(real64), intent(in) :: x
r = a * x * x + b * x + c
end function residual
end program f7
For the equation x**2 + 1e8*x + 1 = 0 the roots are about -1e8 and
-1e-8. The textbook formula returns -7.45e-9 for one root and -1e8 for the other,
an error of a quarter on the small one, and the residual there is 0.25 instead
of zero. The stable form returns -1e-8 with a residual of 1.1e-16, and the product
of the two roots is exactly the constant term. The same rearrangement applies
to the vector length, the law of cosines and any expression that subtracts
nearly equal quantities.
$ gfortran f7.f90 -o f7 $ ./f7 naive root 1: -7.4505806E-09 naive root 2: -1.0000000E+08 residual of root 1: 2.5494194E-01 stable root 1: -1.0000000E+08 stable root 2: -1.0000000E-08 residual of root 2: 1.1102230E-16 product of roots: 1.0000000E+00 c / a: 1.0000000E+00
Next steps
The data types article introduces
real64 and the other kinds whose limits are printed here, and the
literals article shows the kind suffixes that
select them. The intrinsic numeric
functions article covers sqrt, exp and the
rounding functions, the
formatted output article explains the
es and g descriptors used throughout, and the
bit manipulation article supplies the
transfer and bit operations that decode the representation.
Source
Fortran 2018 standard (J3/18-007r1)
Fortran-lang: learning resources
In this article we have read the limits of a real kind with
epsilon, spacing, digits and
range, decoded the sign, exponent and significand of a single
precision value, seen why 0.1 plus 0.2 misses 0.3 and why a small addition
disappears into a large one, compared a naive sum with Kahan summation,
constructed and tested NaN and infinity, read the exception flags and changed
the rounding mode, and rewritten the quadratic formula to avoid catastrophic
cancellation. The previous article,
Fortran bit manipulation, explains the
bit operations used to inspect the representation.
Author
List all Fortran tutorials.