Modern Computer Arithmetic

Richard P. Brent and Paul Zimmermann

Version 0.5.9 of 7 October 2010

 iii

Copyright c 2003-2010 Richard P. Brent and Paul Zimmermann

This electronic version is distributed under the terms and conditions of the
Creative Commons license “Attribution-Noncommercial-No Derivative Works
3.0”. You are free to copy, distribute and transmit this book under the following
conditions:
• Attribution. You must attribute the work in the manner specified by the
author or licensor (but not in any way that suggests that they endorse you or
your use of the work).
• Noncommercial. You may not use this work for commercial purposes.
• No Derivative Works. You may not alter, transform, or build upon this
work.
For any reuse or distribution, you must make clear to others the license terms
of this work. The best way to do this is with a link to the web page below. Any
of the above conditions can be waived if you get permission from the copyright
holder. Nothing in this license impairs or restricts the author’s moral rights.

For more information about the license, visit

http://creativecommons.org/licenses/by-nc-nd/3.0/
 Contents

Contents page iv
Preface ix
Acknowledgements xi
Notation xiii
1 Integer arithmetic 1
1.1 Representation and notations 1
1.2 Addition and subtraction 2
1.3 Multiplication 3
1.3.1 Naive multiplication 4
1.3.2 Karatsuba’s algorithm 5
1.3.3 Toom–Cook multiplication 6
1.3.4 Use of the fast Fourier transform (FFT) 8
1.3.5 Unbalanced multiplication 8
1.3.6 Squaring 11
1.3.7 Multiplication by a constant 13
1.4 Division 14
1.4.1 Naive division 14
1.4.2 Divisor preconditioning 16
1.4.3 Divide and conquer division 18
1.4.4 Newton’s method 21
1.4.5 Exact division 21
1.4.6 Only quotient or remainder wanted 22
1.4.7 Division by a single word 23
1.4.8 Hensel’s division 24
1.5 Roots 25
1.5.1 Square root 25
1.5.2 kth root 27
 Contents v

1.5.3 Exact root 28

1.6 Greatest common divisor 29
1.6.1 Naive GCD 29
1.6.2 Extended GCD 32
1.6.3 Half binary GCD, divide and conquer GCD 33
1.7 Base conversion 37
1.7.1 Quadratic algorithms 37
1.7.2 Subquadratic algorithms 38
1.8 Exercises 39
1.9 Notes and references 44
2 Modular arithmetic and the FFT 47
2.1 Representation 47
2.1.1 Classical representation 47
2.1.2 Montgomery’s form 48
2.1.3 Residue number systems 48
2.1.4 MSB vs LSB algorithms 49
2.1.5 Link with polynomials 49
2.2 Modular addition and subtraction 50
2.3 The Fourier transform 50
2.3.1 Theoretical setting 50
2.3.2 The fast Fourier transform 51
2.3.3 The Schönhage–Strassen algorithm 55
2.4 Modular multiplication 58
2.4.1 Barrett’s algorithm 58
2.4.2 Montgomery’s multiplication 60
2.4.3 McLaughlin’s algorithm 63
2.4.4 Special moduli 65
2.5 Modular division and inversion 65
2.5.1 Several inversions at once 67
2.6 Modular exponentiation 68
2.6.1 Binary exponentiation 70
2.6.2 Exponentiation with a larger base 70
2.6.3 Sliding window and redundant representation 72
2.7 Chinese remainder theorem 73
2.8 Exercises 75
2.9 Notes and references 77
3 Floating-point arithmetic 79
3.1 Representation 79
3.1.1 Radix choice 80
vi Contents

3.1.2 Exponent range 81

3.1.3 Special values 82
3.1.4 Subnormal numbers 82
3.1.5 Encoding 83
3.1.6 Precision: local, global, operation, operand 84
3.1.7 Link to integers 86
3.1.8 Ziv’s algorithm and error analysis 86
3.1.9 Rounding 87
3.1.10 Strategies 90
3.2 Addition, subtraction, comparison 91
3.2.1 Floating-point addition 92
3.2.2 Floating-point subtraction 93
3.3 Multiplication 95
3.3.1 Integer multiplication via complex FFT 98
3.3.2 The middle product 99
3.4 Reciprocal and division 101
3.4.1 Reciprocal 102
3.4.2 Division 106
3.5 Square root 111
3.5.1 Reciprocal square root 112
3.6 Conversion 114
3.6.1 Floating-point output 115
3.6.2 Floating-point input 117
3.7 Exercises 118
3.8 Notes and references 120
4 Elementary and special function evaluation 125
4.1 Introduction 125
4.2 Newton’s method 126
4.2.1 Newton’s method for inverse roots 127
4.2.2 Newton’s method for reciprocals 128
4.2.3 Newton’s method for (reciprocal) square roots 129
4.2.4 Newton’s method for formal power series 129
4.2.5 Newton’s method for functional inverses 130
4.2.6 Higher-order Newton-like methods 131
4.3 Argument reduction 132
4.3.1 Repeated use of a doubling formula 134
4.3.2 Loss of precision 134
4.3.3 Guard digits 135
4.3.4 Doubling versus tripling 136
 Contents vii

