Is there an equivalent of Perl's big float data type in python? If I do factorial(500) in python

Is there an equivalent of Perl's BigFloat data type in Python? The reason I ask is that I want to compute nCk using its definition, n!/k!*(n-k)!.

With Perl's BigFloat data type, computing from the definition works correctly for any n and k. For example

factorial(500)/factorial(10)*factorial(490)

produces the exact answer when n and k are BigFloats.

In Python, both factorial(500) and factorial(10)*factorial(49) give the exact answers using whatever Python uses for ints. Thus it appears that python can do very high precision arithmetic. However the quotient

int(factorial(500)/(factorial(10)*factorial(490))

comes close to the exact answer, but falls a little short?

Is there a way to get an exact answers out of python for expressions like the above?

Solution

Python's int objects can grow as big as necessary (subject only to how much memory is available), so they can be used for calculations involving huge numbers. In Python 2, there were 2 integer types, int and long, with int being used for values that fit into machine integers, but in Python 3 they've been amalgamated into a single int type.

Python doesn't have a built-in BigFloat type, but the standard library does have the decimal module which can do basic arithmetic, including square roots, to any desired precision. If you need to do arbitrary precision mathematics with more functions, please see the excellent 3rd-party library, mpmath.

You can safely use // floor division when computing binomial coefficients, since the terms in the denominator are guaranteed to divide the numerator. Eg

from math import factorial

a = (factorial(500) // factorial(490)) // factorial(10)
print(a)

output

245810588801891098700

However, it may be more efficient to calculate the binomial coefficient with a simple loop, rather than calculating those huge factorials.

def binomial(n, r):
    ''' Binomial coefficients '''
    if not 0 <= r <= n:
        return 0
    p = 1
    r = min(r, n - r)
    for i in range(1, r+1):
        p *= n
        p //= i
        n -= 1
    return p

# Test

print(binomial(500, 10), '\n')

for i in range(10):
    print([binomial(i, j) for j in range(i+1)])

output

245810588801891098700 

[1]
[1, 1]
[1, 2, 1]
[1, 3, 3, 1]
[1, 4, 6, 4, 1]
[1, 5, 10, 10, 5, 1]
[1, 6, 15, 20, 15, 6, 1]
[1, 7, 21, 35, 35, 21, 7, 1]
[1, 8, 28, 56, 70, 56, 28, 8, 1]
[1, 9, 36, 84, 126, 126, 84, 36, 9, 1]