11.001001000011111101101010100010001000 Arithmazium
Home

compute_C_and_inverse()

When we explore underflow in the code to come, we want a starting point tiny but not too close to the underflow threshold. This function returns the value intended to be assigned to C. We'll not forget that it's C that's tiny, makeing 1/C huge. 

def compute_C_and_inverse():
    """Return untility constant C = 1 / B**k, safely above the underflow
    threshold, and its reciprocal.

    C is a useful starting point for investigation of underflow.
    It has the form 1 / B**k with k chosen so that C is farther than
    a factor of 1 / B**PRECISION of underflowing to zero.
    C is a normal number. It's reciprocal is a large power of B.

    Returns:
        C, 1 / C

    Basic 4370-4400
    """

First, we compute d, typical just an ulp of values less than \( 1 \). When the number of digits is fractional, we expand the last digit to full size. Here is the case for 22 bits in hexadecimal.

ONE    100000 5 1/2 hex digits ULP_ONE_PLUS    000004 last hex digit = 0100 in binary ONE_MINUS_ULP    0fffffc just 22 significant bits in 5 1/2 digits ULP_ONE_MINUS    0000004 last hex digit = 0010 in binary # 1 / 16**ceil(5.5)    0000001 finds the LSB of full last hex digit

The value just computed is an ulp of numbers just less than \( 1 \) in six hex digits. 

    d = ULP_OF_ONE_MINUS  # default value for integer precision
    if PRECISION != floor(PRECISION):
        # Compute 1 / B**ceiling(PRECISION)
        d = ONE_OVER_B
        for _ in range(floor(PRECISION)):
            d = d * ONE_OVER_B

The next loop computes value c that is roughly the square root or perhaps the fourth root of the underflow threshold. This is one quick leap, before the loop that follows walks to the underflow threshold in measured steps.

This code exhibits the Paranoia idiom for detecting underflow. The typical step is to compute z = y * y for a tiny value y. The two cases tested are

  1. z + z <= z which arises when z underflows to zero or to an epsilon value \epsilon of nonzero but indeterminate size that behaves well in comparisons
  2. y <= z which arises when underflow ultimately pins at a nonzero value the compares reasonably with itself

Typically, it's the first case that terminates the loop. In this case, c is left pessimistically far from the underfulow threshold, resulting in extra steps in the next loop. 

    # d = 1 / B**k for k = PRECISION or floor(PRECISION)+1
    # Compute 1/d * 1/d**2 * 1/d**4 * ... until the multiplier
    # tapers off (due to underflow). This is the quick route to
    # a very tiny 1 / B**k in ever-inceasing steps.
    y = ONE
    z = d
    while True:
        c = y
        y = z
        z = y * y
        if y <= z or z + z <= z: break

The second loop approaches the underflow threshold step by step by a factor of d. In an IEEE 754 system, z underflows to zero, y is possibly subnormal, and c is a tiny normal number. In arithmetic that just flushes to zero past a certain point, both y and c will be normal. In any case, c is tiny and normal. 

    # Finally, step toward the minimum normal number by multiples of d.
    # c stays one step behind y, so that at termination
    # 1  >>  c  >  B**(min+PRECISION)  >  y  >  B**mnin  >  z (underflows)
    y = c
    z = y * d
    while True:
        c = y
        y = z
        z = y * d
        if y <= z or z + z <= z: break
    return c, ONE / c

Your turn

Exercise: Compute c for the case of six-digit decimal arithmetic in which any value smaller than \( 10 ^ {-99} \) flushes to zero. As a hint, an ulp of a number just less than one is \( 10 ^ {-6} \).

Exercise: If you have some familiarity with IEEE 754 arithmetic, compute c for the single format, with gradual underflow. An ulp of a number just less than one is \( 2 ^ {-24} \).

Home