Reputation: 593

Factorial of a large number in python

Here's my approach to factorials:

def factorial(n):
    '''Returns factorial of n'''
    r = 1
    for i in range(1, n + 1):
        r *= i
    return r

I think it's pretty straightforward, though I guess you could make something more efficient, because it takes ages for large numbers like 100000. My question is, is there? math.factorial() is no good either, it takes roughly the same amount of time.

Upvotes: 11

Answers (9)

Daniel Fischer

Reputation: 183868

Multiplying the numbers in sequence,

r = 1
for i in range(1, n + 1):
    r *= i
return r

creates a large number (as in tens of thousands of bits) very quickly, and then you have a lot of multiplications of one huge number and one small number. Multiplications where at least one of the factors is huge are slow.

You can speed it up considerably by reducing the number of multiplications involving huge numbers, for example

def range_prod(lo,hi):
    if lo+1 < hi:
        mid = (hi+lo)//2
        return range_prod(lo,mid) * range_prod(mid+1,hi)
    if lo == hi:
        return lo
    return lo*hi

def treefactorial(n):
    if n < 2:
        return 1
    return range_prod(1,n)

produces, timing the computation of 100000! % 100019 (I first tried len(str(fun(100000)), but the conversion to string is abominably slow, so that made the difference seem smaller than it is):

$ python factorial.py 
81430
math.factorial took 4.06193709373 seconds
81430
factorial took 3.84716391563 seconds
81430
treefactorial took 0.344486951828 seconds

so a more than 10× speedup for 100000!.

Upvotes: 27

Robotbugs

Reputation: 4387

It's actually unusual to really need the true factorial value n! in many areas of application. Often it's way more realistic to use the natural log of the factorial. I can't think of any applications where the log can't be used as a better alternative, because factorials are most often used to compute values related to probabilities of choosing combinations of things.

A common thing to compute is probabilities that are based on factorials such as choosing the binomial coefficient (n k) = n! / (k!(n-k)!). Given this is a ratio of factorials then log(n k) = log(n!)-log(k!)-log((n-k)!) which is reliably computed using one of the various log factorial approximations. And if you do a lot of probability math it's generally best to do it in the log domain anyway (measuring probability in decibels) because it often involves extremely wide ranges of numbers less than 1, and so math precision will fall apart very quickly using common floating point representations if the log version is not used.

E.T.Jaynes was a famous mathematician and an expert in probability theory and I'd recommend his book "Probability Theory: The Logic of Science" as a very readable source on this topic and Bayesian reasoning and information theory using log probabilities.

Upvotes: 2

Paddy3118

Reputation: 4772

You can use the reduce function rather than explicit looping thus:

>>> from functools import reduce
>>> mul = int.__mul__
>>> len(str(reduce(mul, range(2,100001), 1)))
456574
>>>

In Python 2 you need to use longs: long.__mul__, and len(str(reduce(mul, range(2L,100001L), 1L)))

Upvotes: 1

dstromberg

Reputation: 7167

If you just need an approximation, Ramanujan's factorial approximation is supposed to be more accurate than Stirling's.

If you need (or want) something precise, you might try GMP, the GNU Multiple Precision library. I've used it successfully for primality testing of large numbers in Python.

Upvotes: 3

Peter de Rivaz

Reputation: 33509

Factorials get very large, so it is often better to deal with logarithms of the number.

Many languages have an lgamma library function which computes the natural logarithm of the factorial of n-1.

This means that you can compute the natural logarithm of factorial(n) via lgamma(n+1).

You can divide by log10 to turn this into a base 10 logarithm.

So if you just want the number of digits, then this Python code will give the answer immediately:

from math import *
print ceil(lgamma(100000+1)/log(10))

Upvotes: 11

Martin Drautzburg

Reputation: 5243

The slowdown in caused by a quadradic effect: as n gets larger you have to do more multiplications, but you also have to multiply larger numbers.

Finding a better algorithm won't be easy. You can try to exploit symmetries (as in FFT). It could also well pay to do multiplications in a different order, with intermediate results, such that you end up with multiplying only a few very big numbers at the end, but I haven't thought that to the end. In any case, you will have to find a law to exploit.

Look here for further inspiration.

Upvotes: 0