4.4 Power series 136

4.4.1 Direct power series evaluation 140
4.4.2 Power series with argument reduction 140
4.4.3 Rectangular series splitting 141
4.5 Asymptotic expansions 144
4.6 Continued fractions 150
4.7 Recurrence relations 152
4.7.1 Evaluation of Bessel functions 153
4.7.2 Evaluation of Bernoulli and tangent numbers 154
4.8 Arithmetic-geometric mean 158
4.8.1 Elliptic integrals 158
4.8.2 First AGM algorithm for the logarithm 159
4.8.3 Theta functions 160
4.8.4 Second AGM algorithm for the logarithm 162
4.8.5 The complex AGM 163
4.9 Binary splitting 163
4.9.1 A binary splitting algorithm for sin, cos 166
4.9.2 The bit-burst algorithm 167
4.10 Contour integration 169
4.11 Exercises 171
4.12 Notes and references 179
5 Implementations and pointers 185
5.1 Software tools 185
5.1.1 CLN 185
5.1.2 GNU MP (GMP) 185
5.1.3 MPFQ 186
5.1.4 GNU MPFR 187
5.1.5 Other multiple-precision packages 187
5.1.6 Computational algebra packages 188
5.2 Mailing lists 189
5.2.1 The GMP lists 189
5.2.2 The MPFR list 190
5.3 On-line documents 190
References 191
Index 207
 Preface

This is a book about algorithms for performing arithmetic, and their imple-
mentation on modern computers. We are concerned with software more than
hardware – we do not cover computer architecture or the design of computer
hardware since good books are already available on these topics. Instead, we
focus on algorithms for efficiently performing arithmetic operations such as
addition, multiplication, and division, and their connections to topics such
as modular arithmetic, greatest common divisors, the fast Fourier transform
(FFT), and the computation of special functions.
The algorithms that we present are mainly intended for arbitrary-precision
arithmetic. That is, they are not limited by the computer wordsize of 32 or 64
bits, only by the memory and time available for the computation. We consider
both integer and real (floating-point) computations.
The book is divided into four main chapters, plus one short chapter (essen-
tially an appendix). Chapter 1 covers integer arithmetic. This has, of course,
been considered in many other books and papers. However, there has been
much recent progress, inspired in part by the application to public key cryp-
tography, so most of the published books are now partly out of date or incom-
plete. Our aim is to present the latest developments in a concise manner. At the
same time, we provide a self-contained introduction for the reader who is not
an expert in the field.
Chapter 2 is concerned with modular arithmetic and the FFT, and their appli-
cations to computer arithmetic. We consider different number representations,
fast algorithms for multiplication, division and exponentiation, and the use of
the Chinese remainder theorem (CRT).
Chapter 3 covers floating-point arithmetic. Our concern is with high-
precision floating-point arithmetic, implemented in software if the precision
provided by the hardware (typically IEEE standard 53-bit significand) is
x Preface

inadequate. The algorithms described in this chapter focus on correct round-

ing, extending the IEEE standard to arbitrary precision.
Chapter 4 deals with the computation, to arbitrary precision, of functions
such as sqrt, exp, ln, sin, cos, and more generally functions defined by power
series or continued fractions. Of course, the computation of special functions is
a huge topic so we have had to be selective. In particular, we have concentrated
on methods that are efficient and suitable for arbitrary-precision computations.
The last chapter contains pointers to implementations, useful web sites,
mailing lists, and so on. Finally, at the end there is a one-page Summary of
complexities which should be a useful aide-mémoire.
The chapters are fairly self-contained, so it is possible to read them out of
order. For example, Chapter 4 could be read before Chapters 1–3, and Chap-
ter 5 can be consulted at any time. Some topics, such as Newton’s method,
appear in different guises in several chapters. Cross-references are given where
appropriate.
For details that are omitted, we give pointers in the Notes and references
sections of each chapter, as well as in the bibliography. We have tried, as far
as possible, to keep the main text uncluttered by footnotes and references, so
most references are given in the Notes and references sections.
The book is intended for anyone interested in the design and implementation
of efficient algorithms for computer arithmetic, and more generally efficient
numerical algorithms. We did our best to present algorithms that are ready to
implement in your favorite language, while keeping a high-level description
and not getting too involved in low-level or machine-dependent details. An
alphabetical list of algorithms can be found in the index.
Although the book is not specifically intended as a textbook, it could be
used in a graduate course in mathematics or computer science, and for this
reason, as well as to cover topics that could not be discussed at length in the
text, we have included exercises at the end of each chapter. The exercises vary
considerably in difficulty, from easy to small research projects, but we have
not attempted to assign them a numerical rating. For solutions to the exercises,
please contact the authors.
We welcome comments and corrections. Please send them to either of the
authors.

