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