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 "dat.h"
6 b3f61791 2004-03-21 devnull mpfactorial(ulong n)
10 b3f61791 2004-03-21 devnull unsigned cnt;
11 b3f61791 2004-03-21 devnull int max, mmax;
12 b3f61791 2004-03-21 devnull mpdigit p, pp[2];
13 b3f61791 2004-03-21 devnull mpint *r, *s, *stk[31];
16 b3f61791 2004-03-21 devnull max = mmax = -1;
18 b3f61791 2004-03-21 devnull r = mpnew(0);
19 b3f61791 2004-03-21 devnull for(k=2; k<=n; k++){
20 b3f61791 2004-03-21 devnull pp[0] = 0;
21 b3f61791 2004-03-21 devnull pp[1] = 0;
22 b3f61791 2004-03-21 devnull mpvecdigmuladd(&p, 1, (mpdigit)k, pp);
23 b3f61791 2004-03-21 devnull if(pp[1] == 0) /* !overflow */
24 b3f61791 2004-03-21 devnull p = pp[0];
27 b3f61791 2004-03-21 devnull if((cnt & 1) == 0){
28 b3f61791 2004-03-21 devnull s = stk[max];
29 b3f61791 2004-03-21 devnull mpbits(r, Dbits*(s->top+1+1));
30 b3f61791 2004-03-21 devnull memset(r->p, 0, Dbytes*(s->top+1+1));
31 b3f61791 2004-03-21 devnull mpvecmul(s->p, s->top, &p, 1, r->p);
32 b3f61791 2004-03-21 devnull r->sign = 1;
33 b3f61791 2004-03-21 devnull r->top = s->top+1+1; /* XXX: norm */
34 b3f61791 2004-03-21 devnull mpassign(r, s);
35 b3f61791 2004-03-21 devnull for(i=4; (cnt & (i-1)) == 0; i=i<<1){
36 b3f61791 2004-03-21 devnull mpmul(stk[max], stk[max-1], r);
37 b3f61791 2004-03-21 devnull mpassign(r, stk[max-1]);
42 b3f61791 2004-03-21 devnull if(max > mmax){
44 b3f61791 2004-03-21 devnull if(max > nelem(stk))
46 b3f61791 2004-03-21 devnull stk[max] = mpnew(Dbits);
48 b3f61791 2004-03-21 devnull stk[max]->top = 1;
49 b3f61791 2004-03-21 devnull stk[max]->p[0] = p;
51 b3f61791 2004-03-21 devnull p = (mpdigit)k;
54 b3f61791 2004-03-21 devnull if(max < 0){
55 b3f61791 2004-03-21 devnull mpbits(r, Dbits);
56 b3f61791 2004-03-21 devnull r->top = 1;
57 b3f61791 2004-03-21 devnull r->sign = 1;
58 b3f61791 2004-03-21 devnull r->p[0] = p;
60 b3f61791 2004-03-21 devnull s = stk[max--];
61 b3f61791 2004-03-21 devnull mpbits(r, Dbits*(s->top+1+1));
62 b3f61791 2004-03-21 devnull memset(r->p, 0, Dbytes*(s->top+1+1));
63 b3f61791 2004-03-21 devnull mpvecmul(s->p, s->top, &p, 1, r->p);
64 b3f61791 2004-03-21 devnull r->sign = 1;
65 b3f61791 2004-03-21 devnull r->top = s->top+1+1; /* XXX: norm */
68 b3f61791 2004-03-21 devnull while(max >= 0)
69 b3f61791 2004-03-21 devnull mpmul(r, stk[max--], r);
70 b3f61791 2004-03-21 devnull for(max=mmax; max>=0; max--)
71 b3f61791 2004-03-21 devnull mpfree(stk[max]);
72 b3f61791 2004-03-21 devnull mpnorm(r);
73 b3f61791 2004-03-21 devnull return r;