How to multiply super large number with super small number in python?

Question

I am doing some probability calculation.
In one of my tasks, I need to multiply the number of combinations by choosing 8,000 samples from 10,000 items with 0.8 ** 8000.
The combination number is long long-number, and with numpy I get the result 0.8**8000as 5.2468172239242176864e-776.
But when I try to multiply these two numbers, I got [9] 34845 segmentation fault ipython -i.
How can I do this multiplication?

PS: This is part of my code.

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

Answer 1

, . , , !

, , , , "" . , , .

, , 8000 10000 . , n 10000, r 8000. Exclam (!) , math.factorial python.

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

, - , , 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 . , .

EDIT: , gammaln , - -. - , factorial(n) == gamma(n+1). , comb(n,r) == gamma(n+1)/(gamma(n-r+1)*gamma(r+1)). .

Gamma . .

Answer 3

gmpy2, .

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

Answer 4

By creating a large wim answer, you can also save this number as Fractionby building a list of simple factors, making any cancellations and multiplying all together.

I have included a rather naive implementation for this problem. It returns a fraction in less than a minute as soon as, but if you implement a bit smarter factorization, you can do it even faster.

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)

Looks like 8.686 * 10 ^ 1395