xref: /plan9/sys/src/cmd/unix/drawterm/libmp/mpexp.c (revision 8ccd4a6360d974db7bd7bbd4f37e7018419ea908)
1*8ccd4a63SDavid du Colombier #include "os.h"
2*8ccd4a63SDavid du Colombier #include <mp.h>
3*8ccd4a63SDavid du Colombier #include "dat.h"
4*8ccd4a63SDavid du Colombier 
5*8ccd4a63SDavid du Colombier // res = b**e
6*8ccd4a63SDavid du Colombier //
7*8ccd4a63SDavid du Colombier // knuth, vol 2, pp 398-400
8*8ccd4a63SDavid du Colombier 
9*8ccd4a63SDavid du Colombier enum {
10*8ccd4a63SDavid du Colombier 	Freeb=	0x1,
11*8ccd4a63SDavid du Colombier 	Freee=	0x2,
12*8ccd4a63SDavid du Colombier 	Freem=	0x4,
13*8ccd4a63SDavid du Colombier };
14*8ccd4a63SDavid du Colombier 
15*8ccd4a63SDavid du Colombier //int expdebug;
16*8ccd4a63SDavid du Colombier 
17*8ccd4a63SDavid du Colombier void
mpexp(mpint * b,mpint * e,mpint * m,mpint * res)18*8ccd4a63SDavid du Colombier mpexp(mpint *b, mpint *e, mpint *m, mpint *res)
19*8ccd4a63SDavid du Colombier {
20*8ccd4a63SDavid du Colombier 	mpint *t[2];
21*8ccd4a63SDavid du Colombier 	int tofree;
22*8ccd4a63SDavid du Colombier 	mpdigit d, bit;
23*8ccd4a63SDavid du Colombier 	int i, j;
24*8ccd4a63SDavid du Colombier 
25*8ccd4a63SDavid du Colombier 	t[0] = mpcopy(b);
26*8ccd4a63SDavid du Colombier 	t[1] = res;
27*8ccd4a63SDavid du Colombier 
28*8ccd4a63SDavid du Colombier 	tofree = 0;
29*8ccd4a63SDavid du Colombier 	if(res == b){
30*8ccd4a63SDavid du Colombier 		b = mpcopy(b);
31*8ccd4a63SDavid du Colombier 		tofree |= Freeb;
32*8ccd4a63SDavid du Colombier 	}
33*8ccd4a63SDavid du Colombier 	if(res == e){
34*8ccd4a63SDavid du Colombier 		e = mpcopy(e);
35*8ccd4a63SDavid du Colombier 		tofree |= Freee;
36*8ccd4a63SDavid du Colombier 	}
37*8ccd4a63SDavid du Colombier 	if(res == m){
38*8ccd4a63SDavid du Colombier 		m = mpcopy(m);
39*8ccd4a63SDavid du Colombier 		tofree |= Freem;
40*8ccd4a63SDavid du Colombier 	}
41*8ccd4a63SDavid du Colombier 
42*8ccd4a63SDavid du Colombier 	// skip first bit
43*8ccd4a63SDavid du Colombier 	i = e->top-1;
44*8ccd4a63SDavid du Colombier 	d = e->p[i];
45*8ccd4a63SDavid du Colombier 	for(bit = mpdighi; (bit & d) == 0; bit >>= 1)
46*8ccd4a63SDavid du Colombier 		;
47*8ccd4a63SDavid du Colombier 	bit >>= 1;
48*8ccd4a63SDavid du Colombier 
49*8ccd4a63SDavid du Colombier 	j = 0;
50*8ccd4a63SDavid du Colombier 	for(;;){
51*8ccd4a63SDavid du Colombier 		for(; bit != 0; bit >>= 1){
52*8ccd4a63SDavid du Colombier 			mpmul(t[j], t[j], t[j^1]);
53*8ccd4a63SDavid du Colombier 			if(bit & d)
54*8ccd4a63SDavid du Colombier 				mpmul(t[j^1], b, t[j]);
55*8ccd4a63SDavid du Colombier 			else
56*8ccd4a63SDavid du Colombier 				j ^= 1;
57*8ccd4a63SDavid du Colombier 			if(m != nil && t[j]->top > m->top){
58*8ccd4a63SDavid du Colombier 				mpmod(t[j], m, t[j^1]);
59*8ccd4a63SDavid du Colombier 				j ^= 1;
60*8ccd4a63SDavid du Colombier 			}
61*8ccd4a63SDavid du Colombier 		}
62*8ccd4a63SDavid du Colombier 		if(--i < 0)
63*8ccd4a63SDavid du Colombier 			break;
64*8ccd4a63SDavid du Colombier 		bit = mpdighi;
65*8ccd4a63SDavid du Colombier 		d = e->p[i];
66*8ccd4a63SDavid du Colombier 	}
67*8ccd4a63SDavid du Colombier 	if(m != nil){
68*8ccd4a63SDavid du Colombier 		mpmod(t[j], m, t[j^1]);
69*8ccd4a63SDavid du Colombier 		j ^= 1;
70*8ccd4a63SDavid du Colombier 	}
71*8ccd4a63SDavid du Colombier 	if(t[j] == res){
72*8ccd4a63SDavid du Colombier 		mpfree(t[j^1]);
73*8ccd4a63SDavid du Colombier 	} else {
74*8ccd4a63SDavid du Colombier 		mpassign(t[j], res);
75*8ccd4a63SDavid du Colombier 		mpfree(t[j]);
76*8ccd4a63SDavid du Colombier 	}
77*8ccd4a63SDavid du Colombier 
78*8ccd4a63SDavid du Colombier 	if(tofree){
79*8ccd4a63SDavid du Colombier 		if(tofree & Freeb)
80*8ccd4a63SDavid du Colombier 			mpfree(b);
81*8ccd4a63SDavid du Colombier 		if(tofree & Freee)
82*8ccd4a63SDavid du Colombier 			mpfree(e);
83*8ccd4a63SDavid du Colombier 		if(tofree & Freem)
84*8ccd4a63SDavid du Colombier 			mpfree(m);
85*8ccd4a63SDavid du Colombier 	}
86*8ccd4a63SDavid du Colombier }
87