|
|
| Author |
Message |
Richard J. Fateman
science forum addict
Joined: 04 May 2005
Posts: 81
|
 Posted: Wed Apr 05, 2006 5:23 pm Post subject:
factorial of largish numbers
|
|
|
I invite comments on a paper posted at
http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf
on computing factorial of, say, 20,000 or more.
Programs in lisp and some comparison with Maple and
Mathematica. are included.
I've used (Allegro 7.0) Lisp + GMP for best times, but maybe
you have other ideas.
Thanks
RJF |
|
| Back to top |
|
 |
Alexander Schmolck
science forum beginner
Joined: 05 Apr 2006
Posts: 4
|
| Posted: Wed Apr 05, 2006 6:34 pm Post subject:
Re: factorial of largish numbers
|
|
|
"Richard J. Fateman" <fateman@eecs.berkeley.edu> writes:
The reason for using lisp is somewhat bogus. For example this python program
def k(n, m=1):
if n <= m: return n
else: return k(n,2*m) * k(n-m,2*m)
will compute k(20000) just fine in any non-ancient version of python (it is in
fact several orders of magnitude faster than compiled cmucl which I eventually
^C'ed and seems to be about the same speed as clisp).
'as |
|
| Back to top |
|
|
Richard J. Fateman
science forum addict
Joined: 04 May 2005
Posts: 81
|
| Posted: Wed Apr 05, 2006 6:58 pm Post subject:
Re: factorial of largish numbers
|
|
|
OK, Python now does bignums too. thanks.
CMUCL may benefit if you set the optimization to (speed 3) and
also (declare (fixnum n m)). It is unreasonable to compute factorials
of numbers that are not fixnums. Python and CLISP presumably both
use GMP arithmetic.
RJF
Alexander Schmolck wrote:
| Quote: | "Richard J. Fateman" <fateman@eecs.berkeley.edu> writes:
I invite comments on a paper posted at
http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf
The reason for using lisp is somewhat bogus. For example this python program
def k(n, m=1):
if n <= m: return n
else: return k(n,2*m) * k(n-m,2*m)
will compute k(20000) just fine in any non-ancient version of python (it is in
fact several orders of magnitude faster than compiled cmucl which I eventually
^C'ed and seems to be about the same speed as clisp).
'as |
|
|
| Back to top |
|
|
Alexander Schmolck
science forum beginner
Joined: 05 Apr 2006
Posts: 4
|
| Posted: Wed Apr 05, 2006 8:11 pm Post subject:
Re: factorial of largish numbers
|
|
|
"Richard J. Fateman" <fateman@eecs.berkeley.edu> writes:
| Quote: | OK, Python now does bignums too. thanks.
|
As do ruby, smalltalk, J etc. On the other hand some popular scheme
implementations (e.g. bigloo, chicken and stalin) either don't have bignums at
all or require you to jump through hoops to get them -- so I think it is more
a question of high-levelness vs. focus on speed than algol vs lisp-like
syntax.
| Quote: | Python and CLISP presumably both use GMP arithmetic.
|
Despite the developer overlap for clisp/GMP neither does IIRC; I haven't
checked again though.
'as |
|
| Back to top |
|
|
Jens Axel Søgaard
science forum beginner
Joined: 08 May 2005
Posts: 28
|
| Posted: Wed Apr 05, 2006 8:51 pm Post subject:
Re: factorial of largish numbers
|
|
|
Alexander Schmolck wrote:
Why? Luschny is into factorial progams and uses Java:
<http://www.luschny.de/math/factorial/index.html>. It makes
sense to choose a language with bignums to showcase factorial
algorithms. The ones mentioned in the paper haven't got builtin
bignums and Lisp does.
--
Jens Axel Søgaard |
|
| Back to top |
|
|
Peter Luschny
science forum beginner
Joined: 13 May 2005
Posts: 36
|
| Posted: Wed Apr 05, 2006 9:10 pm Post subject:
Re: factorial of largish numbers
|
|
|
Richard J. Fateman wrote:
It's nice to see you contemplating this seemingly
trivial problem. It is not trivial at all!
And it can be used to study implementations of arithmetic,
computational models and computational complexity in detail.
For example Arnold Schoenhage studied in his book
'Fast Algorithms: a Multitape Turing Machine Implementation'
(with A.Grotefeld and E.Vetter) his routines for integer
arithmetic comparing different implementations of n!.
See chapter 6.5 'How to compute n!'.
The running time of his algorithm is asymptotically bounded
by the same order of growth as the time for integer
multiplication for binary numbers of length log(n!),
improving a result of P.B. Borwein 'On the complexity
of calculating factorials' by a factor of order log log n.
On the other hand Schoenhage et alia wrote:
"As it stands, it is a bit disappointing to see [this
algorithm] merely better by a factor of two - in
contrast to the theoretical improvement by a factor
'of order' log n .."
Having implemented his algorithm and many others
in a completely different computational environment
(namely Java), I made the same observation.
For those who care to look:
http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
Regrettably I have no implementations in Lisp. I do would like to
include some. Perhaps someone in this group could give me a
helping hand and show me the Lisp implementation of the algorithm
which I like best? I put it in an appendix, written in a C-like language.
Regards Peter
Integer Factorial(int n)
{
if(n < 2) return 1;
Integer p = 1, r = 1;
N = 1;
int h = 0, shift = 0, high = 1;
int log2n = Floor(Log2(n));
while(h != n)
{
shift += h;
h = n >> log2n--;
int len = high;
high = (h & 1) == 1 ? h : h - 1;
len = (high - len) / 2;
if(len > 0)
{
p *= Product(len);
r *= p;
}
}
return r << shift;
}
int N;
Integer Product(int n)
{
int m = n / 2;
if( m == 0 ) return new Integer(N += 2);
if( n == 2 ) return new Integer((N += 2)*(N += 2));
return Product(n - m) * Product(m);
} |
|
| Back to top |
|
|
Richard J. Fateman
science forum addict
Joined: 04 May 2005
Posts: 81
|
| Posted: Wed Apr 05, 2006 9:19 pm Post subject:
Re: factorial of largish numbers
|
|
|
Peter Luschny wrote:
| Quote: | Perhaps someone in this group could give me a
helping hand and show me the Lisp implementation of the algorithm
which I like best? I put it in an appendix, written in a C-like language.
........ |
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)))))) |
|
| Back to top |
|
|
Pascal Bourguignon
science forum beginner
Joined: 07 Feb 2006
Posts: 4
|
| Posted: Wed Apr 05, 2006 9:48 pm Post subject:
Re: factorial of largish numbers
|
|
|
"Richard J. Fateman" <fateman@eecs.berkeley.edu> writes:
| Quote: | I invite comments on a paper posted at
http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf
on computing factorial of, say, 20,000 or more.
Programs in lisp and some comparison with Maple and
Mathematica. are included.
I've used (Allegro 7.0) Lisp + GMP for best times, but maybe
you have other ideas.
Thanks
RJF
|
I would write fact5 as:
(defun fact5 (x)
(if (typep x '(integer 0 12))
(aref #(1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800 479001600)
x)
(error "~S is not a valid argument to FACT5" x)))
clisp already uses GMP. On my Athlon 1200 MHz, here is the native
factorial time:
[290]> (time (defparameter *e* (ext:! 20000)))
Real time: 0.070124 sec.
Run time: 0.072005 sec.
Space: 696224 Bytes
*E*
For reference, here's the maxfact time:
[306]> (ext:gc)
1588278
[307]> (time (defparameter *e* (ext:! 20000)))
Real time: 0.069836 sec.
Run time: 0.068004 sec.
Space: 696224 Bytes
*E*
[308]> (mapcar 'compile '(maxfact cloop loopprod))
(MAXFACT CLOOP LOOPPROD)
[309]> (ext:gc)
1587970
[310]> (time (defparameter *m* (maxfact 20000)))
Real time: 0.113327 sec.
Run time: 0.112007 sec.
Space: 3504744 Bytes
GC: 2, GC time: 0.032002 sec.
*M*
[311]> (assert (= *m* *e*))
NIL
[312]>
------------------------------------------------------------------------
;; This, taken from the pdf file:
(defun cloop (n)
(let* ((start (list 1)) (end start))
(dotimes (i (1- n) (setf (cdr end) start))
(setq start (cons (+ 2 i) start)))))
(setf *print-circle* t)
(defun loopprod (l)
"efficient product of all the approximately equal numbers in the loop
trying to keep the sizes of inputs approximately balanced. the
circular loop l is destroyed in the process. "
(cond ((eq (cdr l) l)(car l))
(t (setf (car l) (* (car l) (cadr l)))
(setf (cdr l) (cddr l))
(loopprod (cdr l)))))
(defun maxfact (n)
(let* ((z (if (< n 100) 5 100))
(aloop (cloop z)))
(loop for i from (1+ z) to n
do (setf (car aloop) (* (car aloop) i))
(setf aloop (cdr aloop)))
(loopprod aloop)))
;; gives:
(maxfact 2) --> 120
------------------------------------------------------------------------
Otherwise, it seems to me that in addition to the cloop, it would be
worthwhile to extract the even. For 20000!, we'll have 10000
doubles, so we can easily avoid 10000 bits in the bignum
multiplications, and just shift the last product.
We can even avoid some more bits for the multiples of other powers of two.
So, I'd try to use a cloop of a power of two instead of 5 or 100, and
instead of multiplying i, I'd multiply (ash i -2^p) to the cloop and
increment p to a counter.
Something like:
(defun cloop (n) ; #1=( 1 1 1 ... 1 . #1# )
(let ((result (make-list n :initial-element 1)))
(setf (cdr (last result)) result)))
(defmacro unroll (z aloop pot i &optional (n nil n-p))
;; Should be a macrolet inside maxfact,
;; but it's easier to debug as a defmacro.
(let ((end (gensym)))
`(,(if n-p 'tagbody 'progn)
,@(loop
:for u :from 1 :to z ; =i[z]
:for p = (loop
:for p :from 0
:for x = u :then (truncate x 2)
:while (evenp x)
:finally (return p)) ; i\2^p
:if n-p
:collect `(when (> ,i ,n) (go ,end))
:end
:collect (if (plusp p)
`(setf (car ,aloop) (* (car ,aloop) (ash ,i ,(- p)))
,pot (+ ,pot ,p))
`(setf (car ,aloop) (* (car ,aloop) ,i)))
:collect `(setf ,aloop (cdr ,aloop) ,i (1+ ,i)))
,@(when n-p `(,end)))))
(defun maxfact (n)
(macrolet ((compute (z)
`(let ((aloop (cloop ,z))
(pot 0)
(m (1+ (* ,z (truncate n ,z)))))
;; (print `(m = ,m))
(loop
:with i = 1
:while (< i m)
:do (unroll ,z aloop pot i))
;;(print aloop) (print `(pot = ,pot)) (print `(i = ,i))
(unroll ,z aloop pot m n)
;; (print aloop) (print `(pot = ,pot)) (print `(m = ,m))
(ash (loopprod aloop) pot))))
(if (< n 128)
(compute 4)
(compute 128))))
[296]> (ext:gc)
1597184
[297]> (time (defparameter *e* (ext:! 20000)))
Real time: 0.070717 sec.
Run time: 0.072005 sec.
Space: 696224 Bytes
*E*
[298]> (mapcar 'compile '(maxfact cloop loopprod))
(MAXFACT CLOOP LOOPPROD)
[299]> (ext:gc)
1589146
[300]> (time (defparameter *m* (maxfact 20000)))
Real time: 0.081703 sec.
Run time: 0.080005 sec.
Space: 2667600 Bytes
GC: 1, GC time: 0.016001 sec.
*M*
[301]> (assert (= *m* *e*))
NIL
[316]> (dotimes (i 512) (assert (= (ext:! i) (maxfact i))))
NIL
So I reach almost the internal #+clisp ext:! time, but still have some
bignum consing that the internal ext:! is able to avoid using directly
GMP.
--
__Pascal Bourguignon__ http://www.informatimago.com/
READ THIS BEFORE OPENING PACKAGE: According to certain suggested
versions of the Grand Unified Theory, the primary particles
constituting this product may decay to nothingness within the next
four hundred million years. |
|
| Back to top |
|
|
daly@axiom-developer.org
science forum beginner
Joined: 10 Nov 2005
Posts: 12
|
| Posted: Wed Apr 05, 2006 10:17 pm Post subject:
Re: factorial of largish numbers
|
|
|
I suggest a different algorithm using binary arithmetic.
Call it the "trailing zeros" algorithm.
1 = 1
2 = 10 = 1 * 2^1 ==> 2! = 2
3 = 11 = 3 * 2^0 ==> 3! = 3*1*2^1 = 110 = 6
4 = 100 = 1 * 2^2 ==> 4! = 3 * 2^3 = 11000 = 24
5 = 101 = 5 * 2^2 ==> 5! = 5 * 3 * 2^3 = 1111000 = 120
6 = 110 = 3 * 2^1 ==> 6! = 3 * 5 * 3 * 2^4 = 10 1101 0000 = 720
7 = 111 = 7 * 2^0 ==> 7! = 7 * 3 * 5 * 3 * 2^4 = 1 0011 1011 0000 =
5040
8 = 1000 = 1 * 2^3 ==> 8! = 7 * 3 * 5 * 3 * 2^7 = 1 0011 1011 0000 0000
= 40320
notice that we're using shifting to rid the world of the large trailing
zeros
and that certain computations can be see speedup because we're only
keeping significant factor bits. 7! and 8! obviously are the same
computation.
further speedups are possible with table lookups of values of products
with
certain hex values, rather like special multiplication tables for hex
products
of, say, 7*3, etc.
if you "sieve out" all of the factors that have no trailing zeros you
can get
the "content" of the number without the "fat".
there is also "fat" content with the intermediate zeros (e.g.
100....0001)
which could be removed by dealing with the number first in hex
segments,
then in word segments, then doubleword, quadword, etc.
I'm sure this exists somewhere but cannot find it in a google search.
Tim Daly |
|
| Back to top |
|
|
Christophe Rhodes
science forum beginner
Joined: 05 Apr 2006
Posts: 1
|
| Posted: Wed Apr 05, 2006 10:21 pm Post subject:
Re: factorial of largish numbers
|
|
|
[ note followups ]
Alexander Schmolck <a.schmolck@gmail.com> writes:
| Quote: | will compute k(20000) just fine in any non-ancient version of python
(it is in fact several orders of magnitude faster than compiled
cmucl which I eventually ^C'ed and seems to be about the same speed
as clisp).
|
Be careful of what you are measuring. You don't quote how you made
your observations, but you might be running into a performance problem
in cmucl, where the factorial is in fact computed acceptably fast, but
printing of return value (by the REPL) takes a very long time.
Christophe |
|
| Back to top |
|
|
Peter Luschny
science forum beginner
Joined: 13 May 2005
Posts: 36
|
| Posted: Wed Apr 05, 2006 10:57 pm Post subject:
Re: factorial of largish numbers
|
|
|
| Quote: | Richard J. Fateman schrieb:
Peter Luschny wrote:
Perhaps someone in this group could give me a
helping hand and show me the Lisp implementation of the algorithm
which I like best?
here it is in lisp
|
Superbe! Merci beaucoup! |
|
| Back to top |
|
|
Alexander Schmolck
science forum beginner
Joined: 05 Apr 2006
Posts: 4
|
| Posted: Wed Apr 05, 2006 11:04 pm Post subject:
Re: factorial of largish numbers
|
|
|
Peter Luschny <spamgrube@luschny.de> writes:
| Quote: | Richard J. Fateman schrieb:
Peter Luschny wrote:
Perhaps someone in this group could give me a
helping hand and show me the Lisp implementation of the algorithm
which I like best?
here it is in lisp
|
I think you will also need something like (untested)
(defacro while (condition &body body)
`(loop while ,condition do (progn ,@body)))
'as |
|
| Back to top |
|
|
Alexander Schmolck
science forum beginner
Joined: 05 Apr 2006
Posts: 4
|
| Posted: Wed Apr 05, 2006 11:06 pm Post subject:
Re: factorial of largish numbers
|
|
|
Alexander Schmolck <a.schmolck@gmail.com> writes:
| Quote: | I think you will also need something like (untested)
(defacro while (condition &body body)
(defmacro |
|
|
| Back to top |
|
|
Rob Warnock
science forum beginner
Joined: 05 Apr 2006
Posts: 4
|
| Posted: Wed Apr 05, 2006 11:23 pm Post subject:
Re: factorial of largish numbers
|
|
|
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 |
|
| Back to top |
|
|
Robert B. Israel
science forum Guru
Joined: 24 Mar 2005
Posts: 2151
|
| Posted: Thu Apr 06, 2006 3:45 am Post subject:
Re: factorial of largish numbers
|
|
|
In article <e10udn$2jgi$1@agate.berkeley.edu>,
Richard J. Fateman <fateman@eecs.berkeley.edu> wrote:
| Quote: | I invite comments on a paper posted at
http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf
on computing factorial of, say, 20,000 or more.
Programs in lisp and some comparison with Maple and
Mathematica. are included.
I've used (Allegro 7.0) Lisp + GMP for best times, but maybe
you have other ideas.
Thanks
RJF
|
I believe that what current versions of Maple use is the GMP factorial
function mpz_fac_ui: see <http://www.swox.com/gmp/gmp-man-4.2.pdf>.
It's quite fast: 40000! took only 0.093 seconds on my machine (Pentium
4, 3 GHz running Maple 10.03 Classic under Windows XP). Section
16.7.2 of that document tells about the algorithm used.
Robert Israel israel@math.ubc.ca
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia Vancouver, BC, Canada |
|
| Back to top |
|
|
Google
|
|
| Back to top |
|
|
|
|
The time now is Wed Jan 07, 2009 8:29 pm | All times are GMT
|
|
Bankruptcy | Mortgage Calculator | Mortgages | Debt Consolidation | Credit Cards
|
|
Copyright © 2004-2005 DeniX Solutions SRL
|
|
Other DeniX Solutions sites:
Electronics forum |
Medicine forum |
Unix/Linux blog |
Unix/Linux documentation |
Unix/Linux forums
|
Powered by phpBB © 2001, 2005 phpBB Group
|
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
Forum index
»
Science and Technology
»
Math
»
Symbolic
