I'm skeptical of the need for - and wisdom of - this. Where does it come up? I can't think of any context where this would have been useful, or of any other language or package that does something like this. Long chains of mults are unusual outside of integer combinatorics to begin with.
Closest I can recall came up in the Spambayes project. There we needed to compute the sum of the logs of a potentially large number of probabilities. But log (lack of!) speed proved to be a bottleneck, so I changed it to compute the log of the product of the probabilities. That product was all but certain to underflow to 0.0. But "for real" - no rearrangement of terms would have made any difference to that. Instead, akin to what Mark spelled out, every now & again frexp() was used to push the product bits back into [0.5, 1.0), and the power-of-2 exponent was tracked apart from that.
fsum() is primarily aimed at a very different problem: improving the low-order bits.
How to explain what this change to prod() would do? It may or may not stop spurious overflows or underflows, and the result depends on the order of the multiplicands. But in what way? If we can't/won't define what it does, how can other implementations of Python reproduce CPython's result?
While I don't (yet?) see a real need for it, one thing that could work: compute the product left to right. Full speed - no checks at all. If the result is a 0, a NaN, or an infinity (which is bound to be rare in real life), do it again left to right with the frexp() approach. Then it can be explained: the result is what you'd get from left-to-right multiplication if the hardware had an unbounded exponent field, suffering overflow/underflow at the end only if the true exponent has magnitude too large for the hardware. |