如何在python中将超大数与超小数相乘？

Question

我正在做一些概率计算。
在我的一项任务中，我需要将从 10000 个项目中选择 8000 个样本的组合数乘以 0.8**8000。
组合号是很长的long-number，在numpy的帮助下，我得到0.8**8000的结果为5.2468172239242176864e-776。
但是当我尝试将这两个数字相乘时，我得到了[9] 34845 segmentation fault ipython -i。
那我该怎么做这样的乘法呢？

PS：这是我的一段代码

import numpy
d2 = numpy.float128(0.8) ** 8000
d1 = 165555575235503558460892983752748337696863078099010763950122624527927836980322780662408249953188062227721112100054260160204180655980717428736444016909193193353770953722788106404786520413339850951599929567643032803416164290936680088121145665954509987077953596641237451927908536624592636591471456488142060812180933761408708169972797751139799352908109763166895772281109195968567911923343187466596002627570139321755043803267091330804414889831229832744256038117150720178689066894068507531026417815624234453195871008113238128934831837842040515600131726096039123279876153916504647241693083829553081901075278042326502699324012014817969085443550523855284341221708045253558716789811929298590803855947461554713178815399150688529048306222786951038548880400191620565711291586700534540755526276938422405001345270278335726581375322976014611332999126216550500951669985289322635729053541565465940744524663726205818866513444952048185208697438054246674199211750006230637806394882672053335493831407089830994135058867370833787098758113596190447219426121568324685764151601296948654893782399960327514764114467176417125060133454019708700782282480571935020898204763471121684913190735908414301826140125010936910161942130277906874552721346626800201093026689035996876035329180150478191582393837824731994055511844267891121846403164857127885959745644323971338513739214928092232132691519007718752719466750891748327404893783451436251805894736392433617289459646429204124129760273396235033220480921175386059331059354409267348067375581516003852060360378571075522650956157791058846993826792047806030332676423336065499519953076910418838626376480202828151673161942289092221049283902410699951912366163469099917310239336454637062482599733606299329923589714875696509548029668358723465427602758225427644633549944802010973352599970041918971524450218727345622721744933664742499521140235707102217164259438766026322532351208348119475549696983427008567651685921355966036780080415723688044325099562693124488758728102729947753752228785786200998322978801432511608341549234067324280214361346940194251357867820535466891356019219904248859277399657389914429390105240751239760865282709465029549690591863591028864648910033430400L
print d1 * d2

Answer 1

当用一个非常大的数字乘以一个非常小的数字时，使用浮点数会带来很大的不准确性。在您的情况下，数字的大小会导致溢出错误，因此您遇到的问题不仅仅是不准确！

每当您发现自己处于这种情况时，首先检查是否可以留在整数域中，然后先“按摩”一下数字会很有用。在您的情况下，这是可能的，我将在下面解释。

乘法的一个操作数，非常大的数，是来自 10000 个项目的 8000 个样本。对组合数量使用封闭式方程，其中您的样本大小 n 是 10000，子集大小 r 是 8000。此处的感叹号 (!) 是阶乘，您可以在 python 中的 math.factorial 中找到它。

C(n,r) = n! / r! (n - r)!

另一个操作数0.8 ** 8000 是极小的数字，根据索引法则等于：

8**8000 / 10**8000

所以当我们将这两个数字相乘时，我们想要的答案是：

     10000! * 8**8000
--------------------------
 8000! * 2000! * 10**8000

我们称这个号码为x，然后取两边的对数。在对数域中工作会将乘法转换为加法，将除法转换为减法，从而使事情更易于管理。

from math import log, factorial
numerator = log(factorial(10000)) + 8000*log(8)
denominator = log(factorial(8000)) + log(factorial(2000)) + 8000*log(10)
log_x = numerator - denominator

现在这些数字的数量级可以在 python 中使用。

您会发现log_x 大约等于3214。您现在只需观察exp(log_x) == x 即可找到答案。这是一个非常大但有限的数字。

Answer 2

任意精度整数并不是解决这个问题的真正方法，因为调用 log 会破坏任何精度，所以我会让 scipy.special.gammaln 自己说话（但请参阅下面的编辑):

from math import log, factorial
from scipy.special import gammaln

def comp_integral(n, r, p, q):
    numerator = log(factorial(n)) + r*log(8)
    denominator = log(factorial(r)) + log(factorial(n-r)) + r*log(q)
    return numerator - denominator

def comp_gamma(n, r, p, q):
    comb = gammaln(n+1) - gammaln(n-r+1) - gammaln(r+1)
    expon = r*(log(p) - log(q))
    return comb+expon

In [220]: comp_integral(10000, 8000, 8, 10)
Out[220]: 3214.267963130871

In [221]: comp_gamma(10000, 8000, 8, 10)
Out[221]: 3214.2679631308811

In [222]: %timeit comp_integral(10000, 8000, 8, 10)
10 loops, best of 3: 80.3 ms per loop

In [223]: %timeit comp_gamma(10000, 8000, 8, 10)
100000 loops, best of 3: 11.4 µs per loop

请注意，最多 14 位的输出相同，但 gammaln 版本快了近 8000 倍。如果您要经常这样做，这将很重要。

编辑：gammaln 所做的是计算 gamma 函数的自然对数。伽马函数可以被认为是阶乘的泛化，factorial(n) == gamma(n+1)。所以comb(n,r) == gamma(n+1)/(gamma(n-r+1)*gamma(r+1))。然后把日志变成上面的形式。

Gamma 还具有 fractional inputs and for negative numbers 的值。不过，这并不重要。

Answer 3

我维护gmpy2 库，它可以很容易地做到这一点。

>>> import gmpy2
>>> gmpy2.comb(10000,8000) * gmpy2.mpfr('0.8')**8000
mpfr('8.6863984366232171e+1395')

Answer 4

在 wim 的出色答案的基础上，您还可以将此数字存储为 Fraction，方法是构建一个素因数列表、进行任何取消并将所有内容相乘。

我已经为这个问题包含了一个相当幼稚的实现。它会在不到一分钟的时间内按原样返回一个分数，但如果您实施更智能的因式分解，您肯定可以使其更快。

from collections import Counter
from fractions import Fraction
import gmpy2 as gmpy

def get_factors(n):
    factors = Counter()
    factor = 1
    while n != 1:
        factor = int(gmpy.next_prime(factor))
        while not n % factor:
            n //= factor
            factors[factor] += 1
    return factors

factors = Counter()

# multiply by 10000!
for i in range(10000):
  factors += get_factors(i+1)

# multiply by 8^8000
factors[2] += 3*8000

#divide by 2000!
for i in range(2000):
  factors -= get_factors(i+1)

#divide by 8000!
for i in range(8000):
  factors -= get_factors(i+1)

# divide by 10^8000
factors[2] -= 8000
factors[5] -= 8000

# build Fraction
numer = 1
denom = 1
for f,c in factors.items():
    if c>0:
        numer *= f**c
    elif c<0:
        denom *= f**-c

frac = Fraction(numer, denom)

看起来大约是 8.686*10^1395