math.frexp

[Steven D'Aprano]

I'm experimenting with various implementations of product(values). Here's a
naive version that easily suffers from overflow:

from functools import reduce
from operator import mul

def product_naive(values):
    return reduce(mul, values)


py> product_naive([2, 4, 5])
40
py> product_naive([1e200, 1e200, 1e200])
inf


So I started experimenting with math.frexp, but I'm having trouble
understanding what I'm doing wrong here.


Help on built-in function frexp in module math:

frexp(...)
    frexp(x)

    Return the mantissa and exponent of x, as pair (m, e).
    m is a float and e is an int, such that x = m * 2.**e.
    If x is 0, m and e are both 0.  Else 0.5 <= abs(m) < 1.0.



If x and y are both floats, then given:

m1, e1 = math.frexp(x)
m2, e2 = math.frexp(y)

Then x*y = m1*m2 * 2.0**(e1 + e2). We can test that:


py> from math import frexp
py> x = 2.5
py> y = 3.5
py> x*y
8.75
py> m1, e1 = math.frexp(x)
py> m2, e2 = math.frexp(y)
py> m1*m2 * 2.0**(e1 + e2)
8.75


Looks good to me. So let's try a less naive version of product():


def product_scaled(values):
    scale = 0
    prod = 1.0
    for a in values:
        m1, e1 = math.frexp(a)
        m2, e2 = math.frexp(prod)
        scale += (e1 + e2)
        prod *= (m1*m2)
    return (prod * 2.0**scale)


py> product_scaled([2.5, 3.5])  # expected 8.75
2.734375



What am I doing wrong?




-- 
Steven
“Cheer up,” they said, “things could be worse.” So I cheered up, and sure
enough, things got worse.