find_ulp_of_one_minus() – the iteration

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

Your turn

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.