Subversion Repositories planix.SVN

Rev

Rev 2 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 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 
}