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