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