Как обнаружить / исправить ошибки округления в этой 32-битной реализации умножения FP?

Я пытаюсь понять и реализовать 32-битное умножение с плавающей запятой. Я использую формат FP с одинарной точностью IEEE 754, в соответствии с которым:

Итак, я следую этому пошаговому руководству по алгоритму. Используя bitstream.pack, я могу преобразовать число с плавающей запятой в строку битов в формате IEEE 754, поэтому я делаю это и сравниваю его с результатом, который я получаю, когда пытаюсь вручную реализовать алгоритм.

Однако похоже, что у меня возникают проблемы с округлением как мантиссы, так и экспоненты примерно в 50% случаев. Может ли кто-нибудь увидеть, что я здесь неправильно реализую? Если нет, есть ли способ обнаружить и исправить эти ошибки округления?

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)

И вот мои результаты:

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