Scientific Programming
Editors: Konrad Hinsen, konrad.hinsen@cnrs-orleans.fr | Matthew Turk, matthewturk@gmail.com

Practically Accurate Floating-Point Math
Neil Toronto and Jay McCarthy | Brigham Young University

With the right tools, floating-point code can be debugged like any other code, drastically improving its
accuracy and reliability.

M

ost computing scientists know something about
the dangers of treating floating-point numbers
as if they’re real numbers. We even know somefolk wisdom intended to avoid the dangers.
Test for nearness instead of equality. If the standard libraryhas a function to compute a hypotenuse, use it instead of
x 2 + y 2 . Use (x − y) . (x + y) instead of x 2 − y2. Computesums in order of decreasing magnitude—or is it increasing
magnitude?
Sometimes we apply these rules habitually. Sometimes weapply them after tracing a reproducible error to an obvious
problem. Sometimes we fire them at our code from a shotgun
because it gives obviously wrong answers after many hours ofapparently proper functioning. But most of the time, we can’ttell whether this folk wisdom actually improves anything.
The typical prescription—error analysis—is timeconsuming and requires expertise that most of us can’t seemto find the time to obtain. We’re apparently stuck between
the two extremes of ignorance and certainty. Unfortunately,
the initial and continuing costs of certainty often makeignorance more practical.
In almost every other software development task, testingand debugging is a practical middle ground. We gain knowledge
and confidence about an algorithm’s implementation by its
behavior on representative inputs, using informal arguments
about path coverage and similarity. With the right tools, testingand debugging floating-point code is no different, and can
drastically improve its accuracy and reliability.
These tools and the principles behind them can beused for more than just fixing badly behaved floating-point
code: they can also help us determine whether our researchquestions themselves are buggy. For concreteness, let’s
first concentrate on debugging the code that answers thequestions, and then expand our scope.
80

CISE-16-04-Sciprog.indd 80

Computing in Science & Engineering

You Could Have Invented Floating-Point Math
While easier than error analysis, floating-point debugging
still requires some background knowledge. Therefore, weneed to review what every computing scientist really needs
to know about floating point. This includes
■■
■■
■■

how floats are distributed on the real line;
how error is defined and how to measure it; andhow errors originate and propagate.

For all of these, it helps to know something about how
floating-point numbers are represented.
Floating-point numbers are little more than scien­tific notation encoded as bits. It’s instructive, and notthat difficult, to formalize the rules we’ve all learned
for manipulating numbers in scientific notation bycoding up floating-point math from scratch. Fortunately,
we can ignore all those capricious rules about significant
figures.
Any real number x that can be put in scientific notation
can also be put in the form
x = s . m . 10e,(1)where s is −1 or 1, m is a natural number, and e is an
integer. For example, to put −1.23456789 . 10−10 in thisform, extract the sign, move the decimal point 8 places tothe right, and subtract 8 from the exponent to get −1 .
123456789 . 10−18.
Because computers operate more efficiently on powers of
two than on powers of 10, floating-point numbers represent
real numbers that can be put in the form
x = s . m . 2e,(2)1521-9615/14/$31.00 © 2014 IEEE

Copublished by the IEEE CS and the AIPJuly/August 201428/07/14 7:01 PM

where s ∈ {−1, 1}, m ∈  , and e ∈  , as before. In a float,s is called the sign, m is called the significand (an old-schooland somewhat confusing synonym is mantissa), and e iscalled the exponent.
To code up floats in the Racket programming language(see the sidebar on Racket and its design),1 we first define a
new struct type to represent a float:#lang typed/racket(struct float ([sign : Integer][significand : Natural][exponent : Integer])#:transparent)

The first line of our program, #lang typed/racket, tells
Racket that we’re using the typed dialect. The next few lines
define the struct type for floats, containing three fields
with types Integer, Natural, and Integer. (We could haveused [sign : (U −1 1)] to tell Typed Racket that the sign
can be only −1 or 1, but it would make working with thesign field more difficult.) It’s important to note that Naturaland Integer instances aren’t limited to the system’s word
size or to any other maximum size at all: they’re big integers.
The #:transparent flag tells Racket, among other things,
to not hide a printed float’s fields.
Using x = s . m . 2e as a guide, we define a function toconvert floats to exact rationals:(: rat (−> float Exact-Rational))(define (rat x)(match-define (float s m e) x)(* s m (expt 2 e)))

The first line declares that the function rat takes
arguments of type float and returns Exact-Rational
instances, which are represented internally by reduced
quotients of big integers. The line (match-define (floats m e) x) unpacks argument x’s fields. The next line,(* s m (expt 2 e)) simply encodes the formula s . m . 2 e
in Racket.
If we run the program we have so far in the DrRacket
IDE, an interactions window opens with a friendly welcome
and a “>” prompt. We should try a few test cases.
One way to represent 80 as a float is 1 . 10 . 23. In theinteractions window, we write (float 1 10 3), press Enter,and see> (float 1 10 3)(float 1 10 3)

DrRacket has printed the results of applying the float
function to arguments 1, 10, and 3: a float with fields 1,
www.computer.org/cise			

CISE-16-04-Sciprog.indd 81

The Racket Programming
LanguageRacket’s syntax (see www.racket-lang.org) is similar to
Lisp’s, so expect nested parentheses and (notationprefix). But it goes far beyond Lisp in its ability to define
new languages as libraries. One new language that ships with
Racket is Typed Racket, in which all the code in this article
is written.
Racket and its libraries are open source (Lesser General
Public License, or LGPL). They’ve been under development for
20 years by PLT, a group of programming language researchers
who put correctness first but aim for practical applicability. As aresult, Racket’s libraries are broad and reliable, and its runtime
performance is good. In particular, Racket’s floating point is bothstandards-compliant and fast: within a factor of 3 of C’s speed
when carefully tuned, which Typed Racket accelerates to withina factor of 1.5.
Racket ships with an integrated development environment(IDE) called DrRacket. Originally developed for teaching
programming, it’s unobtrusive and useful for both students andexperts. DrRacket provides programmers with abilities that most
IDEs don’t, such as the ability to manipulate pictures and store
them in programs as literal values, as easily as if they were
simply numbers.

10, and 3. So far, so good. Converting it to an exact rational
works as expected:> (rat (float 1 10 3))
80

Another correct test case is as follows:> (rat (float −1 23 −6))
−

26
.
34

That is, −1 ⋅ 23 ⋅ 2−6 = −

23
.
64

Multiplying float instances is easy, using the propertiesof multiplication and the identity 2e1 ⋅ 2e2 = 2e1+e2 . If
x1 = s1 ⋅ m1 ⋅ 2e1 and x 2 = s 2 ⋅ m2 ⋅ 2e2 , thenx1 ⋅ x 2 = ( s1 ⋅ m1 ⋅ 2e1 ) ⋅ ( s 2 ⋅ m2 ⋅ 2e2 )= ( s1 ⋅ s 2 ) ⋅ (m1 ⋅ m2 ) ⋅ 2e1+e2 .

(3)

	

81

