mpexp.c (1283B)
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 }