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
mpexp(mpint * b,mpint * e,mpint * m,mpint * res)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 }
95