28/07/14 7:01 PM

Scientific Programming(: add (−> float float float));; Returns the exact sum of two floats(define (add x1 x2)(match-define (float s1 m1 e1) x1)(match-define (float s2 m2 e2) x2)(define d (− e1 e2))
(define n(if (< d 0)(+ ( s1 m1) ( s2 m2 (expt 2 (− d))))(+ ( s2 m2) ( s1 m1 (expt 2 d)))))(float (if (< n 0) −1 1) (abs n) (min e1 e2)))

Figure 1. An implementation of exact floating-point addition.

Using the last formula as a guide, we define(: mul (−> float float float))(define (mul x1 x2)(match-define (float s1 m1 e1) x1)(match-define (float s2 m2 e2) x2)(float (* s1 s2) (* m1 m2) (+ e1 e2)))

Checking it against Racket’s exact rational arithmetic,
we get>(* 80 − 23/64)
3−28
4>(rat (mul (float 1 10 3)(float − 1 23 − 6)))3−28
4Addition is worse: the number of bits required to store
the result’s significand depends on the difference in thearguments’ exponents. For example, the following sum has
a 1 + 16 − (−16) = 33-bit significand, computed from two
1-bit significands:> (add (float 1 1 16) (float 1 1 −16))(float 1 4294967297 −16)

Division and irrational functions like square root areeven worse: most results require infinitely many bits. We
clearly need to limit the significand’s range.
An unbounded exponent is also a curse. For example,
it would take billions of terabytes to store just the exponent
field of exp(exp(exp(100))).
At this point, we could review how floats are packed
into bit fields to be stored in memory. But it takes a lot ofwork to turn those details into intuition, which we could getjust as easily by visualizing the floats on the number line. An
exhaustive visualization is infeasible for standard float sizes—
billions of floats can be packed into the smallest standard
size, 32 bits—so let’s take a look at 6-bit floats.
The 6-bit floats we’ll use have one sign bit to represents ∈ {−1, 1}, two significand bits to represent m ∈ [4..7],
and three exponent bits to represent e ∈ [−4..1]. They’re
encoded as bits in a way that gives them the same properties
as standard 32-, 64-, and 128-bit floats. Figure 2 shows theirpositions on the extended real line, except those that encode
the error value not-a-number (NaN).
Because there are finitely many 6-bit floats, they aren’t
infinitely dense like rational numbers or instances of float.
Consequently,■■■■

as expected.
Implementing addition is more involved, requiring
shifting the significand of one float to line up its “decimal
point” with the other float’s, but it’s certainly within thereach of anyone reading this article. Figure 1 gives a Typed
Racket implementation.
Floats on the Number Line
Instances of float are different from typical floating-point
numbers in one crucial way: because a float’s significand andexponent are big integers, there’s no bound on their size. While
this might seem like an improvement, it’s actually kind of a curse.
Multiplication can double the number of bits required tostore the significand. For example, multiplying two floats with
4-bit significands can result in a float with an 8-bit significand:> (mul (float 1 15 0) (float 1 15 0))(float 1 255 0)

82

CISE-16-04-Sciprog.indd 82

■■every float has at least one closest neighbor;floats can be enumerated in their natural order; andevery extended real number is either a float or is ­between
two floats.

As an example of the last consequence, 19 is between 14
and ∞. Also, p is between 3 and 3.5, but is closer to 3, sowith 6-bit floats, p = 3. Many struggling math students, andsome historic lawmaking bodies, would heartily approve.
To keep the numeric labels from overlapping, Figure
2 has to “zoom in” to the floats near zero, twice. Because
there are the same number of floats between every twoneighboring powers of two, floats are quite dense near zeroand quite sparse away from it. This has somewhat surprisingconsequences. For example, every 6-bit float above 22 =4 is an integer, and every 6-bit float above 23 = 8 is even.
Likewise, every 64-bit float above 252 is an integer, andevery 64-bit float above 253 is even.
Floats that can be written s . m . 2e, where (for 6-bit floats)s ∈ {−1, 1}, m ∈ [4..7], and e ∈ [−4..1], are called normal. The
July/August 201428/07/14 7:01 PM

–∞

–14

–12

–10

–8 –7 –6 –5 –4

4

7 8

10

12

14

∞

–4

–3½

–3

–2½

–2 –1¾ –1½ –1¼ –1

1 1¼ 1½ 1¾ 2

2½

3

3½

4

–1

–7/8

–¾

–5/8

–½ –7/16 –3/8 –5/16 –1/4–3/16 –1/8 –1/16 ±0

/

5 8

/

3 4

/

7 8

/

1

subnormal//

1 16 1 8/

3 16 1 4

5

/

6

/

/

/

5 16 3 8 7 16 1 2

subnormal

Figure 2. The positions of all 58 non-not-a-number (NaN), 6-bit floats on the extended real line.

strange range of m helps ensure that the encoding of non-NaN
floats is unique. If we allowed m = 0, then 0 could be represented
by both 1 . 0 . 20 and 1 . 0 . 21, among others.
But if so, 0 can’t be a normal float. One exponent field bit
pattern is thus reserved for specially encoding 0 and the next
three uniformly spaced floats, which are called subnormal. For
standard 64-bit floats, there are 252 −1 positive subnormals,
uniformly spaced in the open interval (0, 2−1,022).
Another exponent field bit pattern is reserved forspecially encoding −∞, ∞, and NaN.
There are two floating-point zeros. It’s occasionally useful
to distinguish their signs, but we’ll regard them as equal.
One fact Figure 2 can’t show is that 6-bit and standard-size
floats are cleverly encoded to make it easy to enumerate them intheir natural order. In particular, the encoding of any non-NaN
float x, interpreted as a signed magnitude integer, or an integer
whose most significant bit determines whether it’s positive or
negative, is the signed index of x on the number line. For example,
the 6-bit encoding of 1 is 0011002, or 12. The encoding of −1 is1011002, which, as a 6-bit signed magnitude integer, is −12.
We’ll define ord(x) to be a float x’s signed ordinal index,
so that ord(1) = 12 and ord(−1) = −12 for 6-bit floats.
Figure 3 gives C implementations of ord and its inverse ord−1
for 64-bit floats. Racket’s 64-bit implementations of ord andord−1 are provided by its math library as flonum−>ordinaland ordinal−>flonum.
We care about the ordinal indexes of floats because
they enable us to measure error and to design effective
randomized testing strategies. For now, let’s focus on
measuring error.
Measuring Floating-Point Error
Functions that return approximations can at best be close tocorrect, so we need a natural way to measure the amount oferror in their outputs. For a float x that approximates a real
number r, the most obvious error measure is the distancewww.computer.org/cise			

CISE-16-04-Sciprog.indd 83

int64_t double_to_ordinal(double x) {double y = fabs(x);int64_t n = *((int64_t*)(&y));return x <= 0.0 ? -n : n;}double ordinal_to_double(int64_t n) {int64_t m = n < 0 ? -n : n;double x = *((double*)(&m));return n < 0 ? -x : x;}

