Fortran complex numbers
last modified November 14, 2026
A complex number is a pair of real numbers written as
a + bi, where i is the imaginary unit with
the property i**2 = -1. The first number, a, is the
real part and the second, b, the imaginary part.
Complex numbers appear wherever something has an amplitude and a phase: an
alternating current, a wave, a rotation in the plane, the coefficients of a
polynomial with no real roots.
Fortran has had a complex type since the first standard, so a
complex value is a built-in type rather than a structure of two reals. This
article constructs complex values, does arithmetic with them, converts them to
and from polar form, calls the complex versions of the mathematical functions,
works with arrays of them, and uses them twice: to solve quadratic equations
and to draw the Mandelbrot set. The last program computes a discrete Fourier
transform, which is what complex arithmetic is famous for.
Constructing and printing complex values
A complex literal is written as two real values in parentheses, separated by a
comma, so (3.0, 4.0) is 3 + 4i. The intrinsic
cmplx(a, b) builds the same value from two expressions of any
numeric type, and an optional third argument selects the kind. The functions
real, aimag, conjg and abs
take a complex number apart: the real part, the imaginary part, the
conjugate (the imaginary part with the opposite sign) and the
modulus (the distance from zero).
program c1
implicit none
complex :: z, w
complex :: from_int
z = (3.0, 4.0)
w = cmplx(2.0, -1.0)
call show('z ', z)
call show('w ', w)
call show('z + w ', z + w)
call show('z - w ', z - w)
call show('z * w ', z * w)
call show('z / w ', z / w)
call show('conjg(z) ', conjg(z))
call show('z**2 ', z**2)
print '(a, f7.2)', 'abs(z) = ', abs(z)
print '(a, f7.2)', 'abs(w) = ', abs(w)
from_int = cmplx(2, 3)
call show('cmplx(2,3)', from_int)
contains
subroutine show(label, value)
character(len=*), intent(in) :: label
complex, intent(in) :: value
character :: sign
if (aimag(value) < 0.0) then
sign = '-'
else
sign = '+'
end if
print '(a, a, f5.1, a, a, f5.1, a)', label, ' = (', real(value), &
' ', sign, abs(aimag(value)), 'i)'
end subroutine show
end program c1
The program prints each value in the familiar (a + bi) form with a
small internal subroutine that chooses the sign, because a formatted output of
the imaginary part alone would show + -1.0i. Addition and
subtraction work component by component; multiplication follows
(a + bi)(c + di) = (ac - bd) + (ad + bc)i, which is why
z * w has a real part of 10 and an imaginary part of 5. Division
multiplies both sides by the conjugate of the divisor, so z / w
is 0.4 + 2.2i, and z**2 is -7 + 24i because the cross terms add.
$ gfortran c1.f90 -o c1 $ ./c1 z = ( 3.0 + 4.0i) w = ( 2.0 - 1.0i) z + w = ( 5.0 + 3.0i) z - w = ( 1.0 + 5.0i) z * w = ( 10.0 + 5.0i) z / w = ( 0.4 + 2.2i) conjg(z) = ( 3.0 - 4.0i) z**2 = ( -7.0 + 24.0i) abs(z) = 5.00 abs(w) = 2.24 cmplx(2,3) = ( 2.0 + 3.0i)
Polar form and Euler's formula
A complex number is a point in the plane, and a point can also be described by
its distance from the origin and its angle. The distance is the
modulus r = abs(z), and the angle is the
argument, computed with atan2(aimag(z), real(z)). In that
form z = r * (cos(theta) + i sin(theta)), and Euler's formula
replaces the bracket with exp(i * theta). Multiplying by
exp(i * theta) rotates a point by theta without
changing its distance, which is the property that makes complex numbers the
natural language of waves and oscillations.
program c2
use iso_fortran_env, only: real64
implicit none
real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64)
complex(real64) :: z, rebuilt, euler, root, point
real(real64) :: r, theta
integer :: k
z = cmplx(3.0_real64, 4.0_real64, real64)
r = abs(z)
theta = atan2(aimag(z), real(z))
print '(a, f9.4)', 'modulus: ', r
print '(a, f9.4, a)', 'argument: ', theta, ' rad'
print '(a, f9.4, a)', ' ', theta * 180.0_real64 / pi, ' deg'
rebuilt = r * cmplx(cos(theta), sin(theta), real64)
call show('rebuilt from polar', rebuilt)
euler = exp(cmplx(0.0_real64, pi, real64))
call show('exp(i*pi) ', euler)
call show('exp(i*pi/2) ', exp(cmplx(0.0_real64, pi / 2, real64)))
call show('sqrt(-1) ', sqrt(cmplx(-1.0_real64, 0.0_real64, real64)))
print *, 'fourth roots of unity:'
do k = 0, 3
root = exp(cmplx(0.0_real64, 2.0_real64 * pi * k / 4.0_real64, real64))
call show(' root', root)
end do
point = (1.0_real64, 0.0_real64)
call show('point ', point)
call show('point * i ', point * (0.0_real64, 1.0_real64))
call show('point * i * i ', point * (0.0_real64, 1.0_real64)**2)
contains
subroutine show(label, value)
character(len=*), intent(in) :: label
complex(real64), intent(in) :: value
print '(a, a, f9.6, a, f9.6, a)', label, ' = (', real(value), &
', ', aimag(value), ')'
end subroutine show
end program c2
The point 3 + 4i lies at distance 5 and at an angle of 0.9273 radians, about
53.13 degrees, and the polar expression rebuilds it to six decimals. Setting
the angle to pi gives -1 with an imaginary part of 1.2e-16, which prints as
0.000000: the exact answer is real, and the residue is the rounding error of
pi and of sin. The fourth roots of unity are 1,
i, -1 and -i, and multiplying the point (1, 0) by i twice turns it into (-1, 0).
$ gfortran c2.f90 -o c2
$ ./c2
modulus: 5.0000
argument: 0.9273 rad
53.1301 deg
rebuilt from polar = ( 3.000000, 4.000000)
exp(i*pi) = (-1.000000, 0.000000)
exp(i*pi/2) = ( 0.000000, 1.000000)
sqrt(-1) = ( 0.000000, 1.000000)
fourth roots of unity:
root = ( 1.000000, 0.000000)
root = ( 0.000000, 1.000000)
root = (-1.000000, 0.000000)
root = (-0.000000, -1.000000)
point = ( 1.000000, 0.000000)
point * i = ( 0.000000, 1.000000)
point * i * i = (-1.000000, 0.000000)
Complex versions of the mathematical functions
Most intrinsic functions accept a complex argument and return the principal
value of the complex function. sqrt returns the root with a
non-negative real part, log the branch whose imaginary part lies
between -pi and pi, and exp,
sin, cos and their relatives follow the same
identities as their real counterparts. The results below are checked against
the closed forms, written with real functions, so the agreement is visible in
the printed digits.
program c3
use iso_fortran_env, only: real64
implicit none
real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64)
complex(real64) :: z, checked
call show('log(-1) ', log(cmplx(-1.0_real64, 0.0_real64, real64)))
call show('log(1+i) ', log(cmplx(1.0_real64, 1.0_real64, real64)))
print '(a, f9.6)', ' 0.5*log(2) = ', 0.5_real64 * log(2.0_real64)
print '(a, f9.6)', ' pi/4 = ', pi / 4.0_real64
call show('sqrt(3+4i) ', sqrt(cmplx(3.0_real64, 4.0_real64, real64)))
call show('sqrt(-4) ', sqrt(cmplx(-4.0_real64, 0.0_real64, real64)))
call show('exp(1+i) ', exp(cmplx(1.0_real64, 1.0_real64, real64)))
print '(a, f9.6)', ' e*cos(1) = ', exp(1.0_real64) * cos(1.0_real64)
print '(a, f9.6)', ' e*sin(1) = ', exp(1.0_real64) * sin(1.0_real64)
z = cmplx(1.0_real64, 1.0_real64, real64)
call show('sin(1+i) ', sin(z))
print '(a, f9.6)', ' sin1*cosh1 = ', sin(1.0_real64) * cosh(1.0_real64)
print '(a, f9.6)', ' cos1*sinh1 = ', cos(1.0_real64) * sinh(1.0_real64)
call show('cos(i) ', cos(cmplx(0.0_real64, 1.0_real64, real64)))
print '(a, f9.6)', ' cosh(1) = ', cosh(1.0_real64)
checked = sin(z)**2 + cos(z)**2
call show('sin^2+cos^2 ', checked)
print '(a, es10.2)', 'distance from 1: ', abs(checked - 1.0_real64)
contains
subroutine show(label, value)
character(len=*), intent(in) :: label
complex(real64), intent(in) :: value
print '(a, a, f10.6, a, f10.6, a)', label, ' = (', real(value), &
', ', aimag(value), ')'
end subroutine show
end program c3
log(-1) is i*pi, because the logarithm of a complex
number is log(r) + i*theta and the negative real axis has
r = 1 and theta = pi. log(1+i) has real
part 0.5*log(2) and angle pi/4, and both match the
printed constants. sqrt(3+4i) is 2 + i, whose square is indeed
3 + 4i, and sin(1+i) is sin(1)cosh(1) +
i cos(1)sinh(1). The identity sin**2 + cos**2 = 1 holds for
complex arguments too, with an error of 2.2e-16.
$ gfortran c3.f90 -o c3 $ ./c3 log(-1) = ( 0.000000, 3.141593) log(1+i) = ( 0.346574, 0.785398) 0.5*log(2) = 0.346574 pi/4 = 0.785398 sqrt(3+4i) = ( 2.000000, 1.000000) sqrt(-4) = ( 0.000000, 2.000000) exp(1+i) = ( 1.468694, 2.287355) e*cos(1) = 1.468694 e*sin(1) = 2.287355 sin(1+i) = ( 1.298458, 0.634964) sin1*cosh1 = 1.298458 cos1*sinh1 = 0.634964 cos(i) = ( 1.543081, -0.000000) cosh(1) = 1.543081 sin^2+cos^2 = ( 1.000000, 0.000000) distance from 1: 2.22E-16
Arrays of complex values and the conjugate
An array of complex numbers is an ordinary array with a complex element type,
so whole-array operations apply: sum, product,
abs, maxloc and the rest work as they do for reals.
The one trap is the inner product. In mathematics the inner product of two
complex vectors conjugates the first one, and Fortran's
dot_product does exactly that, while a plain
sum(a * b) does not. The two differ as soon as the vectors have
imaginary parts.
program c4
use iso_fortran_env, only: real64
implicit none
complex(real64) :: a(5), b(5), mean
complex(real64) :: conjugated, plain
character(len=12) :: label
integer :: i, biggest
a = cmplx([1.0_real64, 2.0_real64, -1.0_real64, 0.5_real64, 3.0_real64], &
[1.0_real64, -1.0_real64, 0.5_real64, 0.0_real64, 2.0_real64], &
real64)
b = cmplx([0.5_real64, 1.0_real64, 2.0_real64, 0.0_real64, 1.0_real64], &
[0.0_real64, 1.0_real64, -2.0_real64, 1.0_real64, -1.0_real64], &
real64)
print *, 'a:'
do i = 1, size(a)
write (label, '(a, i0, a)') ' a(', i, ')'
call show(label, a(i))
end do
print '(a, f8.3)', 'sum of moduli: ', sum(abs(a))
print '(a, i0)', 'largest modulus at index ', maxloc(abs(a))
biggest = maxloc(abs(a), 1)
call show('largest element', a(biggest))
mean = sum(a) / size(a)
call show('mean ', mean)
call show('sum(a) ', sum(a))
call show('product(a) ', product(a))
conjugated = sum(conjg(a) * b)
plain = sum(a * b)
call show('sum(conjg(a)*b) ', conjugated)
call show('sum(a*b) ', plain)
call show('dot_product(a,b) ', dot_product(a, b))
contains
subroutine show(label, value)
character(len=*), intent(in) :: label
complex(real64), intent(in) :: value
print '(a, a, f10.6, a, f10.6, a)', label, ' = (', real(value), &
', ', aimag(value), ')'
end subroutine show
end program c4
The constructors build both arrays from two real arrays, one real part and one
imaginary part per element. The sum of the moduli is 8.874, the largest modulus
belongs to the last element, and the mean is the sum divided by the number of
elements. The last three lines are the instructive part:
sum(conjg(a)*b) and dot_product(a, b) agree at
-0.5 - i, while sum(a*b) is 7.5 + 4i. When a program computes an
inner product by hand, it must conjugate the first factor to match the
intrinsic.
$ gfortran c4.f90 -o c4 $ ./c4 a: a(1) = ( 1.000000, 1.000000) a(2) = ( 2.000000, -1.000000) a(3) = ( -1.000000, 0.500000) a(4) = ( 0.500000, 0.000000) a(5) = ( 3.000000, 2.000000) sum of moduli: 8.874 largest modulus at index 5 largest element = ( 3.000000, 2.000000) mean = ( 1.100000, 0.500000) sum(a) = ( 5.500000, 2.500000) product(a) = ( -5.750000, -2.750000) sum(conjg(a)*b) = ( -0.500000, -1.000000) sum(a*b) = ( 7.500000, 4.000000) dot_product(a,b) = ( -0.500000, -1.000000)
Quadratic equations
The quadratic formula x = (-b +- sqrt(b**2 - 4ac)) / (2a) has no
real answer when the discriminant is negative, and a program that takes a real
square root of a negative value produces a NaN or stops. Complex arithmetic
removes the special case: the discriminant is converted with
cmplx, its complex square root is taken, and the same two
expressions cover real, repeated and complex roots. The roots can then be put
back into the equation to check them.
module quadratics
use iso_fortran_env, only: real64
implicit none
private
public :: solve
contains
subroutine solve(a, b, c, z1, z2, repeated)
real(real64), intent(in) :: a, b, c
complex(real64), intent(out) :: z1, z2
logical, intent(out) :: repeated
real(real64) :: discriminant
complex(real64) :: root
discriminant = b * b - 4.0_real64 * a * c
root = sqrt(cmplx(discriminant, 0.0_real64, real64))
z1 = cmplx(-b, 0.0_real64, real64) + root
z1 = z1 / (2.0_real64 * a)
z2 = cmplx(-b, 0.0_real64, real64) - root
z2 = z2 / (2.0_real64 * a)
repeated = abs(discriminant) < 1.0e-12_real64 * max(1.0_real64, abs(b))
end subroutine solve
end module quadratics
program c5
use iso_fortran_env, only: real64
use quadratics
implicit none
real(real64) :: cases(3, 3)
complex(real64) :: z1, z2, check
logical :: repeated
integer :: i
cases(1, :) = [1.0_real64, -6.0_real64, 13.0_real64]
cases(2, :) = [1.0_real64, 4.0_real64, 4.0_real64]
cases(3, :) = [2.0_real64, 3.0_real64, -2.0_real64]
do i = 1, 3
call solve(cases(i, 1), cases(i, 2), cases(i, 3), z1, z2, repeated)
print '(a, f5.1, a, f5.1, a, f5.1, a)', 'x^2 coefficient ', &
cases(i, 1), ', x coefficient ', cases(i, 2), &
', constant ', cases(i, 3)
call show(' first root ', z1)
call show(' second root', z2)
check = cases(i, 1) * z1**2 + cases(i, 2) * z1 + cases(i, 3)
print '(a, es10.2, a, l1)', ' residual ', abs(check), &
', repeated roots: ', repeated
end do
contains
subroutine show(label, value)
character(len=*), intent(in) :: label
complex(real64), intent(in) :: value
print '(a, a, f9.4, a, f9.4, a)', label, ' = (', real(value), &
', ', aimag(value), ')'
end subroutine show
end program c5
The first equation has a discriminant of -16 and the roots 3 +- 2i, which are conjugates because the coefficients are real. The second has a zero discriminant and both roots are -2, which the tolerance test reports. The third has real roots, 0.5 and -2, and the complex type represents them with a zero imaginary part. The residual of every root is exactly zero at this precision, which is the expected result for these small integers.
$ gfortran c5.f90 -o c5 $ ./c5 x^2 coefficient 1.0, x coefficient -6.0, constant 13.0 first root = ( 3.0000, 2.0000) second root = ( 3.0000, -2.0000) residual 0.00E+00, repeated roots: F x^2 coefficient 1.0, x coefficient 4.0, constant 4.0 first root = ( -2.0000, 0.0000) second root = ( -2.0000, 0.0000) residual 0.00E+00, repeated roots: T x^2 coefficient 2.0, x coefficient 3.0, constant -2.0 first root = ( 0.5000, 0.0000) second root = ( -2.0000, 0.0000) residual 0.00E+00, repeated roots: F
The Mandelbrot set
The Mandelbrot set is the set of complex numbers c for which the
sequence z(0) = 0, z(n+1) = z(n)**2 + c stays
bounded. The test stops when abs(z) exceeds 2, because a value
outside that circle is guaranteed to escape; the number of iterations before
the escape measures how close the point came to the set. The program below
walks a grid of points, maps the iteration count to a character, and prints the
result as text art.
program c6
implicit none
integer, parameter :: width = 60, height = 26, max_iter = 100
character(len=*), parameter :: ramp = ' .:-=+*#%@'
real :: x0, y0, x, y, next_x, dx, dy
integer :: i, j, iter, band
dx = 3.4 / real(width - 1)
dy = 2.4 / real(height - 1)
do j = 1, height
y0 = -1.2 + real(j - 1) * dy
do i = 1, width
x0 = -2.2 + real(i - 1) * dx
x = 0.0
y = 0.0
iter = 0
do while (x * x + y * y <= 4.0 .and. iter < max_iter)
next_x = x * x - y * y + x0
y = 2.0 * x * y + y0
x = next_x
iter = iter + 1
end do
band = int(sqrt(real(iter) / real(max_iter)) * &
real(len(ramp) - 1)) + 1
write (*, '(a)', advance='no') ramp(band:band)
end do
print *
end do
print '(a, i0, a, i0, a, i0)', 'grid ', width, ' by ', height, &
', up to ', max_iter
end program c6
Every character is one complex point, so the picture is a 60 by 26 sampling of the plane between -2.2 and 1.2 on the real axis and -1.2 and 1.2 on the imaginary axis. The square root in the band calculation brightens the points that escape quickly, which would otherwise all fall into the first two characters. The cardioid and the round bulb on its left are the two largest regions of the set, and the fine filaments between them are aliasing: a character grid cannot resolve features smaller than its spacing, but the shape survives.
$ gfortran c6.f90 -o c6
$ ./c6
.................................................
........................::::......................
.......................::::-:-:....................
.......................:::::--:::...................
......................:::::=@@@-:::::................
....................::::::::+@@@*:::::::..............
................::::::=+@-@@@@@@@@+-:=:::.............
.............::::::::::-@@@@@@@@@@@@@@@=:::............
..........::::::::::::#@@@@@@@@@@@@@@@@=:::............
.........:::::::-:::::@@@@@@@@@@@@@@@@@@@@=:............
.......::::::-@+@@@=--@@@@@@@@@@@@@@@@@@@=::............
.....::::::::@@@@@@@@=@@@@@@@@@@@@@@@@@@@+::............
::::::::::-@@@@@@@@@@@@@@@@@@@@@@@@@@@@@%:::............
::::::::::-@@@@@@@@@@@@@@@@@@@@@@@@@@@@@%:::............
.....::::::::@@@@@@@@=@@@@@@@@@@@@@@@@@@@+::............
.......::::::-@+@@@=--@@@@@@@@@@@@@@@@@@@=::............
.........:::::::-:::::@@@@@@@@@@@@@@@@@@@@=:............
..........::::::::::::#@@@@@@@@@@@@@@@@=:::............
.............::::::::::-@@@@@@@@@@@@@@@=:::............
................::::::=+%-@@@@@@@@+-:=:::.............
....................::::::::+@@@*:::::::..............
......................:::::=@@@-:::::................
.......................:::::--:::...................
.......................::::-:-:....................
........................::::......................
.................................................
grid 60 by 26, up to 100
A discrete Fourier transform
The discrete Fourier transform turns a sequence of n
samples into n complex amplitudes, one per frequency. The formula
is X(k) = sum(x(j) * exp(-2*pi*i*j*k/n)), and every term of the sum
is a rotation of the sample by an angle proportional to its position. The
result at index k says how much of the frequency k
cycles per record the signal contains: a real amplitude of A at a
frequency present in the signal appears as A*n/2 at the matching
index and at its mirror image.
program c7
use iso_fortran_env, only: real64
implicit none
real(real64), parameter :: pi = 4.0_real64 * atan(1.0_real64)
integer, parameter :: n = 8
complex(real64) :: samples(n), spectrum(n), turn
real(real64) :: energy_time, energy_freq
integer :: i, k
do i = 1, n
samples(i) = cos(2.0_real64 * pi * real(i - 1) / n) + &
0.5_real64 * sin(4.0_real64 * pi * real(i - 1) / n)
end do
do k = 1, n
spectrum(k) = (0.0_real64, 0.0_real64)
do i = 1, n
turn = exp(cmplx(0.0_real64, -2.0_real64 * pi * &
real((k - 1) * (i - 1)) / n, real64))
spectrum(k) = spectrum(k) + samples(i) * turn
end do
end do
print '(a)', ' k real imag magnitude'
do k = 1, n
print '(i3, 3f12.6)', k - 1, real(spectrum(k)), &
aimag(spectrum(k)), abs(spectrum(k))
end do
energy_time = sum(real(samples, real64)**2)
energy_freq = sum(abs(spectrum)**2) / n
print '(a, f9.4)', 'energy in samples: ', energy_time
print '(a, f9.4)', 'energy in spectrum: ', energy_freq
end program c7
The signal is a cosine of one cycle per record plus half a sine of two cycles,
sampled at eight points. The transform has magnitude 4 at k = 1,
2 at k = 2, and the same magnitudes again at k = 7
and k = 6, which are the mirror images of the two frequencies. The
phases sit in the imaginary parts: the cosine gives a real coefficient and the
sine an imaginary one, with the opposite sign at the mirror index. The energy
check at the end is Parseval's theorem, sum(abs(x)**2) =
sum(abs(X)**2)/n, and both sides give 5.
$ gfortran c7.f90 -o c7 $ ./c7 k real imag magnitude 0 -0.000000 0.000000 0.000000 1 4.000000 -0.000000 4.000000 2 -0.000000 -2.000000 2.000000 3 0.000000 -0.000000 0.000000 4 0.000000 -0.000000 0.000000 5 0.000000 -0.000000 0.000000 6 -0.000000 2.000000 2.000000 7 4.000000 0.000000 4.000000 energy in samples: 5.0000 energy in spectrum: 5.0000
Next steps
The data types article explains the kind
system that complex(real64) uses, and the
intrinsic numeric functions article
covers the real versions of the functions called here. The
arrays article describes the whole-array
operations used on complex arrays, and the
operators article gives the precedence rules
that decide how z**2 + w is grouped. The
formatted output article explains the
f and es descriptors that print the components.
Source
Fortran 2018 standard (J3/18-007r1)
Fortran-lang: learning resources
In this article we have built complex values with literals and
cmplx, printed their components, done the four arithmetic
operations, converted between rectangular and polar form, used Euler's formula
and the fourth roots of unity, called the complex versions of
sqrt, log, exp and sin,
seen how dot_product conjugates its first argument, solved
quadratic equations with complex roots, drawn the Mandelbrot set as text, and
computed a discrete Fourier transform with a Parseval check. The previous
article, Fortran procedure arguments
and callbacks, shows how a routine receives another routine.
Author
List all Fortran tutorials.