The code continues im parallel with the case for ulps
just larger than 1.0, but with a twist.
To investigate values immediately less than 1.0,
we look at the expression
\[ 1 - (1 - \epsilon) \]
where \( \epsilon \) is a tiny value, not necessarily a whole
number of ulps of numbers less than one.
The wrinkle is that we must take care not to lose a digit off the right hand side by consuming a digit in the one's place.
The idiom for computing \( 1 - w \) when \( w < 1 \)is to use two operations \[ 1/2 + (1/2 - w) \] avoiding use of the one's palce in any subexpression. The formula above expands to \[ 1/2 + (1/2 - (1/2 + (1/2 - \epsilon))) \] which you will see in four lines of the code below.
In the loop, y is half a whole number of ulps plus
a small nudge. When subtracted from ONE_HALF the result
rounds to \(1/2 - \delta\), where Greek delta is a whole number of
ulps of values just less than 1.0.
The combined expression ONE_HALF + (ONE_HALF - y) is
exactly \(1 - \delta\).
We then apply the same tecchnique to compute
\[ 1 - (1 - \delta) \]
to isolate \( \delta \), iterating until it is down to
a single ulp.
This binary diagram starts with 3 ulps of values just less
than 1.0. It illustrates the computation of
\( 1 - \delta \)
Next, we subtract from 1.0in two steps to isolate
the difference.
As in the parallel case, the loop terminates when the halved value
falls off to zero when rounded or when u is pinned.
while True:
u = x
y = ONE_HALF * u + THIRTY_TWO * u * u
y = ONE_HALF - y
x = ONE_HALF + y
y = ONE_HALF - x
x = ONE_HALF + y
if u <= x or x <= ZERO: break
return u
Exercise: Complete the binary case
started in the diagram above, assuming rounded arithmetic,
to see that it results in
one ulp of 1.0.
Exercise: The decimal diagram below starts the iteration in 6-digit rounded decimal arithmetic.
Continue the loop to find one ulp of 1.0.
Exercise: Complete the computation in 6-digit chopped octal arithmetic.