Richard Brent and Paul Zimmermann

Canberra and Nancy
[email protected]
[email protected]
 Acknowledgements

We thank the French National Institute for Research in Computer Science and
Control (INRIA), the Australian National University (ANU), and the Aus-
tralian Research Council (ARC), for their support. The book could not have
been written without the contributions of many friends and colleagues, too nu-
merous to mention here, but acknowledged in the text and in the Notes and
references sections at the end of each chapter.
We also thank those who have sent us comments on and corrections to ear-
lier versions of this book: Jörg Arndt, Marco Bodrato, Wolfgang Ehrhardt
(with special thanks), Steven Galbraith, Torbjörn Granlund, Guillaume Han-
rot, Marc Mezzarobba, Jean-Michel Muller, Denis Roegel, Wolfgang Schmid,
Arnold Schönhage, Sidi Mohamed Sedjelmaci, Emmanuel Thomé, and Mark
Wezelenburg. Two anonymous reviewers provided very helpful suggestions.
Jérémie Detrey and Anne Rix helped us in the copy-editing phase.
The Mathematics Genealogy Project (http://www.genealogy.ams.
org/) and Don Knuth’s The Art of Computer Programming [142] were useful
resources for details of entries in the index.
We also thank the authors of the LATEX program, which allowed us to pro-
duce this book, the authors of the gnuplot program, and the authors of the
GNU MP library, which helped us to illustrate several algorithms with concrete
figures.
Finally, we acknowledge the contribution of Erin Brent, who first suggested
writing the book; and thank our wives, Judy-anne and Marie, for their patience
and encouragement.
 Notation

C set of complex numbers

C set of extended complex numbers C ∪ {∞}
N set of natural numbers (nonnegative integers)
N∗ set of positive integers N\{0}
Q set of rational numbers
R set of real numbers
Z set of integers
Z/nZ ring of residues modulo n
Cn set of (real or complex) functions with n continuous derivatives
in the region of interest

ℜ(z) real part of a complex number z

ℑ(z) imaginary part of a complex number z
z̄ conjugate of a complex number z
|z| Euclidean norm of a complex number z,
or absolute value of a scalar z

Bn Bernoulli numbers, n≥0 Bn z n /n! = z/(ez − 1)

Cn scaled Bernoulli numbers, Cn = B2n /(2n)! ,
Cn z 2n = (z/2)/ tanh(z/2)
Tn tangent numbers, Tn z 2n−1 /(2n − 1)! = tan z
n
Hn harmonic number j=1 1/j (0 if n ≤ 0)
n
k binomial coefficient “n choose k” = n!/(k! (n − k)!)
(0 if k < 0 or k > n)
xiv Notation

β “word” base (usually 232 or 264 ) or “radix” (floating-point)

n “precision”: number of base β digits in an integer or in a
floating-point significand, or a free variable
ε “machine precision” β 1−n /2 or (in complexity bounds)
an arbitrarily small positive constant
η smallest positive subnormal number

◦(x), ◦n (x) rounding of real number x in precision n (Definition 3.1)

ulp(x) for a floating-point number x, one unit in the last place

M (n) time to multiply n-bit integers, or polynomials of

degree n − 1, depending on the context
∼ M (n) a function f (n) such that f (n)/M (n) → 1 as n → ∞
(we sometimes lazily omit the “∼” if the meaning is clear)
M (m, n) time to multiply an m-bit integer by an n-bit integer
D(n) time to divide a 2n-bit integer by an n-bit integer,
giving quotient and remainder
D(m, n) time to divide an m-bit integer by an n-bit integer,
giving quotient and remainder

a|b a is a divisor of b, that is b = ka for some k ∈ Z

a = b mod m modular equality, m|(a − b)
q ← a div b assignment of integer quotient to q (0 ≤ a − qb < b)
r ← a mod b assignment of integer remainder to r (0 ≤ r = a − qb < b)
(a, b) greatest common divisor of a and b
a
b or (a|b) Jacobi symbol (b odd and positive)
iff if and only if
i∧j bitwise and of integers i and j,
or logical and of two Boolean expressions
i∨j bitwise or of integers i and j,
or logical or of two Boolean expressions
i⊕j bitwise exclusive-or of integers i and j
i≪k integer i multiplied by 2k
i≫k quotient of division of integer i by 2k
a · b, a × b product of scalars a, b
a∗b cyclic convolution of vectors a, b
ν(n) 2-valuation: largest k such that 2k divides n (ν(0) = ∞)
σ(e) length of the shortest addition chain to compute e
φ(n) Euler’s totient function, #{m : 0 < m ≤ n ∧ (m, n) = 1}
 Notation xv

