1 b3f61791 2004-03-21 devnull #include "os.h"
 2 b3f61791 2004-03-21 devnull #include <mp.h>
 3 b3f61791 2004-03-21 devnull #include <libsec.h>
 4 b3f61791 2004-03-21 devnull #include "dat.h"
 5 b3f61791 2004-03-21 devnull 
 6 b3f61791 2004-03-21 devnull mpint*
 7 b3f61791 2004-03-21 devnull mpfactorial(ulong n)
 8 b3f61791 2004-03-21 devnull {
 9 b3f61791 2004-03-21 devnull 	int i;
10 b3f61791 2004-03-21 devnull 	ulong k;
11 b3f61791 2004-03-21 devnull 	unsigned cnt;
12 b3f61791 2004-03-21 devnull 	int max, mmax;
13 b3f61791 2004-03-21 devnull 	mpdigit p, pp[2];
14 b3f61791 2004-03-21 devnull 	mpint *r, *s, *stk[31];
15 b3f61791 2004-03-21 devnull 
16 b3f61791 2004-03-21 devnull 	cnt = 0;
17 b3f61791 2004-03-21 devnull 	max = mmax = -1;
18 b3f61791 2004-03-21 devnull 	p = 1;
19 b3f61791 2004-03-21 devnull 	r = mpnew(0);
20 b3f61791 2004-03-21 devnull 	for(k=2; k<=n; k++){
21 b3f61791 2004-03-21 devnull 		pp[0] = 0;
22 b3f61791 2004-03-21 devnull 		pp[1] = 0;
23 b3f61791 2004-03-21 devnull 		mpvecdigmuladd(&p, 1, (mpdigit)k, pp);
24 b3f61791 2004-03-21 devnull 		if(pp[1] == 0)	/* !overflow */
25 b3f61791 2004-03-21 devnull 			p = pp[0];
26 b3f61791 2004-03-21 devnull 		else{
27 b3f61791 2004-03-21 devnull 			cnt++;
28 b3f61791 2004-03-21 devnull 			if((cnt & 1) == 0){
29 b3f61791 2004-03-21 devnull 				s = stk[max];
30 b3f61791 2004-03-21 devnull 				mpbits(r, Dbits*(s->top+1+1));
31 b3f61791 2004-03-21 devnull 				memset(r->p, 0, Dbytes*(s->top+1+1));
32 b3f61791 2004-03-21 devnull 				mpvecmul(s->p, s->top, &p, 1, r->p);
33 b3f61791 2004-03-21 devnull 				r->sign = 1;
34 b3f61791 2004-03-21 devnull 				r->top = s->top+1+1;		/* XXX: norm */
35 b3f61791 2004-03-21 devnull 				mpassign(r, s);
36 b3f61791 2004-03-21 devnull 				for(i=4; (cnt & (i-1)) == 0; i=i<<1){
37 b3f61791 2004-03-21 devnull 					mpmul(stk[max], stk[max-1], r);
38 b3f61791 2004-03-21 devnull 					mpassign(r, stk[max-1]);
39 b3f61791 2004-03-21 devnull 					max--;
40 b3f61791 2004-03-21 devnull 				}
41 b3f61791 2004-03-21 devnull 			}else{
42 b3f61791 2004-03-21 devnull 				max++;
43 b3f61791 2004-03-21 devnull 				if(max > mmax){
44 b3f61791 2004-03-21 devnull 					mmax++;
45 b3f61791 2004-03-21 devnull 					if(max > nelem(stk))
46 b3f61791 2004-03-21 devnull 						abort();
47 b3f61791 2004-03-21 devnull 					stk[max] = mpnew(Dbits);
48 b3f61791 2004-03-21 devnull 				}
49 b3f61791 2004-03-21 devnull 				stk[max]->top = 1;
50 b3f61791 2004-03-21 devnull 				stk[max]->p[0] = p;
51 b3f61791 2004-03-21 devnull 			}
52 b3f61791 2004-03-21 devnull 			p = (mpdigit)k;
53 b3f61791 2004-03-21 devnull 		}
54 b3f61791 2004-03-21 devnull 	}
55 b3f61791 2004-03-21 devnull 	if(max < 0){
56 b3f61791 2004-03-21 devnull 		mpbits(r, Dbits);
57 b3f61791 2004-03-21 devnull 		r->top = 1;
58 b3f61791 2004-03-21 devnull 		r->sign = 1;
59 b3f61791 2004-03-21 devnull 		r->p[0] = p;
60 b3f61791 2004-03-21 devnull 	}else{
61 b3f61791 2004-03-21 devnull 		s = stk[max--];
62 b3f61791 2004-03-21 devnull 		mpbits(r, Dbits*(s->top+1+1));
63 b3f61791 2004-03-21 devnull 		memset(r->p, 0, Dbytes*(s->top+1+1));
64 b3f61791 2004-03-21 devnull 		mpvecmul(s->p, s->top, &p, 1, r->p);
65 b3f61791 2004-03-21 devnull 		r->sign = 1;
66 b3f61791 2004-03-21 devnull 		r->top = s->top+1+1;		/* XXX: norm */
67 b3f61791 2004-03-21 devnull 	}
68 b3f61791 2004-03-21 devnull 
69 b3f61791 2004-03-21 devnull 	while(max >= 0)
70 b3f61791 2004-03-21 devnull 		mpmul(r, stk[max--], r);
71 b3f61791 2004-03-21 devnull 	for(max=mmax; max>=0; max--)
72 b3f61791 2004-03-21 devnull 		mpfree(stk[max]);
73 b3f61791 2004-03-21 devnull 	mpnorm(r);
74 b3f61791 2004-03-21 devnull 	return r;
75 b3f61791 2004-03-21 devnull }