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@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
23int 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;
239LBL_RES:mp_clear (&res);
240LBL_MU:mp_clear (&mu);
241LBL_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