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 t[0] = mpcopy(b);
26 t[1] = res;
27
28 tofree = 0;
29 if(res == b){
30 b = mpcopy(b);
31 tofree |= Freeb;
32 }
33 if(res == e){
34 e = mpcopy(e);
35 tofree |= Freee;
36 }
37 if(res == m){
38 m = mpcopy(m);
39 tofree |= Freem;
40 }
41
42 // skip first bit
43 i = e->top-1;
44 d = e->p[i];
45 for(bit = mpdighi; (bit & d) == 0; bit >>= 1)
46 ;
47 bit >>= 1;
48
49 j = 0;
50 for(;;){
51 for(; bit != 0; bit >>= 1){
52 mpmul(t[j], t[j], t[j^1]);
53 if(bit & d)
54 mpmul(t[j^1], b, t[j]);
55 else
56 j ^= 1;
57 if(m != nil && t[j]->top > m->top){
58 mpmod(t[j], m, t[j^1]);
59 j ^= 1;
60 }
61 }
62 if(--i < 0)
63 break;
64 bit = mpdighi;
65 d = e->p[i];
66 }
67 if(m != nil){
68 mpmod(t[j], m, t[j^1]);
69 j ^= 1;
70 }
71 if(t[j] == res){
72 mpfree(t[j^1]);
73 } else {
74 mpassign(t[j], res);
75 mpfree(t[j]);
76 }
77
78 if(tofree){
79 if(tofree & Freeb)
80 mpfree(b);
81 if(tofree & Freee)
82 mpfree(e);
83 if(tofree & Freem)
84 mpfree(m);
85 }
86 }
87