Elementary Functions. Algorithms and Implementation

286 Pages·2006·1.304 MB·English

by  Jean-Michel Muller

Jean-Michel Muller
Elementary Functions
Algorithms and Implementation
Second Edition
Birkha¨user
Boston • Basel • Berlin
Jean-MichelMuller
CNRS-LaboratoireLIP
EcoleNormaleSuperieuredeLyon
46alleed’Italie
69364LyonCedex07
France
CoverdesignbyJosephSherman.
AMSSubjectClassifications:26A09,33Bxx,90Cxx,65D15,65K05,68Wxx,65Yxx
ACMSubjectClassifications:B.2.4,G.1.0.,G.1.2,G.4
LibraryofCongressCataloging-in-PublicationData
Muller,J.M.(Jean-Michel),1961-
Elementaryfunctions:algorithmsandimplementation/Jean-MichelMuller.–2nded.
p.cm.
Includesbibliographicalreferencesandindex.
ISBN0-8176-4372-9(alk.paper)
1.Functions–Dataprocessing2.Algorithms.I.Title.
QA331.M8662005
(cid:1)
518.1–dc22 2005048094
ISBN-100-8176-4372-9 eISBN0-8176-4408-3 Printedonacid-freepaper.
ISBN-13978-0-8176-4372-0
(cid:2)c2006Birkha¨userBoston,2ndedition
(cid:2)c1997Birkha¨userBoston,1stEdition
Allrightsreserved.Thisworkmaynotbetranslatedorcopiedinwholeorinpartwithoutthewrittenpermissionofthe
publisher(Birkha¨userBoston,c/oSpringerScience+BusinessMediaInc.,233SpringStreet,NewYork,NY10013,
USA),exceptforbriefexcerptsinconnectionwithreviewsorscholarlyanalysis.Useinconnectionwithanyformof
informationstorageandretrieval,electronicadaptation,computersoftware,orbysimilarordissimilarmethodology
nowknownorhereafterdevelopedisforbidden.
Theuseinthispublicationoftradenames,trademarks,servicemarksandsimilarterms,eveniftheyarenotidentified
assuch,isnottobetakenasanexpressionofopinionastowhetherornottheyaresubjecttoproprietaryrights.
PrintedintheUnitedStatesofAmerica. (LAP/HP)
987654321
www.birkhauser.com
Contents
ListofFigures xi
ListofTables xv
PrefacetotheSecondEdition xix
PrefacetotheFirstEdition xxi
1 Introduction 1
2 SomeBasicThingsAboutComputerArithmetic 9
2.1 Floating-PointArithmetic . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Floating-pointformats . . . . . . . . . . . . . . . . . . . . . 9
2.1.2 Roundingmodes . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.3 Subnormalnumbersandexceptions . . . . . . . . . . . . . 13
2.1.4 ULPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.1.5 Fusedmultiply-addoperations . . . . . . . . . . . . . . . . 15
2.1.6 Testingyourcomputationalenvironment . . . . . . . . . . 16
2.1.7 Floating-pointarithmeticandproofs . . . . . . . . . . . . . 17
2.1.8 Mapleprogramsthatcomputedouble-precision
approximations . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 RedundantNumberSystems . . . . . . . . . . . . . . . . . . . . . 19
2.2.1 Signed-digitnumbersystems . . . . . . . . . . . . . . . . . 19
2.2.2 Radix-2redundantnumbersystems . . . . . . . . . . . . . 21
I AlgorithmsBasedonPolynomialApproximationand/or
TableLookup,Multiple-PrecisionEvaluationofFunctions 25
3 PolynomialorRationalApproximations 27
3.1 LeastSquaresPolynomialApproximations . . . . . . . . . . . . . 28
3.1.1 Legendrepolynomials . . . . . . . . . . . . . . . . . . . . . 29
3.1.2 Chebyshevpolynomials . . . . . . . . . . . . . . . . . . . . 29
3.1.3 Jacobipolynomials . . . . . . . . . . . . . . . . . . . . . . . 31
vi Contents
3.1.4 Laguerrepolynomials . . . . . . . . . . . . . . . . . . . . . 31
3.1.5 Usingtheseorthogonalpolynomialsinanyinterval . . . . 31
3.2 LeastMaximumPolynomialApproximations . . . . . . . . . . . . 32
3.3 SomeExamples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.4 SpeedofConvergence . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.5 Remez’sAlgorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.6 RationalApproximations. . . . . . . . . . . . . . . . . . . . . . . . 46
3.7 ActualComputationofApproximations . . . . . . . . . . . . . . . 50
3.7.1 Getting“general”approximations . . . . . . . . . . . . . . 50
3.7.2 Gettingapproximationswithspecialconstraints . . . . . . 51
3.8 AlgorithmsandArchitecturesfortheEvaluationofPolynomials . 54
3.8.1 TheE-method . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.8.2 Estrin’smethod . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.9 EvaluationErrorAssumingHorner’sSchemeisUsed . . . . . . . 59
3.9.1 Evaluationusingfloating-pointadditionsand
multiplications . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.9.2 Evaluationusingfusedmultiply-accumulate
instructions . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.10 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4 Table-BasedMethods 67
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2 Table-DrivenAlgorithms . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.1 Tang’salgorithmforexp(x)inIEEEfloating-pointarithmetic 71
4.2.2 ln(x)on[1,2] . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.2.3 sin(x)on[0,π/4] . . . . . . . . . . . . . . . . . . . . . . . . 73
4.3 Gal’sAccurateTablesMethod . . . . . . . . . . . . . . . . . . . . . 73
4.4 TableMethodsRequiringSpecializedHardware . . . . . . . . . . 77
4.4.1 WongandGoto’salgorithmforcomputing
logarithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.4.2 WongandGoto’salgorithmforcomputing
exponentials . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.4.3 Ercegovacetal.’salgorithm . . . . . . . . . . . . . . . . . . 82
4.4.4 Bipartiteandmultipartitemethods . . . . . . . . . . . . . . 83
4.4.5 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5 Multiple-PrecisionEvaluationofFunctions 89
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2 JustaFewWordsonMultiple-PrecisionMultiplication . . . . . . 90
5.2.1 Karatsuba’smethod . . . . . . . . . . . . . . . . . . . . . . 91
5.2.2 FFT-basedmethods . . . . . . . . . . . . . . . . . . . . . . . 92
5.3 Multiple-PrecisionDivisionandSquare-Root . . . . . . . . . . . . 92
5.3.1 Newton–Raphsoniteration . . . . . . . . . . . . . . . . . . 92
Contents vii
5.4 AlgorithmsBasedontheEvaluationof
PowerSeries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.5 TheArithmetic-Geometric(AGM)Mean . . . . . . . . . . . . . . . 95
5.5.1 PresentationoftheAGM. . . . . . . . . . . . . . . . . . . . 95
5.5.2 ComputinglogarithmswiththeAGM . . . . . . . . . . . . 95
5.5.3 ComputingexponentialswiththeAGM . . . . . . . . . . . 98
5.5.4 Veryfastcomputationoftrigonometricfunctions . . . . . 98
II Shift-and-AddAlgorithms 101
6 IntroductiontoShift-and-AddAlgorithms 103
6.1 TheRestoringandNonrestoringAlgorithms . . . . . . . . . . . . 105
6.2 SimpleAlgorithmsforExponentialsandLogarithms. . . . . . . . 109
6.2.1 Therestoringalgorithmforexponentials . . . . . . . . . . 109
6.2.2 Therestoringalgorithmforlogarithms . . . . . . . . . . . 111
6.3 FasterShift-and-AddAlgorithms . . . . . . . . . . . . . . . . . . . 113
6.3.1 Fastercomputationofexponentials . . . . . . . . . . . . . 113
6.3.2 Fastercomputationoflogarithms . . . . . . . . . . . . . . 119
6.4 Baker’sPredictiveAlgorithm . . . . . . . . . . . . . . . . . . . . . 122
6.5 BibliographicNotes . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7 TheCORDICAlgorithm 133
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
7.2 TheConventionalCORDICIteration . . . . . . . . . . . . . . . . . 134
7.3 ScaleFactorCompensation . . . . . . . . . . . . . . . . . . . . . . 139
7.4 CORDICWithRedundantNumberSystemsandaVariableFactor 141
7.4.1 Signed-digitimplementation . . . . . . . . . . . . . . . . . 142
7.4.2 Carry-saveimplementation . . . . . . . . . . . . . . . . . . 143
7.4.3 Thevariablescalefactorproblem . . . . . . . . . . . . . . . 143
7.5 TheDoubleRotationMethod . . . . . . . . . . . . . . . . . . . . . 144
7.6 TheBranchingCORDICAlgorithm . . . . . . . . . . . . . . . . . . 146
7.7 TheDifferentialCORDICAlgorithm . . . . . . . . . . . . . . . . . 150
7.8 Computationofcos−1 andsin−1 Using CORDIC . . . . . . . . . . 153
7.9 VariationsonCORDIC . . . . . . . . . . . . . . . . . . . . . . . . . 156
8 SomeOtherShift-and-AddAlgorithms 157
8.1 High-RadixAlgorithms . . . . . . . . . . . . . . . . . . . . . . . . 157
8.1.1 Ercegovac’sradix-16algorithms . . . . . . . . . . . . . . . 157
8.2 TheBKMAlgorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 162
8.2.1 TheBKMiteration . . . . . . . . . . . . . . . . . . . . . . . 162
8.2.2 Computationoftheexponentialfunction(E-mode) . . . . 162
8.2.3 Computationofthelogarithmfunction(L-mode) . . . . . 166
viii Contents
8.2.4 Applicationtothecomputationofelementary
functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
III RangeReduction,FinalRoundingandExceptions 171
9 RangeReduction 173
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
9.2 CodyandWaite’sMethodforRangeReduction. . . . . . . . . . . 177
9.3 FindingWorstCasesforRangeReduction? . . . . . . . . . . . . . 179
9.3.1 Afewbasicnotionsoncontinuedfractions . . . . . . . . . 179
9.3.2 Findingworstcasesusingcontinuedfractions . . . . . . . 180
9.4 ThePayneandHanekReductionAlgorithm. . . . . . . . . . . . . 184
9.5 TheModularRangeReductionAlgorithm . . . . . . . . . . . . . . 187
9.5.1 Fixed-pointreduction . . . . . . . . . . . . . . . . . . . . . 188
9.5.2 Floating-pointreduction . . . . . . . . . . . . . . . . . . . . 190
9.5.3 Architecturesformodularreduction . . . . . . . . . . . . . 190
9.6 AlternateMethods . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
10 FinalRounding 193
10.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193
10.2 Monotonicity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
10.3 CorrectRounding:PresentationoftheProblem . . . . . . . . . . . 195
10.4 SomeExperiments . . . . . . . . . . . . . . . . . . . . . . . . . . . 198
10.5 A“Probabilistic”ApproachtotheProblem . . . . . . . . . . . . . 198
10.6 UpperBoundsonm . . . . . . . . . . . . . . . . . . . . . . . . . . 202
10.7 ObtainedWorstCasesforDouble-Precision . . . . . . . . . . . . . 203
10.7.1 Specialinputvalues . . . . . . . . . . . . . . . . . . . . . . 203
10.7.2 Lefe`vre’sexperiment . . . . . . . . . . . . . . . . . . . . . . 203
11 Miscellaneous 217
11.1 Exceptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217
11.1.1 NaNs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 218
11.1.2 Exactresults . . . . . . . . . . . . . . . . . . . . . . . . . . . 218
11.2 Notesonxy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
11.3 SpecialFunctions,FunctionsofComplexNumbers . . . . . . . . . 222
12 ExamplesofImplementation 225
12.1 Example1:TheCyrixFastMathProcessor . . . . . . . . . . . . . . 225
12.2 TheINTELFunctionsDesignedfortheItaniumProcessor . . . . . 226
12.2.1 Sineandcosine . . . . . . . . . . . . . . . . . . . . . . . . . 227
12.2.2 Arctangent . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
12.3 TheLIBULTIMLibrary . . . . . . . . . . . . . . . . . . . . . . . . . 229
12.4 TheCRLIBMLibrary . . . . . . . . . . . . . . . . . . . . . . . . . . 229
12.4.1 Computationofsin(x)orcos(x)(quickphase) . . . . . . . 230
Contents ix
12.4.2 Computationofln(x). . . . . . . . . . . . . . . . . . . . . . 230
12.5 SUN’sLIBMCRLibrary . . . . . . . . . . . . . . . . . . . . . . . . 231
12.6 TheHP-UXCompilerfortheItaniumProcessor . . . . . . . . . . 231
Bibliography 233
Index 261

List of Figures
2.1 Different possible roundings of a real number x in a radix-r
floating-pointsystem.Inthisexample,x>0 . . . . . . . . . . . . 12
2.2 Above is the set of the nonnegative, normalized floating-point
numbers(assumingradix2and2-bitmantissas).Inthatset,a−b
isnotexactlyrepresentable,andthefloating-pointcomputation
ofa−bwillreturn0intheroundtonearest,roundto0orround
−∞
to rounding modes. Below, the same set with subnormal
numbers.Now,a−bisexactlyrepresentable,andtheproperties
a(cid:2)=banda(cid:3)b(cid:2)=0(wherea(cid:3)bdenotesthecomputedvalueof
a−b)becomeequivalent . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Computation of 153120+112616 usingAvizienis’ algorithm in
radixr =10witha=6 . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Afull adder (FA) cell. From three bits x, y, and z, it computes
twobitstandusuchthatx+y+z =2t+u . . . . . . . . . . . . 22
2.5 A carry-save adder (bottom), compared to a carry-propagate
adder(top) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.6 APPM cell. From three bits x, y, and z, it computes two bits t
andusuchthatx+y−z =2t−u . . . . . . . . . . . . . . . . . . 23
2.7 Aborrow-saveadder . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.8 A structure for adding a borrow-save number and a
nonredundantnumber(bottom),comparedtoacarry-propagate
subtractor(top) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1 GraphofthepolynomialT7(x) . . . . . . . . . . . . . . . . . . . . 30
3.2 Theexp(−x2)functionanditsdegree-3minimaxapproximation
on the interval [0,3] (dashed line). There are five values where
the maximum approximation error is reached with alternate
signs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.3 The difference between exp(−x2) and its degree-3 minimax
approximationontheinterval[0,3] . . . . . . . . . . . . . . . . . 34