Blob
1 #include "os.h"2 #include <mp.h>3 #include "dat.h"5 // res = b**e6 //7 // knuth, vol 2, pp 398-4009 enum {10 Freeb= 0x1,11 Freee= 0x2,12 Freem= 0x4,13 };15 //int expdebug;17 void18 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;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");33 t[0] = mpcopy(b);34 t[1] = res;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 }50 // skip first bit51 i = e->top-1;52 d = e->p[i];53 for(bit = mpdighi; (bit & d) == 0; bit >>= 1)54 ;55 bit >>= 1;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 else64 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 }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 }