Cool!  So looks like you could also address an accuracy (not out-of-range) thing the frexp() method also does as well as possible: loosen the definition of "underflow" to include losing bits to subnormal products. For example, with the inputs

>>> xs = [1e-200, 1e-121, 1e308]
>>> product(xs)
9.98012604599318e-14
>>> _.hex()
'0x1.c17704f58189cp-44'
>>> math.prod(xs) # same thing
9.98012604599318e-14
>>> prod2(xs) # a tiny implementation of the frexp() method
1e-13
>>> _.hex()
'0x1.c25c268497682p-44'

product/math.prod get only a handful of the leading result bits right, because

>>> 1e-200 * 1e-121
1e-321
>>> _.hex()
'0x0.00000000000cap-1022'

loses all but a handful of the trailing bits due to being in the subnormal range.  Changing the order can repair that:

>>> 1e-200 * 1e308 * 1e-121
1e-13

IOW, so far as "underflow" goes, we _start_ losing bits when the product becomes subnormal, well before it reaches 0.