Figure 3. C implementations of ord(x) and ord –1 (x) for 64-bit
floats. This code should work correctly on any system that stores
floats with the same endianness as integers.

MPFR

T

he Multiple-Precision Float with correct Rounding (MPFR;
www.mpfr.org) library is an open source project whose goal isto provide well-defined, standards-compliant, efficient, multi­
precision floating-point functions. Its contributors take great
pains to prove that the algorithms it uses always produce correctlyrounded results.

between them, or |x − r|. This is called absolute error, and it’susually the wrong measure.
To see why, we’ll first define fl(r) to be the nearest float
to a real number r. In case of a tie, fl(r) is the nearest float
with an even significand. If r is too large, fl(r) = ±∞.
	83

28/07/14 7:02 PM

Scientific ProgrammingIEEE 754-2008 ComplianceT

he first floating-point standard, IEEE 754-1985, was introduced
to stem the tide of proliferating floating-point formats and implementations. It requires all basic functions to be correctly rounded—essentially, to be as accurate as possible. Consumer chip makers,
who didn’t have to support legacy floating-point units, embraced thestandard with uncharacteristic zeal. The 2008 update addressed ambiguities and added a list of recommended basic functions.
As computing scientists, we should seek and emphaticallysupport IEEE 754-2008 compliance. It makes our computations
more accurate and our scientific results repeatable.
It seems most modern desktop hardware is mostly compliant,and system libraries are catching up. Unfortunately, many
compilers still aggressively optimize floating-point code in a waythat breaks compliance.
The following short Racket program tests hardware and systemlibrary compliance with IEEE 754-2008:#lang racket(require math/utils)(test-floating-point 10000)
This program tests Racket’s basic functions, as well as additional
math library functions whose implementations are very sensitive
to noncompliance, on thousands of deterministic and 10,000
random arguments, against MPFR. An automated process runs
a similar program on a compliant system every time a Racket
developer uploads changes to ensure Racket remains compliant.
According to reports from Racket users who have run the
program on their own systems, it appears that■■ arithmetic and square root are almost universally correctlyrounded;■■ on many systems, trigonometric functions are correctly rounded

for arguments near zero, but return values with a few ulps error
for very large arguments; and■■ on a small handful of systems, exp(x) or trigonometric functions are accurate near zero but have error up to hundreds ofulps when x is large.
IEEE 754-2008 doesn’t specify which floating-point functions
are implemented in hardware, so it’s difficult to make a generalstatement about whether the errors are caused by hardware
floating-point units or system libraries.

With 64-bit floats, |fl(exp(100)) − exp(100)| ≈ 1.61 .

1027, which is a ludicrous amount of error for the closest

approximation of exp(100)—it’s about 1.5 times thediameter of our galaxy in micrometers. The problem here is
84

CISE-16-04-Sciprog.indd 84

the sparseness of 64-bit floats near exp(100): the floats oneither side of it are 292 ≈ 4.95 . 1027 apart!
The standard way to measure floating-point error
relativizes absolute error by dividing it by the distance
between neighboring floats near r:
error( x , r ) =x − r .ulp( fl(r ))(4)Here, ulp(y) is y’s unit in last place, which is defined as themagnitude of the lowest-order bit in y’s significand, or equivalently
the distance from y to the next float from y away from zero. (If thenext float from y is an infinity, ulp(y) is the distance to the next
float toward zero.) It’s easy to compute ulp(y) using ord and ord–1.
There’s a nice intuitive interpretation of error(x, r):it corresponds closely with the number of floating-point
numbers from x to r. For example, with 6-bit floats,error(8,14) =8 − 146= = 3.(5)ulp(14) 2This interpretation admits fractional numbers of floats aswell. For example,
error( fl(exp(100)),exp(100)) ≈1.61⋅1027≈ 0.325.(6)292That is, fl(exp(100)) and exp(100) are about 0.325 floats
apart. Generally, if fl(r) is finite, then error( fl(r ), r ) ≤ 0.5,
meaning that the nearest float to a real number r is no morethan half a float, or half an ulp, away.
Racket’s math library provides ulp and error as flulp
and flulp-error. (Racket’s 64-bit floating-point function
names usually start with fl.) To use flulp-error to
debug floating-point functions, we need correct values for
its second argument. Racket’s built-in exact rationals canrepresent correct values, at least for the results of arithmetic.
For example,> (require math)> (flulp-error (/ 1.0 7.0) (/ 1 7))
0.2857142857142857

For irrational functions such as exp, we can appro­
ximate correct values closely with high-precision floats.
Racket’s math library provides such floats by wrapping the
C library MPFR,2 the most accurate multiprecision floatingpoint library in existence (see the “MPFR” sidebar). For
example,> (flulp-error (flexp 100.0)(bigfloat−>real	  (bfexp (bf 100.0))))0.32516258740803655
July/August 201428/07/14 7:02 PM

> (require plot/typed)> (plot (function (fun-error flexp bfexp) 0 2)#:x-label "x"#:y-label "exp(x) error")0.4> (plot3d (contour-intervals3d(fun2d-error fl— bf—)
—1 1 —1 1))

0.5
0.4

exp(x) error

0.3
0.3

0.2
0.1
0
1

0.2

0.5
1

0

0.5

xis

ya

0.1

0

–0.5

0
0

0.5

1
x

1.5

–0.5

2

(a)

–1 –1

xis

xa

(b)

Figure 4. Error plots for basic functions that are correctly rounded in accordance with IEEE 754-2008. (a) Exponentiation and(b) subtraction.

The expression (bf 100.0) returns an MPFR floatingpoint number, or a bigfloat, which by default has a ­128-bit
significand. Then bfexp computes exp(100) with highprecision, and bigfloat−>real converts the result to an
exact rational.
We’ll frequently use bigfloat functions to estimate error
in the outputs of 64-bit floating-point functions, so it’s
helpful to abstract away the details. We define(: fun-error (−> (−> Flonum Flonum)(−> Bigfloat Bigfloat)(−> Real Flonum)))(define ((fun-error f f*) r)(define x (fl r))(flulp-error (f x)(bigfloat−>real (f* (bf x)))))

Given a 64-bit floating-point function f and a bigfloat
function f*, fun-error produces a function that accepts a real
number r and returns an estimate of the error in (f (fl r)),in ulps. To use it to measure the error in flexp, for example,> ((fun-error flexp bfexp) 100)0.32516258740803655

which is the same as before.
www.computer.org/cise			

CISE-16-04-Sciprog.indd 85

Now that we can get flexp’s error as a function, it’seasy to plot. Figure 4a is a transcript from DrRacket’sinteractions window, plotting (fun-error flexp bfexp).
(In Racket, plots are values that are defined in a way thatcauses DrRacket to print them as graphics.) It appears thaton this system, for arguments between 0 and 2, (fun-error
flexp bfexp) never returns a value greater than 0.5. If so,
flexp is accurate within half an ulp—as accurate as it canpossibly be.
According to the IEEE 754-2008 floating-point standard,3 such accuracy is expected. Section 4.3 says
Except where stated otherwise, every operation shall be
performed as if it first produced an intermediate result
correct to infinite precision and with unbounded range,
and then rounded ...

