[New-bugs-announce] [issue41458] Avoid overflow/underflow in math.prod()

[Raymond Hettinger]

New submission from Raymond Hettinger <raymond.hettinger at gmail.com>:

For float inputs, the math.prod() function could be made smarter about avoiding overflow() and underflow().  That would also improve commutativity as well.

Other tools like math.hypot() already take measures to avoid overflow/underflow and to improve commutativity.  This makes them nicer to use than näive implementations.

The recipe that I've been using for years is shown below.  It certainly isn't the best way, but it is O(n) and always avoids overflow and underflow when possible.  It has made for a nice primitive when implementing other functions that need to be as robust as possible.  For example, in the quadratic formula, the √(b²-4ac) factors to b√(1-4ac/b²) and the rightmost term gets implemented in Python as product([4.0, a, c, 1.0/b, 1.0/b]).

>>> from math import prod, fabs
>>> from collections import deque
>>> from itertools import permutations

>>> def product(seq):
    s = deque()
    for x in seq:
        s.appendleft(x) if fabs(x) < 1.0 else s.append(x)
    while len(s) > 1:
        x = s.popleft() * s.pop()
        s.appendleft(x) if fabs(x) < 1.0 else s.append(x)
    return s[0] if s else 1.0

>>> data = [2e300, 2e200, 0.5e-150, 0.5e-175]
>>> for factors in permutations(data):
	print(product(factors), prod(factors), sep='\t')

	
1e+175	inf
1.0000000000000001e+175	inf
1e+175	inf
1e+175	1e+175
1.0000000000000001e+175	inf
1.0000000000000001e+175	1.0000000000000001e+175
1.0000000000000001e+175	inf
1e+175	inf
1.0000000000000001e+175	inf
1.0000000000000001e+175	1.0000000000000001e+175
1e+175	inf
1e+175	1e+175
1e+175	inf
1e+175	1e+175
1.0000000000000001e+175	inf
1.0000000000000001e+175	1.0000000000000001e+175
1e+175	0.0
1.0000000000000001e+175	0.0
1.0000000000000001e+175	inf
1.0000000000000001e+175	1.0000000000000001e+175
1e+175	inf
1e+175	1e+175
1.0000000000000001e+175	0.0
1e+175	0.0

For math.prod(), I think a better implementation would be to run normally until an underflow or overflow is detected, then back up a step and switch-over to pairing low and high magnitude values.  Since this is fallback code, it would only affect speed to the extent that we test for overflow or underflow at every step.  Given how branch prediction works, the extra tests might even be free or at least very cheap.

The math module functions usually go the extra mile to be more robust (and often faster) than a näive implementation.  These are the primary reasons we teach people to prefer sqrt() over x**2, log1p(x) over log(1+x), prod(seq) over reduce(mul, seq, 1.0), log2(x) over log(x, 2), fsum(seq) over sum(seq), and hypot(x,y) over sqrt(x**2 + y**2).  In each case, the library function is some combination of more robust, more accurate, more commutative, and/or faster than a user can easily create for themselves.

----------
components: Library (Lib)
messages: 374693
nosy: mark.dickinson, pablogsal, rhettinger, tim.peters
priority: normal
severity: normal
status: open
title: Avoid overflow/underflow in math.prod()
type: behavior
versions: Python 3.10

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