mpfactorial.c (1347B)
1 #include "os.h" 2 #include <mp.h> 3 #include "dat.h" 4 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]; 14 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 } 67 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 }