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