mult_exp.py

# this was generated by ChatGPT and I'm not sure if it is correct.
# It looks OK but errors might be subtle.
#
# the wiki states:
# Mantissa: (m4 >= 128 ? -1 : +1) * ((m4 | 0x80) >> 8 + m3 >> 16 + m2 >> 24 + m1 >> 32)
# as a C-language like expression, with "x >> y" as "float multiply x by 2↑(-y)" (right-bit-shift operation)

import struct
import numpy as np

def cbm_float_to_python_float(hex_bytes):
    # Ensure we have exactly 5 bytes
    if len(hex_bytes) != 5:
        raise ValueError("CBM float requires exactly 5 bytes")

    if hex_bytes[0] == 0:
        return 0

    # Extract the exponent and bias it by subtracting 128
    exponent = hex_bytes[0] - 128

    # Extract the mantissa as a 4-byte integer
    mantissa = (hex_bytes[1] << 24) | (hex_bytes[2] << 16) | (hex_bytes[3] << 8) | hex_bytes[4]

    # The sign bit is in the MSB of the mantissa; 0 means positive, 1 means negative
    sign = -1 if (mantissa & 0x80000000) else 1

    mantissa |= 0x80000000 # normalize mantissa to be in [0.5;1)

    #print(bin(mantissa), hex(mantissa))

    # Convert mantissa to a float by dividing by 2**31 (to match the fixed point range)
    float_value = sign * (mantissa / (2**32)) * (2 ** exponent)


    return float_value


def cbm_mult(hex1, hex2):
    return cbm_float_to_python_float(hex1) * cbm_float_to_python_float(hex2)


def _ieee_float_32(sign, exponent, mantissa):
    return (
        ((-1) ** sign) *
        (2**(exponent - 127)) *
        (1+sum([((mantissa >> (23-i)) & 1) * (1/2**i) for i in range(0, 24)]))
    )


def python_float_to_cbm_float(f):
    f32 = struct.unpack('!I', struct.pack('!f', f))[0]
    sign = (f32 & 0x80000000) >> 31
    exponent = (f32 & 0x7f800000) >> 23
    mantissa = f32 & 0x007fffff

    extracted_float = _ieee_float_32(
        sign,
        exponent,
        mantissa,
        )

    #print(sign, exponent, mantissa)
    assert abs(extracted_float - f) < 1e-3, extracted_float

    mantissa_bytes = [
        (mantissa >> 16) & 0xff,
        (mantissa >>  8) & 0xff,
        (mantissa >>  0) & 0xff,
    ]

    #print(f'mantissa: {mantissa:024b}')
    #print([f'{n:08b}' for n in mantissa_bytes])

    mantissa_cbm = (
        sign << 31 |
        mantissa_bytes[0] << 24 |
        mantissa_bytes[1] << 16 |
        mantissa_bytes[2] << 8
    )

    #print(bin(mantissa_cbm))

    hex_bytes = [
        exponent + 2, # TODO check rebiasing
        mantissa_cbm >> 24 & 0xff,
        mantissa_cbm >> 16 & 0xff,
        mantissa_cbm >>  8 & 0xff,
        mantissa_cbm >>  0 & 0xff,
    ]

    return hex_bytes


z = cbm_float_to_python_float([0x00, 0x00, 0x00, 0x00, 0x00])
assert z == 0.0, z

z = cbm_float_to_python_float([0x80, 0x00, 0x00, 0x00, 0x00])
assert z == 0.5, z

z = cbm_float_to_python_float([0x81, 0x00, 0x00, 0x00, 0x00])
assert z == 1.0, z

z = cbm_float_to_python_float([0x81, 0x80, 0x00, 0x00, 0x00])
assert z == -1.0, z

z = cbm_float_to_python_float([0x98, 0x35, 0x44, 0x7a, 0x00])
assert z == 11879546.0, z


def approx_mult_naive_py(f1, f2):
    int_cast = lambda x: struct.unpack('!I', struct.pack('!f', x))[0]
    i1 = int_cast(f1)
    i2 = int_cast(f2)
    bias = 0x3f76d000
    mult = i1 + i2 - bias
    return struct.unpack('!f', struct.pack('!I', mult))[0]


def approx_mult_py(f1, f2):  # handles underflow
    int_cast = lambda x: struct.unpack('!I', struct.pack('!f', x))[0]
    i1 = int_cast(f1)
    i2 = int_cast(f2)

    # exponent bias: 0b01111110 (=126)
    # mantissa bias: 0b11101101101000000000000
    bias = 0x3f76d000

    mult = (i1 & 0x7fffffff) + (i2 & 0x7fffffff)
    if mult <= bias:
        mult = 0
    else:
        mult -= bias

    mult |= (i1 ^ i2) & 0x80000000
    return struct.unpack('!f', struct.pack('!I', mult))[0]


