reduce vs. for loop

[Tim Hochberg]

Oops, that should have been:

def fact4(n, step=1):
    from sys import maxint
    mult = 1
    result = 1L
    subresult = 1
    max_sub_result = maxint / n
    for mult in range(1,n+1,step):
        subresult *= mult
        if subresult > max_sub_result:
            result *= subresult
            subresult = 1
    return result * subresult


def fact5(n):
    if n < 1024:
        return fact4(n)
    result = fact4(n, 2) * 2**(n/2) * fact5(n/2)
    return result


Where fact5 is faster for large enough n. These aren't quite as fast as I
claimed before -- I should have been suspicous when fact4 sped up so much
when I changed from a while to a for loop. I trusted that things were still
working 'cause I had some code in their to check the results, but it
depended on assert and I ran those final tests with '-O'. Doh! I beg the
court for mercy -- it was late.

As Fernado notes, these sort of computations are often better done using
approximate methods, but that's not nearly as entertaining.

-tim

"Fernando Pérez" wrote:
> 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.)
> > 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.