factorial of largish numbers

Klueless · science forum beginner Joined: 30 May 2005 Posts: 19

"Richard J. Fateman" <fateman@eecs.berkeley.edu> wrote in message news:e10udn$2jgi$1@agate.berkeley.edu...

Richard Fateman · science forum Guru Wannabe Joined: 03 May 2005 Posts: 181

I looked at the algorithm in gmp, and timed calling it directly from Lisp.
It is
somewhat faster than running a good factorial algorithm in Lisp (There is no
reason for it to be much faster... it is still doing the same big
multiplies, and
doesn't really have a better algorithm).

As I recall it makes some number of buckets and tries to accumulate
products in all of them of about equal size. (Like the version in Maxima.)
This is not too bad, but presumably not as good as Peter Luschny's best,
which requires interleaving a prime number sieve with the calculation,
similar
to what Tim Daly suggests.

I tried running the version (in Lisp)of Lucshny's posted earlier, but
modified to run
in Lisp + GMP. It is good, but not as good as some others, probably
because it requires creating lots of GMPZ numbers to return from the
"prod" routine, instead of storing those numbers destructively on top
of each other.

RJF

"Robert Israel" <israel@math.ubc.ca> wrote in message
news:e122tl$ite$1@nntp.itservices.ubc.ca...

Peter Luschny · science forum beginner Joined: 13 May 2005 Posts: 36

Rob Warnock wrote:

Frank Buss · science forum beginner Joined: 06 Apr 2006 Posts: 3

Rob Warnock wrote:

Peter Luschny · science forum beginner Joined: 13 May 2005 Posts: 36

Frank Buss wrote:

Mark Lawton · science forum beginner Joined: 22 Mar 2006 Posts: 11

Hi Richard,

Not really a serious comment - just interesting to note that the
following Maxima function computes 20,000! correctly to 9 digits:

fact(x) := (y: sum(float(log(n)),n,1,x),
y: float(y/log(10)),
print(10^(y-fix(y)), "e", fix(y)),"")$

If I set floating point precision to 100, and rewrite the function as a
bigfloat, it computes it correctly to 92 digits:

fact(x) := (y: sum(bfloat(log(n)),n,1,x),
y: bfloat(y/log(10)),
print(10^(y-fix(y)), "e", fix(y)),"")$

So if you don't need your factorial's full precision, you can use this
method. I remembered this trick of obtaining factorials by summing the
logs from the early 1980s, when the only tool I had was my TI-58C
programmable calculator.

Richard J. Fateman wrote:

Mark Lawton · science forum beginner Joined: 22 Mar 2006 Posts: 11

I missed the obvious point that there is no need to use logs when using
a CAS - sorry.

fact(x) := (y: 1,
for i: 2 thru x do
y: bfloat(y * i),
y)$

Using the above Maxima function with floating point precision set to
100, the result of 20,000! is produced in 0.5 seconds (circa 2 ghz PC)
and all 100 digits produced are correct.

Mark Lawton wrote:

Richard Fateman · science forum Guru Wannabe Joined: 03 May 2005 Posts: 181

Thanks for all the comments!

The sparse recursive factorial algorithm as Peter L wrote in Java
does a lot of allocation of numbers (as bignums) when they
are really just small numbers. A more efficient way of doing
could be to allocate two bignums for p and r and overwrite them. The
"API" offered by GMP allows you to do this overwriting, and the Lisp
version I provide allows access to these destructive versions.
(The Python implementation, according to at least one note I
found by Googling, did not provide such access.) Anyway, here
is a Lisp version for sparse recursive, with a modification of
the "prod" routine to take another argument, ans. If the old prod(n)
computed a value C(n), the new gprod(n,ans) produces a result
like ans:=C(n)*ans, destroying the old value of ans.
The ordering and size of multiplications is perhaps not as good
in this re-working. I believe it is just an unrolling of the recurrence
f(n,m):=f(n,2*m)*f(n-m,2*m).

This sprfact is not quite as fast as the "k" program, presumably because
it doesn't do as good a job of keeping numbers equally small, and
accumulating
their product "gently". It uses almost no storage except inside the GMP
allocator.

This is the program (it won't run unless you have my gmp library). If
you are familiar with GMP, the programs mpz_.... are named the same as
the GMP entries. I've also changed the naming and the while, and labels to
be
fully ANSI common lisp. I've also prefixed the fixnum operators with
cl:: just in case this might make it compile better. I've overloaded +, *,
etc
with GMP generic functions.

(defun gsprfact(n);;works. in :gmp package
(declare (fixnum n) (optimize (speed 3)(safety 0)) )
(let ((p (into 1)) (r (into 1)) ;initialize by converting 1 into GMPZ
numbers.
(NN 1) (log2n (floor (log n 2)))
(h 0) (shift 0) (high 1) (len 0))
(declare (fixnum log2n h shift high len NN))
(labels ((gprod(n ans)
(declare (fixnum n NN))
(let ((m (ash n -1)))
(declare (fixnum m))
(cond ((cl::= m 0)
(mpz_mul_ui ans ans
(cl::incf NN 2)))
((cl::= n 2)
(mpz_mul_ui ans ans
(cl::* (cl::incf NN 2)
(cl::incf NN 2))))
(t (gprod (cl::- n m) (gprod m ans))))
ans)))
(cond ((< n 0)(error "bad arg ~s to factorial" n))
((< n 2) (gmp::into 1))
(t
(loop while (cl::/= h n) do
(cl::incf shift h)
(setq h (ash n (cl::- log2n)))
(cl::decf log2n)
(setq len high)
(setq high (if (oddp h) h (cl::1- h)))
(setq len (ash (cl::- high len) -1))
(cond ((cl::> len 0)
(gprod len (gmpz-z p));change p
(mpz_mul (gmpz-z r)(gmpz-z r)(gmpz-z p)) )))
(mpz_mul_2exp (gmpz-z r) (gmpz-z r) shift);; arithmetic shift
r)))))

Here's another version, ANSI CL, but less hacking.

(defun sprfact(n)
(declare (fixnum n))
(let ((p 1) (r 1) (NN 1) (log2n (floor (log n 2)))
(h 0) (shift 0) (high 1) (len 0))
(declare (fixnum log2n h shift high len NN))
(labels ((prod(n)
(declare (fixnum n))
(let ((m (ash n -1)))
(cond ((= m 0) (incf NN 2))
((= n 2) (* (incf NN 2)(incf NN 2)))
(t (* (prod (- n m)) (prod m)))))))

(cond ((< n 0)(error "bad arg ~s to factorial" n))
((< n 2) 1)
(t
(loop while (/= h n) do
(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))))))

Geoffrey Summerhayes · science forum beginner Joined: 06 Apr 2006 Posts: 1

<daly@axiom-developer.org> wrote in message
news:1144275433.752025.9500@u72g2000cwu.googlegroups.com...

Richard Fateman · science forum Guru Wannabe Joined: 03 May 2005 Posts: 181

The algorithms given by Peter Luschny cover a lot of ground.

I think Tim's observation is basically that if k(n,1) == n!, then

k(n,m) := if even(n) then 2^n*k(n/2,m)*k(n-m,m) // rule for (2*n)!
else n*k(n-1,m) // where n-1 will be even.

Note that multiplication by 2^n for integer n can be written as an
arithmetic shift in languages that permit it. In lisp, (ash x n) .

There are similar rules for other factorials of numbers
with other divisors, e.g. if n is divisible by 3.

Multiplying x by 3 can be done by (+ x (ash x 1)) though mult by
powers of 3 are not as fast as binary shifts.

There is a potential for speedup by observing that the answer can always be
factored,
but then if you want to multiply it through to get a single bignum, you
have to figure out the right sequence of multiplies to save time.
Interleaving a prime number sieve with the computation can win, but
only by about a factor of 2, according to Luschny.

For GMP-based systems you should try to get equal-sized inputs to multiply.
For systems which do not use 'advanced' bignum multiply, different
heuristics apply: for them it may be faster to have inputs that are
smallnum X bignum.

There is perhaps a simpler way of expressing some of the
idea in that Geoff's version, eg.

;; (h2 n) is factorial of n, but faster.
(defun hf2(n);; n is positive integer
(if (< n 2) 1
(if (evenp n)
(let ((q (ash n -1)))
(ash (* (hf2 q)(f (- n 1) 2)) q))
(* n (hf2 (1- n))))))

(defun f (n m) (if (<= n 1) 1 (* n (f (- n m) m)))) ;; could be hacked, too.

Richard Fateman · science forum Guru Wannabe Joined: 03 May 2005 Posts: 181

Peter, Frank (and others who seem to have been caught by this
problem)...

Here is a more lisp-like version of the split-recursive factorial

(defun fact (n)
(let ((shift 0))
(labels
((kg1 (n m)
(cond ((and (evenp n)(> m 1))
(incf shift (ash n -1))
(kg1 (ash n -1) (ash m -1)))
((<= n m) n)
(t (* (kg1 n (ash m 1))
(kg1 (- n m)(ash m 1)))))))
(ash (kg1 n 1) shift) )))

It runs at the same speed, pretty much, as the unrolled program posted
earlier. I recommend that
you set the optimization flags to (speed 3)(safety 0) and declare that
n,m,shift are fixnums.

To make it run competitively in Allegro CL, I changed it to use the GMP
library,
but then had to change it more so it used a small collection of preallocated
GMP numbers
instead of madly allocating small GMP numbers just to toss them out after
one or
two arithmetic operations. A lisp system using GMP "natively" might be able
to
use the above program. Discussion at
http://www.cs.berkeley.edu/~fateman/papers/factorial.pdf
has been updated.
Thanx
RJF

"Peter Luschny" <spamgrube@luschny.de> wrote in message
news:e12qh7$noh$1@online.de...

Peter Luschny · science forum beginner Joined: 13 May 2005 Posts: 36

Richard Fateman wrote:

Google