def approx_mult_cbm_v2(h1, h2):
    from numpy import uint8, uint16

    cbm_bytes = uint8([0, 0, 0, 0, 0])
    h1 = uint8(h1)
    h2 = uint8(h2)
    overflow = 0

    for i in [4, 3, 2]:
        tmp = h1[i].astype(uint16) + h2[i] + overflow
        overflow = tmp >> 8
        cbm_bytes[i] = tmp

    # this byte has the sign bit
    i = 1
    tmp = (h1[i] & 0x7f) + (h2[i] & 0x7f) + overflow
    overflow = tmp >> 7
    cbm_bytes[i] = (tmp & 0x7f) | ((h1[i] ^ h2[i]) & 0x80)

    # this byte has the exponent
    i = 0
    tmp = h1[i].astype(uint16) + h2[i] + overflow
    if tmp < 128:
        print('AWHHHHH, {:016b}'.format(tmp))
        tmp = 0
    else:
        print('{:016b}'.format(tmp))
        tmp = tmp - 129

    cbm_bytes[i] = tmp

    return [int(n) for n in cbm_bytes]


def approx_mult_cbm_v1(h1, h2):
    from numpy import uint8, uint16

    bias = 0x3f76d00000

    # CBM uses -128 exponent offset instead of -127.
    # in the original we have
    #   exponent = x - 0x3f - 127
    #            = x - 126 - 127
    #            = x -
    #
    # TODO maybe 0x7d is better instead of 0x81?
    bias = 0x8176d00000
    bias_bytes = uint8([
        (bias >> 32) & 0xff,
        (bias >> 24) & 0xff,
        (bias >> 16) & 0xff,
        (bias >>  8) & 0xff,
        (bias >>  0) & 0xff,
    ])
    cbm_bytes = uint8([0, 0, 0 ,0, 0])

    h1 = uint8(h1)
    h2 = uint8(h2)

    print(' h1:', [f'{n:08b}' for n in h1])
    print(' h2:', [f'{n:08b}' for n in h2])
    print('1+2:', [f'{(n+m):08b}' for n, m in zip(h1, h2)])
    print('bia:', [f'{n:08b}' for n in bias_bytes])

    # mantissa
    i = 4
    tmp = h1[i].astype(uint16) + h2[i] - bias_bytes[i]
    overflow = (tmp >> 8) & 1
    cbm_bytes[i] = tmp.astype(uint8)
    i = 3
    tmp = h1[i].astype(uint16) + h2[i] - bias_bytes[i] + overflow
    overflow = (tmp >> 8) & 1
    cbm_bytes[i] = tmp.astype(uint8)
    i = 2
    tmp = h1[i].astype(uint16) + h2[i] + overflow
    overflow = (tmp >> 8) & 1

    if overflow:
        tmp -= bias_bytes[i]
    else:
        tmp = np.uint8(0)

    cbm_bytes[i] = tmp.astype(uint8)
    i = 1
    # this byte includes the sign bit and needs special treatment
    # FIXME underflow is not correctly handled here IMO
    tmp = (h1[i] & 0x7f).astype(uint16) + (h2[i] & 0x7f) + overflow
    overflow = int((tmp & 0x80) > 0)

    if overflow:
        tmp -= bias_bytes[i]
    else:
        tmp = np.uint8(0)

    cbm_bytes[i] = (tmp & 0x7f) | ((h1[i] ^ h2[i]) & 0x80)
    cbm_bytes[i] = tmp.astype(uint8)
    # exponent
    i = 0
    tmp = h1[i].astype(uint16) + h2[i] + overflow
    overflow = (tmp >> 8) & 1  # TODO use?
    cbm_bytes[i] = tmp.astype(uint8)

    # if we see a overflow bit we know that the byte can be subtracted
    # with the bias
    if overflow:
        cbm_bytes[i] -= bias_bytes[i]
    else:
        cbm_bytes[i] = 0

    overflow = (cbm_bytes[i] >> 8) & 1  # TODO use?


    print('cbm:',[f'{n:08b}' for n in cbm_bytes])

    return [int(n) for n in cbm_bytes]


approx_mult_cbm = approx_mult_cbm_v2


def print_error(f1, f2, threshold=0.1):
    error = abs(f1-f2)/(f2+1e-4)
    end = "OK" if abs(error) < threshold else "\033[1m\033[31mFAILED\033[0m"
    end += "\n"
    print(f"result: {f1:.5f}, ref: {f2:.5f} error: {abs(f1-f2)/(f2+1e-4):.3f}: ", end=end)


