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