Modern Computer Arithmetic Richard P. Brent and Paul Zimmermann Version 0.2 Copyright c 2008 Richard P. Brent and Paul Zimmermann (cid:13) 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. Youmaynotusethisworkforcommercialpurposes. • No Derivative Works. You may not alter, transform, or build upon • this work. For any reuse or distribution, you must makeclear 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/ 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—wedo notcovercomputerarchitectureorthedesignofcomputer 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. Thealgorithmsthatwepresentaremainlyintendedforarbitrary-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 fourmainchapters, plus anappendix. Chapter1 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 cryptography, so most of the published books are now partly out of date or incomplete. Our aim has been to present the latest developments in a concise manner. Chapter 2 is concerned with the FFT and modular arithmetic, and their applications to computer arithmetic. We consider different number represen- tations, 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 64-bit arithmetic) is in- adequate. Thealgorithmsdescribedinthischapterfocusoncorrect rounding, extending the IEEE standard to arbitrary precision. Chapter 4 deals with the computation, to arbitrary precision, of functions 3 4 Modern Computer Arithmetic, version 0.2 of June 26, 2008 such as sqrt, exp, ln, sin, cos, and more generally functions defined by power series or continued fractions. We also consider the computation of certain constants, such as π and (Euler’s constant) γ. Of course, the computation of special functions is a huge topic so we have had to be selective. In particular, wehaveconcentratedonmethodsthatareefficientandsuitableforarbitrary- precision computations. For details that are omitted we give pointers in the Notes and References sections of each chapter, and in the bibliography. Finally, the Appendix contains pointers to implementations, useful web sites, mailing lists, and so on. The book is intended for anyone interested in the design and implemen- tation 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. 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. For solutions to the exercises, please contact the authors. WethanktheFrenchNationalInstituteforResearchinComputerScience and Control (INRIA), the Australian National University (ANU), and the Australian Research Council (ARC), for their support. The book could not have been written without the contributions of many friends and colleagues, toonumeroustomentionhere, but acknowledgedinthetextandinthe Notes and References sections at the end of each chapter. Finally, we thank Erin Brent, who first suggested writing the book, and our wives Judy-anne and Marie, for their patience and encouragement. This is a preliminary version — there are still a few sections to be com- pleted. We welcome comments and corrections. Please send them to either of the authors. Richard Brent and Paul Zimmermann [email protected] [email protected] Canberra and Nancy, June 2008 Contents Preface 3 Notation 9 1 Integer Arithmetic 11 1.1 Representation and Notations . . . . . . . . . . . . . . . . . . 11 1.2 Addition and Subtraction . . . . . . . . . . . . . . . . . . . . 12 1.3 Multiplication . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3.1 Naive Multiplication . . . . . . . . . . . . . . . . . . . 14 1.3.2 Karatsuba’s Algorithm . . . . . . . . . . . . . . . . . . 15 1.3.3 Toom-Cook Multiplication . . . . . . . . . . . . . . . . 16 1.3.4 Fast Fourier Transform (FFT) . . . . . . . . . . . . . . 18 1.3.5 Unbalanced Multiplication . . . . . . . . . . . . . . . . 18 1.3.6 Squaring . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.3.7 Multiplication by a Constant . . . . . . . . . . . . . . 20 1.4 Division . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.4.1 Naive Division . . . . . . . . . . . . . . . . . . . . . . . 22 1.4.2 Divisor Preconditioning . . . . . . . . . . . . . . . . . 24 1.4.3 Divide and Conquer Division . . . . . . . . . . . . . . 25 1.4.4 Newton’s Method . . . . . . . . . . . . . . . . . . . . . 28 1.4.5 Exact Division . . . . . . . . . . . . . . . . . . . . . . 28 1.4.6 Only Quotient or Remainder Wanted . . . . . . . . . . 29 1.4.7 Division by a Constant . . . . . . . . . . . . . . . . . . 30 1.4.8 Hensel’s Division . . . . . . . . . . . . . . . . . . . . . 31 1.5 Roots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 1.5.1 Square Root . . . . . . . . . . . . . . . . . . . . . . . . 32 1.5.2 k-th Root . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.5.3 Exact Root . . . . . . . . . . . . . . . . . . . . . . . . 35 5 6 Modern Computer Arithmetic, version 0.2 of June 26, 2008 1.6 Greatest Common Divisor . . . . . . . . . . . . . . . . . . . . 36 1.6.1 Naive GCD . . . . . . . . . . . . . . . . . . . . . . . . 37 1.6.2 Extended GCD . . . . . . . . . . . . . . . . . . . . . . 39 1.6.3 Half GCD, Divide and Conquer GCD . . . . . . . . . . 40 1.7 Base Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . 43 1.7.1 Quadratic Algorithms . . . . . . . . . . . . . . . . . . 43 1.7.2 Subquadratic Algorithms . . . . . . . . . . . . . . . . . 43 1.8 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 1.9 Notes and References . . . . . . . . . . . . . . . . . . . . . . . 49 2 The FFT and Modular Arithmetic 51 2.1 Representation . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.1.1 Classical Representation . . . . . . . . . . . . . . . . . 51 2.1.2 Montgomery’s Form . . . . . . . . . . . . . . . . . . . 52 2.1.3 Residue Number Systems . . . . . . . . . . . . . . . . 52 2.1.4 MSB vs LSB Algorithms . . . . . . . . . . . . . . . . . 53 2.1.5 Link with Polynomials . . . . . . . . . . . . . . . . . . 53 2.2 Addition and Subtraction . . . . . . . . . . . . . . . . . . . . 54 2.3 Multiplication . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.3.1 Barrett’s Algorithm . . . . . . . . . . . . . . . . . . . . 55 2.3.2 Montgomery’s Multiplication . . . . . . . . . . . . . . 56 2.3.3 Mihailescu’s Algorithm . . . . . . . . . . . . . . . . . . 60 2.3.4 Special Moduli . . . . . . . . . . . . . . . . . . . . . . 61 2.3.5 Fast Multiplication Over GF(2)[x] . . . . . . . . . . . . 62 2.4 Division and Inversion . . . . . . . . . . . . . . . . . . . . . . 68 2.4.1 Several Inversions at Once . . . . . . . . . . . . . . . . 70 2.5 Exponentiation . . . . . . . . . . . . . . . . . . . . . . . . . . 72 2.5.1 Binary Exponentiation . . . . . . . . . . . . . . . . . . 73 2.5.2 Base 2k Exponentiation . . . . . . . . . . . . . . . . . 74 2.5.3 Sliding Window and Redundant Representation . . . . 75 2.6 Chinese Remainder Theorem . . . . . . . . . . . . . . . . . . . 76 2.7 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 2.8 Notes and References . . . . . . . . . . . . . . . . . . . . . . . 78 3 Floating-Point Arithmetic 81 3.1 Representation . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.1.1 Radix Choice . . . . . . . . . . . . . . . . . . . . . . . 82 3.1.2 Exponent Range . . . . . . . . . . . . . . . . . . . . . 83 Modern Computer Arithmetic, 0.0 7 § 3.1.3 Special Values . . . . . . . . . . . . . . . . . . . . . . . 84 3.1.4 Subnormal Numbers . . . . . . . . . . . . . . . . . . . 84 3.1.5 Encoding . . . . . . . . . . . . . . . . . . . . . . . . . 85 3.1.6 Precision: Local, Global, Operation, Operand . . . . . 87 3.1.7 Link to Integers . . . . . . . . . . . . . . . . . . . . . . 88 3.1.8 Ziv’s Algorithm and Error Analysis . . . . . . . . . . . 88 3.1.9 Rounding . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.1.10 Strategies . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.2 Addition, Subtraction, Comparison . . . . . . . . . . . . . . . 94 3.2.1 Floating-Point Addition . . . . . . . . . . . . . . . . . 95 3.2.2 Floating-Point Subtraction . . . . . . . . . . . . . . . . 96 3.3 Multiplication . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 3.3.1 Integer Multiplication via Complex FFT . . . . . . . . 102 3.3.2 The Middle Product . . . . . . . . . . . . . . . . . . . 103 3.4 Reciprocal and Division . . . . . . . . . . . . . . . . . . . . . 105 3.4.1 Reciprocal . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.4.2 Division . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.5 Square Root . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 3.5.1 Reciprocal Square Root . . . . . . . . . . . . . . . . . 115 3.6 Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 3.6.1 Floating-Point Output . . . . . . . . . . . . . . . . . . 118 3.6.2 Floating-Point Input . . . . . . . . . . . . . . . . . . . 121 3.7 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 3.8 Notes and References . . . . . . . . . . . . . . . . . . . . . . . 124 4 Newton’s Method and Function Evaluation 127 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.2 Newton’s Method . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.2.1 Newton’s Method for Inverse Roots . . . . . . . . . . . 130 4.2.2 Newton’s Method for Reciprocals . . . . . . . . . . . . 130 4.2.3 Newton’s Method for (Reciprocal) Square Roots . . . . 131 4.2.4 Newton’s Method for Formal Power Series . . . . . . . 132 4.2.5 Newton’s Method for Functional Inverses . . . . . . . . 133 4.2.6 Higher Order Newton-like Methods . . . . . . . . . . . 134 4.3 Argument Reduction . . . . . . . . . . . . . . . . . . . . . . . 135 4.3.1 Repeated Use of Doubling Formulæ . . . . . . . . . . . 136 4.3.2 Loss of Precision . . . . . . . . . . . . . . . . . . . . . 136 4.3.3 Guard Digits . . . . . . . . . . . . . . . . . . . . . . . 137 8 Modern Computer Arithmetic, version 0.2 of June 26, 2008 4.3.4 Doubling versus Tripling Formula . . . . . . . . . . . . 137 4.4 Power Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 4.4.1 Direct Power Series Evaluation . . . . . . . . . . . . . 142 4.4.2 Power Series With Argument Reduction . . . . . . . . 142 4.4.3 The Rectangular Series Splitting . . . . . . . . . . . . 143 4.5 Asymptotic Expansions . . . . . . . . . . . . . . . . . . . . . . 147 4.6 Continued Fractions . . . . . . . . . . . . . . . . . . . . . . . 149 4.7 Recurrence Relations . . . . . . . . . . . . . . . . . . . . . . . 150 4.8 Arithmetic-Geometric Mean . . . . . . . . . . . . . . . . . . . 150 4.8.1 Elliptic Integrals . . . . . . . . . . . . . . . . . . . . . 151 4.8.2 First AGM Algorithm for the Logarithm . . . . . . . . 152 4.8.3 Theta Functions . . . . . . . . . . . . . . . . . . . . . . 153 4.8.4 Second AGM Algorithm for the Logarithm . . . . . . . 155 4.8.5 The Complex AGM . . . . . . . . . . . . . . . . . . . . 157 4.9 Binary Splitting . . . . . . . . . . . . . . . . . . . . . . . . . . 158 4.9.1 A Binary Splitting Algorithm for sin,cos . . . . . . . . 160 4.9.2 The Bit-Burst Algorithm . . . . . . . . . . . . . . . . . 161 4.10 Contour Integration . . . . . . . . . . . . . . . . . . . . . . . . 164 4.11 Other Special Functions . . . . . . . . . . . . . . . . . . . . . 164 4.12 Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 4.13 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 4.14 Notes and References . . . . . . . . . . . . . . . . . . . . . . . 167 5 Appendix: Implementations and Pointers 169 5.1 Software Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 5.1.1 CLN . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 5.1.2 GNU MP . . . . . . . . . . . . . . . . . . . . . . . . . 169 5.1.3 MPFR . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 5.1.4 ZEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 5.2 Mailing Lists . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 5.2.1 The BNIS Mailing List . . . . . . . . . . . . . . . . . . 170 5.2.2 The GMP Lists . . . . . . . . . . . . . . . . . . . . . . 170 5.3 On-Line Documents . . . . . . . . . . . . . . . . . . . . . . . . 171 Bibliography 173 Index 186 Notation C set of complex numbers N set of natural numbers (nonnegative integers) 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 the complex number z z Euclidean norm of the complex number z | | β “word” base (usually 232 or 264) n number of base β digits in integer or in floating-point significand, depending on the context ε “machine precision” 1β1−n 2 η smallest positive subnormal number (x) rounding of real number x ◦ 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(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 D(m,n) time to divide an m-bit integer by an n-bit integer 9 10 Modern Computer Arithmetic, version 0.2 of June 26, 2008 sign(n) +1 if n > 0, 1 if n< 0, and 0 if n =0 − r :=a modb integer remainder (0 r < b) ≤ q := adiv b integer quotient (0 a qb< b) ≤ − i j bitwise and of integers i and j, ∧ or logical and 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 ≫ ν(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 { ≤ ∧ } deg(A) for a polynomial A, the degree of A ord(A) for a power series A = a zj, ord(A) = min j : a = 0 j j { j 6 } (note the special case ord(0)= + ) P ∞ log, ln natural logarithm log , lg base-2 logarithm 2 nbits(n) lg(n) +1 if n > 0, 0 if n = 0 ⌊ ⌋ a t[a,b] or [a,b]t vector b (cid:18) (cid:19) f(n) = O(g(n)) c,n such that f(n) cg(n) for all n n 0 0 ∃ | |≤ ≥ f(n) = Θ(g(n)) f(n)= O(g(n)) and g(n) = O(f(n)) a b [a,b;c,d] 2 2 matrix × c d (cid:18) (cid:19) xxx.yyy a number xxx.yyy written in base ρ; ρ for example, the decimal number 3.25 is 11.01 in binary 2