1 #include <tommath.h> 2 #ifdef BN_S_MP_EXPTMOD_C 3 /* LibTomMath, multiple-precision integer library -- Tom St Denis 4 * 5 * LibTomMath is a library that provides multiple-precision 6 * integer arithmetic as well as number theoretic functionality. 7 * 8 * The library was designed directly after the MPI library by 9 * Michael Fromberger but has been written from scratch with 10 * additional optimizations in place. 11 * 12 * The library is free for all purposes without any express 13 * guarantee it works. 14 * 15 * Tom St Denis, tomstdenis (at) gmail.com, http://math.libtomcrypt.com 16 */ 17 #ifdef MP_LOW_MEM 18 #define TAB_SIZE 32 19 #else 20 #define TAB_SIZE 256 21 #endif 22 23 int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode) 24 { 25 mp_int M[TAB_SIZE], res, mu; 26 mp_digit buf; 27 int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize; 28 int (*redux)(mp_int*,mp_int*,mp_int*); 29 30 /* find window size */ 31 x = mp_count_bits (X); 32 if (x <= 7) { 33 winsize = 2; 34 } else if (x <= 36) { 35 winsize = 3; 36 } else if (x <= 140) { 37 winsize = 4; 38 } else if (x <= 450) { 39 winsize = 5; 40 } else if (x <= 1303) { 41 winsize = 6; 42 } else if (x <= 3529) { 43 winsize = 7; 44 } else { 45 winsize = 8; 46 } 47 48 #ifdef MP_LOW_MEM 49 if (winsize > 5) { 50 winsize = 5; 51 } 52 #endif 53 54 /* init M array */ 55 /* init first cell */ 56 if ((err = mp_init(&M[1])) != MP_OKAY) { 57 return err; 58 } 59 60 /* now init the second half of the array */ 61 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 62 if ((err = mp_init(&M[x])) != MP_OKAY) { 63 for (y = 1<<(winsize-1); y < x; y++) { 64 mp_clear (&M[y]); 65 } 66 mp_clear(&M[1]); 67 return err; 68 } 69 } 70 71 /* create mu, used for Barrett reduction */ 72 if ((err = mp_init (&mu)) != MP_OKAY) { 73 goto LBL_M; 74 } 75 76 if (redmode == 0) { 77 if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) { 78 goto LBL_MU; 79 } 80 redux = mp_reduce; 81 } else { 82 if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) { 83 goto LBL_MU; 84 } 85 redux = mp_reduce_2k_l; 86 } 87 88 /* create M table 89 * 90 * The M table contains powers of the base, 91 * e.g. M[x] = G**x mod P 92 * 93 * The first half of the table is not 94 * computed though accept for M[0] and M[1] 95 */ 96 if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) { 97 goto LBL_MU; 98 } 99 100 /* compute the value at M[1<<(winsize-1)] by squaring 101 * M[1] (winsize-1) times 102 */ 103 if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) { 104 goto LBL_MU; 105 } 106 107 for (x = 0; x < (winsize - 1); x++) { 108 /* square it */ 109 if ((err = mp_sqr (&M[1 << (winsize - 1)], 110 &M[1 << (winsize - 1)])) != MP_OKAY) { 111 goto LBL_MU; 112 } 113 114 /* reduce modulo P */ 115 if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) { 116 goto LBL_MU; 117 } 118 } 119 120 /* create upper table, that is M[x] = M[x-1] * M[1] (mod P) 121 * for x = (2**(winsize - 1) + 1) to (2**winsize - 1) 122 */ 123 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) { 124 if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) { 125 goto LBL_MU; 126 } 127 if ((err = redux (&M[x], P, &mu)) != MP_OKAY) { 128 goto LBL_MU; 129 } 130 } 131 132 /* setup result */ 133 if ((err = mp_init (&res)) != MP_OKAY) { 134 goto LBL_MU; 135 } 136 mp_set (&res, 1); 137 138 /* set initial mode and bit cnt */ 139 mode = 0; 140 bitcnt = 1; 141 buf = 0; 142 digidx = X->used - 1; 143 bitcpy = 0; 144 bitbuf = 0; 145 146 for (;;) { 147 /* grab next digit as required */ 148 if (--bitcnt == 0) { 149 /* if digidx == -1 we are out of digits */ 150 if (digidx == -1) { 151 break; 152 } 153 /* read next digit and reset the bitcnt */ 154 buf = X->dp[digidx--]; 155 bitcnt = (int) DIGIT_BIT; 156 } 157 158 /* grab the next msb from the exponent */ 159 y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1; 160 buf <<= (mp_digit)1; 161 162 /* if the bit is zero and mode == 0 then we ignore it 163 * These represent the leading zero bits before the first 1 bit 164 * in the exponent. Technically this opt is not required but it 165 * does lower the # of trivial squaring/reductions used 166 */ 167 if (mode == 0 && y == 0) { 168 continue; 169 } 170 171 /* if the bit is zero and mode == 1 then we square */ 172 if (mode == 1 && y == 0) { 173 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 174 goto LBL_RES; 175 } 176 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 177 goto LBL_RES; 178 } 179 continue; 180 } 181 182 /* else we add it to the window */ 183 bitbuf |= (y << (winsize - ++bitcpy)); 184 mode = 2; 185 186 if (bitcpy == winsize) { 187 /* ok window is filled so square as required and multiply */ 188 /* square first */ 189 for (x = 0; x < winsize; x++) { 190 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 191 goto LBL_RES; 192 } 193 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 194 goto LBL_RES; 195 } 196 } 197 198 /* then multiply */ 199 if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) { 200 goto LBL_RES; 201 } 202 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 203 goto LBL_RES; 204 } 205 206 /* empty window and reset */ 207 bitcpy = 0; 208 bitbuf = 0; 209 mode = 1; 210 } 211 } 212 213 /* if bits remain then square/multiply */ 214 if (mode == 2 && bitcpy > 0) { 215 /* square then multiply if the bit is set */ 216 for (x = 0; x < bitcpy; x++) { 217 if ((err = mp_sqr (&res, &res)) != MP_OKAY) { 218 goto LBL_RES; 219 } 220 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 221 goto LBL_RES; 222 } 223 224 bitbuf <<= 1; 225 if ((bitbuf & (1 << winsize)) != 0) { 226 /* then multiply */ 227 if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) { 228 goto LBL_RES; 229 } 230 if ((err = redux (&res, P, &mu)) != MP_OKAY) { 231 goto LBL_RES; 232 } 233 } 234 } 235 } 236 237 mp_exch (&res, Y); 238 err = MP_OKAY; 239 LBL_RES:mp_clear (&res); 240 LBL_MU:mp_clear (&mu); 241 LBL_M: 242 mp_clear(&M[1]); 243 for (x = 1<<(winsize-1); x < (1 << winsize); x++) { 244 mp_clear (&M[x]); 245 } 246 return err; 247 } 248 #endif 249 250 /* $Source: /cvs/libtom/libtommath/bn_s_mp_exptmod.c,v $ */ 251 /* $Revision: 1.4 $ */ 252 /* $Date: 2006/03/31 14:18:44 $ */ 253
