ref: 1018d7a7091fd5e0433f501a3474cdc3fec322cc
dir: /libmp/port/mpfactorial.c/
#include "os.h" #include <mp.h> #include <libsec.h> #include "dat.h" mpint* mpfactorial(ulong n) { int i; ulong k; unsigned cnt; int max, mmax; mpdigit p, pp[2]; mpint *r, *s, *stk[31]; cnt = 0; max = mmax = -1; p = 1; r = mpnew(0); for(k=2; k<=n; k++){ pp[0] = 0; pp[1] = 0; mpvecdigmuladd(&p, 1, (mpdigit)k, pp); if(pp[1] == 0) /* !overflow */ p = pp[0]; else{ cnt++; if((cnt & 1) == 0){ s = stk[max]; mpbits(r, Dbits*(s->top+1+1)); memset(r->p, 0, Dbytes*(s->top+1+1)); mpvecmul(s->p, s->top, &p, 1, r->p); r->sign = 1; r->top = s->top+1+1; /* XXX: norm */ mpassign(r, s); for(i=4; (cnt & (i-1)) == 0; i=i<<1){ mpmul(stk[max], stk[max-1], r); mpassign(r, stk[max-1]); max--; } }else{ max++; if(max > mmax){ mmax++; if(max >= nelem(stk)){ while(--max >= 0) mpfree(stk[max]); mpfree(r); sysfatal("mpfactorial: stack overflow"); } stk[max] = mpnew(Dbits); } stk[max]->top = 1; stk[max]->p[0] = p; } p = (mpdigit)k; } } if(max < 0){ mpbits(r, Dbits); r->top = 1; r->sign = 1; r->p[0] = p; }else{ s = stk[max--]; mpbits(r, Dbits*(s->top+1+1)); memset(r->p, 0, Dbytes*(s->top+1+1)); mpvecmul(s->p, s->top, &p, 1, r->p); r->sign = 1; r->top = s->top+1+1; /* XXX: norm */ } while(max >= 0) mpmul(r, stk[max--], r); for(max=mmax; max>=0; max--) mpfree(stk[max]); mpnorm(r); return r; }