deg(A) for a polynomial A, the degree of A

ord(A) for a power series A = j aj z j ,
ord(A) = min{j : aj = 0} (ord(0) = +∞)

exp(x) or ex exponential function

ln(x) natural logarithm
logb (x) base-b logarithm ln(x)/ ln(b)
lg(x) base-2 logarithm ln(x)/ ln(2) = log2 (x)
log(x) logarithm to any fixed base
logk (x) (log x)k

⌈x⌉ ceiling function, min{n ∈ Z : n ≥ x}

⌊x⌋ floor function, max{n ∈ Z : n ≤ x}
⌊x⌉ nearest integer function, ⌊x + 1/2⌋

sign(n) +1 if n > 0, −1 if n < 0, and 0 if n = 0

nbits(n) ⌊lg(n)⌋ + 1 if n > 0, 0 if n = 0

[a, b] closed interval {x ∈ R : a ≤ x ≤ b} (empty if a > b)

(a, b) open interval {x ∈ R : a < x < b} (empty if a ≥ b)
[a, b), (a, b] half-open intervals, a ≤ x < b, a < x ≤ b respectively

a
t
[a, b] or [a, b]t column vector
b
a b
[a, b; c, d] 2 × 2 matrix
c d
aj element of the (forward) Fourier transform of vector a
aj element of the backward Fourier transform of vector a

f (n) = O(g(n)) ∃c, n0 such that |f (n)| ≤ cg(n) for all n ≥ n0

f (n) = Ω(g(n)) ∃c > 0, n0 such that |f (n)| ≥ cg(n) for all n ≥ n0
f (n) = Θ(g(n)) f (n) = O(g(n)) and g(n) = O(f (n))
f (n) ∼ g(n) f (n)/g(n) → 1 as n → ∞
f (n) = o(g(n)) f (n)/g(n) → 0 as n → ∞
f (n) ≪ g(n) f (n) = O(g(n))
f (n) ≫ g(n) g(n) ≪ f (n)
n n
f (x) ∼ 0 aj /xj f (x) − 0 aj /xj = o(1/xn ) as x → +∞

123 456 789 123456789 (for large integers, we may use a space after
every third digit)
xvi Notation

xxx.yyyρ a number xxx.yyy written in base ρ;

for example, the decimal number 3.25 is 11.012 in binary
a c e
b+ d+ f + ··· continued fraction a/(b + c/(d + e/(f + · · · )))

a b
|A| determinant of a matrix A, e.g. = ad − bc
c d
b
PV a
f (x) dx Cauchy principal value integral, defined by a limit
if f has a singularity in (a, b)

s || t concatenation of strings s and t

⊲ <text> comment in an algorithm
end of a proof
 1
Integer arithmetic

In this chapter, our main topic is integer arithmetic. However, we

shall see that many algorithms for polynomial arithmetic are sim-
ilar to the corresponding algorithms for integer arithmetic, but
simpler due to the lack of carries in polynomial arithmetic. Con-
sider for example addition: the sum of two polynomials of degree
n always has degree at most n, whereas the sum of two n-digit in-
tegers may have n + 1 digits. Thus, we often describe algorithms
for polynomials as an aid to understanding the corresponding
algorithms for integers.

1.1 Representation and notations

We consider in this chapter algorithms working on integers. We distinguish
between the logical – or mathematical – representation of an integer, and its
physical representation on a computer. Our algorithms are intended for “large”
integers – they are not restricted to integers that can be represented in a single
computer word.
Several physical representations are possible. We consider here only the
most common one, namely a dense representation in a fixed base. Choose an
integral base β > 1. (In case of ambiguity, β will be called the internal base.)
A positive integer A is represented by the length n and the digits ai of its base
β expansion
A = an−1 β n−1 + · · · + a1 β + a0 ,
where 0 ≤ ai ≤ β − 1, and an−1 is sometimes assumed to be non-zero.
Since the base β is usually fixed in a given program, only the length n and
the integers (ai )0≤i<n need to be stored. Some common choices for β are
232 on a 32-bit computer, or 264 on a 64-bit machine; other possible choices
2 Integer arithmetic

are respectively 109 and 1019 for a decimal representation, or 253 when using
double-precision floating-point registers. Most algorithms given in this chapter
work in any base; the exceptions are explicitly mentioned.
We assume that the sign is stored separately from the absolute value. This
is known as the “sign-magnitude” representation. Zero is an important special
case; to simplify the algorithms we assume that n = 0 if A = 0, and we usually
assume that this case is treated separately.
Except when explicitly mentioned, we assume that all operations are off-line,
i.e. all inputs (resp. outputs) are completely known at the beginning (resp. end)
of the algorithm. Different models include lazy and relaxed algorithms, and
are discussed in the Notes and references (§1.9).

