## 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.

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

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 {
::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: [1 ]).

``` 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.

 Category Concept Arts and Crafts of Tcl-Tk Programming Category Mathematics