tmpexp.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tmpexp.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 }