1.2 Addition and subtraction

As an explanatory example, here is an algorithm for integer addition. In the
algorithm, d is a carry bit.
Our algorithms are given in a language that mixes mathematical notation
and syntax similar to that found in many high-level computer languages. It
should be straightforward to translate into a language such as C. Note that
“:=” indicates a definition, and “←” indicates assignment. Line numbers are
included if we need to refer to individual lines in the description or analysis of
the algorithm.

Algorithm 1.1 IntegerAddition

n−1 n−1
Input: A = 0 ai β i , B = 0 bi β i , carry-in 0 ≤ din ≤ 1
n−1
Output: C := 0 ci β i and 0 ≤ d ≤ 1 such that A + B + din = dβ n + C
1: d ← din
2: for i from 0 to n − 1 do
3: s ← ai + bi + d
4: (d, ci ) ← (s div β, s mod β)
5: return C, d.

Let T be the number of different values taken by the data type representing
the coefficients ai , bi . (Clearly, β ≤ T , but equality does not necessarily hold,
for example β = 109 and T = 232 .) At step 3, the value of s can be as
large as 2β − 1, which is not representable if β = T . Several workarounds
are possible: either use a machine instruction that gives the possible carry of
ai + bi , or use the fact that, if a carry occurs in ai + bi , then the computed
 1.3 Multiplication 3

sum – if performed modulo T – equals t := ai + bi − T < ai ; thus, comparing

t and ai will determine if a carry occurred. A third solution is to keep a bit in
reserve, taking β ≤ T /2.
The subtraction code is very similar. Step 3 simply becomes s ← ai −bi +d,
where d ∈ {−1, 0} is the borrow of the subtraction, and −β ≤ s < β. The
other steps are unchanged, with the invariant A − B + din = dβ n + C.
We use the arithmetic complexity model, where cost is measured by the
number of machine instructions performed, or equivalently (up to a constant
factor) the time on a single processor.
Addition and subtraction of n-word integers cost O(n), which is negligible
compared to the multiplication cost. However, it is worth trying to reduce the
constant factor implicit in this O(n) cost. We shall see in §1.3 that “fast” mul-
tiplication algorithms are obtained by replacing multiplications by additions
(usually more additions than the multiplications that they replace). Thus, the
faster the additions are, the smaller will be the thresholds for changing over to
the “fast” algorithms.

1.3 Multiplication
A nice application of large integer multiplication is the Kronecker–Schönhage
trick, also called segmentation or substitution by some authors. Assume we
want to multiply two polynomials, A(x) and B(x), with non-negative integer
coefficients (see Exercise 1.1 for negative coefficients). Assume both polyno-
mials have degree less than n, and the coefficients are bounded by ρ. Now take
a power X = β k > nρ2 of the base β, and multiply the integers a = A(X) and
b = B(X) obtained by evaluating A and B at x = X. If C(x) = A(x)B(x) =
ci xi , we clearly have C(X) = ci X i . Now since the ci are bounded by
nρ < X, the coefficients ci can be retrieved by simply “reading” blocks of k
2

words in C(X). Assume for example that we want to compute

(6x5 + 6x4 + 4x3 + 9x2 + x + 3)(7x4 + x3 + 2x2 + x + 7),
with degree less than n = 6, and coefficients bounded by ρ = 9. We can take
X = 103 > nρ2 , and perform the integer multiplication
6 006 004 009 001 003 × 7 001 002 001 007
= 42 048 046 085 072 086 042 070 010 021,
from which we can read off the product
42x9 + 48x8 + 46x7 + 85x6 + 72x5 + 86x4 + 42x3 + 70x2 + 10x + 21.
4 Integer arithmetic

Conversely, suppose we want to multiply two integers a = 0≤i<n ai β i

and b = 0≤j<n bj β j . Multiply the polynomials A(x) = 0≤i<n ai xi and
B(x) = 0≤j<n bj xj , obtaining a polynomial C(x), then evaluate C(x) at
x = β to obtain ab. Note that the coefficients of C(x) may be larger than β, in
fact they may be up to about nβ 2 . For example, with a = 123, b = 456, and
β = 10, we obtain A(x) = x2 + 2x + 3, B(x) = 4x2 + 5x + 6, with product
C(x) = 4x4 + 13x3 + 28x2 + 27x + 18, and C(10) = 56088. These examples
demonstrate the analogy between operations on polynomials and integers, and
also show the limits of the analogy.
A common and very useful notation is to let M (n) denote the time to mul-
tiply n-bit integers, or polynomials of degree n − 1, depending on the context.
In the polynomial case, we assume that the cost of multiplying coefficients is
constant; this is known as the arithmetic complexity model, whereas the bit
complexity model also takes into account the cost of multiplying coefficients,
and thus their bit-size.

1.3.1 Naive multiplication

