How to detect/fix rounding errors in this 32-bit FP multiplication implementation?

Question

I'm trying to understand and implement 32-bit floating point multiplication. I'm using the IEEE 754 single-precision FP format whereby:

So I'm following this walkthrough of the algorithm. Using bitstream.pack I can convert a float to a string of bits in IEEE 754 format, so I'm doing that and comparing it to the result I get when I try to manually implement the algorithm.

However it looks like I'm getting rounding problems on both the mantissa and the exponent about 50% of the time. Can anyone see what I'm implementing incorrectly here? If not, is there a way to detect and correct these rounding errors?

import bitstring, random 

successes = 0
failures = 0

def ieee754(float):
    b = bitstring.pack('>f', float)
    sbit, wbits, pbits = b[:1], b[1:12], b[12:]
    return sbit.bin + wbits.bin + pbits.bin

def extract_parts(bits):
    m = bits[9:]     # 23 bits 0-22 are mantissa
    e = bits[1:1+8]  # 8 bits 23-31 are exponent
    s = bits[:1]     # bit 32 is sign
    return s, e, m 

tests = [(-18.0, 9.5), (134.0625, -2.25), (-14.5, -0.375), (7.5, 15.5), (1.2, 23.1), (-0.5, -0.2)]

for a, b in tests:
    #a = random.uniform(-10, 10)
    #b = random.uniform(-10, 10)

    a_s, a_e, a_m = extract_parts(ieee754(a))
    b_s, b_e, b_m = extract_parts(ieee754(b))

    # sign is exclusive-or of signs
    s = '1' if int(a_s) != int(b_s) else '0'

    # exponent is sum of exponents minus 127 'bias'
    e = int(a_e, 2) + int(b_e, 2) - 127

    # mantissa is product of mantissas with a 1 added as their MSB
    # then we ignore the MSB of the result
    m = '{0:023b}'.format(int('1' + a_m, 2) * int('1' + b_m, 2))[1:24]

    # convert to binary strings for comparison
    e = '{0:08b}'.format(e)

    print("Calculated:\t", (s, e, m));
    print("Expected:\t", extract_parts(ieee754(a*b)))

    if((s, e, m) == extract_parts(ieee754(a*b))):
        print("Success with", a, b); successes += 1
    else:
        print("Failure with", a, b); failures += 1

print("Successes", successes, "Failures", failures)

And these are my results:

Calculated:  ('1', '10000110', '01010110000000000000000')
Expected:    ('1', '10000110', '01010110000000000000000')
Success with -18.0 9.5

Calculated:  ('1', '10000111', '00101101101001000000000')
Expected:    ('1', '10000111', '00101101101001000000000')
Success with 134.0625 -2.25

Calculated:  ('0', '10000000', '01011100000000000000000')
Expected:    ('0', '10000001', '01011100000000000000000')
Failure with -14.5 -0.375

Calculated:  ('0', '10000100', '11010001000000000000000')
Expected:    ('0', '10000101', '11010001000000000000000')
Failure with 7.5 15.5

Calculated:  ('0', '10000011', '10111011100001010010000')
Expected:    ('0', '10000011', '10111011100001010001111')
Failure with 1.2 23.1

Calculated:  ('0', '01111011', '10011001100110011001101')
Expected:    ('0', '01111011', '10011001100110011001101')
Success with -0.5 -0.2

Answer 1

Two (three?) issues that I see.

1. If the product of the significands is greater than or equal to 2, then the exponent is wrong. This explains the off-by-one errors in the exponents.

2. Instead of truncating the product of the significands, you need to apply round-to-even logic. This explains the off-by-one errors in the significands.

3. You're not handling subnormals, infinities, or NaN correctly (but you probably know that).