Programming

that's one of the aspects of the current implementation

I didn't have the patience to do it myself bit wanted to see just how complex it would get:

fp32_mul() {
    local a=$1 b=$2
    local sa=$(( (a >> 31) & 1 ))
    local sb=$(( (b >> 31) & 1 ))
    local sign=$((sa ^ sb))

    local ea=$(( (a >> 23) & 0xff ))
    local eb=$(( (b >> 23) & 0xff ))
    local fa=$(( a & 0x7fffff ))
    local fb=$(( b & 0x7fffff ))

    # NaN / infinity / zero handling
    if (( ea == 255 )); then
        if (( fa != 0 )); then
            printf '%08x\n' $((0x7fc00000))
            return
        fi
        if (( eb == 0 && fb == 0 )); then
            printf '%08x\n' $((0x7fc00000))   # inf * 0 = NaN
            return
        fi
        printf '%08x\n' $(((sign << 31) | 0x7f800000))
        return
    fi

    if (( eb == 255 )); then
        if (( fb != 0 )); then
            printf '%08x\n' $((0x7fc00000))
            return
        fi
        if (( ea == 0 && fa == 0 )); then
            printf '%08x\n' $((0x7fc00000))
            return
        fi
        printf '%08x\n' $(((sign << 31) | 0x7f800000))
        return
    fi

    if (( ea == 0 && fa == 0 || eb == 0 && fb == 0 )); then
        printf '%08x\n' $((sign << 31))
        return
    fi

    # Convert subnormals to a normalized significand/exponent.
    # m is a 24-bit significand for normals.
    local ma mb
    if (( ea == 0 )); then
        ma=$fa
        ea=1
        while (( (ma & 0x800000) == 0 )); do
            ma=$((ma << 1))
            ((ea--))
        done
    else
        ma=$((fa | 0x800000))
    fi

    if (( eb == 0 )); then
        mb=$fb
        eb=1
        while (( (mb & 0x800000) == 0 )); do
            mb=$((mb << 1))
            ((eb--))
        done
    else
        mb=$((fb | 0x800000))
    fi

    # Multiply the two 24-bit significands.
    # Product is up to 48 bits.
    local p=$((ma * mb))
    local e=$((ea + eb - 127))

    # Normalize product.
    #
    # ma*mb has binary point after bit 46.  If bit 47 is set,
    # product is [2,4), otherwise [1,2).
    local shift
    if (( p & 0x800000000000 )); then
        shift=24
        ((e++))
    else
        shift=23
    fi

    # Extract 23 fraction bits plus guard/round/sticky information.
    local frac=$(( (p >> shift) & 0x7fffff ))
    local guard=$(( (p >> (shift - 1)) & 1 ))
    local round=$(( (p >> (shift - 2)) & 1 ))
    local sticky=0

    if (( shift >= 3 )); then
        local mask=$(( (1 << (shift - 2)) - 1 ))
        (( (p & mask) != 0 )) && sticky=1
    fi

    # Round-to-nearest, ties-to-even.
    if (( guard && (round || sticky || (frac & 1)) )); then
        ((frac++))
        if (( frac == 0x800000 )); then
            frac=0
            ((e++))
        fi
    fi

    # Overflow -> infinity.
    if (( e >= 255 )); then
        printf '%08x\n' $(((sign << 31) | 0x7f800000))
        return
    fi

    # Normal result.
    if (( e > 0 )); then
        printf '%08x\n' $(((sign << 31) | (e << 23) | frac))
        return
    fi

    # Underflow into the subnormal range.
    #
    # At this point the normalized significand represented by
    # (1.frac) must be shifted right by 1-e positions.
    local mant=$((0x800000 | frac))
    local rshift=$((1 - e))
    local lost=0
    local halfway=0
    local low=0

    if (( rshift >= 25 )); then
        # Everything rounds to zero (unless the exact value is
        # sufficiently close, which it cannot be here).
        mant=0
    else
        low=$((mant & ((1 << rshift) - 1)))
        mant=$((mant >> rshift))

        halfway=$((1 << (rshift - 1)))

        if (( low > halfway || (low == halfway && (mant & 1)) )); then
            ((mant++))
        fi
    fi

    # Rounding a subnormal can produce the smallest normal.
    if (( mant >= 0x800000 )); then
        printf '%08x\n' $(((sign << 31) | (1 << 23)))
    else
        printf '%08x\n' $(((sign << 31) | mant))
    fi
}