Algorithm 1.2 BasecaseMultiply

m−1 n−1
Input: A = 0 a i β i , B = 0 bj β j
m+n−1
Output: C = AB := 0 ck β k
1: C ← A · b0
2: for j from 1 to n − 1 do
3: C ← C + β j (A · bj )
4: return C.

Theorem 1.1 Algorithm BasecaseMultiply computes the product AB

correctly, and uses Θ(mn) word operations.

The multiplication by β j at step 3 is trivial with the chosen dense representa-

tion; it simply requires shifting by j words towards the most significant words.
The main operation in Algorithm BasecaseMultiply is the computation of
A · bj and its accumulation into C at step 3. Since all fast algorithms rely on
multiplication, the most important operation to optimize in multiple-precision
software is thus the multiplication of an array of m words by one word, with
accumulation of the result in another array of m + 1 words.
 1.3 Multiplication 5

We sometimes call Algorithm BasecaseMultiply schoolbook multiplication

since it is close to the “long multiplication” algorithm that used to be taught at
school.
Since multiplication with accumulation usually makes extensive use of the
pipeline, it is best to give it arrays that are as long as possible, which means
that A rather than B should be the operand of larger size (i.e. m ≥ n).

1.3.2 Karatsuba’s algorithm

Karatsuba’s algorithm is a “divide and conquer” algorithm for multiplication
of integers (or polynomials). The idea is to reduce a multiplication of length n
to three multiplications of length n/2, plus some overhead that costs O(n).
In the following, n0 ≥ 2 denotes the threshold between naive multiplica-
tion and Karatsuba’s algorithm, which is used for n0 -word and larger inputs.
The optimal “Karatsuba threshold” n0 can vary from about ten to about 100
words, depending on the processor and on the relative cost of multiplication
and addition (see Exercise 1.6).

Algorithm 1.3 KaratsubaMultiply

n−1 n−1
Input: A = 0 ai β i , B = 0 bj β j
2n−1
Output: C = AB := 0 ck β k
if n < n0 then return BasecaseMultiply(A, B)
k ← ⌈n/2⌉
(A0 , B0 ) := (A, B) mod β k , (A1 , B1 ) := (A, B) div β k
sA ← sign(A0 − A1 ), sB ← sign(B0 − B1 )
C0 ← KaratsubaMultiply(A0 , B0 )
C1 ← KaratsubaMultiply(A1 , B1 )
C2 ← KaratsubaMultiply(|A0 − A1 |, |B0 − B1 |)
return C := C0 + (C0 + C1 − sA sB C2 )β k + C1 β 2k .

Theorem 1.2 Algorithm KaratsubaMultiply computes the product AB

correctly, using K(n) = O(nα ) word multiplications, with α = lg 3 ≈ 1.585.

Proof. Since sA |A0 − A1 | = A0 − A1 and sB |B0 − B1 | = B0 − B1 , we

have sA sB |A0 − A1 ||B0 − B1 | = (A0 − A1 )(B0 − B1 ), and thus C =
A0 B0 +(A0 B1 + A1 B0 )β k + A1 B1 β 2k .
Since A0 , B0 , |A0 − A1 | and |B0 − B1 | have (at most) ⌈n/2⌉ words, and A1
and B1 have (at most) ⌊n/2⌋ words, the number K(n) of word multiplications
6 Integer arithmetic

satisfies the recurrence K(n) = n2 for n < n0 , and K(n) = 2K(⌈n/2⌉) +

K(⌊n/2⌋) for n ≥ n0 . Assume 2ℓ−1 n0 < n ≤ 2ℓ n0 with ℓ ≥ 1. Then K(n)
is the sum of three K(j) values with j ≤ 2ℓ−1 n0 , so at most 3ℓ K(j) with
j ≤ n0 . Thus, K(n) ≤ 3ℓ max(K(n0 ), (n0 − 1)2 ), which gives K(n) ≤ Cnα
with C = 31−lg(n0 ) max(K(n0 ), (n0 − 1)2 ).
Different variants of Karatsuba’s algorithm exist; the variant presented here
is known as the subtractive version. Another classical one is the additive ver-
sion, which uses A0 +A1 and B0 +B1 instead of |A0 −A1 | and |B0 −B1 |. How-
ever, the subtractive version is more convenient for integer arithmetic, since it
avoids the possible carries in A0 + A1 and B0 + B1 , which require either an
extra word in these sums, or extra additions.
The efficiency of an implementation of Karatsuba’s algorithm depends heav-
ily on memory usage. It is important to avoid allocating memory for the inter-
mediate results |A0 − A1 |, |B0 − B1 |, C0 , C1 , and C2 at each step (although
modern compilers are quite good at optimizing code and removing unneces-
sary memory references). One possible solution is to allow a large temporary
storage of m words, used both for the intermediate results and for the recur-
sive calls. It can be shown that an auxiliary space of m = 2n words – or even
m = O(log n) – is sufficient (see Exercises 1.7 and 1.8).
Since the product C2 is used only once, it may be faster to have auxiliary
routines KaratsubaAddmul and KaratsubaSubmul that accumulate their re-
sults, calling themselves recursively, together with KaratsubaMultiply (see
Exercise 1.10).
The version presented here uses ∼ 4n additions (or subtractions): 2 × (n/2)
to compute |A0 − A1 | and |B0 − B1 |, then n to add C0 and C1 , again n to
add or subtract C2 , and n to add (C0 + C1 − sA sB C2 )β k to C0 + C1 β 2k . An
improved scheme uses only ∼ 7n/2 additions (see Exercise 1.9).
When considered as algorithms on polynomials, most fast multiplication
algorithms can be viewed as evaluation/interpolation algorithms. Karatsuba’s
algorithm regards the inputs as polynomials A0 +A1 x and B0 +B1 x evaluated
at x = β k ; since their product C(x) is of degree 2, Lagrange’s interpolation
theorem says that it is sufficient to evaluate C(x) at three points. The subtrac-
tive version evaluates1 C(x) at x = 0, −1, ∞, whereas the additive version
uses x = 0, +1, ∞.

