find_ulp_of_one_plus() – the iteration

So far, we have computed a value x with one or more ulps of 1.0. Next, we exploit a powerful numerical idiom of Paranoia:

When searching for the smallest value with a given property, start with a handy value at least as big and iterate downward.

The iteration exhibits another useful technique: the creation of nonzero low-order digits to encourage an upward round, if there's nearly a half-way case. The binary diagram below shows how u is halved, but with some added nonzero low-order bits. Since u is some ulps of \(1.0\), u * u is some ulps of u. The factor THIRTY_TWO shifts those bits well back into u in radix 2, 8, or 16. The diagram continues the binary tale.

The expression x + ONE rounds to 1.0 plus some ulps of one. Subtracting ONE is exact, leaving just the ulps. Successively halving the value walks it down to a single ulp.

    while True:
        u = x
        x = ONE_HALF * u + THIRTY_TWO * u * u
        x = x + ONE
        x = x - ONE
        if u <= x or x <= ZERO: break
    return u

Your turn

Exercise: Complete the binary case started in the diagram above 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.