More precisely, if f is a floating-point implementationof f, and x is the exact floating-point representation of areal value x, then (f x) must be the exact floating-point
representation of fl( f (x)). Operations that are computed
this accurately, as mandated by IEEE 754-2008, are called
correctly rounded.
Let’s test another function. First, we’ll define an abstraction
similar to fun-error but for  × →  functions:	85

28/07/14 7:02 PM

Scientific Programming> (plot (function(lambda (r)(define x (fl r))(flulp-error (fllog (flnext x))(bigfloat—>real(bflog (bf x)))))
0 2)#:x-label "x"#:y-label "log(x+ulp(x)) error")writers to adhere to IEEE 754-2008 (see the related sidebar).
Unfortunately, we’re not so lucky.
Error Propagation
The most disappointing and difficult fact about floatingpoint error is that it isn’t compositional: low error in aformula’s subterms doesn’t guarantee low error in the formula
itself. One consequence is that correctly rounded functions
are guaranteed to have low error only when their arguments
are exact. We’ll use this consequence to demonstrate theproblem.
Suppose we give fllog an argument with just one ulp
error. We can do that using flnext, which returns the next
64-bit float toward positive infinity. We then measure theerror of its log:

500

log(x+ulp(x)) error

400

300

> (flulp-error (fllog (flnext 1.1))(bigfloat−>real(bflog (bf 1.1))))
14.572897331748424

200

100

0

0.5

1

1.5

2

x
Figure 5. How log(x) magnifies error in x.

(: fun2d-error(−> (−> Flonum Flonum Flonum)(−> Bigfloat Bigfloat Bigfloat)(−> Real Real Flonum)))(define ((fun2d-error f f*) rx ry)(define x (fl rx))
(define y (fl ry))
(flulp-error (f x y) (bigfloat−>real(f* (bf x) (bf y)))))The plot of (fun2d-error fl− bf−) in Figure 4b
suggests subtraction is correctly rounded on this system. The
trough in the center demonstrates an additional, useful fact:
when x and y are within a factor of two of each other, or
1 2 ≤ x /y ≤ 2 , then x − y is not only correctly rounded but
is exact.
It would be nice if low error in basic functions’ outputs
meant that every floating-point formula had low error.
If that were true, this article could end now with no morethan an impassioned plea to hardware makers and compiler
86

CISE-16-04-Sciprog.indd 86This isn’t terrible. 64-bit floats have 53-bit significands.
Being 14.6 ulps off means that the log2(14.6) ≈ 3.87
low-order significand bits are bad, so the result’s 53-bit
significand still has about 49 bits precision.
But really, how bad can it possibly get?
Unfortunately, there’s effectively no bound on how bad
it can get:> (flulp-error (fllog (flnext 1.0001))(bigfloat−>real(bflog (bf 1.0001))))
16381.87167545197> (flulp-error (fllog (flnext 1.00000000001))(bigfloat−>real(bflog (bf 1.00000000001))))137438953470.5061

In the last example, we have only 16 bits precision. For
perspective, this is just a little more precise than 22 7 ≈ π.
This sensitivity to argument error seems to get worse as
the argument approaches 1. Let’s plot how sensitive log is
to error in arguments between 0 and 2 to find out if that’strue. This time, we’ll plot a lambda, or an anonymous
function, instead of using a function like fun-error to create
a function to plot.
Figure 5 contains the transcript from DrRacket’sinteractions window. The spike at x = 1 looks perfectly
dreadful, but it’s actually worse than it looks. If we were
to zoom into it by plotting at ever-narrower intervals that
contain 1, we’d find that the spike goes all the way up to
infinity. There’s no limit to how much log can magnify its
argument’s error.
July/August 201428/07/14 7:02 PM

There’s a classic joke pertinent to our situation:

Log-One-Pea? Exp-Em-One?

[Doctor arrives to find patient bending joint backwards,
biting armpit, or other grotesque action.]
Patient: Doctor, it hurts when I do this.
Doctor: Then don’t do that!

I

We can’t demand that fllog be more accurate if it’s alreadyas accurate as it can possibly be. We can’t apply fllog only to
exact arguments because the outputs of other floating-point
functions are almost always approximate. All we have left is,“Then don’t do that!” Don’t apply fllog to an approximate
result that’s near 1.
Is this normal? Unfortunately, yes. Most common
real-valued functions, even when their floating-point
implementations are correctly rounded, have subdomains on
which they magnify argument error in their outputs. These
places are often called ill-conditioned, but that term doesn’t
fill anyone with the right amount of dread. Instead, we’ll callthose subdomains the badlands.
High error in floating-point code isn’t caused byrounding exact results to the nearest float. It’s caused bywandering into the badlands with a rounded result.
Avoiding the Badlands
The main task in debugging floating-point math is figuring
out when we’ve wandered into the badlands with anapproximation and then rewriting parts of our code so we
don’t do that.
It’s therefore essential to know where common
functions’ badlands are. Figures 6 and 7 have plots andformulas for the badlands of the most common floatingpoint functions. To derive them, we’ve made the definitionof badlands more precise: they’re the parts of the domain inwhich error is magnified about 4 or more times, without
considering additional rounding error. In other words, thebadlands reduce precision by two bits or more, as well as
add rounding error.
Perhaps surprisingly, some common functions don’t
have badlands, such as multiplication and arctangent. We
can confidently apply them to approximations.
Figure 6 identifies the badlands for x − y as1
≤ x /y ≤ 2 .(7)2
In other words, the badlands are where x and y arewithin a factor of two of each other. This is the same
subdomain on which floating-point subtraction is exact
when x and y are exact, so whether subtraction in thebadlands is prudent or deadly depends only on exactness.
The universe is clearly messing with us, but at least it’smade this important subdomain of  ×  easier toremember.
www.computer.org/cise			

CISE-16-04-Sciprog.indd 87

t might seem like a stretch to say log1p(x) = log(1 + x) and
expm1(x) = exp(x) − 1 are common functions. They shouldbe, because they so often provide an escape from log’s andexp’s badlands. If the language you use for scientific computing
doesn’t have them, submit a bug report.

If f :  →  is differentiable, its badlands, as we’ve
defined the term, can be estimated by solvingx ⋅ f ′( x )≥ 4 (8)f (x )
for x, where f ′ is the first derivative of f. The left-handside is called f ’s condition number at x. There are moregeneral condition number formulas for multivariate
functions, but they’re more complicated than we need to
get into.
We say the badlands can be “estimated” because thecondition number formula doesn’t account for any floatingpoint details. But the fact that such a simple formula canestimate how much f magnifies floating-point error is
telling: it means that f ’s badlands are primarily a propertyof f itself, and much less a property of its floating-point
implementation.
We can use the condition number formula to getsome intuition about how to avoid the badlands. Whenthe denominator f (x) is near zero and the numerator
isn’t, the condition number is large. We should therefore
avoid zero crossings that aren’t at the origin. A canonicalexample is log(1) = 0. Another trait to avoid isexponential growth away from the origin, which makes
the numerator large.
Figures 6 and 7 plot badlands that can be estimated using
condition numbers, which estimate how error is magnified
only for normal floats. Floating-point implementations of somefunctions, particularly x/y, x , and log(x), additionally magnify
error without bound when their arguments are approximate andsubnormal.
A sometimes convenient, but complicating fact isthat high error isn’t compositional, either: high error in aformula’s subterms doesn’t guarantee high error in theformula itself. For example, applying fllog far enough fromzero reduces error:> (flulp-error (fllog (flstep 1e10 40))
(bigfloat−>real(bflog (bf 1e10))))2.11101355742103
	

87

28/07/14 7:02 PM

Scientific Programming

8

Operation

Badlands

x+y

1

x–y

1

x·y

None

x /y

x or y subnormalx8log(x + ulp(x)) error

log1p(x + ulp(x)) error

2 ≤ –x / y ≥ 2
2 ≤ x /y ≤ 2

6

6

4

4

x subnormallog(x)

exp – 1 4 ≤ x ≤ exp

exp(x)|x | ≥ 4

log1p(x)

x ≤ –

15

1

4 or x subnormal

exp(–1/4) exp(1/4)
exp(–1/2)
exp(1/2)
0

16

x ≥ 4

expm1(x)

2

2

0

0.5

1

2

1.5

2.5

(a)

–1/2
–0.6

–0.4

–0.2

0

(b)

8

(x + ulp(x)) – (y – ulp(y)) error–15/16
0
–1
–0.8

8

exp(x + ulp(x)) error

8
6
4
2
0
1

6

6

0.5

4

4

2

2

expm1(x + ulp(x)) error

0
–0.5
–1–1

0

–0.5

(c)

0.5

1
0

–4

–2

0

2

4

(d)

0
–4

–2

0

2

4

6

8

(e)

Figure 6. Formulas and plots of the badlands for common 64-bit floating-point functions.

Here, (flstep 1e10 40) is equivalent to 40 iterations offlnext, starting with 1e10. Computing the log has reduced
error from 40 ulps to 2 ulps.
It’s occasionally useful to reformulate floating-point
math to reduce error. We’ll let error reduction happen by
accident, and use testing to find out where it doesn’t make up
for wandering into the badlands.
Floating-Point Test-Debug Cycle
A good test-debug cycle for a floating-point function f proceeds as follows (also see the sidebar, “Testing and Debugging Support”). First, write a straightforward version f. Next,
write a version f* with higher precision. Then repeat the following steps as necessary:
1. Test f against f* to find high-error subdomains, using
only exact arguments to f.
2. Reformulate f to reduce error.
3. Make the same changes to f*, if necessary.
We’ll demonstrate this test-debug cycle on three functions,which will illustrate the answers to the following questions.
88

CISE-16-04-Sciprog.indd 88



Q. Are you seriously advocating writing two versions of every
floating-point function?

A. Yes. Admittedly, that answer isn’t very motivating, sohere’s a counterquestion: How certain do you want to bethat your good, bad, or unexpected results aren’t an artifact
of floating-point error?
Q. What kind of numbers should f* use?

A. For arithmetic functions, use exact rationals if yourlanguage or libraries provide them. Otherwise, multiprecision
floats with 128-bit significands seem to work well for testing
most 64-bit floating-point functions.
Q. Why test f only with exact arguments?

A. How f magnifies error when given approximate arguments
is a property of the function it implements. We can’t doanything about it.
Q. Why make corresponding changes to f*?

A. Unless f* is written using exact rational arithmetic, if
f wanders into the badlands with an approximation, so
does f*.
July/August 201428/07/14 7:02 PM

Operation

Badlands

tan(x)x ≥ arctan( 7)sin(x)x ≥cos(x)xarctan(x)Nonearcsin(x)x ≥arccos(x)

5
6

8

π

≥ π or1
3

π ≤ x ≤

2
3

π

63

6

6

4

4

–sqrt(63/64)
1
2

or x ≤ –255
256–sqrt(15/16)
0–1–0.975 –0.95|0.95

0.975

1

arctan(sqrt(7))

111–sqrt(63/64)
0–1–0.995

3
3
3

2/3 π

3

0.9

1

4
1/3 π

5/6 π

2

3

2
1

1

–5/6 π–5/6 π

-1
-1
-1

(c)

–0.99|0.8

cos(θ + ulp(θ)) error

4

2/3 π

5/6 π

–arctan(sqrt(7))

cos(1/2)(b)sin(θ + ulp(θ)) error

444–sqrt(255/256)

sqrt(63/64)
sqrt(15/16)(a)tan(θ + ulp(θ)) error

arccos(x + ulp(x)) error

2

2

64

x ≥ cos8

arcsin(x + ulp(x)) error–2/3 π–1/3 π–2/3 π(d)(e)

Figure 7. Formulas and plots of the badlands for trigonometric 64-bit floating-point functions.

Q. Why not just use f* instead of f?

A. It’s almost certainly hundreds to thousands of times
slower. Perhaps more importantly, the initial version of f*,though more accurate than the initial version of f, might not
be accurate enough.
One question without an easy answer is “What kind of
accuracy should I aim for?” This depends strongly on how
f will be used. Racket’s math library generally aims for lessthan a few ulps error on the entire domain. But library codeis used in many different, unpredictable ways, so such highstandards make sense. Your project might tolerate 10,000
ulps error, or even millions, and almost certainly doesn’ttouch every floating-point function’s entire domain.
Whatever your project is, we recommend reducing thesubdomain on which error is unbounded, if not eliminating
unbounded error altogether.
Debugging the Geometric Distribution
We’ll start debugging with a real task from creating
Racket’s math library. The Geometric distribution’s
www.computer.org/cise			

CISE-16-04-Sciprog.indd 89

Testing and Debugging Support

For testing and debugging floating-point code as described inthis article, a language or library needs■■ 64-bit implementations of ord, ord –1, and ulp;■■ multiprecision floats (MPFR recommended);■■ for each kind of high-precision number, implementations of fl(correctly rounded) and error;■■ plotting on any subinterval of (−∞, ∞) without additional error,

underflow, or overflow; and■■ random float and 64-bit integer generation.

Other high-precision numbers, such as exact rationals anddouble-doubles, are a plus.

	89

28/07/14 7:02 PM

Scientific Programming

6

6x108

4

4x108

2

2x108

0
1

1x10–8
8x10–90.8
0.6

6x10–9

0.8

xis

0.4

ya

xis

ya

1
0.6
0.4

0.2
0 0

0.2

1
0.8

4x10–90.6
0.4

2x10–9

xis

xa

0.2
0 0

(a)

xis

xa

(b)

40
1.5

30

1

20

0.5

10

0
1

0
8x10–240.8
0.6

6x10–24

0.8
0.6
0.4

0.2
0 0

0.2

x

xis

0.4

ya

xis

ya

1

0.6

2x10–24

s
axi

(c)

1
0.8

4x10–24
0.4
0 0

0.2

xis

xa

(d)

Figure 8. Plots of (fun2d-error f f*) at two debugging stages, in different domains. (a) Initial implementation, entire domain.
(b) Initial implementation, approaching y = 0. (c) Implementation using log1p, entire domain. (d) Implementation using log1p, without
updating the bigfloat implementation to use log1p, near y = 2 · 10 −24.

implementation is partly defined in terms of f : (0, 1) ×
(0, 1) → , defined byf (x, y ) =log(1− y ).(9)log( x )Following the test-debug procedure, we first define twostraightforward implementations: f, for 64-bit floats, andf*, for arbitrary-precision bigfloats:(: f (−> Flonum Flonum Flonum))(define (f x y)(/ (fllog (− 1.0 y))(fllog x)))90

CISE-16-04-Sciprog.indd 90

(: f* (−> Bigfloat Bigfloat Bigfloat))(define (f* x y)(bf/ (bflog (bf− (bf 1.0) y))(bflog x)))We start testing f by plotting (fun2d-error f f*) on
the entire domain. Figure 8a shows the results, which aren’tterrible so far. We apparently have error less than 4 ulps fory > 0.1 or so. The plot only suggests that error increases as yapproaches zero.
But really, how bad could it possibly get?
There’s effectively no bound on how bad it can get.
Figure 8b plots the error for y ≤ 10−8. Instead of the
maximum error at just under 8 ulps, in this plot, error is justJuly/August 201428/07/14 7:02 PM

under 8 . 108 ulps. Further plots would show that it only
gets worse closer to y = 0. Ugh.
Is this normal? Unfortunately, yes. It’s utterly typicalfor a straightforward floating-point implementation of areal function to have low error on most of its domain andunbounded error on a small subdomain.
There must be some badlands-wandering going on.
Knowing that f’s error is high when y is near 0.0, it’s not
hard to find the problem: if y is near 0.0, then (− 1.0 y)
is approximate and near 1.0, which is in log’s badlands, so(fllog (− 1.0 y)) has high error.
What about the rest of f’s definition? The outer division
doesn’t magnify error when its arguments are normal.
(We’ll deal with subnormal arguments only if testing turns
up problems.) For a different reason, (fllog x) should be
fine: x is assumed exact because it’s a function argument, so
it doesn’t matter if x is in log’s badlands.
To fix (fllog (− 1.0 y)), we need a log-like function
with different badlands, such as log1p(x) = log(1 + x),
whose badlands are shown in Figure 6. Replacing log(1 − y)in the definition of the real function f with log1p(−y)
wouldn’t change f. We could similarly replace (fllog(− 1.0 y)) in the definition of the floating-point function
f with (fllog1p (− y)), but to change f in a way thatchanges only its floating-point error.
Using fllog1p will definitely reduce error: y is exact andnegation is exact, so even if (− y) is in log1p’s badlands, it
won’t matter. So we change f’s definition to(: f (−> Flonum Flonum Flonum))(define (f x y)(/ (fllog1p (− y))(fllog x)))and plot its error. Figure 8c shows the result. It looks great—error is apparently bounded by 2 ulps—but we’re not done.
We still need to update f*. If (bf− (bf 1.0) y) is in
log’s badlands, then (bflog (bf− (bf 1.0) y)) at least
quadruples error, even with bigfloats. Because the amountthat error can be magnified is unbounded, at some point,
log’s badlands could make f* less accurate than f! In thiscase, it happens near 2e−24, as shown in Figure 8d.
Instead of updating f*, we could increase the numberof bits in the bigfloats that f* computes with. In this case,it takes bigfloats with 1,074-bit significands to make f*more accurate than f again! Clearly, using more precision
isn’t always the best approach. Changing f* to avoid log’s
badlands the same way f does is more reliable and resultsin faster high-precision code, which means faster and morereliable testing.
Whether by using more significand bits or changing
f*, we would find through more testing that our updated
implementation of f has less than 2 ulps error everywhere.
www.computer.org/cise			

CISE-16-04-Sciprog.indd 91

Debugging Log of Quotient
The function
h(x, y) = log(x/y)(10)comes up often when dealing with logarithms, and has anannoyingly large subdomain for which a straightforward
implementation has high error. Accurate implementation
requires a common technique we didn’t use in the last example:
splitting up the domain. For simplicity, we’ll restrict h to thefirst quadrant; that is, h : [0, ∞) × (0, ∞) → .
Figure 9a plots the error in the straightforward 64-bit
float implementation h against the bigfloat implementation
h*. The ridge in the middle suggests that error gets up to atleast 30 ulps near the diagonal x = y.
But really, how bad could it possibly get?
This time, we won’t dignify that silly question with ananswer.
The problem can’t be division, which just adds rounding
error. The problem must be applying floating-point log to therounded quotient. According to Figure 6, we can avoid log’s
badlands by not applying it if exp ( −1 4) ≤ x /y ≤ exp (1 4),
but we’ll need a reformulation on that subdomain.
An obvious reformulation to try is log(x/y) = log(x) −log(y), but this is doomed to fail in the same way. The error
plot looks similar, but has a fatter diagonal ridge. Because(fllog x) and (fllog y) are approximate, subtracting
them results in high error when they’re within a factor of
two, which happens when x and y are near the diagonal.
Again, we’ll turn to log1p. We need to computelog( x /y ) = log(1 + x /y − 1)
= log 1p( x /y − 1)
 x − y . (11)= log 1p  y 
At this point, we could find the conditions under which(x − y)/y is in log1p’s badlands and determine whether theyintersect exp (−1 4) ≤ x /y ≤ exp (1 4) . Instead, we’ll justcode it up:(: h (−> Flonum Flonum Flonum))(define (h x y)
(define q (/ x y))(cond [(<= (flexp −0.25) q (flexp 0.25))(fllog1p (/ (− x y) y))][else
(fllog q)]))Then we’ll test it, as shown in Figure 9b. We apparently
have error bounded by about 2 ulps. Of course, we should
test it more thoroughly than we can demonstrate in thisarticle.
	

91

28/07/14 7:02 PM

Scientific Programming

25
20
15
10
5
0
10

1.5
1
0.5
0
10
8

8
6

6

8

4

xis

xis

ya

ya

10
6
4

2
0 0

2

4

6
4

2

is

x ax

(a)

10
8

0 0

2

is

x ax

(b)

Figure 9. Plots of (fun2d-error h h*) before and after debugging. (a) Initial implementation and (b) debugged implementation.

We’re not done debugging h yet, but we stress that, at
this point, it’s good enough for most purposes. It’s also much
more accurate than it was before, with error bounded byabout 2 ulps near the diagonal instead of unbounded. It’s still
fairly simple.
But suppose we’re writing library code, or h’s callers
need it to be extremely solid, or someone using h hasreported there’s occasionally a problem with it. Let’ssee how robust we can make h using just testing anddebugging.
The first problem we’ll find is that it’s hard to find aproblem, at least by generating plots. We need a way tocreate unusual arguments. For this, we turn to random
sampling. One effective way to generate unusual com­binations of floats is to generate uniformly random ordinalindexes and convert them to floats using ord–1. A Typed
Racket generator that chooses floats this way between a andb (excluding b) is(: random-flonum (−> Flonum Flonum Flonum));; Returns a random flonum in [a,b)(define (random-flonum a b)
(ordinal−>flonum(random-integer (flonum−>ordinal a)(flonum−>ordinal b))))

We can now use (random-flonum −max.0 +inf.0)
to sample from all finite 64-bit floats, where −max.0 isdefined by the math library as the largest negative float(approximately −1.8·10308) and +inf.0 is the next float
92

CISE-16-04-Sciprog.indd 92

after +max.0. This tends to return very large and verysmall numbers because the exponent field is uniformly
distributed.
It takes only a few random samples to find an x and y forwhich ((fun2d-error h h*) x y) is infinite. We’ll definethem temporarily in the interactions window to keep fromhaving to retype them:> (define x 5.295588382145396e+229)> (define y 7.503806844818842e−122)

We get infinite error because (h x y) is infinite but (h*x y) isn’t:> (h x y)+inf.0> (bigfloat−>flonum(h* (bf x) (bf y)))807.858831263043

Manually computing (h x y) reveals what happened:> (/ x y)+inf.0> (fllog +inf.0)+inf.0The exact quotient of x and y is too large tobe approximated by a float. We say that (/ x y)
overflows.
July/August 201428/07/14 7:02 PM

Similarly, when an exact result is too close to zero, it
underflows to 0.0 or −0.0. That happens, for instance, if weswap x and y:> (/ y x)
0.0

Fortunately, the reformulation log(x/y) = log(x) − log(y)seems to work with low error for these arguments, and would
keep h from prematurely overflowing or underflowing:> (− (fllog x) (fllog y))807.8588312630429

Randomized testing finds another, less frequent problem:
when (/ x y) is subnormal, (fllog (/ x y)) has high error.
Using (− (fllog x) (fllog y)) may solve this problem
as well because it doesn’t use division. More testing should
reveal whether (fllog x) and (fllog y) are in subtraction’s
badlands.
The definition of h is now(: h (−> Flonum Flonum Flonum))(define (h x y)
(define q (/ x y))(cond [(<= (flexp −0.25) q (flexp 0.25))(fllog1p (/ (− x y) y))][(or (<= q +max-subnormal.0)(= q +inf.0))(− (fllog x) (fllog y))][else
(fllog q)]))where +max-subnormal.0 is the largest subnormal 64-bit
float, or the float just before 2−1022. Racket takes abouta minute to apply this version of h to 1 million random
argument pairs and compare its outputs to the high-precision outputs of h*. The largest reported error is less than
2.1 ulps.
When updating h* to match h, it’s not necessary to test
whether q overflows or underflows. MPFR’s multiprecision
floats can have much larger exponents than 64-bit floats, sothey don’t overflow or underflow as easily.
Debugging 3D Dot Product
What if we can’t find an algebraic trick to get out of thebadlands? When this happens, sometimes the only goodanswer is to use more precision. Of course, we can use
the high-precision implementations we use for testing,but they’re usually quite slow. Fortunately, on systems
with correctly rounded arithmetic—which is the lowest
possible requirement we can think of for a scientificcomputing platform—there are faster high-precision
options.
www.computer.org/cise			

CISE-16-04-Sciprog.indd 93

Extended-Precision Floats

Some languages provide access to the hardware’s 80-bit
floats, which system libraries use to implement 64-bit­functions such as pow(x, y) = x y accurately in software. Whilean ­attractive alternative to other higher-precision numbers,
80-bit floats ­usually don’t seem to have enough precision to
make up for wandering into the badlands.
Some compilers silently replace intermediate 64-bit resultswith 80-bit results, which isn’t compliant with IEEE 754-2008.
In practical terms, these compilers break important error-reducing
64-bit algorithms such as Kahan summation and double-double
arithmetic, and ruin portability and scientific repeatability. If
the language you use for scientific computing does this, submit
a bug report.

We’ll debug the dot product to demonstrate some
practically accurate, efficient options. We’ll use the definition
dot(x1, y1, z1, x 2 , y2 , z 2)(12)		= x1 . x 2 + y1 · y2 + z1 . z 2.
This might seem perfectly innocuous, but accurate
implementation is beastly. The first problem is that if
the products x1 . x 2, y1 . y2 and z1 . z 2 don’t all have thesame sign, the implementation can wander into addition’s
badlands. The second problem is premature overflow. For
example, even if fl(x1 . x 2 + y1 . y2) < ∞, we still might
get fl(x1 . x 2) = ∞ as an intermediate result, if y1 . y2 isabout as large as x1 . x 2 and has the opposite sign.
As before, we start by coding up straightforward implementations dot and dot*. For dot*, we can use exact rational arithmetic instead of bigfloats. If we do, we won’t haveto update dot* to avoid any badlands.
We can’t plot error on a domain with six dimensions, so
we’ll have to settle for random sampling. To avoid dealing
with overflow at first, we’ll sample between −1 and 1.
We’ll have our generator flip a coin to determine whether
to sample uniformly from the floats in [−1, 1] to represent
typical arguments or to sample uniformly from the encodingsof floats in [−1, 1] to represent atypical arguments:(: random-dot-flonum (−> Flonum))(define (random-dot-flonum)(if (< (random) 0.5)
(random-flonum −1.0 (flnext 1.0))(* 2.0 (− (random) 0.5))))Using random-dot-flonum to generate arguments for
dot and dot*, we occasionally get errors in the thousands ofulps. By coding up the addition badlands formula in Figure 6,
	93

28/07/14 7:02 PM

Scientific Programming

About That Folk WisdomLooking to use some floating-point folk wisdom in debugging
your project? We now have the tools to evaluate it.

Test nearness instead of equality. Wise, if done right. Floatingpoint epsilon is defined as ε = ulp(1), the distance between 1
and the next float. If x is a normal float, then |x| . ε is between
ulp(x) and 2 . ulp(x). If |x − y| ≤ |x| . (4 . ε), then y is within 4
to 8 ulps of x.
Use the library’s hypotenuse function. Often wise. A libraryfunction should avoid premature overflow from squaring, and may
be correctly rounded.
Use (x − y) . (x + y) instead of x 2 − y 2. Wise. Half thetime, x 2 − y 2 wanders into subtraction’s badlands with two
approximations. But (x − y) . (x + y) multiplies approximations,which doesn’t magnify their error.
Compute sums in some specified order or grouping. Not usually wise.
No ordering or grouping strategy is best in all cases. Use Kahansummation (http://en.wikipedia.org/wiki/Kahan_summation_
algorithm), double-doubles, or a specialized summation algorithm.1

Reference
1. J. Demmel and Y. Hida, “Accurate and Efficient Floating Point
Summation,” SIAM J. Scientific Computing, vol. 25, no. 4,
2003, pp. 1214–1248.

we start writing a version of dot that avoids wandering into
addition’s badlands with a rounded result:(define (dot x1 y1 z1 x2 y2 z2)(define x (* x1 x2))
(define y (* y1 y2))(cond [(<= 0.5 (− (/ x y)) 2.0)
<compute-dot-somehow>][else(define d (+ x y))(define z (* z1 z2))
(if (<= 0.5 (− (/ d z)) −2.0)<compute-dot-somehow>(+ d z))]))

We could replace <compute-dot-somehow> with codethat converts the arguments to exact rationals, applies dot*to them, and converts the exact result to a float. This works,
and testing the updated dot suggests that its error is boundedby about 2.5 ulps. This is great, but it’s very slow: about 40
times slower than the initial dot on random arguments. If we
94

CISE-16-04-Sciprog.indd 94



try using bigfloats instead of exact rationals, we’ll find thatthe resulting code is hardly faster.
If we could somehow get exact multiplication out of 64bit floats, we’d be set. It turns out that there’s a way to dothat, using two 64-bit floats to represent one real number.
These are often called double-doubles. Without going into
details, a double-double is a pair of 64-bit floats, that eachhave a 53-bit significand, whose exact sum has 107 bits
precision. Operations on double-doubles work correctly only
on systems with correctly rounded arithmetic.
(Of course, 53 + 53 = 106. The second float’s sign bit
carries information equivalent to one significand bit.)
Racket’s math library exports functions that operate on
double-doubles. Those that only return double-doubles havethe suffix /error (read “with error”). For example,> (fl*/error 1.1 2.1)
2.3100000000000005-2.1316282072803005e-16

The exact sum (not the floating-point sum) of these
returned float values is the exact product of 1.1 and 2.1.
Alternatively, the first float value is an approximation, andthe second is its signed rounding error.
Math library functions that accept double-doubles asarguments have the prefix fl2; for example, fl2+ accepts
four arguments and returns two results. We can use fl2+
and fl*/error to define a high-precision, but reasonably
fast, version of dot, as follows:(define (slowish-dot x1 y1 z1 x2 y2 z2)(define-values (xhi xlo) (fl*/error x1 x2))
(define-values (yhi ylo) (fl*/error y1 y2))(define-values (zhi zlo) (fl*/error z1 z2))
(define-values (dhi dlo) (fl2+ xhi xlo yhi ylo))(define-values (ehi elo) (fl2+ dhi dlo zhi zlo))
ehi)

Replacing <compute-dot-somehow> with (slowish-dot
x1 y1 z1 x2 y2 z2) in the definition of dot yields animplementation that’s only about four times slower than theinitial implementation on random arguments, with error
again bounded by about 2.5 ulps.
On the system we use for testing, the updated dot
computes 25 million dot products per second instead of theinitial implementation’s 100 million per second. This mightseem like a bad deal, but there are other considerations. First,
in scientific computing, we should worry more about ourprograms being faithful to our models than about computation
speed. Second, so much other code usually surrounds floatingpoint computations—for example, code that shuffles dataaround, indexes arrays, or decides what to compute—that
taking even a fourfold speed hit hardly matters.
July/August 201428/07/14 7:34 PM

If we’re willing to tolerate more error,
it’s possible to get good performanceusing the exact-rational dot*. We can giveourselves a speed versus accuracy knob toturn by implementing
1≤ −x / y ≤ t (13)t

Build a parallel
high-precision
implementationfrom the bottom
up. Rememberthat every small
improvement canhelp a lot.

instead of the original addition badlands
formula. Setting the “knob” t to 2 yields
the original formula. With t = 1.1, using
the exact-rational dot* as a fallback makes
dot about nine times slower than the initial
implementation and bounded by about 12 ulps error. Witht = 1.01, dot with exact-rational fallback is only four times
slower but is bounded by about 300 ulps error, according torandomized testing.
Our current implementation might prematurely overflow when given arguments outside of about [−10150, 10150 ].
We’ll leave debugging premature overflow as an exercise.
A simple solution punts to dot* when dot is given large­arguments or returns an infinity. A more complicated solution scales the arguments by powers of two, which is exact
as long as it doesn’t overflow or underflow.

It’s daunting to think of the millions of lines of floatingpoint code that have never been tested and debugged, simply for lack of tools and a bit of know-how.
It’s probably more daunting to consider testing anddebugging your own project, which is certainly muchlarger than the three small examples we’ve shown here. It
can be done, and testing and debugging scales reasonably
well. Start small and work toward larger functions. Movetroublesome formulas into their own functions to debug
them on their own. Document their badlands and avoid
them. Build a parallel high-precision implementation fromthe bottom up. Remember that every small improvement
can help a lot.
Also remember to step back and look at the big picture:
you should test and debug your research question, too.
Floating-point data is an approximation of real data, withup to half an ulp error. This data is often an approximation
itself, with potentially much more error, of some real-life
quantity. Your functions and algorithms have badlands. If
the data is in the badlands and the computed answers arenumbers, these two approximations could make all of theirdigits nonsense.
Fortunately, testing a research question can be done
similarly to plotting log’s error sensitivity in Figure 5. Repeatedly
run the code that answers the question on randomly perturbed,
representative data. Are the answers’ errors millions of times
www.computer.org/cise			

CISE-16-04-Sciprog.indd 95

larger than the perturbations? Is some ofthe data nonzero but the answers near zero?
Is some of the data far from zero and in anarea of exponential growth? If the answer
to any of these questions is “yes,” your data
might be in the badlands. If so, either yourdata and code need more precision, in whichcase you’ll be very glad to have a parallel highprecision implementation, or you need to ask
a different question.
We’ll end with a dirty secret: almost
everything we presented and demonstrated
is error analysis in disguise. The condition
number formula, for example, is straightout of any numerical methods textbook. But we replaced
the burden of proof with testing, which reduces the burden
on you.

References
1. M. Flatt, “Reference: Racket,” tech. report PLT-TR-2010-1,
PLT Inc., 2010; http://racket-lang.org/tr1.
2. L. Fousse et al., “MPFR: A Multiple-Precision Binary
Floating-Point Library with Correct Rounding,” ACM
Trans. Mathematical Software, vol. 33, no. 2, 2007, pp.
13:1–13:15.
3. “IEEE Standard for Floating-Point Arithmetic,” IEEE Std
754-2008, Aug. 2008, pp. 1–70.
Neil Toronto recently graduated with a PhD from Brigham
Young University, and is a postdoctoral researcher at Universityof Maryland, College Park. His research focuses on programming
language support for mathematical computation, and currently
emphasizes Bayesian modeling and inference. Contact him at neil.
toronto@gmail.com.
Jay McCarthy is a visiting assistant professor of computer science at
Vassar College and formerly an assistant professor at Brigham Young
University. His research focuses on the nexus of systems, verification,and security from the perspective of applied programming language
techniques. McCarthy has a BS in mathematics from the University of Massachusetts, Lowell, and a PhD in computer science from
Brown University. Contact him at jay.mccarthy@gmail.com.

Selected articles and columns from IEEE Computer Society
publications are also available for free at http://ComputingNow.computer.org.
	

95

28/07/14 7:02 PM