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