Message375052
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. |
|
Date |
User |
Action |
Args |
2020-08-08 17:30:17 | tim.peters | set | recipients:
+ tim.peters, rhettinger, mark.dickinson, veky, pablogsal, Jeffrey.Kintscher |
2020-08-08 17:30:17 | tim.peters | set | messageid: <1596907817.78.0.829600653701.issue41458@roundup.psfhosted.org> |
2020-08-08 17:30:17 | tim.peters | link | issue41458 messages |
2020-08-08 17:30:17 | tim.peters | create | |
|