| rsc | b3f6179 | 2004-03-21 14:06:38 +0000 | [diff] [blame] | 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 | } |
