Blob
1 #include "os.h"2 #include <mp.h>3 #include <libsec.h>4 #include "dat.h"6 mpint*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];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 }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 }