1.3.3 Toom–Cook multiplication

Karatsuba’s idea readily generalizes to what is known as Toom–Cook r-way
multiplication. Write the inputs as a0 +· · ·+ar−1 xr−1 and b0 +· · ·+br−1 xr−1 ,
1 Evaluating C(x) at ∞ means computing the product A1 B1 of the leading coefficients.
 1.3 Multiplication 7

with x = β k , and k = ⌈n/r⌉. Since their product C(x) is of degree 2r − 2,

it suffices to evaluate it at 2r − 1 distinct points to be able to recover C(x),
and in particular C(β k ). If r is chosen√optimally, Toom–Cook multiplication
of n-word numbers takes time n1+O(1/ log n) .
Most references, when describing subquadratic multiplication algorithms,
only describe Karatsuba and FFT-based algorithms. Nevertheless, the Toom–
Cook algorithm is quite interesting in practice.
Toom–Cook r-way reduces one n-word product to 2r − 1 products of about
n/r words, thus costs O(nν ) with ν = log(2r − 1)/ log r. However, the con-
stant hidden by the big-O notation depends strongly on the evaluation and
interpolation formulæ, which in turn depend on the chosen points. One possi-
bility is to take −(r − 1), . . . , −1, 0, 1, . . . , (r − 1) as evaluation points.
The case r = 2 corresponds to Karatsuba’s algorithm (§1.3.2). The case
r = 3 is known as Toom–Cook 3-way, sometimes simply called “the Toom–
Cook algorithm”. Algorithm ToomCook3 uses the evaluation points 0, 1, −1,
2, ∞, and tries to optimize the evaluation and interpolation formulæ.

Algorithm 1.4 ToomCook3

Input: two integers 0 ≤ A, B < β n
Output: AB := c0 + c1 β k + c2 β 2k + c3 β 3k + c4 β 4k with k = ⌈n/3⌉
Require: a threshold n1 ≥ 3
1: if n < n1 then return KaratsubaMultiply(A, B)
2: write A = a0 + a1 x + a2 x2 , B = b0 + b1 x + b2 x2 with x = β k .
3: v0 ← ToomCook3(a0 , b0 )
4: v1 ← ToomCook3(a02 +a1 , b02 +b1 ) where a02 ← a0 +a2 , b02 ← b0 +b2
5: v−1 ← ToomCook3(a02 − a1 , b02 − b1 )
6: v2 ← ToomCook3(a0 + 2a1 + 4a2 , b0 + 2b1 + 4b2 )
7: v∞ ← ToomCook3(a2 , b2 )
8: t1 ← (3v0 + 2v−1 + v2 )/6 − 2v∞ , t2 ← (v1 + v−1 )/2
9: c0 ← v0 , c1 ← v1 − t1 , c2 ← t2 − v0 − v∞ , c3 ← t1 − t2 , c4 ← v∞ .

The divisions at step 8 are exact; if β is a power of two, the division by 6

can be done using a division by 2 – which consists of a single shift – followed
by a division by 3 (see §1.4.7).
Toom–Cook r-way has to invert a (2r − 1) × (2r − 1) Vandermonde matrix
with parameters the evaluation points; if we choose consecutive integer points,
the determinant of that matrix contains all primes up to 2r − 2. This proves
that division by (a multiple of) 3 can not be avoided for Toom–Cook 3-way
with consecutive integer points. See Exercise 1.14 for a generalization of this
result.
8 Integer arithmetic

1.3.4 Use of the fast Fourier transform (FFT)

