Handling very very small numbers in Python

Question

I need to multiply about 1e6 numbers that are of the order of 0.01. The expected result is of order 1e-100000000. Obviously the typical floating-point arithmetic cannot handle this.

Doing some research on the web I found the decimal library which seems to fix this problem. However it appears to have limitations that make it useless for my needs:

>>> Decimal('.01')**Decimal('1e5') # Seems to handle this
Decimal('1E-200000')
>>> Decimal('.01')**Decimal('1e5')*Decimal('1E200000') # Yeah! It works!
Decimal('1')
>>> Decimal('.01')**Decimal('1e6') # This result is strange...
Decimal('0E-1000026')
>>> Decimal('.01')**Decimal('1e6')*Decimal('0E1000026') # Wrong result
Decimal('0')

Does anyone know any solution to this?

You could always maintain the exponent yourself in a separate int and keep the product scaled to be between .5 and 1.0. There are a lot of possibilities, this is just one example. See `math.frexp()` for example. Maybe taking logs is the easiest thing, it depends somewhat on what you intend to do with the result. — President James K. Polk, Nov 09 '19 at 22:57

ingvar · Accepted Answer · 2019-11-09T23:06:02.810

Your result is incorrect because decimal has precision too (decimal is fixed point math), so you get underflow issue here too:

Decimal('.01')**Decimal('1e6')

Decimal('0E-1000026')

But:

getcontext().prec = 1000000000   # sets precision to 1000000000
Decimal('.01')**Decimal('1e6')

Decimal('1E-2000000')

You can fix your issue by manualy setting precision as in example above or manually calculate powers, for example:

Decimal('.01')**Decimal('1e6')

can be converted to

Decimal('1e-2') ** Decimal('1e6')

and later to

1 ** ((-2) ** 1e6) = 1 ** (-2000000)

Decimal module documentation

score 3 · Answer 2 · answered Nov 09 '19 at 23:41

Why not use logarithms?

You want to compute:

RESULT  = x1 * x2 * x3 * x4 ... * xn

Represent that as:

ln(RESULT) = ln(x1) + ln(x2) + ln(x3) + ln(x4) ... + ln(xn)

Very small positive numbers store nicely into floats if you store their natural logarithm:

ln(0.000001) ≈ -13.81551

Instead of storing the numbers themselves, store the log of the values.

Suppose you add ln(0.0000011) to itself 10^6 times. You get approximately -13815510.558. Less precision is lost on that as a float than 0.000001^(10^6)

Whatever number you get in the end, you know that your result is just the number e raised to that power. For example, RESULT = e^-13815510.558

You can use the code below:

import math

class TinyNum:
    def __init__(self, other=None, *, pow=None):
        """
        x = TinyNum(0.0000912922)
        x = TinyNum("0.12345")     # strings are okay too
        x = TinyNum(pow = -110)    # e^-110
        y = TinyNum(x)             # copy constructor
        """
        if other:
            if isinstance(other, type(self)):
                self._power = other._power
            else:
                self._power = math.log(float(str(other)))
        else: # other == None
            self._power = float(str(pow))

    def __str__(self):
        return "e^"+str(self._power)

    def __mul__(lhs, rhs):
        rhs = type(lhs)(rhs)
        return type(lhs)(pow=lhs._power + rhs._power)

    def __rmul__(rhs, lhs):
        lhs = type(rhs)(lhs)
        return type(rhs)(pow=lhs._power + rhs._power)

    def __imul__(total, margin):
        total._power = total._power + type(total)(margin)._power


lyst = [
    0.00841369,
    0.004766949,
    0.003188046,
    0.002140916,
    0.004780032
]

sneaky_lyst = map(TinyNum, lyst)

print(math.prod(sneaky_lyst))

The message printed to the console is:

e^-27.36212057035477

Handling very very small numbers in Python

2 Answers2