diff options
author | rsc <devnull@localhost> | 2004-03-21 14:06:38 +0000 |
---|---|---|
committer | rsc <devnull@localhost> | 2004-03-21 14:06:38 +0000 |
commit | b3f61791f1e9095ce8ae9c6d6415b4ee94e2f7eb (patch) | |
tree | deaa85427ae73a045ddef39395bd286f9d3dc2ff /src/libmp/port/mpfactorial.c | |
parent | 498bb22174aa2c76493e8c67b92949271131ebfb (diff) | |
download | plan9port-b3f61791f1e9095ce8ae9c6d6415b4ee94e2f7eb.tar.gz plan9port-b3f61791f1e9095ce8ae9c6d6415b4ee94e2f7eb.tar.bz2 plan9port-b3f61791f1e9095ce8ae9c6d6415b4ee94e2f7eb.zip |
Add libmp.
Diffstat (limited to 'src/libmp/port/mpfactorial.c')
-rw-r--r-- | src/libmp/port/mpfactorial.c | 75 |
1 files changed, 75 insertions, 0 deletions
diff --git a/src/libmp/port/mpfactorial.c b/src/libmp/port/mpfactorial.c new file mode 100644 index 00000000..01079c44 --- /dev/null +++ b/src/libmp/port/mpfactorial.c @@ -0,0 +1,75 @@ +#include "os.h" +#include <mp.h> +#include <libsec.h> +#include "dat.h" + +mpint* +mpfactorial(ulong n) +{ + int i; + ulong k; + unsigned cnt; + int max, mmax; + mpdigit p, pp[2]; + mpint *r, *s, *stk[31]; + + cnt = 0; + max = mmax = -1; + p = 1; + r = mpnew(0); + for(k=2; k<=n; k++){ + pp[0] = 0; + pp[1] = 0; + mpvecdigmuladd(&p, 1, (mpdigit)k, pp); + if(pp[1] == 0) /* !overflow */ + p = pp[0]; + else{ + cnt++; + if((cnt & 1) == 0){ + s = stk[max]; + mpbits(r, Dbits*(s->top+1+1)); + memset(r->p, 0, Dbytes*(s->top+1+1)); + mpvecmul(s->p, s->top, &p, 1, r->p); + r->sign = 1; + r->top = s->top+1+1; /* XXX: norm */ + mpassign(r, s); + for(i=4; (cnt & (i-1)) == 0; i=i<<1){ + mpmul(stk[max], stk[max-1], r); + mpassign(r, stk[max-1]); + max--; + } + }else{ + max++; + if(max > mmax){ + mmax++; + if(max > nelem(stk)) + abort(); + stk[max] = mpnew(Dbits); + } + stk[max]->top = 1; + stk[max]->p[0] = p; + } + p = (mpdigit)k; + } + } + if(max < 0){ + mpbits(r, Dbits); + r->top = 1; + r->sign = 1; + r->p[0] = p; + }else{ + s = stk[max--]; + mpbits(r, Dbits*(s->top+1+1)); + memset(r->p, 0, Dbytes*(s->top+1+1)); + mpvecmul(s->p, s->top, &p, 1, r->p); + r->sign = 1; + r->top = s->top+1+1; /* XXX: norm */ + } + + while(max >= 0) + mpmul(r, stk[max--], r); + for(max=mmax; max>=0; max--) + mpfree(stk[max]); + mpnorm(r); + return r; +} |