Subject: Re: factorial of largish numbers
From: rpw3@rpw3.org (Rob Warnock)
Date: Wed, 05 Apr 2006 18:23:38 -0500
Newsgroups: sci.math.symbolic,comp.lang.lisp
Message-ID: <ULWdnVUQaeznzKnZ4p2dnA@speakeasy.net>
Richard J. Fateman <fateman@eecs.berkeley.edu> wrote:
+---------------
| Peter Luschny wrote:
|  > Perhaps someone in this group could give me a
| > helping hand and show me the Lisp implementation of the algorithm
| 
| here it is in lisp
| 
| (defun sprfact(n)
|    (cond ((< n 0)(error "bad arg ~s to factorial" n))
| 	((< n 2) 1)
| 	(t
| 	 (let ((p 1) (r 1) (N 1) (log2n (floor (log n 2)))
| 	       (h 0) (shift 0) (high 1) (len 0))
| 	   (declare (special N))
| 	   (while (/= h n)
| 	     (incf shift h)
| 	     (setf h (ash n  (- log2n)))
| 	     (decf log2n)
| 	     (setf len high)
| 	     (setf high (if (oddp h) h (1- h)))
| 	     (setf len (ash (- high len) -1))
| 	     (cond ((> len 0)
| 		    (setf p (* p (prod len)))
| 		    (setf r (* r p)))))
| 	   (ash r shift)))))
| 
| (defun prod(n)
|    (declare (special N))
|    (let ((m (ash n -1)))
|      (cond ((= m 0) (incf N 2))  ((= n 2) (* (incf N 2)(incf N 2)))
| 	  (t (* (prod (- n m)) (prod m))))))
+---------------

This will not run in Common Lisp as as written, because of the
WHILE, and even when that is corrected will give incorrect answers
[e.g., (SPRFACT 5) => 3] because the default READTABLE-CASE in
Common Lisp is :UPCASE. If you put this line at the top of your
program:

    (setf (readtable-case *readtable*) :invert)

and then define WHILE as follows [yeah, yeah, I know, using DO
would be safer]:

    (defmacro while (bool &body body)
      `(loop while ,bool do ,@body))

and only *THEN* type [or paste] the above definitions for SPRFACT/PROD,
then it will work and give correct answers.

But on CMUCL it's more than twice as slow [0.21f0 s versus 0.1f0 s]
than the previous K/FACT version mentioned in this thread. [And, yes,
I did compile it.]  I suspect that if you rewrite it without using
special variables it might improve [e.g., perhaps by making PROD be
an FLET inside the LET binding within SPRFACT].


-Rob

-----
Rob Warnock			<rpw3@rpw3.org>
627 26th Avenue			<URL:http://rpw3.org/>
San Mateo, CA 94403		(650)572-2607