[issue41458] Avoid overflow/underflow in math.prod()

[Tim Peters]

Tim Peters <tim at python.org> added the comment:

Well, that can't work:  the most likely result for a long input is 0.0 (try it!). frexp() forces the mantissa into range [0.5, 1.0).  Multiply N of those, and the result _can_ be as small as 2**-N. So, as in Mark's code, every thousand times (2**-1000 is nearing the subnormal boundary) frexp() is used again to force the product-so-far back into range. That's straightforward when going "left to right".

With fancier reduction schemes, "it varies". Aiming for "obviously correct" rather than for maximal cleverness ;-) , here I'll associate each partial product with an integer e such that it's guaranteed (in the absence of infinities, NaNs, zeroes), abs(partial_product) >= 2^^-e. Then quick integer arithmetic can be used in advance to guard against a partial product underflowing:

    def frpair(seq):
        from math import frexp, ldexp
        stack = []
        exp = 0
        for i, x in enumerate(seq, start=1):
            m, e = frexp(x)
            exp += e
            stack += [(m, 1)]
            while not i&1:
                i >>= 1
                (x, e1), (y, e2) = stack[-2:]
                esum = e1 + e2
                if esum >= 1000:
                    x, e = frexp(x)
                    exp += e
                    y, e = frexp(y)
                    exp += e
                    esum = 2
                stack[-2:] = [(x * y, esum)]
        total = 1.0
        totale = 0
        while stack:
            x, e = stack.pop()
            totale += e
            if totale >= 1000:
               total, e = frexp(total)
               exp += e
               x, e = frexp(x)
               exp += e
               totale = 2
            total *= x
        return ldexp(total, exp)

But I see no obvious improvement in accuracy over "left to right" for the uniformly random test cases I've tried.

----------

_______________________________________
Python tracker <report at bugs.python.org>
<https://bugs.python.org/issue41458>
_______________________________________