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