[Tutor] factorial

[Greg Chapman]

I'm leaving this pretty much intact since it wasn't originally posted to the
list.

>You may want to try the following. fact(n) and fact2(n) will compute
>1000000!. Fact3(n) uses
>De Moivre's approximation also called Stirling's approximation and is good
>only for about 15 places.
>
>The output for 15! is shown.
>>>>fact(15)
>1307674368000L
>>>>fact2(15)
>1307674368000
>n! for  n    =   15
>Actual  digits   13
>Formula digits   13
>Seconds to calc  0.0
>>>>fact3(15)
>1.307674368e+012
>>>>
>
>
>Python has great power. The integer precision is fantastic. Any way to get
>float to be arbitrary?
>I wish the ^ operator could replace ** for exponentiation AND to  be able
to
>edit script and go right
>back in and run without having to close the dos >>> box and re-enter and
>then from factorial import
>* . should be able to  from factorial import * without closing the >>>
>session.

you should be able to do: >>>reload(factorial)

>Have fun
>
>
>#                           Program factorial.py
>#                  Cino Hilliard hillcino368@hotmail.com
># Compute n! using python and the ConText editor which supports Python
># Fact2 tests the accuracy of my formula for the number of digits in n!
># Fact3 tests the accuracy of my modification of De moivre's Formula
>#                                                   1                3
>#          N                  ( -N + 1/(12*N) - 1/(360*N) +. . .)
># N! ~ N * sqrt(2*PI*N) * e
>#
># Perm calculates the permutations of n things taken r at a time
># Comb calculates the combinations of n things taken r at a time
># digits = int((n+.5)*log(n)/log(10) - n*log(exp(1.0))/log(10)+.3991)+1
># usage: >>> from factorial import *  >>> fact2(number)
>#
>from math import *
>def fact(n):
>    f = j = 1
>    while j <= n:
>            f*=j
>            j+=1
>    return f
># digits = int((n+.5)*log(n)/log(10) - n*log(exp(1.0))/log(10)+.3991)+1
># Test this fantastic formula for digits in n!
>def fact2(n):
>    import time
>    t1=time.time()
>    x = fact(n)
>    t2=time.time()
>    print x
>    d = int((n+.5)*log(n)/log(10) - n*log(exp(1.0))/log(10)+.3991)+1
>    print "n! for  n    =  ",n
>    print "Actual  digits  ",len(str(x))
>    print "Formula digits  ",d
>    print "Seconds to calc ",t2 - t1
>
># Test a modification of De moivre's formula for n!
>def fact3(n):
>    series=n
>    nsqr=n*n
>    n2pi =n*2*3.1415926535897932
>    sqrt2pi = sqrt(n2pi)
>    ptn = .5/n
>    tmp = 6
>    ptn /= tmp
>    series -= ptn
>    tmp=30
>    ptn /= nsqr
>    ptn /= tmp
>    series += ptn
>    tmp= 3.49999997364305
>    ptn /= nsqr
>    ptn /= tmp
>    series -= ptn
>    tmp = 1.3333346851
>    ptn /= nsqr
>    ptn /= tmp
>    series += ptn
>    tmp = exp(series)
>    stirtmp = (n**n)/tmp
>    stirling = sqrt2pi*stirtmp
>    print stirling
>
>def comb(n,r):
>        return fact(n)/fact(n-r)/fact(r)
>
>def perm(n,r):
>        return fact(n)/fact(n-r)
>
>
>2^3 + 13^3 + 33^3 + 43^3 = 49^3

I assumed that there would be some advantage to using the approximation, so
I modified the fact2() function to run and time both the fact() function and
the fact3() function:

def stopwatch(func, p):           #I just thought this would be a neat
utility function
    from time import clock        #clock() gives better resolution than
time()
    t1 = clock()
    x  = func(p)
    t2 = clock()
    return (x, t2 - t1)

def fact2(n):
    a = stopwatch(fact, n)
    b = stopwatch(fact3, n)
    print '%e' % a[0]            #I changed this so I could play with bigger
numbers
    print '%e' % b[0]
    d = int((n+.5)*log(n)/log(10) - n*log(exp(1.0))/log(10)+.3991)+1
    print "n! for  n    =         ",n
    print "Actual  digits         ",len(str(a[0]))
    print "Formula digits         ",d
    print "Using fact() time was  ",a[1],"sec"
    print "Using fact3(), time was",b[1],"sec"

>>> fact2(100)
9.332622e+157
9.332622e+157
n! for  n    =          100
Actual  digits          158
Formula digits          158
Using fact() time was   0.000508724584506 sec
Using fact3(), time was 0.00019527648783 sec

but if I switch from scientific notation and force b[0] to long(b[0]), then
we can see the lack of precision:

>>> fact2(30)
265252859812191058636308480000000
265252859812190167497676244910080
n! for  n    =          30
Actual  digits          33
Formula digits          33
Using fact() time was   0.000237181313423 sec
Using fact3(), time was 0.00014834308331 sec

so there is definitely a performance benefit, if you're doing a lot of large
factorial calculations and only need the 14 significant digits. so I
wondered if there would be any application for that.  I added a function
comb3() which is the same as comb() except that it calls fact3() instead of
fact().

>>> x = factorial.comb(100, 5)
>>> print x
75287520
>>> y = factorial.comb3(100, 5)
>>> print int(y)
75287520

and it looks like for some problems, approximation does come close enough.

Interesting stuff.  Thanks Cino.

greg