xref: /plan9/sys/src/libmp/port/mpfactorial.c (revision f7db61556a577f91350f05658e9c0724969b02c3)
1da51d93aSDavid du Colombier #include "os.h"
2da51d93aSDavid du Colombier #include <mp.h>
3da51d93aSDavid du Colombier #include <libsec.h>
4da51d93aSDavid du Colombier #include "dat.h"
5da51d93aSDavid du Colombier 
6da51d93aSDavid du Colombier mpint*
mpfactorial(ulong n)7da51d93aSDavid du Colombier mpfactorial(ulong n)
8da51d93aSDavid du Colombier {
9da51d93aSDavid du Colombier 	int i;
10da51d93aSDavid du Colombier 	ulong k;
11da51d93aSDavid du Colombier 	unsigned cnt;
12da51d93aSDavid du Colombier 	int max, mmax;
13da51d93aSDavid du Colombier 	mpdigit p, pp[2];
14da51d93aSDavid du Colombier 	mpint *r, *s, *stk[31];
15da51d93aSDavid du Colombier 
16da51d93aSDavid du Colombier 	cnt = 0;
17da51d93aSDavid du Colombier 	max = mmax = -1;
18da51d93aSDavid du Colombier 	p = 1;
19da51d93aSDavid du Colombier 	r = mpnew(0);
20da51d93aSDavid du Colombier 	for(k=2; k<=n; k++){
21da51d93aSDavid du Colombier 		pp[0] = 0;
22da51d93aSDavid du Colombier 		pp[1] = 0;
23da51d93aSDavid du Colombier 		mpvecdigmuladd(&p, 1, (mpdigit)k, pp);
24da51d93aSDavid du Colombier 		if(pp[1] == 0)	/* !overflow */
25da51d93aSDavid du Colombier 			p = pp[0];
26da51d93aSDavid du Colombier 		else{
27da51d93aSDavid du Colombier 			cnt++;
28da51d93aSDavid du Colombier 			if((cnt & 1) == 0){
29da51d93aSDavid du Colombier 				s = stk[max];
30da51d93aSDavid du Colombier 				mpbits(r, Dbits*(s->top+1+1));
31da51d93aSDavid du Colombier 				memset(r->p, 0, Dbytes*(s->top+1+1));
32da51d93aSDavid du Colombier 				mpvecmul(s->p, s->top, &p, 1, r->p);
33da51d93aSDavid du Colombier 				r->sign = 1;
34da51d93aSDavid du Colombier 				r->top = s->top+1+1;		/* XXX: norm */
35da51d93aSDavid du Colombier 				mpassign(r, s);
36da51d93aSDavid du Colombier 				for(i=4; (cnt & (i-1)) == 0; i=i<<1){
37da51d93aSDavid du Colombier 					mpmul(stk[max], stk[max-1], r);
38da51d93aSDavid du Colombier 					mpassign(r, stk[max-1]);
39da51d93aSDavid du Colombier 					max--;
40da51d93aSDavid du Colombier 				}
41da51d93aSDavid du Colombier 			}else{
42da51d93aSDavid du Colombier 				max++;
43da51d93aSDavid du Colombier 				if(max > mmax){
44da51d93aSDavid du Colombier 					mmax++;
45*f7db6155SDavid du Colombier 					if(max >= nelem(stk))
46da51d93aSDavid du Colombier 						abort();
47da51d93aSDavid du Colombier 					stk[max] = mpnew(Dbits);
48da51d93aSDavid du Colombier 				}
49da51d93aSDavid du Colombier 				stk[max]->top = 1;
50da51d93aSDavid du Colombier 				stk[max]->p[0] = p;
51da51d93aSDavid du Colombier 			}
52da51d93aSDavid du Colombier 			p = (mpdigit)k;
53da51d93aSDavid du Colombier 		}
54da51d93aSDavid du Colombier 	}
55da51d93aSDavid du Colombier 	if(max < 0){
56da51d93aSDavid du Colombier 		mpbits(r, Dbits);
57da51d93aSDavid du Colombier 		r->top = 1;
58da51d93aSDavid du Colombier 		r->sign = 1;
59da51d93aSDavid du Colombier 		r->p[0] = p;
60da51d93aSDavid du Colombier 	}else{
61da51d93aSDavid du Colombier 		s = stk[max--];
62da51d93aSDavid du Colombier 		mpbits(r, Dbits*(s->top+1+1));
63da51d93aSDavid du Colombier 		memset(r->p, 0, Dbytes*(s->top+1+1));
64da51d93aSDavid du Colombier 		mpvecmul(s->p, s->top, &p, 1, r->p);
65da51d93aSDavid du Colombier 		r->sign = 1;
66da51d93aSDavid du Colombier 		r->top = s->top+1+1;		/* XXX: norm */
67da51d93aSDavid du Colombier 	}
68da51d93aSDavid du Colombier 
69da51d93aSDavid du Colombier 	while(max >= 0)
70da51d93aSDavid du Colombier 		mpmul(r, stk[max--], r);
71da51d93aSDavid du Colombier 	for(max=mmax; max>=0; max--)
72da51d93aSDavid du Colombier 		mpfree(stk[max]);
73da51d93aSDavid du Colombier 	mpnorm(r);
74da51d93aSDavid du Colombier 	return r;
75da51d93aSDavid du Colombier }
76