def test_conversion():
    convert = cbm_float_to_python_float
    is_close = lambda x, y, tol=1e-4: abs(x-y) < tol

    assert is_close(convert(python_float_to_cbm_float(0.5)), 0.5)
    assert is_close(convert(python_float_to_cbm_float(1.0)), 1.0)
    assert is_close(convert(python_float_to_cbm_float(2.0)), 2.0)
    assert is_close(convert(python_float_to_cbm_float(0.25)), 0.25)
    assert is_close(convert(python_float_to_cbm_float(0.33)), 0.33)
    assert is_close(convert(python_float_to_cbm_float(0.1234)), 0.1234)
    assert is_close(convert(python_float_to_cbm_float(0.12345678)), 0.1234578)
    assert is_close(convert(python_float_to_cbm_float(-0.12345678)), -0.12345678)
    assert is_close(convert(python_float_to_cbm_float(11879546.0)), 11879546.0)


def test_mult_python():
    print_error(approx_mult_py(1, 2), 2)
    print_error(approx_mult_py(1.5, 2), 3)
    print_error(approx_mult_py(1.2345, 2), 2.469)
    print_error(approx_mult_py(1.2345, 1.2345), 1.5239)
    print_error(approx_mult_py(0.5, 0), 0)

    print('wide range test')
    for x in np.linspace(-30, 30, 20):
        print(f'{x} * {x} ?= {x*x}')
        print_error(approx_mult_py(x, x), x * x)



if __name__ == "__main__":
    #test_conversion()
    #test_mult_python()

    """
    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(1),
            python_float_to_cbm_float(2),
        )
    ), 2)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(1.5),
            python_float_to_cbm_float(2),
        )
    ), 3)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(1.2345),
            python_float_to_cbm_float(2),
        )
    ), 2.469)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(-1.2345),
            python_float_to_cbm_float(2),
        )
    ), -2.469)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(1.2345),
            python_float_to_cbm_float(1.2345),
        )
    ), 1.5239)

    """

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(1.2345),
            python_float_to_cbm_float(0),
        )
    ), 0)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(0),
            python_float_to_cbm_float(0),
        )
    ), 0)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(0.5),
            python_float_to_cbm_float(0),
        )
    ), 0)

    print('wide range test')
    for x in np.linspace(-30, 30, 100):
        """
        print(f'{x} * {x} ?= {x*x}')
        print_error(cbm_float_to_python_float(
            approx_mult_cbm(
                python_float_to_cbm_float(x),
                python_float_to_cbm_float(x),
            )
        ), x*x)
        """
        print_error(cbm_float_to_python_float(
            approx_mult_cbm(
                python_float_to_cbm_float(3.14159),
                python_float_to_cbm_float(x),
            )
        ), 3.14159*x)

    print('rand test')
    np.random.seed(42)
    for x in np.random.uniform(size=100):
        print_error(cbm_float_to_python_float(
            approx_mult_cbm(
                python_float_to_cbm_float(3.14159),
                python_float_to_cbm_float(x),
            )
        ), 3.14159*x)

    print([hex(n) for n in python_float_to_cbm_float(3.1415)])
    print([hex(n) for n in python_float_to_cbm_float(-10)])

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(3.1415),
            python_float_to_cbm_float(-10),
        )), 3.1415*-10)


    def example_numbers_to_asm(f1, f2, addr=0xc480):
        h1 = python_float_to_cbm_float(f1)
        h2 = python_float_to_cbm_float(f2)

        expected = cbm_float_to_python_float(
            approx_mult_cbm(
                python_float_to_cbm_float(f1),
                python_float_to_cbm_float(f2),
            ))

        print(f"; code for multiplying {f1} and {f2}")
        print(f"; inspect ${addr+5+5:04x} for result")
        print(f"; expected for approx mult: {expected}")
        for i, hi in enumerate(h1):
            print(f"lda #${hi:02x}")
            print(f"sta ${addr+i:04x}")
        print(f"+float_to_fac1 ${addr:04x}")
        for i, hi in enumerate(h2):
            print(f"lda #${hi:02x}")
            print(f"sta ${addr+5+i:04x}")
        print(f"+fmult ${addr+5:04x}")
        print(f"+movmf ${addr+5+5:04x}")


    example_numbers_to_asm(3.1415, 1000, addr=0xc480)

    print_error(cbm_float_to_python_float(
        approx_mult_cbm(
            python_float_to_cbm_float(100 - 0.5),
            python_float_to_cbm_float(1),
        )) + 20, 119.5)