Version 9 of factorial

Richard Suchenwirth 2004-02-08 - Factorial (n!) is a popular function with super-exponential growth. Mathematically put,

   0! = 1
   n! = n (n-1)! if n >0, else undefined

In Tcl, we can have it pretty similarly:

 proc fact n {expr {$n<2? 1: $n * [fact [incr n -1]]}}

But this very soon crosses the limits of integers, giving wrong results.

Weekend reading in an older math book showed me the Stirling approximation to n! for large n (at Tcl's precisions, "large" is > 20 ...), so I built that in:

 proc fact n {expr {
     $n<2? 1: 
     $n>20? pow($n,$n)*exp(-$n)*sqrt(2*acos(-1)*$n):
            wide($n)*[fact [incr n -1]]}
 }

Just in case somebody needs approximated large factorials... But for n>143 we reach the domain limit of floating point numbers. In fact, as Additional math functions points out, the float limit is at n>170, so an intermediate result in the Stirling formula must have busted at 144. For such few values it is most efficient to just look them up in a pre-built table, as Tcllib's math::factorial does.

Incidentally, playing with that, I noticed a bug: as table lookup is guarded with

 if {$x <= [llength $factorialList]} {return [lindex $factorialList $x]}

that function has a blind spot at x=171: it returns "" (list index overrun, should test with < ;-) instead of throwing the error, which it does for x>171...

FPX I have a story to tell on the subject of factorials. In 1993, I was in competition with a friend to come up with the largest factorials. It started when he called me in despair, he was in a bet to show that you could compute 10000!, but his implementation was crapping up. The day after, I mailed him a disk with the full number. The fun part is that he considered C "too slow," yet I kept beating his assembly code by using better algorithms, e.g., by doing all computations "base 10000" (four digits per word), using 32 bit multiplications -- four times faster than using base 100 (two digits per byte). I keep telling this story to warn about the dangers of premature optimization: first get your algorithms right, and only then concentrate on optimizing hot spots.

Eventually, I extended my program to work over a network, and had 70 workstations grinding away for 48 hours on the computation of 1000000!. Never was CPU time more pointlessly wasted.

So all of the above is peanuts. Try this:

 proc fact {n} {
    set e 0; set m 0
    for {set i 2} {$i <= $n} {incr i} {
        set l [expr {log10(double($i))}]
        incr e [expr {int(floor($l))}]
        if {[set m [expr {$m+(fmod($l,1))}]]>=1} {
            incr e; set m [expr {$m-1}]
        }
    }
    return "[expr {pow(10,$m)}]e$e"
 }

 % fact 200
 7.88657867365e374

 % fact 10000
 2.84625968092e35659

Playing Joy of course involves factorial too: see that page for an older, Polish take of mine, where the generic recursion construct primrec can be wrapped as

 interp alias {} factorial {} primrec 1 *

and, soon on the Wiki in RPN again, the reverse Polish form

 rpn /factorial {1 /* primrec} def

It may look strange, but it's valid Tcl...

Factorial puzzle: From the Tcl chatroom on 2004-08-02:

<antirez_> Many many years ago (at least 12), a computer magazine here launched a price. The goal was to calculate the last non-zero digit of 10000! but the average computer was a 286. Of course the solution was to only care about the latest part of the number but many tried to compute the whole number hee

suchenwi Many many years ago (at least 12), a computer magazine here launched a price. The goal was to calculate the last non-zero digit of 10000! but the average computer was a 286. Of course the solution was to only care about the latest part of the number but many tried to compute the whole number hee

 proc lastnonzero max {
        set prod 1
        for {set i 2} {$i<=$max} {incr i} {
           set prod [string index [string trimright [expr {$i*$prod}] 0] end]
        }
        set prod
 }
 % lastnonzero 10000
 8

Took 933 msec on my 200 MHz box.

<antirez_> suchenwi: haha, good answer!

suchenwi Well, I suppose Tcl's "string think" makes it pretty easy to implement your spec... I tested correctness with some small factorials, so I assume by total induction that 8 is the right answer.

<antirez_> Yeah, your implementation is very smart. And 8 is, I computed the real number with the C version of my lib

suchenwi :D