Fortran random numbers
last modified November 14, 2026
A computer cannot produce a truly random number, because its next state is
determined by its previous one. What it produces instead is a
pseudorandom sequence: numbers that pass statistical tests for
randomness but repeat exactly when the generator starts from the same
seed. That repeatability is a feature, not a defect, since a simulation
that cannot be reproduced cannot be debugged. Fortran provides one generator
through the intrinsic random_number and controls its state through
random_seed.
This article draws uniform values and checks their distribution, seeds the generator twice to show that the sequence repeats, converts the values to dice rolls, estimates pi by Monte Carlo sampling, shuffles an array without bias, generates normally distributed numbers with the Box-Muller transform, and simulates the birthday problem. Every program sets the seed explicitly, so every output below can be reproduced.
Uniform values between zero and one
call random_number(x) assigns a value from the interval from 0
inclusive to 1 exclusive, and the same call with an array fills every element.
The values are uniformly distributed, so each tenth of the interval should
receive about a tenth of a large sample, and the mean should approach 0.5. The
program draws a hundred thousand values, prints the smallest and largest, and
counts how many fall into each of ten equal bins.
program r1
implicit none
integer, parameter :: samples = 100000, bins = 10
real :: values(samples), x
integer :: counts(bins), i, bin
call seed_with(12345)
call random_number(x)
print '(a, f8.6)', 'one value: ', x
call random_number(values)
print '(a, f9.5)', 'mean of 100000: ', sum(values) / samples
print '(a, f9.6)', 'smallest: ', minval(values)
print '(a, f9.6)', 'largest: ', maxval(values)
counts = 0
do i = 1, samples
bin = int(values(i) * bins) + 1
if (bin > bins) bin = bins
counts(bin) = counts(bin) + 1
end do
print '(a)', 'histogram of ten bins:'
do i = 1, bins
print '(a, f3.1, a, f3.1, a, i6)', ' ', real(i - 1) / bins, &
' to ', real(i) / bins, ': ', counts(i)
end do
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r1
The mean is 0.50087, close to one half, and the extreme values are 0.000004 and 0.999999, which confirms that both ends of the interval are reachable but the upper one is never included. The ten bins hold between 9861 and 10146 values each, a spread of about one and a half percent, which is the expected size of the fluctuation for ten thousand samples per bin. The histogram carries no trend from the first bin to the last, which is the visible sign of a uniform distribution.
$ gfortran r1.f90 -o r1 $ ./r1 one value: 0.465394 mean of 100000: 0.50087 smallest: 0.000004 largest: 0.999999 histogram of ten bins: 0.0 to 0.1: 10088 0.1 to 0.2: 10017 0.2 to 0.3: 9967 0.3 to 0.4: 9885 0.4 to 0.5: 9861 0.5 to 0.6: 9953 0.6 to 0.7: 9940 0.7 to 0.8: 10046 0.8 to 0.9: 10097 0.9 to 1.0: 10146
Seeds and reproducible runs
The generator keeps its state in an array of integers. random_seed(size=n)
reports how many, random_seed(get=state) copies the current state
out, and random_seed(put=state) restores it. Without an explicit
put, gfortran seeds the generator from the operating system when
the program starts, so two runs differ; with a fixed seed the whole sequence is
fixed, and saving the state in the middle of a run allows the program to return
to that point later.
program r2
implicit none
real :: first(3), again(3), other(3)
integer :: n
integer, allocatable :: state(:)
call random_seed(size=n)
print '(a, i0)', 'integers in the seed: ', n
call seed_with(7)
call random_number(first)
call seed_with(7)
call random_number(again)
call seed_with(8)
call random_number(other)
print '(a, 3f9.6)', 'seed 7: ', first
print '(a, 3f9.6)', 'seed 7: ', again
print '(a, 3f9.6)', 'seed 8: ', other
print '(a, es10.2)', 'difference after reseeding: ', &
maxval(abs(first - again))
print '(a, es10.2)', 'difference for seed 8: ', &
maxval(abs(first - other))
call seed_with(7)
call random_number(first)
allocate (state(n))
call random_seed(get=state)
call random_number(again)
call random_seed(put=state)
call random_number(other)
print '(a, es10.2)', 'difference after restoring: ', &
maxval(abs(again - other))
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r2
This build uses eight integers of seed. Seeding twice with 7 produces exactly the same three values, and seeding with 8 produces a different stream, which the differences printed in scientific notation confirm: zero for the repeated seed and 0.176 for the changed one. The last three lines are the more useful trick: the state is saved after one draw, two more draws are made, the saved state is put back, and the next values match the ones drawn before the restore. A program can therefore checkpoint a long simulation and resume it exactly.
$ gfortran r2.f90 -o r2 $ ./r2 integers in the seed: 8 seed 7: 0.470820 0.116794 0.135278 seed 7: 0.470820 0.116794 0.135278 seed 8: 0.470819 0.116792 0.311059 difference after reseeding: 0.00E+00 difference for seed 8: 1.76E-01 difference after restoring: 0.00E+00
Integers in a range
A die throw needs one of six values with equal probability. Multiplying a uniform value by six and taking the integer part gives 0 to 5, and adding one shifts the range to 1 to 6; the multiplication also scales the equal widths, so each face keeps the same probability. Rounding instead of truncating would give the two end faces half the width of the middle ones, because the rounding boundaries fall at the edges of the interval. The mean of a fair die is 3.5, whatever the number of rolls.
program r3
implicit none
integer, parameter :: rolls = 60000
integer :: counts(6), i, face
real :: u, total
call seed_with(2024)
counts = 0
do i = 1, rolls
call random_number(u)
face = int(u * 6.0) + 1
counts(face) = counts(face) + 1
end do
total = 0.0
do face = 1, 6
print '(a, i0, a, i6, a, f6.2, a, a)', 'face ', face, ': ', &
counts(face), ' ', 100.0 * counts(face) / rolls, '% ', &
repeat('#', counts(face) / 500)
total = total + real(face) * counts(face)
end do
print '(a, f7.4)', 'mean of the rolls: ', total / rolls
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r3
Sixty thousand rolls give face counts between 9812 and 10231, the largest about two and a half percent away from the expected ten thousand, and the mean is 3.5079 instead of 3.5. The bar length is the count divided by 500, so a full bar would need fifty thousand rolls, and the differences between the bars are the sampling noise. Increasing the number of rolls narrows the spread of the counts but never removes it; the width falls with the square root of the number of rolls.
$ gfortran r3.f90 -o r3 $ ./r3 face 1: 9861 16.43% ################### face 2: 9812 16.35% ################### face 3: 10165 16.94% #################### face 4: 10231 17.05% #################### face 5: 10017 16.69% #################### face 6: 9914 16.52% ################### mean of the rolls: 3.5079
Monte Carlo estimation of pi
A point chosen uniformly in the unit square falls inside the quarter circle of
radius one with probability pi/4. Counting the points inside and
multiplying the fraction by four therefore estimates pi, and the estimate
improves as the square root of the number of points. This is the simplest
example of a Monte Carlo method: replace an integral by the average of
a function over random samples.
program r4
use iso_fortran_env, only: real64
implicit none
integer, parameter :: points = 1000000
integer :: inside, i
real(real64) :: x, y, estimate, pi
call seed_with(31415)
pi = acos(-1.0_real64)
inside = 0
do i = 1, points
call random_number(x)
call random_number(y)
if (x * x + y * y <= 1.0_real64) inside = inside + 1
select case (i)
case (1000, 10000, 100000, 1000000)
estimate = 4.0_real64 * inside / i
print '(a, i8, a, f9.6, a, f9.6)', 'after ', i, &
' points: ', estimate, ', error ', abs(estimate - pi)
end select
end do
print '(a, f10.7)', 'the value of pi: ', pi
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r4
The estimate after a thousand points is 3.168, after ten thousand 3.1408, and after a million 3.139604, an error of 0.002. The checkpoints do not improve monotonically, because each one is a different prefix of the same sequence and the error of a random estimate fluctuates; what falls is the typical size of the error, roughly by a factor of three for every tenfold increase in the sample. One million points buy about three correct decimals, which is why Monte Carlo is used for integrals in many dimensions, where the error does not depend on the dimension at all.
$ gfortran r4.f90 -o r4 $ ./r4 after 1000 points: 3.168000, error 0.026407 after 10000 points: 3.140800, error 0.000793 after 100000 points: 3.144000, error 0.002407 after 1000000 points: 3.139604, error 0.001989 the value of pi: 3.1415927
Shuffling without bias
A shuffle must produce each of the n! orderings with the same
probability. The Fisher-Yates algorithm achieves that with one pass:
for every position from the end down to the second, choose a position from the
part of the array not yet fixed and swap the two. Choosing the second position
from the whole array instead of the remaining part would bias the result
towards some orderings, which is the mistake the algorithm avoids. The program
shuffles ten cards and then checks that every value still appears exactly once.
program r5
implicit none
integer, parameter :: n = 10
integer :: deck(n), i, j, swap
logical :: seen(n), complete
real :: u
call seed_with(99)
deck = [(i, i = 1, n)]
do i = n, 2, -1
call random_number(u)
j = int(u * real(i)) + 1
swap = deck(i)
deck(i) = deck(j)
deck(j) = swap
end do
print '(a, 10(i3))', 'shuffled deck: ', deck
seen = .false.
complete = .true.
do i = 1, n
if (seen(deck(i))) complete = .false.
seen(deck(i)) = .true.
end do
print '(a, l1)', 'every card appears once: ', complete
print '(a, i0)', 'cards in their old place: ', count(deck == [(i, i = 1, n)])
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r5
The deck comes out as 7, 1, 6, 8, 9, 4, 3, 10, 2, 5, and the check confirms that no card was lost or duplicated. The count of cards that stayed in their original position is zero, which is a little unusual for ten cards, where one fixed point is the most likely outcome, but entirely possible. The same algorithm shuffles a playlist, deals a card game or randomises the order of experiments, and it runs in time proportional to the number of elements.
$ gfortran r5.f90 -o r5 $ ./r5 shuffled deck: 7 1 6 8 9 4 3 10 2 5 every card appears once: T cards in their old place: 0
Normally distributed numbers
A uniform value contains no more information than its position in the interval, and many natural quantities need the bell-shaped normal distribution instead. The Box-Muller transform turns two independent uniform values into two independent normal ones: take the square root of minus twice the logarithm of the first, and multiply it by the cosine of two pi times the second. The result has mean zero and standard deviation one, so about 68 percent of the values lie within one unit of zero, 95 percent within two and 99.7 percent within three.
program r6
use iso_fortran_env, only: real64
implicit none
integer, parameter :: samples = 100000
integer :: i, within1, within2, within3
real(real64) :: u1, u2, z, total, total2, mean, sd, pi
pi = acos(-1.0_real64)
call seed_with(2718)
total = 0.0_real64
total2 = 0.0_real64
within1 = 0
within2 = 0
within3 = 0
do i = 1, samples
call random_number(u1)
call random_number(u2)
z = sqrt(-2.0_real64 * log(u1)) * cos(2.0_real64 * pi * u2)
total = total + z
total2 = total2 + z * z
if (abs(z) <= 1.0_real64) within1 = within1 + 1
if (abs(z) <= 2.0_real64) within2 = within2 + 1
if (abs(z) <= 3.0_real64) within3 = within3 + 1
end do
mean = total / samples
sd = sqrt(total2 / samples - mean**2)
print '(a, f9.5)', 'mean: ', mean
print '(a, f9.5)', 'standard dev: ', sd
print '(a, f7.4)', 'within one sigma:', real(within1) / samples
print '(a, f7.4)', 'within two: ', real(within2) / samples
print '(a, f7.4)', 'within three: ', real(within3) / samples
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r6
A hundred thousand values give a mean of -0.003 and a standard deviation of 0.99982, and the three fractions are 0.6834, 0.9549 and 0.9972, each within a thousandth of the theoretical value. The transformation needs the logarithm, which is why the first uniform value must not be zero; the generator returns values from zero inclusive, so a program that wants to be strictly safe rejects a zero and draws again, even though a value of exactly zero is extremely unlikely.
$ gfortran r6.f90 -o r6 $ ./r6 mean: -0.00309 standard dev: 0.99982 within one sigma: 0.6834 within two: 0.9549 within three: 0.9972
The birthday problem
The birthday problem asks how many people are needed before two of them share
a birthday with probability greater than one half. The answer, 23, surprises
most people, because they compare each person with one chosen person instead of
comparing every pair. A simulation draws the birthdays of a group, stops as
soon as a date repeats, and repeats the experiment ten thousand times; the
fraction of trials with a match estimates the probability, and the exact value
follows from multiplying the probabilities that the first k
birthdays are all different.
program r7
use iso_fortran_env, only: real64
implicit none
integer, parameter :: trials = 10000, days = 365
integer :: sizes(5), size_index, trial, person, day, hits
logical :: seen(days), shared
real :: u
real(real64) :: exact
sizes = [5, 10, 23, 30, 50]
call seed_with(2026)
print '(a)', 'people simulated exact'
do size_index = 1, size(sizes)
hits = 0
do trial = 1, trials
seen = .false.
shared = .false.
do person = 1, sizes(size_index)
call random_number(u)
day = int(u * real(days)) + 1
if (seen(day)) then
shared = .true.
exit
end if
seen(day) = .true.
end do
if (shared) hits = hits + 1
end do
exact = 1.0_real64
do person = 0, sizes(size_index) - 1
exact = exact * real(days - person, real64) / real(days, real64)
end do
print '(i6, f11.4, f8.4)', sizes(size_index), &
real(hits) / trials, 1.0_real64 - exact
end do
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program r7
Groups of 5, 10, 23, 30 and 50 people share a birthday in 2.57, 11.65, 51.03, 70.01 and 97.03 percent of the trials, against exact values of 2.71, 11.69, 50.73, 70.63 and 97.04 percent. The simulation and the formula agree to about a percentage point at the small sizes, where ten thousand trials are not enough to resolve a probability of a few percent, and to a tenth of a point at 50. The same loop answers the question for any number of days, which makes it a cheap way to check a probability argument before trusting it.
$ gfortran r7.f90 -o r7
$ ./r7
people simulated exact
5 0.0257 0.0271
10 0.1165 0.1169
23 0.5103 0.5073
30 0.7001 0.7063
50 0.9703 0.9704
Next steps
The loops article explains the counted and
conditional loops that drive every simulation here, and the
arrays article covers the array arguments of
random_number and the swaps of the shuffle. The
floating-point article describes the
precision of the estimates and the logarithm that the Box-Muller transform
needs, the numeric
functions article covers sqrt, log and
acos, and the big integers
article counts the exact factorials that a simulation can only estimate.
Source
Fortran 2018 standard (J3/18-007r1)
Fortran-lang: learning resources
In this article we have drawn uniform values and checked them with a mean and a
histogram, seeded the generator and restored a saved state, rolled dice with
int(u*6) + 1, estimated pi by Monte Carlo sampling, shuffled an
array with the Fisher-Yates algorithm and verified the permutation, generated
normal values with the Box-Muller transform and checked the sigma fractions,
and simulated the birthday problem against its exact formula. The previous
article, Fortran big integers, computes
exactly the numbers that a simulation can only sample.
Author
List all Fortran tutorials.