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... Or see Functional programming for how one can do it with

 interp alias {} factorial {} o product iota1

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:

 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

See also Factorial using event loop

So what would the code look like making use of bignum in pure Tcl? - RS: Like so, for instance:

 package require math ;# bignum is in more recent Tcllib
 proc fac x {
    set prod [math::bignum::fromstr 1]
    while {$x>1} {
      set prod [math::bignum::mul $prod [math::bignum::fromstr $x]]
      incr x -1
    }
    math::bignum::tostr $prod
 }

Testing:

 % fac 100
 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Lars H: The bignum module in infix provides the ! operation (postfix) for factorials and the !! operation for semifactorials. It uses the following recursive procedure (whose first step is to expand n! as n!!*(n-1)!!) to make all multiplications of two roughly equal-sized numbers.

 proc xfactorial {k n} {
    if {[::math::bignum::le $n 1]} then {
       return 1
    } elseif {[::math::bignum::le $n $k]} then {
       return $n
    } else {
       set kk [::math::bignum::add $k $k]
       ::math::bignum::mul [
          xfactorial $kk $n
       ] [
          xfactorial $kk [::math::bignum::sub $n $k]
       ]
    }
 }

Ordinary factorials become

 proc factorial {n} {math::bignum::tostr [xfactorial 1 [math::bignum::fromstr $n]]}

whereas the semifactorial is

 proc semifactorial {n} {
    math::bignum::tostr [xfactorial [math::bignum::fromstr 2] [math::bignum::fromstr $n]]
 }

This seems to work without problems in Tcl 8.5a6 (obviously: [L1 ]).

 for {set x 1; set y 0} {$x<=10000} {incr x} {set y [expr {$x * $y}]}
 set y

Conversion to string for printing the number takes about 20 times as long as the calculation itself. Obviously, I won't be pasting the result here.