Most subquadratic multiplication algorithms can be seen as evaluation-inter-
polation algorithms. They mainly differ in the number of evaluation points, and
the values of those points. However, the evaluation and interpolation formulæ
become intricate in Toom–Cook r-way for large r, since they involve O(r2 )
scalar operations. The fast Fourier transform (FFT) is a way to perform evalu-
ation and interpolation efficiently for some special points (roots of unity) and
special values of r. This explains why multiplication algorithms with the best
known asymptotic complexity are based on the FFT.
There are different flavours of FFT multiplication, depending on the ring
where the operations are performed. The Schönhage–Strassen algorithm, with
a complexity of O(n log n log log n), works in the ring Z/(2n + 1)Z. Since it
is based on modular computations, we describe it in Chapter 2.
Other commonly used algorithms work with floating-point complex num-
bers. A drawback is that, due to the inexact nature of floating-point computa-
tions, a careful error analysis is required to guarantee the correctness of the im-
plementation, assuming an underlying arithmetic with rigorous error bounds.
See Theorem 3.6 in Chapter 3.
We say that multiplication is in the FFT range if n is large and the multi-
plication algorithm satisfies M (2n) ∼ 2M (n). For example, this is true if the
Schönhage–Strassen multiplication algorithm is used, but not if the classical
algorithm or Karatsuba’s algorithm is used.

1.3.5 Unbalanced multiplication

The subquadratic algorithms considered so far (Karatsuba and Toom–Cook)
work with equal-size operands. How do we efficiently multiply integers of dif-
ferent sizes with a subquadratic algorithm? This case is important in practice,
but is rarely considered in the literature. Assume the larger operand has size
m, and the smaller has size n ≤ m, and denote by M (m, n) the corresponding
multiplication cost.
If evaluation-interpolation algorithms are used, the cost depends mainly on
the size of the result, i.e. m + n, so we have M (m, n) ≤ M ((m + n)/2), at
least approximately. We can do better than M ((m+n)/2) if n is much smaller
than m, for example M (m, 1) = O(m).
When m is an exact multiple of n, say m = kn, a trivial strategy is to cut the
larger operand into k pieces, giving M (kn, n) = kM (n) + O(kn). However,
this is not always the best strategy, see Exercise 1.16.
 1.3 Multiplication 9

When m is not an exact multiple of n, several strategies are possible:

• split the two operands into an equal number of pieces of unequal sizes;
• or split the two operands into different numbers of pieces.
Each strategy has advantages and disadvantages. We discuss each in turn.

First strategy: equal number of pieces of unequal sizes

Consider for example Karatsuba multiplication, and let K(m, n) be the num-
ber of word-products for an m × n product. Take for example m = 5, n = 3.
A natural idea is to pad the smaller operand to the size of the larger one. How-
ever, there are several ways to perform this padding, as shown in the following
figure, where the “Karatsuba cut” is represented by a double column:

a4 a3 a2 a1 a0 a4 a3 a2 a1 a0 a4 a3 a2 a1 a0
b2 b1 b0 b2 b1 b0 b2 b1 b0
A×B A × (βB) A × (β 2 B)

The left variant leads to two products of size 3, i.e. 2K(3, 3), the middle one to
K(2, 1)+K(3, 2)+K(3, 3), and the right one to K(2, 2)+K(3, 1)+K(3, 3),
which give respectively 14, 15, 13 word-products.
However, whenever m/2 ≤ n ≤ m, any such “padding variant” will re-
quire K(⌈m/2⌉, ⌈m/2⌉) for the product of the differences (or sums) of the
low and high parts from the operands, due to a “wrap-around” effect when
subtracting the parts from the smaller operand; this will ultimately lead to a
cost similar to that of an m×m product. The “odd–even scheme” of Algorithm
OddEvenKaratsuba (see also Exercise 1.13) avoids this wrap-around. Here is
an example of this algorithm for m = 3 and n = 2. Take A = a2 x2 + a1 x + a0
and B = b1 x + b0 . This yields A0 = a2 x + a0 , A1 = a1 , B0 = b0 , B1 = b1 ;
thus, C0 = (a2 x + a0 )b0 , C1 = (a2 x + a0 + a1 )(b0 + b1 ), C2 = a1 b1 .

Algorithm 1.5 OddEvenKaratsuba

m−1 n−1
Input: A = 0 a i x i , B = 0 bj x j , m ≥ n ≥ 1
Output: A · B
m−1
if n = 1 then return 0 a i b0 x i
write A = A0 (x ) + xA1 (x2 ), B = B0 (x2 ) + xB1 (x2 )
2

C0 ← OddEvenKaratsuba(A0 , B0 )
C1 ← OddEvenKaratsuba(A0 + A1 , B0 + B1 )
C2 ← OddEvenKaratsuba(A1 , B1 )
return C0 (x2 ) + x(C1 − C0 − C2 )(x2 ) + x2 C2 (x2 ).