reduce vs. for loop

[Fernando Pérez]

Tim Hochberg wrote:

> Try this out! For n=100, it's only twice as fast, but try it for n=1000, or
> n=2000 and it really shines.... (I think it even gives the right answer,
> which is a plus.)

Well, 

In [22]: def fact4(n):
   ....:     from sys import maxint
   ....:     mult = 1
   ....:     result = 1L
   ....:     subresult = 1
   ....:     max_sub_result = maxint / n
   ....:     for mult in range(n+1):
   ....:         subresult *= mult
   ....:         if subresult > max_sub_result:
   ....:             result *= subresult
   ....:             subresult = 1
   ....:     return result * subresult
   ....:

In [23]: fact4 10
-------> fact4 (10)
Out[23]: 0L

That doesn't look right to me.

It's hard to say what the op was after, but in general when you need large n 
factorials it's rare to truly need exact results. If you can live with an 
approximate answer (good to 10 sig figs or so), then a constant-time solution 
is vastly preferrable. Note that this is far from optimal, it's a 3 minute 
port of the standard routine from NumRecipes. But it will do in a pinch:

from Numeric import *
#-----------------------------------------------------------------------------
def gammln(x):
    """Return the natural log of gamma(x)"""

    cof = array([76.18009172947146,-86.50532032941677,
                 24.01409824083091,-1.231739572450155,
                 0.1208650973866179e-2,-0.5395239384953e-5],Float)

    y = arange(x+1,x+7,typecode = Float)
    tmp = x+5.5
    tmp -= (x+0.5)*log(tmp)
    ser = 1.000000000190015 + sum(cof/y)
    return -tmp+log(2.5066282746310005*ser/x)

def fact5(n):
    return long(round(exp(gammln(n+1))))
#-----------------------------------------------------------------------------

This solution obviously has two limitations:

- not exact, but quite good. Look at a numerical comparison with an exact (but 
super-slow) solution:

def fact6(n):
  return  reduce(operator.mul, range(1,n+1), 1L )

In [27]: fact5 50
-------> fact5 (50)
Out[27]: 30414093203509286326338827883471693198716826006464232959556714496L

In [28]: fact6 50
-------> fact6 (50)
Out[28]: 30414093201713378043612608166064768844377641568960512000000000000L

In [29]: 1.0*(_28-_27)/_27  # force float div
Out[29]: -5.9048555901676581e-11

So that's pretty good numerically.

- can't go outside the float range.

In most 'typical' cases, neither of these is a problem (since you can write 
most large n problems in terms of logs anyway, which is numerically much 
saner).

Cheers,

f.