在Python中处理非常小的数字

8

我需要对大约 1e6 个数进行乘法运算,这些数的数量级大约为 0.01。预期结果的数量级为 1e-100000000。显然,典型的浮点运算无法处理这种情况。

通过在网上进行一些研究,我发现了似乎可以解决这个问题的十进制库(decimal library)。然而,它似乎有一些限制,使它对我的需求无用:

>>> 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')

有人知道这个问题的解决方案吗?

4
你可以在一个单独的整数中维护指数,并将乘积缩放到0.5到1.0之间。有很多可能性,这只是其中一个例子。例如,请参见math.frexp()。取对数可能是最简单的方法,具体取决于你打算如何使用结果。 - President James K. Polk
2个回答

5

你的结果不正确是因为十进制数也有精度(十进制数是固定点数学),所以你在这里也会遇到下溢问题:

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

十进制数('0E-1000026')

但是:

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

十进制('1E-2000000')

您可以通过手动设置精度来解决问题,就像上面的示例一样,或者手动计算指数,例如:

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

可以转换为

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

并且后来到

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

高精度计算模块文档


3

为什么不使用对数?

你想要计算:

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

将其表示为:

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

如果将非常小的正数存储为其自然对数,它们会很好地存储到浮点数中:

ln(0.000001) ≈ -13.81551

不要存储数字本身,而是存储这些值的对数。

假设您将ln(0.0000011)加起来10^6次。 您将得到大约-13815510.558。 作为float,与0.000001^(10^6)相比,更少的精度会丢失。

无论您最终得到什么数字,您都知道您的结果只是该数字的指数。
例如,RESULT = e^-13815510.558

您可以使用以下代码:

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))

控制台打印的消息是:
e^-27.36212057035477

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接