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