This innocent-looking function features two splendid numerical idioms of Paranoia, one in this section, one in the next.
The task is simple enough. In a six-digit decimal system, the value of
the an ulp of 1.0 is just \(10 ^ {-5}\).
The value is easily obtained by a look at the specification of the number system. In the IEEE 754 single format with 24 significant bits, an ulp of one is \(2 ^ {-23}\).
Actually computing an ulp of 1.0 without knowledge of the radix
or precision is another matter. Our strategy is to compute
\(1 / 3\) or \(2 / 3\) of an ulp exactly and then triple
the value. This guarantees we have exactly one or
two ulps of 1.0. Poor arithmetic
may give us several extra ulps, but that's no problem.
In the second half of the function we'll walk down to exactly
one ulp.
We start with the value \(4 / 3\). We know the decimal expression \(1.3333\cdots\) and in a moment we'll see that the value in binary is \(1.010101\cdots _{2}\). The value must round in any radix not a multiple of \(3\).
What matters is that the rounding error is one or two thirds of an ulp
of 1.0, perhaps plus one or more further ulps of
1.0 if the division
is especially sloppy. To see this, look at the decimal number line:
The tick marks signify the (highly magnified)
representable numbers in each decade.
Their spacing defines ulps of the leftward bounding
power of \(10\), which is \(1.0\) here.
The rounding error
in 4.0 / 3.0 is at least \(0.3333\cdots\)
or \(0.6666\cdots\)
ulp of 1.0, depending on which representable
value is computed.
Any addtional error is in whole ulps.
In decimal, the error is \(-1 / 3, -4 / 3, -7 / 3, \ldots\)
or \(+2 / 3, +5 / 3, +8 / 3, \ldots\)
ulps of 1.0.
That is, there's a guaranteed one or two thirds ulp in the rounding error.
The same logic applies to binary, octal, and hex arithmetic.
After 4.0 / 3.0 is rounded, the rest of the setup is done
in what should be exact arithmetic. The one or two thirds
of an ulp error are baked into the following exact calculations.
Subtracting 1.0 is exact, leaving the variable
third with the value
\(1/3 + \epsilon\).
This remarkable value captures two interesting but unrepresentable
numbers in one machine number. We know that \(\epsilon\) is
exactly \(1 / 3\) of an ulp of 1.0, perhaps with
additional ulps.
The first numerical trick featured in this function
is the use of artful cancellation
to isolate \(\epsilon\) and triple
it, giving us our ulp of 1.0, exactly.
The value third lies in the range of numbers
just less than 1.0, where the numbers are B
times as dense as just above 1.0. Our next operations
will be exact. The code computes the value
x given exactly by
\[ | (1/2 - (1/3 + \epsilon) \times 2) - (1/3 + \epsilon) | \]
which simplifies to \( | 3 \epsilon |\).
We have isolated and tripled \(\epsilon\).
The value x is a
small number of ulps of 1.0, exactly.
The diagram below shows the case of binary arithmetic
with an even number of bits.
Note the fun fact that in this case, unlike decimal,
the nereast value to \(4 / 3\) is gotten
by rounding up \(1 / 3\) of an ulp.
The lowest-order bits of values are not shown as the subtractions proceed.
They're all zero, by design, in this isolation of the ulps of 1.0.
This section gives some optional mathematical elaboration. We rewrite \(4/3\) to expand it in binary. \[ 4/3 = {1 \over {1 - {1/4}}} = 1 + 1/4 + 1/16 + 1/64 + 1/256 + \cdots \] which is the binary expression \(1.010101... _{2}\) above.
Similary, we rearrange \(4/3\) to the form \[ 4/3 = -2 + 3 \times 10/9 = - 2 + 3 \times {1 \over {1 - {1/10}}} = 1 + 3/10 + 3/100 + \cdots \] which gives the familiar expansion \(1.3333\cdots\) in decimal.
While we're here, we can observe that \(3 \times 0.3333\cdots = 0.9999\cdots\) has the value \(1\), which can cause comfusion on decimal calculators with finite precision. In the spirit of the previous examples, we write \(1 = 9/10 \times 10/9\) and expand \[ 10/9 = {1 \over {1 - 1/10}} = 1.1111\cdots \] So \(1 = 9/10 \times 1.1111\cdots = 0.9999\cdots\) as expected from \(3 \times 1/3\).
The code chooses the larger of the computation and the input guess, trusting the iteration to drive down to just one ulp.
def find_ulp_of_one_plus(guess):
"""Compute u such at 1+u is the next larger value than 1."""
x = FOUR / THREE # rounds in ulps of 1
third = x - ONE
sixth = ONE_HALF - third
x = sixth + sixth
x = fabs(x - third) # sums two rounding errors
if x < guess:
x = guess
# Now x = (unknown no.) ulps of 1+.
The code continues in the next section.
Exercise: Replay the binary case with odd
precision to see that
it also proudces an ulp of 1.0.
Exercise: The diagram below starts the computation in 6-digit decimal arithmetic.
Can you complete the computation of x?
Exercise: This diagram
illustrates the computation of x
in 6-digit octal arithmetic. We assume chopped arithmetic,
as found on Burroughs computers of the 1960s.
Now rerun the computation in 7-digit octal.