← Back to branch summary

~ubuntu-branches/ubuntu/vivid/wpasupplicant/vivid

« back to all changes in this revision

Viewing changes to src/tls/libtommath.c

  • Committer: Bazaar Package Importer
  • Author(s): Kel Modderman
  • Date: 2008-03-12 20:03:04 UTC
  • mfrom: (1.1.10 upstream)
  • mto: This revision was merged to the branch mainline in revision 4.
  • Revision ID: james.westby@ubuntu.com-20080312200304-4331y9wj46pdd34z
Tags: 0.6.3-1
http://bugs.debian.org/466910
* New upstream release.
* Drop patches applied upstream:
  - debian/patches/30_wpa_gui_qt4_eventhistoryui_rework.patch
  - debian/patches/31_wpa_gui_qt4_eventhistory_always_scrollbar.patch
  - debian/patches/32_wpa_gui_qt4_eventhistory_scroll_with_events.patch
  - debian/patches/40_dbus_ssid_data.patch
* Tidy up the clean target of debian/rules. Now that the madwifi headers are
  handled differently we no longer need to do any cleanup.
* Fix formatting error in debian/ifupdown/wpa_action.8 to make lintian
  quieter.
* Add patch to fix formatting errors in manpages build from sgml source. Use
  <emphasis> tags to hightlight keywords instead of surrounding them in
  strong quotes.
  - debian/patches/41_manpage_format_fixes.patch
* wpasupplicant binary package no longer suggests pcscd, guessnet, iproute
  or wireless-tools, nor does it recommend dhcp3-client. These are not
  needed.
* Add debian/patches/10_silence_siocsiwauth_icotl_failure.patch to disable
  ioctl failure messages that occur under normal conditions.
* Cherry pick two upstream git commits concerning the dbus interface:
  - debian/patches/11_avoid_dbus_version_namespace.patch
  - debian/patches/12_fix_potential_use_after_free.patch
* Add debian/patches/42_manpage_explain_available_drivers.patch to explain
  that not all of the driver backends are available in the provided
  wpa_supplicant binary, and that the canonical list of supported driver
  backends can be retrieved from the wpa_supplicant -h (help) output.
  (Closes: #466910)
* Add debian/patches/20_wpa_gui_qt4_disable_link_prl.patch to remove
  link_prl CONFIG compile flag added by qmake-qt4 >= 4.3.4-2 to avoid excess
  linking.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/tls/libtommath.c
 
1
/*
 
2
 * Minimal code for RSA support from LibTomMath 0.3.9
 
3
 * http://math.libtomcrypt.com/
 
4
 * http://math.libtomcrypt.com/files/ltm-0.39.tar.bz2
 
5
 * This library was released in public domain by Tom St Denis.
 
6
 *
 
7
 * The combination in this file is not using many of the optimized algorithms
 
8
 * (e.g., Montgomery reduction) and is considerable slower than the LibTomMath
 
9
 * with its default of SC_RSA_1 settins. The main purpose of having this
 
10
 * version here is to make it easier to build bignum.c wrapper without having
 
11
 * to install and build an external library. However, it is likely worth the
 
12
 * effort to use the full library with SC_RSA_1 instead of this minimized copy.
 
13
 * Including the optimized algorithms may increase the size requirements by
 
14
 * 15 kB or so (measured with x86 build).
 
15
 *
 
16
 * If CONFIG_INTERNAL_LIBTOMMATH is defined, bignum.c includes this
 
17
 * libtommath.c file instead of using the external LibTomMath library.
 
18
 */
 
19
 
 
20
#ifndef CHAR_BIT
 
21
#define CHAR_BIT 8
 
22
#endif
 
23
 
 
24
#define BN_MP_INVMOD_C
 
25
#define BN_S_MP_EXPTMOD_C /* Note: #undef in tommath_superclass.h; this would
 
26
                           * require BN_MP_EXPTMOD_FAST_C instead */
 
27
#define BN_S_MP_MUL_DIGS_C
 
28
#define BN_MP_INVMOD_SLOW_C
 
29
#define BN_S_MP_SQR_C
 
30
#define BN_S_MP_MUL_HIGH_DIGS_C /* Note: #undef in tommath_superclass.h; this
 
31
                                 * would require other than mp_reduce */
 
32
 
 
33
 
 
34
/* from tommath.h */
 
35
 
 
36
#ifndef MIN
 
37
   #define MIN(x,y) ((x)<(y)?(x):(y))
 
38
#endif
 
39
 
 
40
#ifndef MAX
 
41
   #define MAX(x,y) ((x)>(y)?(x):(y))
 
42
#endif
 
43
 
 
44
#define  OPT_CAST(x)
 
45
 
 
46
typedef unsigned long mp_digit;
 
47
typedef u64 mp_word;
 
48
 
 
49
#define DIGIT_BIT          28
 
50
#define MP_28BIT
 
51
 
 
52
 
 
53
#define XMALLOC  os_malloc
 
54
#define XFREE    os_free
 
55
#define XREALLOC os_realloc
 
56
 
 
57
 
 
58
#define MP_MASK          ((((mp_digit)1)<<((mp_digit)DIGIT_BIT))-((mp_digit)1))
 
59
 
 
60
#define MP_LT        -1   /* less than */
 
61
#define MP_EQ         0   /* equal to */
 
62
#define MP_GT         1   /* greater than */
 
63
 
 
64
#define MP_ZPOS       0   /* positive integer */
 
65
#define MP_NEG        1   /* negative */
 
66
 
 
67
#define MP_OKAY       0   /* ok result */
 
68
#define MP_MEM        -2  /* out of mem */
 
69
#define MP_VAL        -3  /* invalid input */
 
70
 
 
71
#define MP_YES        1   /* yes response */
 
72
#define MP_NO         0   /* no response */
 
73
 
 
74
typedef int           mp_err;
 
75
 
 
76
/* define this to use lower memory usage routines (exptmods mostly) */
 
77
#define MP_LOW_MEM
 
78
 
 
79
/* default precision */
 
80
#ifndef MP_PREC
 
81
   #ifndef MP_LOW_MEM
 
82
      #define MP_PREC                 32     /* default digits of precision */
 
83
   #else
 
84
      #define MP_PREC                 8      /* default digits of precision */
 
85
   #endif   
 
86
#endif
 
87
 
 
88
/* size of comba arrays, should be at least 2 * 2**(BITS_PER_WORD - BITS_PER_DIGIT*2) */
 
89
#define MP_WARRAY               (1 << (sizeof(mp_word) * CHAR_BIT - 2 * DIGIT_BIT + 1))
 
90
 
 
91
/* the infamous mp_int structure */
 
92
typedef struct  {
 
93
    int used, alloc, sign;
 
94
    mp_digit *dp;
 
95
} mp_int;
 
96
 
 
97
 
 
98
/* ---> Basic Manipulations <--- */
 
99
#define mp_iszero(a) (((a)->used == 0) ? MP_YES : MP_NO)
 
100
#define mp_iseven(a) (((a)->used > 0 && (((a)->dp[0] & 1) == 0)) ? MP_YES : MP_NO)
 
101
#define mp_isodd(a)  (((a)->used > 0 && (((a)->dp[0] & 1) == 1)) ? MP_YES : MP_NO)
 
102
 
 
103
 
 
104
/* prototypes for copied functions */
 
105
#define s_mp_mul(a, b, c) s_mp_mul_digs(a, b, c, (a)->used + (b)->used + 1)
 
106
static int s_mp_exptmod(mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode);
 
107
static int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs);
 
108
static int s_mp_sqr(mp_int * a, mp_int * b);
 
109
static int s_mp_mul_high_digs(mp_int * a, mp_int * b, mp_int * c, int digs);
 
110
 
 
111
static int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs);
 
112
 
 
113
static int mp_init_multi(mp_int *mp, ...);
 
114
static void mp_clear_multi(mp_int *mp, ...);
 
115
static int mp_lshd(mp_int * a, int b);
 
116
static void mp_set(mp_int * a, mp_digit b);
 
117
static void mp_clamp(mp_int * a);
 
118
static void mp_exch(mp_int * a, mp_int * b);
 
119
static void mp_rshd(mp_int * a, int b);
 
120
static void mp_zero(mp_int * a);
 
121
static int mp_mod_2d(mp_int * a, int b, mp_int * c);
 
122
static int mp_div_2d(mp_int * a, int b, mp_int * c, mp_int * d);
 
123
static int mp_init_copy(mp_int * a, mp_int * b);
 
124
static int mp_mul_2d(mp_int * a, int b, mp_int * c);
 
125
static int mp_div_2(mp_int * a, mp_int * b);
 
126
static int mp_copy(mp_int * a, mp_int * b);
 
127
static int mp_count_bits(mp_int * a);
 
128
static int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d);
 
129
static int mp_mod(mp_int * a, mp_int * b, mp_int * c);
 
130
static int mp_grow(mp_int * a, int size);
 
131
static int mp_cmp_mag(mp_int * a, mp_int * b);
 
132
static int mp_invmod(mp_int * a, mp_int * b, mp_int * c);
 
133
static int mp_abs(mp_int * a, mp_int * b);
 
134
static int mp_invmod_slow(mp_int * a, mp_int * b, mp_int * c);
 
135
static int mp_sqr(mp_int * a, mp_int * b);
 
136
static int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d);
 
137
static int mp_reduce_2k_setup_l(mp_int *a, mp_int *d);
 
138
static int mp_2expt(mp_int * a, int b);
 
139
static int mp_reduce_setup(mp_int * a, mp_int * b);
 
140
static int mp_reduce(mp_int * x, mp_int * m, mp_int * mu);
 
141
static int mp_init_size(mp_int * a, int size);
 
142
 
 
143
 
 
144
 
 
145
/* functions from bn_<func name>.c */
 
146
 
 
147
 
 
148
/* reverse an array, used for radix code */
 
149
static void bn_reverse (unsigned char *s, int len)
 
150
{
 
151
  int     ix, iy;
 
152
  unsigned char t;
 
153
 
 
154
  ix = 0;
 
155
  iy = len - 1;
 
156
  while (ix < iy) {
 
157
    t     = s[ix];
 
158
    s[ix] = s[iy];
 
159
    s[iy] = t;
 
160
    ++ix;
 
161
    --iy;
 
162
  }
 
163
}
 
164
 
 
165
 
 
166
/* low level addition, based on HAC pp.594, Algorithm 14.7 */
 
167
static int s_mp_add (mp_int * a, mp_int * b, mp_int * c)
 
168
{
 
169
  mp_int *x;
 
170
  int     olduse, res, min, max;
 
171
 
 
172
  /* find sizes, we let |a| <= |b| which means we have to sort
 
173
   * them.  "x" will point to the input with the most digits
 
174
   */
 
175
  if (a->used > b->used) {
 
176
    min = b->used;
 
177
    max = a->used;
 
178
    x = a;
 
179
  } else {
 
180
    min = a->used;
 
181
    max = b->used;
 
182
    x = b;
 
183
  }
 
184
 
 
185
  /* init result */
 
186
  if (c->alloc < max + 1) {
 
187
    if ((res = mp_grow (c, max + 1)) != MP_OKAY) {
 
188
      return res;
 
189
    }
 
190
  }
 
191
 
 
192
  /* get old used digit count and set new one */
 
193
  olduse = c->used;
 
194
  c->used = max + 1;
 
195
 
 
196
  {
 
197
    register mp_digit u, *tmpa, *tmpb, *tmpc;
 
198
    register int i;
 
199
 
 
200
    /* alias for digit pointers */
 
201
 
 
202
    /* first input */
 
203
    tmpa = a->dp;
 
204
 
 
205
    /* second input */
 
206
    tmpb = b->dp;
 
207
 
 
208
    /* destination */
 
209
    tmpc = c->dp;
 
210
 
 
211
    /* zero the carry */
 
212
    u = 0;
 
213
    for (i = 0; i < min; i++) {
 
214
      /* Compute the sum at one digit, T[i] = A[i] + B[i] + U */
 
215
      *tmpc = *tmpa++ + *tmpb++ + u;
 
216
 
 
217
      /* U = carry bit of T[i] */
 
218
      u = *tmpc >> ((mp_digit)DIGIT_BIT);
 
219
 
 
220
      /* take away carry bit from T[i] */
 
221
      *tmpc++ &= MP_MASK;
 
222
    }
 
223
 
 
224
    /* now copy higher words if any, that is in A+B 
 
225
     * if A or B has more digits add those in 
 
226
     */
 
227
    if (min != max) {
 
228
      for (; i < max; i++) {
 
229
        /* T[i] = X[i] + U */
 
230
        *tmpc = x->dp[i] + u;
 
231
 
 
232
        /* U = carry bit of T[i] */
 
233
        u = *tmpc >> ((mp_digit)DIGIT_BIT);
 
234
 
 
235
        /* take away carry bit from T[i] */
 
236
        *tmpc++ &= MP_MASK;
 
237
      }
 
238
    }
 
239
 
 
240
    /* add carry */
 
241
    *tmpc++ = u;
 
242
 
 
243
    /* clear digits above oldused */
 
244
    for (i = c->used; i < olduse; i++) {
 
245
      *tmpc++ = 0;
 
246
    }
 
247
  }
 
248
 
 
249
  mp_clamp (c);
 
250
  return MP_OKAY;
 
251
}
 
252
 
 
253
 
 
254
/* low level subtraction (assumes |a| > |b|), HAC pp.595 Algorithm 14.9 */
 
255
static int s_mp_sub (mp_int * a, mp_int * b, mp_int * c)
 
256
{
 
257
  int     olduse, res, min, max;
 
258
 
 
259
  /* find sizes */
 
260
  min = b->used;
 
261
  max = a->used;
 
262
 
 
263
  /* init result */
 
264
  if (c->alloc < max) {
 
265
    if ((res = mp_grow (c, max)) != MP_OKAY) {
 
266
      return res;
 
267
    }
 
268
  }
 
269
  olduse = c->used;
 
270
  c->used = max;
 
271
 
 
272
  {
 
273
    register mp_digit u, *tmpa, *tmpb, *tmpc;
 
274
    register int i;
 
275
 
 
276
    /* alias for digit pointers */
 
277
    tmpa = a->dp;
 
278
    tmpb = b->dp;
 
279
    tmpc = c->dp;
 
280
 
 
281
    /* set carry to zero */
 
282
    u = 0;
 
283
    for (i = 0; i < min; i++) {
 
284
      /* T[i] = A[i] - B[i] - U */
 
285
      *tmpc = *tmpa++ - *tmpb++ - u;
 
286
 
 
287
      /* U = carry bit of T[i]
 
288
       * Note this saves performing an AND operation since
 
289
       * if a carry does occur it will propagate all the way to the
 
290
       * MSB.  As a result a single shift is enough to get the carry
 
291
       */
 
292
      u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
 
293
 
 
294
      /* Clear carry from T[i] */
 
295
      *tmpc++ &= MP_MASK;
 
296
    }
 
297
 
 
298
    /* now copy higher words if any, e.g. if A has more digits than B  */
 
299
    for (; i < max; i++) {
 
300
      /* T[i] = A[i] - U */
 
301
      *tmpc = *tmpa++ - u;
 
302
 
 
303
      /* U = carry bit of T[i] */
 
304
      u = *tmpc >> ((mp_digit)(CHAR_BIT * sizeof (mp_digit) - 1));
 
305
 
 
306
      /* Clear carry from T[i] */
 
307
      *tmpc++ &= MP_MASK;
 
308
    }
 
309
 
 
310
    /* clear digits above used (since we may not have grown result above) */
 
311
    for (i = c->used; i < olduse; i++) {
 
312
      *tmpc++ = 0;
 
313
    }
 
314
  }
 
315
 
 
316
  mp_clamp (c);
 
317
  return MP_OKAY;
 
318
}
 
319
 
 
320
 
 
321
/* init a new mp_int */
 
322
static int mp_init (mp_int * a)
 
323
{
 
324
  int i;
 
325
 
 
326
  /* allocate memory required and clear it */
 
327
  a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * MP_PREC);
 
328
  if (a->dp == NULL) {
 
329
    return MP_MEM;
 
330
  }
 
331
 
 
332
  /* set the digits to zero */
 
333
  for (i = 0; i < MP_PREC; i++) {
 
334
      a->dp[i] = 0;
 
335
  }
 
336
 
 
337
  /* set the used to zero, allocated digits to the default precision
 
338
   * and sign to positive */
 
339
  a->used  = 0;
 
340
  a->alloc = MP_PREC;
 
341
  a->sign  = MP_ZPOS;
 
342
 
 
343
  return MP_OKAY;
 
344
}
 
345
 
 
346
 
 
347
/* clear one (frees)  */
 
348
static void mp_clear (mp_int * a)
 
349
{
 
350
  int i;
 
351
 
 
352
  /* only do anything if a hasn't been freed previously */
 
353
  if (a->dp != NULL) {
 
354
    /* first zero the digits */
 
355
    for (i = 0; i < a->used; i++) {
 
356
        a->dp[i] = 0;
 
357
    }
 
358
 
 
359
    /* free ram */
 
360
    XFREE(a->dp);
 
361
 
 
362
    /* reset members to make debugging easier */
 
363
    a->dp    = NULL;
 
364
    a->alloc = a->used = 0;
 
365
    a->sign  = MP_ZPOS;
 
366
  }
 
367
}
 
368
 
 
369
 
 
370
/* high level addition (handles signs) */
 
371
static int mp_add (mp_int * a, mp_int * b, mp_int * c)
 
372
{
 
373
  int     sa, sb, res;
 
374
 
 
375
  /* get sign of both inputs */
 
376
  sa = a->sign;
 
377
  sb = b->sign;
 
378
 
 
379
  /* handle two cases, not four */
 
380
  if (sa == sb) {
 
381
    /* both positive or both negative */
 
382
    /* add their magnitudes, copy the sign */
 
383
    c->sign = sa;
 
384
    res = s_mp_add (a, b, c);
 
385
  } else {
 
386
    /* one positive, the other negative */
 
387
    /* subtract the one with the greater magnitude from */
 
388
    /* the one of the lesser magnitude.  The result gets */
 
389
    /* the sign of the one with the greater magnitude. */
 
390
    if (mp_cmp_mag (a, b) == MP_LT) {
 
391
      c->sign = sb;
 
392
      res = s_mp_sub (b, a, c);
 
393
    } else {
 
394
      c->sign = sa;
 
395
      res = s_mp_sub (a, b, c);
 
396
    }
 
397
  }
 
398
  return res;
 
399
}
 
400
 
 
401
 
 
402
/* high level subtraction (handles signs) */
 
403
static int mp_sub (mp_int * a, mp_int * b, mp_int * c)
 
404
{
 
405
  int     sa, sb, res;
 
406
 
 
407
  sa = a->sign;
 
408
  sb = b->sign;
 
409
 
 
410
  if (sa != sb) {
 
411
    /* subtract a negative from a positive, OR */
 
412
    /* subtract a positive from a negative. */
 
413
    /* In either case, ADD their magnitudes, */
 
414
    /* and use the sign of the first number. */
 
415
    c->sign = sa;
 
416
    res = s_mp_add (a, b, c);
 
417
  } else {
 
418
    /* subtract a positive from a positive, OR */
 
419
    /* subtract a negative from a negative. */
 
420
    /* First, take the difference between their */
 
421
    /* magnitudes, then... */
 
422
    if (mp_cmp_mag (a, b) != MP_LT) {
 
423
      /* Copy the sign from the first */
 
424
      c->sign = sa;
 
425
      /* The first has a larger or equal magnitude */
 
426
      res = s_mp_sub (a, b, c);
 
427
    } else {
 
428
      /* The result has the *opposite* sign from */
 
429
      /* the first number. */
 
430
      c->sign = (sa == MP_ZPOS) ? MP_NEG : MP_ZPOS;
 
431
      /* The second has a larger magnitude */
 
432
      res = s_mp_sub (b, a, c);
 
433
    }
 
434
  }
 
435
  return res;
 
436
}
 
437
 
 
438
 
 
439
/* high level multiplication (handles sign) */
 
440
static int mp_mul (mp_int * a, mp_int * b, mp_int * c)
 
441
{
 
442
  int     res, neg;
 
443
  neg = (a->sign == b->sign) ? MP_ZPOS : MP_NEG;
 
444
 
 
445
  /* use Toom-Cook? */
 
446
#ifdef BN_MP_TOOM_MUL_C
 
447
  if (MIN (a->used, b->used) >= TOOM_MUL_CUTOFF) {
 
448
    res = mp_toom_mul(a, b, c);
 
449
  } else 
 
450
#endif
 
451
#ifdef BN_MP_KARATSUBA_MUL_C
 
452
  /* use Karatsuba? */
 
453
  if (MIN (a->used, b->used) >= KARATSUBA_MUL_CUTOFF) {
 
454
    res = mp_karatsuba_mul (a, b, c);
 
455
  } else 
 
456
#endif
 
457
  {
 
458
    /* can we use the fast multiplier?
 
459
     *
 
460
     * The fast multiplier can be used if the output will 
 
461
     * have less than MP_WARRAY digits and the number of 
 
462
     * digits won't affect carry propagation
 
463
     */
 
464
#ifdef BN_FAST_S_MP_MUL_DIGS_C
 
465
    int     digs = a->used + b->used + 1;
 
466
 
 
467
    if ((digs < MP_WARRAY) &&
 
468
        MIN(a->used, b->used) <= 
 
469
        (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
 
470
      res = fast_s_mp_mul_digs (a, b, c, digs);
 
471
    } else 
 
472
#endif
 
473
#ifdef BN_S_MP_MUL_DIGS_C
 
474
      res = s_mp_mul (a, b, c); /* uses s_mp_mul_digs */
 
475
#else
 
476
#error mp_mul could fail
 
477
      res = MP_VAL;
 
478
#endif
 
479
 
 
480
  }
 
481
  c->sign = (c->used > 0) ? neg : MP_ZPOS;
 
482
  return res;
 
483
}
 
484
 
 
485
 
 
486
/* d = a * b (mod c) */
 
487
static int mp_mulmod (mp_int * a, mp_int * b, mp_int * c, mp_int * d)
 
488
{
 
489
  int     res;
 
490
  mp_int  t;
 
491
 
 
492
  if ((res = mp_init (&t)) != MP_OKAY) {
 
493
    return res;
 
494
  }
 
495
 
 
496
  if ((res = mp_mul (a, b, &t)) != MP_OKAY) {
 
497
    mp_clear (&t);
 
498
    return res;
 
499
  }
 
500
  res = mp_mod (&t, c, d);
 
501
  mp_clear (&t);
 
502
  return res;
 
503
}
 
504
 
 
505
 
 
506
/* c = a mod b, 0 <= c < b */
 
507
static int mp_mod (mp_int * a, mp_int * b, mp_int * c)
 
508
{
 
509
  mp_int  t;
 
510
  int     res;
 
511
 
 
512
  if ((res = mp_init (&t)) != MP_OKAY) {
 
513
    return res;
 
514
  }
 
515
 
 
516
  if ((res = mp_div (a, b, NULL, &t)) != MP_OKAY) {
 
517
    mp_clear (&t);
 
518
    return res;
 
519
  }
 
520
 
 
521
  if (t.sign != b->sign) {
 
522
    res = mp_add (b, &t, c);
 
523
  } else {
 
524
    res = MP_OKAY;
 
525
    mp_exch (&t, c);
 
526
  }
 
527
 
 
528
  mp_clear (&t);
 
529
  return res;
 
530
}
 
531
 
 
532
 
 
533
/* this is a shell function that calls either the normal or Montgomery
 
534
 * exptmod functions.  Originally the call to the montgomery code was
 
535
 * embedded in the normal function but that wasted alot of stack space
 
536
 * for nothing (since 99% of the time the Montgomery code would be called)
 
537
 */
 
538
static int mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y)
 
539
{
 
540
  int dr;
 
541
 
 
542
  /* modulus P must be positive */
 
543
  if (P->sign == MP_NEG) {
 
544
     return MP_VAL;
 
545
  }
 
546
 
 
547
  /* if exponent X is negative we have to recurse */
 
548
  if (X->sign == MP_NEG) {
 
549
#ifdef BN_MP_INVMOD_C
 
550
     mp_int tmpG, tmpX;
 
551
     int err;
 
552
 
 
553
     /* first compute 1/G mod P */
 
554
     if ((err = mp_init(&tmpG)) != MP_OKAY) {
 
555
        return err;
 
556
     }
 
557
     if ((err = mp_invmod(G, P, &tmpG)) != MP_OKAY) {
 
558
        mp_clear(&tmpG);
 
559
        return err;
 
560
     }
 
561
 
 
562
     /* now get |X| */
 
563
     if ((err = mp_init(&tmpX)) != MP_OKAY) {
 
564
        mp_clear(&tmpG);
 
565
        return err;
 
566
     }
 
567
     if ((err = mp_abs(X, &tmpX)) != MP_OKAY) {
 
568
        mp_clear_multi(&tmpG, &tmpX, NULL);
 
569
        return err;
 
570
     }
 
571
 
 
572
     /* and now compute (1/G)**|X| instead of G**X [X < 0] */
 
573
     err = mp_exptmod(&tmpG, &tmpX, P, Y);
 
574
     mp_clear_multi(&tmpG, &tmpX, NULL);
 
575
     return err;
 
576
#else 
 
577
#error mp_exptmod would always fail
 
578
     /* no invmod */
 
579
     return MP_VAL;
 
580
#endif
 
581
  }
 
582
 
 
583
/* modified diminished radix reduction */
 
584
#if defined(BN_MP_REDUCE_IS_2K_L_C) && defined(BN_MP_REDUCE_2K_L_C) && defined(BN_S_MP_EXPTMOD_C)
 
585
  if (mp_reduce_is_2k_l(P) == MP_YES) {
 
586
     return s_mp_exptmod(G, X, P, Y, 1);
 
587
  }
 
588
#endif
 
589
 
 
590
#ifdef BN_MP_DR_IS_MODULUS_C
 
591
  /* is it a DR modulus? */
 
592
  dr = mp_dr_is_modulus(P);
 
593
#else
 
594
  /* default to no */
 
595
  dr = 0;
 
596
#endif
 
597
 
 
598
#ifdef BN_MP_REDUCE_IS_2K_C
 
599
  /* if not, is it a unrestricted DR modulus? */
 
600
  if (dr == 0) {
 
601
     dr = mp_reduce_is_2k(P) << 1;
 
602
  }
 
603
#endif
 
604
    
 
605
  /* if the modulus is odd or dr != 0 use the montgomery method */
 
606
#ifdef BN_MP_EXPTMOD_FAST_C
 
607
  if (mp_isodd (P) == 1 || dr !=  0) {
 
608
    return mp_exptmod_fast (G, X, P, Y, dr);
 
609
  } else {
 
610
#endif
 
611
#ifdef BN_S_MP_EXPTMOD_C
 
612
    /* otherwise use the generic Barrett reduction technique */
 
613
    return s_mp_exptmod (G, X, P, Y, 0);
 
614
#else
 
615
#error mp_exptmod could fail
 
616
    /* no exptmod for evens */
 
617
    return MP_VAL;
 
618
#endif
 
619
#ifdef BN_MP_EXPTMOD_FAST_C
 
620
  }
 
621
#endif
 
622
}
 
623
 
 
624
 
 
625
/* compare two ints (signed)*/
 
626
static int mp_cmp (mp_int * a, mp_int * b)
 
627
{
 
628
  /* compare based on sign */
 
629
  if (a->sign != b->sign) {
 
630
     if (a->sign == MP_NEG) {
 
631
        return MP_LT;
 
632
     } else {
 
633
        return MP_GT;
 
634
     }
 
635
  }
 
636
  
 
637
  /* compare digits */
 
638
  if (a->sign == MP_NEG) {
 
639
     /* if negative compare opposite direction */
 
640
     return mp_cmp_mag(b, a);
 
641
  } else {
 
642
     return mp_cmp_mag(a, b);
 
643
  }
 
644
}
 
645
 
 
646
 
 
647
/* compare a digit */
 
648
static int mp_cmp_d(mp_int * a, mp_digit b)
 
649
{
 
650
  /* compare based on sign */
 
651
  if (a->sign == MP_NEG) {
 
652
    return MP_LT;
 
653
  }
 
654
 
 
655
  /* compare based on magnitude */
 
656
  if (a->used > 1) {
 
657
    return MP_GT;
 
658
  }
 
659
 
 
660
  /* compare the only digit of a to b */
 
661
  if (a->dp[0] > b) {
 
662
    return MP_GT;
 
663
  } else if (a->dp[0] < b) {
 
664
    return MP_LT;
 
665
  } else {
 
666
    return MP_EQ;
 
667
  }
 
668
}
 
669
 
 
670
 
 
671
/* hac 14.61, pp608 */
 
672
static int mp_invmod (mp_int * a, mp_int * b, mp_int * c)
 
673
{
 
674
  /* b cannot be negative */
 
675
  if (b->sign == MP_NEG || mp_iszero(b) == 1) {
 
676
    return MP_VAL;
 
677
  }
 
678
 
 
679
#ifdef BN_FAST_MP_INVMOD_C
 
680
  /* if the modulus is odd we can use a faster routine instead */
 
681
  if (mp_isodd (b) == 1) {
 
682
    return fast_mp_invmod (a, b, c);
 
683
  }
 
684
#endif
 
685
 
 
686
#ifdef BN_MP_INVMOD_SLOW_C
 
687
  return mp_invmod_slow(a, b, c);
 
688
#endif
 
689
 
 
690
#ifndef BN_FAST_MP_INVMOD_C
 
691
#ifndef BN_MP_INVMOD_SLOW_C
 
692
#error mp_invmod would always fail
 
693
#endif
 
694
#endif
 
695
  return MP_VAL;
 
696
}
 
697
 
 
698
 
 
699
/* get the size for an unsigned equivalent */
 
700
static int mp_unsigned_bin_size (mp_int * a)
 
701
{
 
702
  int     size = mp_count_bits (a);
 
703
  return (size / 8 + ((size & 7) != 0 ? 1 : 0));
 
704
}
 
705
 
 
706
 
 
707
/* hac 14.61, pp608 */
 
708
static int mp_invmod_slow (mp_int * a, mp_int * b, mp_int * c)
 
709
{
 
710
  mp_int  x, y, u, v, A, B, C, D;
 
711
  int     res;
 
712
 
 
713
  /* b cannot be negative */
 
714
  if (b->sign == MP_NEG || mp_iszero(b) == 1) {
 
715
    return MP_VAL;
 
716
  }
 
717
 
 
718
  /* init temps */
 
719
  if ((res = mp_init_multi(&x, &y, &u, &v, 
 
720
                           &A, &B, &C, &D, NULL)) != MP_OKAY) {
 
721
     return res;
 
722
  }
 
723
 
 
724
  /* x = a, y = b */
 
725
  if ((res = mp_mod(a, b, &x)) != MP_OKAY) {
 
726
      goto LBL_ERR;
 
727
  }
 
728
  if ((res = mp_copy (b, &y)) != MP_OKAY) {
 
729
    goto LBL_ERR;
 
730
  }
 
731
 
 
732
  /* 2. [modified] if x,y are both even then return an error! */
 
733
  if (mp_iseven (&x) == 1 && mp_iseven (&y) == 1) {
 
734
    res = MP_VAL;
 
735
    goto LBL_ERR;
 
736
  }
 
737
 
 
738
  /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
 
739
  if ((res = mp_copy (&x, &u)) != MP_OKAY) {
 
740
    goto LBL_ERR;
 
741
  }
 
742
  if ((res = mp_copy (&y, &v)) != MP_OKAY) {
 
743
    goto LBL_ERR;
 
744
  }
 
745
  mp_set (&A, 1);
 
746
  mp_set (&D, 1);
 
747
 
 
748
top:
 
749
  /* 4.  while u is even do */
 
750
  while (mp_iseven (&u) == 1) {
 
751
    /* 4.1 u = u/2 */
 
752
    if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
 
753
      goto LBL_ERR;
 
754
    }
 
755
    /* 4.2 if A or B is odd then */
 
756
    if (mp_isodd (&A) == 1 || mp_isodd (&B) == 1) {
 
757
      /* A = (A+y)/2, B = (B-x)/2 */
 
758
      if ((res = mp_add (&A, &y, &A)) != MP_OKAY) {
 
759
         goto LBL_ERR;
 
760
      }
 
761
      if ((res = mp_sub (&B, &x, &B)) != MP_OKAY) {
 
762
         goto LBL_ERR;
 
763
      }
 
764
    }
 
765
    /* A = A/2, B = B/2 */
 
766
    if ((res = mp_div_2 (&A, &A)) != MP_OKAY) {
 
767
      goto LBL_ERR;
 
768
    }
 
769
    if ((res = mp_div_2 (&B, &B)) != MP_OKAY) {
 
770
      goto LBL_ERR;
 
771
    }
 
772
  }
 
773
 
 
774
  /* 5.  while v is even do */
 
775
  while (mp_iseven (&v) == 1) {
 
776
    /* 5.1 v = v/2 */
 
777
    if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
 
778
      goto LBL_ERR;
 
779
    }
 
780
    /* 5.2 if C or D is odd then */
 
781
    if (mp_isodd (&C) == 1 || mp_isodd (&D) == 1) {
 
782
      /* C = (C+y)/2, D = (D-x)/2 */
 
783
      if ((res = mp_add (&C, &y, &C)) != MP_OKAY) {
 
784
         goto LBL_ERR;
 
785
      }
 
786
      if ((res = mp_sub (&D, &x, &D)) != MP_OKAY) {
 
787
         goto LBL_ERR;
 
788
      }
 
789
    }
 
790
    /* C = C/2, D = D/2 */
 
791
    if ((res = mp_div_2 (&C, &C)) != MP_OKAY) {
 
792
      goto LBL_ERR;
 
793
    }
 
794
    if ((res = mp_div_2 (&D, &D)) != MP_OKAY) {
 
795
      goto LBL_ERR;
 
796
    }
 
797
  }
 
798
 
 
799
  /* 6.  if u >= v then */
 
800
  if (mp_cmp (&u, &v) != MP_LT) {
 
801
    /* u = u - v, A = A - C, B = B - D */
 
802
    if ((res = mp_sub (&u, &v, &u)) != MP_OKAY) {
 
803
      goto LBL_ERR;
 
804
    }
 
805
 
 
806
    if ((res = mp_sub (&A, &C, &A)) != MP_OKAY) {
 
807
      goto LBL_ERR;
 
808
    }
 
809
 
 
810
    if ((res = mp_sub (&B, &D, &B)) != MP_OKAY) {
 
811
      goto LBL_ERR;
 
812
    }
 
813
  } else {
 
814
    /* v - v - u, C = C - A, D = D - B */
 
815
    if ((res = mp_sub (&v, &u, &v)) != MP_OKAY) {
 
816
      goto LBL_ERR;
 
817
    }
 
818
 
 
819
    if ((res = mp_sub (&C, &A, &C)) != MP_OKAY) {
 
820
      goto LBL_ERR;
 
821
    }
 
822
 
 
823
    if ((res = mp_sub (&D, &B, &D)) != MP_OKAY) {
 
824
      goto LBL_ERR;
 
825
    }
 
826
  }
 
827
 
 
828
  /* if not zero goto step 4 */
 
829
  if (mp_iszero (&u) == 0)
 
830
    goto top;
 
831
 
 
832
  /* now a = C, b = D, gcd == g*v */
 
833
 
 
834
  /* if v != 1 then there is no inverse */
 
835
  if (mp_cmp_d (&v, 1) != MP_EQ) {
 
836
    res = MP_VAL;
 
837
    goto LBL_ERR;
 
838
  }
 
839
 
 
840
  /* if its too low */
 
841
  while (mp_cmp_d(&C, 0) == MP_LT) {
 
842
      if ((res = mp_add(&C, b, &C)) != MP_OKAY) {
 
843
         goto LBL_ERR;
 
844
      }
 
845
  }
 
846
  
 
847
  /* too big */
 
848
  while (mp_cmp_mag(&C, b) != MP_LT) {
 
849
      if ((res = mp_sub(&C, b, &C)) != MP_OKAY) {
 
850
         goto LBL_ERR;
 
851
      }
 
852
  }
 
853
  
 
854
  /* C is now the inverse */
 
855
  mp_exch (&C, c);
 
856
  res = MP_OKAY;
 
857
LBL_ERR:mp_clear_multi (&x, &y, &u, &v, &A, &B, &C, &D, NULL);
 
858
  return res;
 
859
}
 
860
 
 
861
 
 
862
/* compare maginitude of two ints (unsigned) */
 
863
static int mp_cmp_mag (mp_int * a, mp_int * b)
 
864
{
 
865
  int     n;
 
866
  mp_digit *tmpa, *tmpb;
 
867
 
 
868
  /* compare based on # of non-zero digits */
 
869
  if (a->used > b->used) {
 
870
    return MP_GT;
 
871
  }
 
872
  
 
873
  if (a->used < b->used) {
 
874
    return MP_LT;
 
875
  }
 
876
 
 
877
  /* alias for a */
 
878
  tmpa = a->dp + (a->used - 1);
 
879
 
 
880
  /* alias for b */
 
881
  tmpb = b->dp + (a->used - 1);
 
882
 
 
883
  /* compare based on digits  */
 
884
  for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
 
885
    if (*tmpa > *tmpb) {
 
886
      return MP_GT;
 
887
    }
 
888
 
 
889
    if (*tmpa < *tmpb) {
 
890
      return MP_LT;
 
891
    }
 
892
  }
 
893
  return MP_EQ;
 
894
}
 
895
 
 
896
 
 
897
/* reads a unsigned char array, assumes the msb is stored first [big endian] */
 
898
static int mp_read_unsigned_bin (mp_int * a, const unsigned char *b, int c)
 
899
{
 
900
  int     res;
 
901
 
 
902
  /* make sure there are at least two digits */
 
903
  if (a->alloc < 2) {
 
904
     if ((res = mp_grow(a, 2)) != MP_OKAY) {
 
905
        return res;
 
906
     }
 
907
  }
 
908
 
 
909
  /* zero the int */
 
910
  mp_zero (a);
 
911
 
 
912
  /* read the bytes in */
 
913
  while (c-- > 0) {
 
914
    if ((res = mp_mul_2d (a, 8, a)) != MP_OKAY) {
 
915
      return res;
 
916
    }
 
917
 
 
918
#ifndef MP_8BIT
 
919
      a->dp[0] |= *b++;
 
920
      a->used += 1;
 
921
#else
 
922
      a->dp[0] = (*b & MP_MASK);
 
923
      a->dp[1] |= ((*b++ >> 7U) & 1);
 
924
      a->used += 2;
 
925
#endif
 
926
  }
 
927
  mp_clamp (a);
 
928
  return MP_OKAY;
 
929
}
 
930
 
 
931
 
 
932
/* store in unsigned [big endian] format */
 
933
static int mp_to_unsigned_bin (mp_int * a, unsigned char *b)
 
934
{
 
935
  int     x, res;
 
936
  mp_int  t;
 
937
 
 
938
  if ((res = mp_init_copy (&t, a)) != MP_OKAY) {
 
939
    return res;
 
940
  }
 
941
 
 
942
  x = 0;
 
943
  while (mp_iszero (&t) == 0) {
 
944
#ifndef MP_8BIT
 
945
      b[x++] = (unsigned char) (t.dp[0] & 255);
 
946
#else
 
947
      b[x++] = (unsigned char) (t.dp[0] | ((t.dp[1] & 0x01) << 7));
 
948
#endif
 
949
    if ((res = mp_div_2d (&t, 8, &t, NULL)) != MP_OKAY) {
 
950
      mp_clear (&t);
 
951
      return res;
 
952
    }
 
953
  }
 
954
  bn_reverse (b, x);
 
955
  mp_clear (&t);
 
956
  return MP_OKAY;
 
957
}
 
958
 
 
959
 
 
960
/* shift right by a certain bit count (store quotient in c, optional remainder in d) */
 
961
static int mp_div_2d (mp_int * a, int b, mp_int * c, mp_int * d)
 
962
{
 
963
  mp_digit D, r, rr;
 
964
  int     x, res;
 
965
  mp_int  t;
 
966
 
 
967
 
 
968
  /* if the shift count is <= 0 then we do no work */
 
969
  if (b <= 0) {
 
970
    res = mp_copy (a, c);
 
971
    if (d != NULL) {
 
972
      mp_zero (d);
 
973
    }
 
974
    return res;
 
975
  }
 
976
 
 
977
  if ((res = mp_init (&t)) != MP_OKAY) {
 
978
    return res;
 
979
  }
 
980
 
 
981
  /* get the remainder */
 
982
  if (d != NULL) {
 
983
    if ((res = mp_mod_2d (a, b, &t)) != MP_OKAY) {
 
984
      mp_clear (&t);
 
985
      return res;
 
986
    }
 
987
  }
 
988
 
 
989
  /* copy */
 
990
  if ((res = mp_copy (a, c)) != MP_OKAY) {
 
991
    mp_clear (&t);
 
992
    return res;
 
993
  }
 
994
 
 
995
  /* shift by as many digits in the bit count */
 
996
  if (b >= (int)DIGIT_BIT) {
 
997
    mp_rshd (c, b / DIGIT_BIT);
 
998
  }
 
999
 
 
1000
  /* shift any bit count < DIGIT_BIT */
 
1001
  D = (mp_digit) (b % DIGIT_BIT);
 
1002
  if (D != 0) {
 
1003
    register mp_digit *tmpc, mask, shift;
 
1004
 
 
1005
    /* mask */
 
1006
    mask = (((mp_digit)1) << D) - 1;
 
1007
 
 
1008
    /* shift for lsb */
 
1009
    shift = DIGIT_BIT - D;
 
1010
 
 
1011
    /* alias */
 
1012
    tmpc = c->dp + (c->used - 1);
 
1013
 
 
1014
    /* carry */
 
1015
    r = 0;
 
1016
    for (x = c->used - 1; x >= 0; x--) {
 
1017
      /* get the lower  bits of this word in a temp */
 
1018
      rr = *tmpc & mask;
 
1019
 
 
1020
      /* shift the current word and mix in the carry bits from the previous word */
 
1021
      *tmpc = (*tmpc >> D) | (r << shift);
 
1022
      --tmpc;
 
1023
 
 
1024
      /* set the carry to the carry bits of the current word found above */
 
1025
      r = rr;
 
1026
    }
 
1027
  }
 
1028
  mp_clamp (c);
 
1029
  if (d != NULL) {
 
1030
    mp_exch (&t, d);
 
1031
  }
 
1032
  mp_clear (&t);
 
1033
  return MP_OKAY;
 
1034
}
 
1035
 
 
1036
 
 
1037
static int mp_init_copy (mp_int * a, mp_int * b)
 
1038
{
 
1039
  int     res;
 
1040
 
 
1041
  if ((res = mp_init (a)) != MP_OKAY) {
 
1042
    return res;
 
1043
  }
 
1044
  return mp_copy (b, a);
 
1045
}
 
1046
 
 
1047
 
 
1048
/* set to zero */
 
1049
static void mp_zero (mp_int * a)
 
1050
{
 
1051
  int       n;
 
1052
  mp_digit *tmp;
 
1053
 
 
1054
  a->sign = MP_ZPOS;
 
1055
  a->used = 0;
 
1056
 
 
1057
  tmp = a->dp;
 
1058
  for (n = 0; n < a->alloc; n++) {
 
1059
     *tmp++ = 0;
 
1060
  }
 
1061
}
 
1062
 
 
1063
 
 
1064
/* copy, b = a */
 
1065
static int mp_copy (mp_int * a, mp_int * b)
 
1066
{
 
1067
  int     res, n;
 
1068
 
 
1069
  /* if dst == src do nothing */
 
1070
  if (a == b) {
 
1071
    return MP_OKAY;
 
1072
  }
 
1073
 
 
1074
  /* grow dest */
 
1075
  if (b->alloc < a->used) {
 
1076
     if ((res = mp_grow (b, a->used)) != MP_OKAY) {
 
1077
        return res;
 
1078
     }
 
1079
  }
 
1080
 
 
1081
  /* zero b and copy the parameters over */
 
1082
  {
 
1083
    register mp_digit *tmpa, *tmpb;
 
1084
 
 
1085
    /* pointer aliases */
 
1086
 
 
1087
    /* source */
 
1088
    tmpa = a->dp;
 
1089
 
 
1090
    /* destination */
 
1091
    tmpb = b->dp;
 
1092
 
 
1093
    /* copy all the digits */
 
1094
    for (n = 0; n < a->used; n++) {
 
1095
      *tmpb++ = *tmpa++;
 
1096
    }
 
1097
 
 
1098
    /* clear high digits */
 
1099
    for (; n < b->used; n++) {
 
1100
      *tmpb++ = 0;
 
1101
    }
 
1102
  }
 
1103
 
 
1104
  /* copy used count and sign */
 
1105
  b->used = a->used;
 
1106
  b->sign = a->sign;
 
1107
  return MP_OKAY;
 
1108
}
 
1109
 
 
1110
 
 
1111
/* shift right a certain amount of digits */
 
1112
static void mp_rshd (mp_int * a, int b)
 
1113
{
 
1114
  int     x;
 
1115
 
 
1116
  /* if b <= 0 then ignore it */
 
1117
  if (b <= 0) {
 
1118
    return;
 
1119
  }
 
1120
 
 
1121
  /* if b > used then simply zero it and return */
 
1122
  if (a->used <= b) {
 
1123
    mp_zero (a);
 
1124
    return;
 
1125
  }
 
1126
 
 
1127
  {
 
1128
    register mp_digit *bottom, *top;
 
1129
 
 
1130
    /* shift the digits down */
 
1131
 
 
1132
    /* bottom */
 
1133
    bottom = a->dp;
 
1134
 
 
1135
    /* top [offset into digits] */
 
1136
    top = a->dp + b;
 
1137
 
 
1138
    /* this is implemented as a sliding window where 
 
1139
     * the window is b-digits long and digits from 
 
1140
     * the top of the window are copied to the bottom
 
1141
     *
 
1142
     * e.g.
 
1143
 
 
1144
     b-2 | b-1 | b0 | b1 | b2 | ... | bb |   ---->
 
1145
                 /\                   |      ---->
 
1146
                  \-------------------/      ---->
 
1147
     */
 
1148
    for (x = 0; x < (a->used - b); x++) {
 
1149
      *bottom++ = *top++;
 
1150
    }
 
1151
 
 
1152
    /* zero the top digits */
 
1153
    for (; x < a->used; x++) {
 
1154
      *bottom++ = 0;
 
1155
    }
 
1156
  }
 
1157
  
 
1158
  /* remove excess digits */
 
1159
  a->used -= b;
 
1160
}
 
1161
 
 
1162
 
 
1163
/* swap the elements of two integers, for cases where you can't simply swap the 
 
1164
 * mp_int pointers around
 
1165
 */
 
1166
static void mp_exch (mp_int * a, mp_int * b)
 
1167
{
 
1168
  mp_int  t;
 
1169
 
 
1170
  t  = *a;
 
1171
  *a = *b;
 
1172
  *b = t;
 
1173
}
 
1174
 
 
1175
 
 
1176
/* trim unused digits 
 
1177
 *
 
1178
 * This is used to ensure that leading zero digits are
 
1179
 * trimed and the leading "used" digit will be non-zero
 
1180
 * Typically very fast.  Also fixes the sign if there
 
1181
 * are no more leading digits
 
1182
 */
 
1183
static void mp_clamp (mp_int * a)
 
1184
{
 
1185
  /* decrease used while the most significant digit is
 
1186
   * zero.
 
1187
   */
 
1188
  while (a->used > 0 && a->dp[a->used - 1] == 0) {
 
1189
    --(a->used);
 
1190
  }
 
1191
 
 
1192
  /* reset the sign flag if used == 0 */
 
1193
  if (a->used == 0) {
 
1194
    a->sign = MP_ZPOS;
 
1195
  }
 
1196
}
 
1197
 
 
1198
 
 
1199
/* grow as required */
 
1200
static int mp_grow (mp_int * a, int size)
 
1201
{
 
1202
  int     i;
 
1203
  mp_digit *tmp;
 
1204
 
 
1205
  /* if the alloc size is smaller alloc more ram */
 
1206
  if (a->alloc < size) {
 
1207
    /* ensure there are always at least MP_PREC digits extra on top */
 
1208
    size += (MP_PREC * 2) - (size % MP_PREC);
 
1209
 
 
1210
    /* reallocate the array a->dp
 
1211
     *
 
1212
     * We store the return in a temporary variable
 
1213
     * in case the operation failed we don't want
 
1214
     * to overwrite the dp member of a.
 
1215
     */
 
1216
    tmp = OPT_CAST(mp_digit) XREALLOC (a->dp, sizeof (mp_digit) * size);
 
1217
    if (tmp == NULL) {
 
1218
      /* reallocation failed but "a" is still valid [can be freed] */
 
1219
      return MP_MEM;
 
1220
    }
 
1221
 
 
1222
    /* reallocation succeeded so set a->dp */
 
1223
    a->dp = tmp;
 
1224
 
 
1225
    /* zero excess digits */
 
1226
    i        = a->alloc;
 
1227
    a->alloc = size;
 
1228
    for (; i < a->alloc; i++) {
 
1229
      a->dp[i] = 0;
 
1230
    }
 
1231
  }
 
1232
  return MP_OKAY;
 
1233
}
 
1234
 
 
1235
 
 
1236
/* b = |a| 
 
1237
 *
 
1238
 * Simple function copies the input and fixes the sign to positive
 
1239
 */
 
1240
static int mp_abs (mp_int * a, mp_int * b)
 
1241
{
 
1242
  int     res;
 
1243
 
 
1244
  /* copy a to b */
 
1245
  if (a != b) {
 
1246
     if ((res = mp_copy (a, b)) != MP_OKAY) {
 
1247
       return res;
 
1248
     }
 
1249
  }
 
1250
 
 
1251
  /* force the sign of b to positive */
 
1252
  b->sign = MP_ZPOS;
 
1253
 
 
1254
  return MP_OKAY;
 
1255
}
 
1256
 
 
1257
 
 
1258
/* set to a digit */
 
1259
static void mp_set (mp_int * a, mp_digit b)
 
1260
{
 
1261
  mp_zero (a);
 
1262
  a->dp[0] = b & MP_MASK;
 
1263
  a->used  = (a->dp[0] != 0) ? 1 : 0;
 
1264
}
 
1265
 
 
1266
 
 
1267
/* b = a/2 */
 
1268
static int mp_div_2(mp_int * a, mp_int * b)
 
1269
{
 
1270
  int     x, res, oldused;
 
1271
 
 
1272
  /* copy */
 
1273
  if (b->alloc < a->used) {
 
1274
    if ((res = mp_grow (b, a->used)) != MP_OKAY) {
 
1275
      return res;
 
1276
    }
 
1277
  }
 
1278
 
 
1279
  oldused = b->used;
 
1280
  b->used = a->used;
 
1281
  {
 
1282
    register mp_digit r, rr, *tmpa, *tmpb;
 
1283
 
 
1284
    /* source alias */
 
1285
    tmpa = a->dp + b->used - 1;
 
1286
 
 
1287
    /* dest alias */
 
1288
    tmpb = b->dp + b->used - 1;
 
1289
 
 
1290
    /* carry */
 
1291
    r = 0;
 
1292
    for (x = b->used - 1; x >= 0; x--) {
 
1293
      /* get the carry for the next iteration */
 
1294
      rr = *tmpa & 1;
 
1295
 
 
1296
      /* shift the current digit, add in carry and store */
 
1297
      *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
 
1298
 
 
1299
      /* forward carry to next iteration */
 
1300
      r = rr;
 
1301
    }
 
1302
 
 
1303
    /* zero excess digits */
 
1304
    tmpb = b->dp + b->used;
 
1305
    for (x = b->used; x < oldused; x++) {
 
1306
      *tmpb++ = 0;
 
1307
    }
 
1308
  }
 
1309
  b->sign = a->sign;
 
1310
  mp_clamp (b);
 
1311
  return MP_OKAY;
 
1312
}
 
1313
 
 
1314
 
 
1315
/* shift left by a certain bit count */
 
1316
static int mp_mul_2d (mp_int * a, int b, mp_int * c)
 
1317
{
 
1318
  mp_digit d;
 
1319
  int      res;
 
1320
 
 
1321
  /* copy */
 
1322
  if (a != c) {
 
1323
     if ((res = mp_copy (a, c)) != MP_OKAY) {
 
1324
       return res;
 
1325
     }
 
1326
  }
 
1327
 
 
1328
  if (c->alloc < (int)(c->used + b/DIGIT_BIT + 1)) {
 
1329
     if ((res = mp_grow (c, c->used + b / DIGIT_BIT + 1)) != MP_OKAY) {
 
1330
       return res;
 
1331
     }
 
1332
  }
 
1333
 
 
1334
  /* shift by as many digits in the bit count */
 
1335
  if (b >= (int)DIGIT_BIT) {
 
1336
    if ((res = mp_lshd (c, b / DIGIT_BIT)) != MP_OKAY) {
 
1337
      return res;
 
1338
    }
 
1339
  }
 
1340
 
 
1341
  /* shift any bit count < DIGIT_BIT */
 
1342
  d = (mp_digit) (b % DIGIT_BIT);
 
1343
  if (d != 0) {
 
1344
    register mp_digit *tmpc, shift, mask, r, rr;
 
1345
    register int x;
 
1346
 
 
1347
    /* bitmask for carries */
 
1348
    mask = (((mp_digit)1) << d) - 1;
 
1349
 
 
1350
    /* shift for msbs */
 
1351
    shift = DIGIT_BIT - d;
 
1352
 
 
1353
    /* alias */
 
1354
    tmpc = c->dp;
 
1355
 
 
1356
    /* carry */
 
1357
    r    = 0;
 
1358
    for (x = 0; x < c->used; x++) {
 
1359
      /* get the higher bits of the current word */
 
1360
      rr = (*tmpc >> shift) & mask;
 
1361
 
 
1362
      /* shift the current word and OR in the carry */
 
1363
      *tmpc = ((*tmpc << d) | r) & MP_MASK;
 
1364
      ++tmpc;
 
1365
 
 
1366
      /* set the carry to the carry bits of the current word */
 
1367
      r = rr;
 
1368
    }
 
1369
    
 
1370
    /* set final carry */
 
1371
    if (r != 0) {
 
1372
       c->dp[(c->used)++] = r;
 
1373
    }
 
1374
  }
 
1375
  mp_clamp (c);
 
1376
  return MP_OKAY;
 
1377
}
 
1378
 
 
1379
 
 
1380
static int mp_init_multi(mp_int *mp, ...) 
 
1381
{
 
1382
    mp_err res = MP_OKAY;      /* Assume ok until proven otherwise */
 
1383
    int n = 0;                 /* Number of ok inits */
 
1384
    mp_int* cur_arg = mp;
 
1385
    va_list args;
 
1386
 
 
1387
    va_start(args, mp);        /* init args to next argument from caller */
 
1388
    while (cur_arg != NULL) {
 
1389
        if (mp_init(cur_arg) != MP_OKAY) {
 
1390
            /* Oops - error! Back-track and mp_clear what we already
 
1391
               succeeded in init-ing, then return error.
 
1392
            */
 
1393
            va_list clean_args;
 
1394
            
 
1395
            /* end the current list */
 
1396
            va_end(args);
 
1397
            
 
1398
            /* now start cleaning up */            
 
1399
            cur_arg = mp;
 
1400
            va_start(clean_args, mp);
 
1401
            while (n--) {
 
1402
                mp_clear(cur_arg);
 
1403
                cur_arg = va_arg(clean_args, mp_int*);
 
1404
            }
 
1405
            va_end(clean_args);
 
1406
            res = MP_MEM;
 
1407
            break;
 
1408
        }
 
1409
        n++;
 
1410
        cur_arg = va_arg(args, mp_int*);
 
1411
    }
 
1412
    va_end(args);
 
1413
    return res;                /* Assumed ok, if error flagged above. */
 
1414
}
 
1415
 
 
1416
 
 
1417
static void mp_clear_multi(mp_int *mp, ...) 
 
1418
{
 
1419
    mp_int* next_mp = mp;
 
1420
    va_list args;
 
1421
    va_start(args, mp);
 
1422
    while (next_mp != NULL) {
 
1423
        mp_clear(next_mp);
 
1424
        next_mp = va_arg(args, mp_int*);
 
1425
    }
 
1426
    va_end(args);
 
1427
}
 
1428
 
 
1429
 
 
1430
/* shift left a certain amount of digits */
 
1431
static int mp_lshd (mp_int * a, int b)
 
1432
{
 
1433
  int     x, res;
 
1434
 
 
1435
  /* if its less than zero return */
 
1436
  if (b <= 0) {
 
1437
    return MP_OKAY;
 
1438
  }
 
1439
 
 
1440
  /* grow to fit the new digits */
 
1441
  if (a->alloc < a->used + b) {
 
1442
     if ((res = mp_grow (a, a->used + b)) != MP_OKAY) {
 
1443
       return res;
 
1444
     }
 
1445
  }
 
1446
 
 
1447
  {
 
1448
    register mp_digit *top, *bottom;
 
1449
 
 
1450
    /* increment the used by the shift amount then copy upwards */
 
1451
    a->used += b;
 
1452
 
 
1453
    /* top */
 
1454
    top = a->dp + a->used - 1;
 
1455
 
 
1456
    /* base */
 
1457
    bottom = a->dp + a->used - 1 - b;
 
1458
 
 
1459
    /* much like mp_rshd this is implemented using a sliding window
 
1460
     * except the window goes the otherway around.  Copying from
 
1461
     * the bottom to the top.  see bn_mp_rshd.c for more info.
 
1462
     */
 
1463
    for (x = a->used - 1; x >= b; x--) {
 
1464
      *top-- = *bottom--;
 
1465
    }
 
1466
 
 
1467
    /* zero the lower digits */
 
1468
    top = a->dp;
 
1469
    for (x = 0; x < b; x++) {
 
1470
      *top++ = 0;
 
1471
    }
 
1472
  }
 
1473
  return MP_OKAY;
 
1474
}
 
1475
 
 
1476
 
 
1477
/* returns the number of bits in an int */
 
1478
static int mp_count_bits (mp_int * a)
 
1479
{
 
1480
  int     r;
 
1481
  mp_digit q;
 
1482
 
 
1483
  /* shortcut */
 
1484
  if (a->used == 0) {
 
1485
    return 0;
 
1486
  }
 
1487
 
 
1488
  /* get number of digits and add that */
 
1489
  r = (a->used - 1) * DIGIT_BIT;
 
1490
  
 
1491
  /* take the last digit and count the bits in it */
 
1492
  q = a->dp[a->used - 1];
 
1493
  while (q > ((mp_digit) 0)) {
 
1494
    ++r;
 
1495
    q >>= ((mp_digit) 1);
 
1496
  }
 
1497
  return r;
 
1498
}
 
1499
 
 
1500
 
 
1501
/* calc a value mod 2**b */
 
1502
static int mp_mod_2d (mp_int * a, int b, mp_int * c)
 
1503
{
 
1504
  int     x, res;
 
1505
 
 
1506
  /* if b is <= 0 then zero the int */
 
1507
  if (b <= 0) {
 
1508
    mp_zero (c);
 
1509
    return MP_OKAY;
 
1510
  }
 
1511
 
 
1512
  /* if the modulus is larger than the value than return */
 
1513
  if (b >= (int) (a->used * DIGIT_BIT)) {
 
1514
    res = mp_copy (a, c);
 
1515
    return res;
 
1516
  }
 
1517
 
 
1518
  /* copy */
 
1519
  if ((res = mp_copy (a, c)) != MP_OKAY) {
 
1520
    return res;
 
1521
  }
 
1522
 
 
1523
  /* zero digits above the last digit of the modulus */
 
1524
  for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++) {
 
1525
    c->dp[x] = 0;
 
1526
  }
 
1527
  /* clear the digit that is not completely outside/inside the modulus */
 
1528
  c->dp[b / DIGIT_BIT] &=
 
1529
    (mp_digit) ((((mp_digit) 1) << (((mp_digit) b) % DIGIT_BIT)) - ((mp_digit) 1));
 
1530
  mp_clamp (c);
 
1531
  return MP_OKAY;
 
1532
}
 
1533
 
 
1534
 
 
1535
/* slower bit-bang division... also smaller */
 
1536
static int mp_div(mp_int * a, mp_int * b, mp_int * c, mp_int * d)
 
1537
{
 
1538
   mp_int ta, tb, tq, q;
 
1539
   int    res, n, n2;
 
1540
 
 
1541
  /* is divisor zero ? */
 
1542
  if (mp_iszero (b) == 1) {
 
1543
    return MP_VAL;
 
1544
  }
 
1545
 
 
1546
  /* if a < b then q=0, r = a */
 
1547
  if (mp_cmp_mag (a, b) == MP_LT) {
 
1548
    if (d != NULL) {
 
1549
      res = mp_copy (a, d);
 
1550
    } else {
 
1551
      res = MP_OKAY;
 
1552
    }
 
1553
    if (c != NULL) {
 
1554
      mp_zero (c);
 
1555
    }
 
1556
    return res;
 
1557
  }
 
1558
        
 
1559
  /* init our temps */
 
1560
  if ((res = mp_init_multi(&ta, &tb, &tq, &q, NULL) != MP_OKAY)) {
 
1561
     return res;
 
1562
  }
 
1563
 
 
1564
 
 
1565
  mp_set(&tq, 1);
 
1566
  n = mp_count_bits(a) - mp_count_bits(b);
 
1567
  if (((res = mp_abs(a, &ta)) != MP_OKAY) ||
 
1568
      ((res = mp_abs(b, &tb)) != MP_OKAY) || 
 
1569
      ((res = mp_mul_2d(&tb, n, &tb)) != MP_OKAY) ||
 
1570
      ((res = mp_mul_2d(&tq, n, &tq)) != MP_OKAY)) {
 
1571
      goto LBL_ERR;
 
1572
  }
 
1573
 
 
1574
  while (n-- >= 0) {
 
1575
     if (mp_cmp(&tb, &ta) != MP_GT) {
 
1576
        if (((res = mp_sub(&ta, &tb, &ta)) != MP_OKAY) ||
 
1577
            ((res = mp_add(&q, &tq, &q)) != MP_OKAY)) {
 
1578
           goto LBL_ERR;
 
1579
        }
 
1580
     }
 
1581
     if (((res = mp_div_2d(&tb, 1, &tb, NULL)) != MP_OKAY) ||
 
1582
         ((res = mp_div_2d(&tq, 1, &tq, NULL)) != MP_OKAY)) {
 
1583
           goto LBL_ERR;
 
1584
     }
 
1585
  }
 
1586
 
 
1587
  /* now q == quotient and ta == remainder */
 
1588
  n  = a->sign;
 
1589
  n2 = (a->sign == b->sign ? MP_ZPOS : MP_NEG);
 
1590
  if (c != NULL) {
 
1591
     mp_exch(c, &q);
 
1592
     c->sign  = (mp_iszero(c) == MP_YES) ? MP_ZPOS : n2;
 
1593
  }
 
1594
  if (d != NULL) {
 
1595
     mp_exch(d, &ta);
 
1596
     d->sign = (mp_iszero(d) == MP_YES) ? MP_ZPOS : n;
 
1597
  }
 
1598
LBL_ERR:
 
1599
   mp_clear_multi(&ta, &tb, &tq, &q, NULL);
 
1600
   return res;
 
1601
}
 
1602
 
 
1603
 
 
1604
#ifdef MP_LOW_MEM
 
1605
   #define TAB_SIZE 32
 
1606
#else
 
1607
   #define TAB_SIZE 256
 
1608
#endif
 
1609
 
 
1610
static int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
 
1611
{
 
1612
  mp_int  M[TAB_SIZE], res, mu;
 
1613
  mp_digit buf;
 
1614
  int     err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
 
1615
  int (*redux)(mp_int*,mp_int*,mp_int*);
 
1616
 
 
1617
  /* find window size */
 
1618
  x = mp_count_bits (X);
 
1619
  if (x <= 7) {
 
1620
    winsize = 2;
 
1621
  } else if (x <= 36) {
 
1622
    winsize = 3;
 
1623
  } else if (x <= 140) {
 
1624
    winsize = 4;
 
1625
  } else if (x <= 450) {
 
1626
    winsize = 5;
 
1627
  } else if (x <= 1303) {
 
1628
    winsize = 6;
 
1629
  } else if (x <= 3529) {
 
1630
    winsize = 7;
 
1631
  } else {
 
1632
    winsize = 8;
 
1633
  }
 
1634
 
 
1635
#ifdef MP_LOW_MEM
 
1636
    if (winsize > 5) {
 
1637
       winsize = 5;
 
1638
    }
 
1639
#endif
 
1640
 
 
1641
  /* init M array */
 
1642
  /* init first cell */
 
1643
  if ((err = mp_init(&M[1])) != MP_OKAY) {
 
1644
     return err; 
 
1645
  }
 
1646
 
 
1647
  /* now init the second half of the array */
 
1648
  for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
 
1649
    if ((err = mp_init(&M[x])) != MP_OKAY) {
 
1650
      for (y = 1<<(winsize-1); y < x; y++) {
 
1651
        mp_clear (&M[y]);
 
1652
      }
 
1653
      mp_clear(&M[1]);
 
1654
      return err;
 
1655
    }
 
1656
  }
 
1657
 
 
1658
  /* create mu, used for Barrett reduction */
 
1659
  if ((err = mp_init (&mu)) != MP_OKAY) {
 
1660
    goto LBL_M;
 
1661
  }
 
1662
  
 
1663
  if (redmode == 0) {
 
1664
     if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
 
1665
        goto LBL_MU;
 
1666
     }
 
1667
     redux = mp_reduce;
 
1668
  } else {
 
1669
     if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
 
1670
        goto LBL_MU;
 
1671
     }
 
1672
     redux = mp_reduce_2k_l;
 
1673
  }    
 
1674
 
 
1675
  /* create M table
 
1676
   *
 
1677
   * The M table contains powers of the base, 
 
1678
   * e.g. M[x] = G**x mod P
 
1679
   *
 
1680
   * The first half of the table is not 
 
1681
   * computed though accept for M[0] and M[1]
 
1682
   */
 
1683
  if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
 
1684
    goto LBL_MU;
 
1685
  }
 
1686
 
 
1687
  /* compute the value at M[1<<(winsize-1)] by squaring 
 
1688
   * M[1] (winsize-1) times 
 
1689
   */
 
1690
  if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
 
1691
    goto LBL_MU;
 
1692
  }
 
1693
 
 
1694
  for (x = 0; x < (winsize - 1); x++) {
 
1695
    /* square it */
 
1696
    if ((err = mp_sqr (&M[1 << (winsize - 1)], 
 
1697
                       &M[1 << (winsize - 1)])) != MP_OKAY) {
 
1698
      goto LBL_MU;
 
1699
    }
 
1700
 
 
1701
    /* reduce modulo P */
 
1702
    if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
 
1703
      goto LBL_MU;
 
1704
    }
 
1705
  }
 
1706
 
 
1707
  /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
 
1708
   * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
 
1709
   */
 
1710
  for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
 
1711
    if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
 
1712
      goto LBL_MU;
 
1713
    }
 
1714
    if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
 
1715
      goto LBL_MU;
 
1716
    }
 
1717
  }
 
1718
 
 
1719
  /* setup result */
 
1720
  if ((err = mp_init (&res)) != MP_OKAY) {
 
1721
    goto LBL_MU;
 
1722
  }
 
1723
  mp_set (&res, 1);
 
1724
 
 
1725
  /* set initial mode and bit cnt */
 
1726
  mode   = 0;
 
1727
  bitcnt = 1;
 
1728
  buf    = 0;
 
1729
  digidx = X->used - 1;
 
1730
  bitcpy = 0;
 
1731
  bitbuf = 0;
 
1732
 
 
1733
  for (;;) {
 
1734
    /* grab next digit as required */
 
1735
    if (--bitcnt == 0) {
 
1736
      /* if digidx == -1 we are out of digits */
 
1737
      if (digidx == -1) {
 
1738
        break;
 
1739
      }
 
1740
      /* read next digit and reset the bitcnt */
 
1741
      buf    = X->dp[digidx--];
 
1742
      bitcnt = (int) DIGIT_BIT;
 
1743
    }
 
1744
 
 
1745
    /* grab the next msb from the exponent */
 
1746
    y     = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
 
1747
    buf <<= (mp_digit)1;
 
1748
 
 
1749
    /* if the bit is zero and mode == 0 then we ignore it
 
1750
     * These represent the leading zero bits before the first 1 bit
 
1751
     * in the exponent.  Technically this opt is not required but it
 
1752
     * does lower the # of trivial squaring/reductions used
 
1753
     */
 
1754
    if (mode == 0 && y == 0) {
 
1755
      continue;
 
1756
    }
 
1757
 
 
1758
    /* if the bit is zero and mode == 1 then we square */
 
1759
    if (mode == 1 && y == 0) {
 
1760
      if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
 
1761
        goto LBL_RES;
 
1762
      }
 
1763
      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 
1764
        goto LBL_RES;
 
1765
      }
 
1766
      continue;
 
1767
    }
 
1768
 
 
1769
    /* else we add it to the window */
 
1770
    bitbuf |= (y << (winsize - ++bitcpy));
 
1771
    mode    = 2;
 
1772
 
 
1773
    if (bitcpy == winsize) {
 
1774
      /* ok window is filled so square as required and multiply  */
 
1775
      /* square first */
 
1776
      for (x = 0; x < winsize; x++) {
 
1777
        if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
 
1778
          goto LBL_RES;
 
1779
        }
 
1780
        if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 
1781
          goto LBL_RES;
 
1782
        }
 
1783
      }
 
1784
 
 
1785
      /* then multiply */
 
1786
      if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
 
1787
        goto LBL_RES;
 
1788
      }
 
1789
      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 
1790
        goto LBL_RES;
 
1791
      }
 
1792
 
 
1793
      /* empty window and reset */
 
1794
      bitcpy = 0;
 
1795
      bitbuf = 0;
 
1796
      mode   = 1;
 
1797
    }
 
1798
  }
 
1799
 
 
1800
  /* if bits remain then square/multiply */
 
1801
  if (mode == 2 && bitcpy > 0) {
 
1802
    /* square then multiply if the bit is set */
 
1803
    for (x = 0; x < bitcpy; x++) {
 
1804
      if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
 
1805
        goto LBL_RES;
 
1806
      }
 
1807
      if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 
1808
        goto LBL_RES;
 
1809
      }
 
1810
 
 
1811
      bitbuf <<= 1;
 
1812
      if ((bitbuf & (1 << winsize)) != 0) {
 
1813
        /* then multiply */
 
1814
        if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
 
1815
          goto LBL_RES;
 
1816
        }
 
1817
        if ((err = redux (&res, P, &mu)) != MP_OKAY) {
 
1818
          goto LBL_RES;
 
1819
        }
 
1820
      }
 
1821
    }
 
1822
  }
 
1823
 
 
1824
  mp_exch (&res, Y);
 
1825
  err = MP_OKAY;
 
1826
LBL_RES:mp_clear (&res);
 
1827
LBL_MU:mp_clear (&mu);
 
1828
LBL_M:
 
1829
  mp_clear(&M[1]);
 
1830
  for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
 
1831
    mp_clear (&M[x]);
 
1832
  }
 
1833
  return err;
 
1834
}
 
1835
 
 
1836
 
 
1837
/* computes b = a*a */
 
1838
static int mp_sqr (mp_int * a, mp_int * b)
 
1839
{
 
1840
  int     res;
 
1841
 
 
1842
#ifdef BN_MP_TOOM_SQR_C
 
1843
  /* use Toom-Cook? */
 
1844
  if (a->used >= TOOM_SQR_CUTOFF) {
 
1845
    res = mp_toom_sqr(a, b);
 
1846
  /* Karatsuba? */
 
1847
  } else 
 
1848
#endif
 
1849
#ifdef BN_MP_KARATSUBA_SQR_C
 
1850
if (a->used >= KARATSUBA_SQR_CUTOFF) {
 
1851
    res = mp_karatsuba_sqr (a, b);
 
1852
  } else 
 
1853
#endif
 
1854
  {
 
1855
#ifdef BN_FAST_S_MP_SQR_C
 
1856
    /* can we use the fast comba multiplier? */
 
1857
    if ((a->used * 2 + 1) < MP_WARRAY && 
 
1858
         a->used < 
 
1859
         (1 << (sizeof(mp_word) * CHAR_BIT - 2*DIGIT_BIT - 1))) {
 
1860
      res = fast_s_mp_sqr (a, b);
 
1861
    } else
 
1862
#endif
 
1863
#ifdef BN_S_MP_SQR_C
 
1864
      res = s_mp_sqr (a, b);
 
1865
#else
 
1866
#error mp_sqr could fail
 
1867
      res = MP_VAL;
 
1868
#endif
 
1869
  }
 
1870
  b->sign = MP_ZPOS;
 
1871
  return res;
 
1872
}
 
1873
 
 
1874
 
 
1875
/* reduces a modulo n where n is of the form 2**p - d 
 
1876
   This differs from reduce_2k since "d" can be larger
 
1877
   than a single digit.
 
1878
*/
 
1879
static int mp_reduce_2k_l(mp_int *a, mp_int *n, mp_int *d)
 
1880
{
 
1881
   mp_int q;
 
1882
   int    p, res;
 
1883
   
 
1884
   if ((res = mp_init(&q)) != MP_OKAY) {
 
1885
      return res;
 
1886
   }
 
1887
   
 
1888
   p = mp_count_bits(n);    
 
1889
top:
 
1890
   /* q = a/2**p, a = a mod 2**p */
 
1891
   if ((res = mp_div_2d(a, p, &q, a)) != MP_OKAY) {
 
1892
      goto ERR;
 
1893
   }
 
1894
   
 
1895
   /* q = q * d */
 
1896
   if ((res = mp_mul(&q, d, &q)) != MP_OKAY) { 
 
1897
      goto ERR;
 
1898
   }
 
1899
   
 
1900
   /* a = a + q */
 
1901
   if ((res = s_mp_add(a, &q, a)) != MP_OKAY) {
 
1902
      goto ERR;
 
1903
   }
 
1904
   
 
1905
   if (mp_cmp_mag(a, n) != MP_LT) {
 
1906
      s_mp_sub(a, n, a);
 
1907
      goto top;
 
1908
   }
 
1909
   
 
1910
ERR:
 
1911
   mp_clear(&q);
 
1912
   return res;
 
1913
}
 
1914
 
 
1915
 
 
1916
/* determines the setup value */
 
1917
static int mp_reduce_2k_setup_l(mp_int *a, mp_int *d)
 
1918
{
 
1919
   int    res;
 
1920
   mp_int tmp;
 
1921
   
 
1922
   if ((res = mp_init(&tmp)) != MP_OKAY) {
 
1923
      return res;
 
1924
   }
 
1925
   
 
1926
   if ((res = mp_2expt(&tmp, mp_count_bits(a))) != MP_OKAY) {
 
1927
      goto ERR;
 
1928
   }
 
1929
   
 
1930
   if ((res = s_mp_sub(&tmp, a, d)) != MP_OKAY) {
 
1931
      goto ERR;
 
1932
   }
 
1933
   
 
1934
ERR:
 
1935
   mp_clear(&tmp);
 
1936
   return res;
 
1937
}
 
1938
 
 
1939
 
 
1940
/* computes a = 2**b 
 
1941
 *
 
1942
 * Simple algorithm which zeroes the int, grows it then just sets one bit
 
1943
 * as required.
 
1944
 */
 
1945
static int mp_2expt (mp_int * a, int b)
 
1946
{
 
1947
  int     res;
 
1948
 
 
1949
  /* zero a as per default */
 
1950
  mp_zero (a);
 
1951
 
 
1952
  /* grow a to accomodate the single bit */
 
1953
  if ((res = mp_grow (a, b / DIGIT_BIT + 1)) != MP_OKAY) {
 
1954
    return res;
 
1955
  }
 
1956
 
 
1957
  /* set the used count of where the bit will go */
 
1958
  a->used = b / DIGIT_BIT + 1;
 
1959
 
 
1960
  /* put the single bit in its place */
 
1961
  a->dp[b / DIGIT_BIT] = ((mp_digit)1) << (b % DIGIT_BIT);
 
1962
 
 
1963
  return MP_OKAY;
 
1964
}
 
1965
 
 
1966
 
 
1967
/* pre-calculate the value required for Barrett reduction
 
1968
 * For a given modulus "b" it calulates the value required in "a"
 
1969
 */
 
1970
static int mp_reduce_setup (mp_int * a, mp_int * b)
 
1971
{
 
1972
  int     res;
 
1973
  
 
1974
  if ((res = mp_2expt (a, b->used * 2 * DIGIT_BIT)) != MP_OKAY) {
 
1975
    return res;
 
1976
  }
 
1977
  return mp_div (a, b, a, NULL);
 
1978
}
 
1979
 
 
1980
 
 
1981
/* reduces x mod m, assumes 0 < x < m**2, mu is 
 
1982
 * precomputed via mp_reduce_setup.
 
1983
 * From HAC pp.604 Algorithm 14.42
 
1984
 */
 
1985
static int mp_reduce (mp_int * x, mp_int * m, mp_int * mu)
 
1986
{
 
1987
  mp_int  q;
 
1988
  int     res, um = m->used;
 
1989
 
 
1990
  /* q = x */
 
1991
  if ((res = mp_init_copy (&q, x)) != MP_OKAY) {
 
1992
    return res;
 
1993
  }
 
1994
 
 
1995
  /* q1 = x / b**(k-1)  */
 
1996
  mp_rshd (&q, um - 1);         
 
1997
 
 
1998
  /* according to HAC this optimization is ok */
 
1999
  if (((unsigned long) um) > (((mp_digit)1) << (DIGIT_BIT - 1))) {
 
2000
    if ((res = mp_mul (&q, mu, &q)) != MP_OKAY) {
 
2001
      goto CLEANUP;
 
2002
    }
 
2003
  } else {
 
2004
#ifdef BN_S_MP_MUL_HIGH_DIGS_C
 
2005
    if ((res = s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
 
2006
      goto CLEANUP;
 
2007
    }
 
2008
#elif defined(BN_FAST_S_MP_MUL_HIGH_DIGS_C)
 
2009
    if ((res = fast_s_mp_mul_high_digs (&q, mu, &q, um)) != MP_OKAY) {
 
2010
      goto CLEANUP;
 
2011
    }
 
2012
#else 
 
2013
    { 
 
2014
#error mp_reduce would always fail
 
2015
      res = MP_VAL;
 
2016
      goto CLEANUP;
 
2017
    }
 
2018
#endif
 
2019
  }
 
2020
 
 
2021
  /* q3 = q2 / b**(k+1) */
 
2022
  mp_rshd (&q, um + 1);         
 
2023
 
 
2024
  /* x = x mod b**(k+1), quick (no division) */
 
2025
  if ((res = mp_mod_2d (x, DIGIT_BIT * (um + 1), x)) != MP_OKAY) {
 
2026
    goto CLEANUP;
 
2027
  }
 
2028
 
 
2029
  /* q = q * m mod b**(k+1), quick (no division) */
 
2030
  if ((res = s_mp_mul_digs (&q, m, &q, um + 1)) != MP_OKAY) {
 
2031
    goto CLEANUP;
 
2032
  }
 
2033
 
 
2034
  /* x = x - q */
 
2035
  if ((res = mp_sub (x, &q, x)) != MP_OKAY) {
 
2036
    goto CLEANUP;
 
2037
  }
 
2038
 
 
2039
  /* If x < 0, add b**(k+1) to it */
 
2040
  if (mp_cmp_d (x, 0) == MP_LT) {
 
2041
    mp_set (&q, 1);
 
2042
    if ((res = mp_lshd (&q, um + 1)) != MP_OKAY) {
 
2043
      goto CLEANUP;
 
2044
    }
 
2045
    if ((res = mp_add (x, &q, x)) != MP_OKAY) {
 
2046
      goto CLEANUP;
 
2047
    }
 
2048
  }
 
2049
 
 
2050
  /* Back off if it's too big */
 
2051
  while (mp_cmp (x, m) != MP_LT) {
 
2052
    if ((res = s_mp_sub (x, m, x)) != MP_OKAY) {
 
2053
      goto CLEANUP;
 
2054
    }
 
2055
  }
 
2056
  
 
2057
CLEANUP:
 
2058
  mp_clear (&q);
 
2059
 
 
2060
  return res;
 
2061
}
 
2062
 
 
2063
 
 
2064
/* multiplies |a| * |b| and only computes upto digs digits of result
 
2065
 * HAC pp. 595, Algorithm 14.12  Modified so you can control how 
 
2066
 * many digits of output are created.
 
2067
 */
 
2068
static int s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
 
2069
{
 
2070
  mp_int  t;
 
2071
  int     res, pa, pb, ix, iy;
 
2072
  mp_digit u;
 
2073
  mp_word r;
 
2074
  mp_digit tmpx, *tmpt, *tmpy;
 
2075
 
 
2076
  /* can we use the fast multiplier? */
 
2077
  if (((digs) < MP_WARRAY) &&
 
2078
      MIN (a->used, b->used) < 
 
2079
          (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
 
2080
    return fast_s_mp_mul_digs (a, b, c, digs);
 
2081
  }
 
2082
 
 
2083
  if ((res = mp_init_size (&t, digs)) != MP_OKAY) {
 
2084
    return res;
 
2085
  }
 
2086
  t.used = digs;
 
2087
 
 
2088
  /* compute the digits of the product directly */
 
2089
  pa = a->used;
 
2090
  for (ix = 0; ix < pa; ix++) {
 
2091
    /* set the carry to zero */
 
2092
    u = 0;
 
2093
 
 
2094
    /* limit ourselves to making digs digits of output */
 
2095
    pb = MIN (b->used, digs - ix);
 
2096
 
 
2097
    /* setup some aliases */
 
2098
    /* copy of the digit from a used within the nested loop */
 
2099
    tmpx = a->dp[ix];
 
2100
    
 
2101
    /* an alias for the destination shifted ix places */
 
2102
    tmpt = t.dp + ix;
 
2103
    
 
2104
    /* an alias for the digits of b */
 
2105
    tmpy = b->dp;
 
2106
 
 
2107
    /* compute the columns of the output and propagate the carry */
 
2108
    for (iy = 0; iy < pb; iy++) {
 
2109
      /* compute the column as a mp_word */
 
2110
      r       = ((mp_word)*tmpt) +
 
2111
                ((mp_word)tmpx) * ((mp_word)*tmpy++) +
 
2112
                ((mp_word) u);
 
2113
 
 
2114
      /* the new column is the lower part of the result */
 
2115
      *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
 
2116
 
 
2117
      /* get the carry word from the result */
 
2118
      u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
 
2119
    }
 
2120
    /* set carry if it is placed below digs */
 
2121
    if (ix + iy < digs) {
 
2122
      *tmpt = u;
 
2123
    }
 
2124
  }
 
2125
 
 
2126
  mp_clamp (&t);
 
2127
  mp_exch (&t, c);
 
2128
 
 
2129
  mp_clear (&t);
 
2130
  return MP_OKAY;
 
2131
}
 
2132
 
 
2133
 
 
2134
/* Fast (comba) multiplier
 
2135
 *
 
2136
 * This is the fast column-array [comba] multiplier.  It is 
 
2137
 * designed to compute the columns of the product first 
 
2138
 * then handle the carries afterwards.  This has the effect 
 
2139
 * of making the nested loops that compute the columns very
 
2140
 * simple and schedulable on super-scalar processors.
 
2141
 *
 
2142
 * This has been modified to produce a variable number of 
 
2143
 * digits of output so if say only a half-product is required 
 
2144
 * you don't have to compute the upper half (a feature 
 
2145
 * required for fast Barrett reduction).
 
2146
 *
 
2147
 * Based on Algorithm 14.12 on pp.595 of HAC.
 
2148
 *
 
2149
 */
 
2150
static int fast_s_mp_mul_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
 
2151
{
 
2152
  int     olduse, res, pa, ix, iz;
 
2153
  mp_digit W[MP_WARRAY];
 
2154
  register mp_word  _W;
 
2155
 
 
2156
  /* grow the destination as required */
 
2157
  if (c->alloc < digs) {
 
2158
    if ((res = mp_grow (c, digs)) != MP_OKAY) {
 
2159
      return res;
 
2160
    }
 
2161
  }
 
2162
 
 
2163
  /* number of output digits to produce */
 
2164
  pa = MIN(digs, a->used + b->used);
 
2165
 
 
2166
  /* clear the carry */
 
2167
  _W = 0;
 
2168
  for (ix = 0; ix < pa; ix++) { 
 
2169
      int      tx, ty;
 
2170
      int      iy;
 
2171
      mp_digit *tmpx, *tmpy;
 
2172
 
 
2173
      /* get offsets into the two bignums */
 
2174
      ty = MIN(b->used-1, ix);
 
2175
      tx = ix - ty;
 
2176
 
 
2177
      /* setup temp aliases */
 
2178
      tmpx = a->dp + tx;
 
2179
      tmpy = b->dp + ty;
 
2180
 
 
2181
      /* this is the number of times the loop will iterrate, essentially 
 
2182
         while (tx++ < a->used && ty-- >= 0) { ... }
 
2183
       */
 
2184
      iy = MIN(a->used-tx, ty+1);
 
2185
 
 
2186
      /* execute loop */
 
2187
      for (iz = 0; iz < iy; ++iz) {
 
2188
         _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
 
2189
 
 
2190
      }
 
2191
 
 
2192
      /* store term */
 
2193
      W[ix] = ((mp_digit)_W) & MP_MASK;
 
2194
 
 
2195
      /* make next carry */
 
2196
      _W = _W >> ((mp_word)DIGIT_BIT);
 
2197
 }
 
2198
 
 
2199
  /* setup dest */
 
2200
  olduse  = c->used;
 
2201
  c->used = pa;
 
2202
 
 
2203
  {
 
2204
    register mp_digit *tmpc;
 
2205
    tmpc = c->dp;
 
2206
    for (ix = 0; ix < pa+1; ix++) {
 
2207
      /* now extract the previous digit [below the carry] */
 
2208
      *tmpc++ = W[ix];
 
2209
    }
 
2210
 
 
2211
    /* clear unused digits [that existed in the old copy of c] */
 
2212
    for (; ix < olduse; ix++) {
 
2213
      *tmpc++ = 0;
 
2214
    }
 
2215
  }
 
2216
  mp_clamp (c);
 
2217
  return MP_OKAY;
 
2218
}
 
2219
 
 
2220
 
 
2221
/* init an mp_init for a given size */
 
2222
static int mp_init_size (mp_int * a, int size)
 
2223
{
 
2224
  int x;
 
2225
 
 
2226
  /* pad size so there are always extra digits */
 
2227
  size += (MP_PREC * 2) - (size % MP_PREC);     
 
2228
  
 
2229
  /* alloc mem */
 
2230
  a->dp = OPT_CAST(mp_digit) XMALLOC (sizeof (mp_digit) * size);
 
2231
  if (a->dp == NULL) {
 
2232
    return MP_MEM;
 
2233
  }
 
2234
 
 
2235
  /* set the members */
 
2236
  a->used  = 0;
 
2237
  a->alloc = size;
 
2238
  a->sign  = MP_ZPOS;
 
2239
 
 
2240
  /* zero the digits */
 
2241
  for (x = 0; x < size; x++) {
 
2242
      a->dp[x] = 0;
 
2243
  }
 
2244
 
 
2245
  return MP_OKAY;
 
2246
}
 
2247
 
 
2248
 
 
2249
/* low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16 */
 
2250
static int s_mp_sqr (mp_int * a, mp_int * b)
 
2251
{
 
2252
  mp_int  t;
 
2253
  int     res, ix, iy, pa;
 
2254
  mp_word r;
 
2255
  mp_digit u, tmpx, *tmpt;
 
2256
 
 
2257
  pa = a->used;
 
2258
  if ((res = mp_init_size (&t, 2*pa + 1)) != MP_OKAY) {
 
2259
    return res;
 
2260
  }
 
2261
 
 
2262
  /* default used is maximum possible size */
 
2263
  t.used = 2*pa + 1;
 
2264
 
 
2265
  for (ix = 0; ix < pa; ix++) {
 
2266
    /* first calculate the digit at 2*ix */
 
2267
    /* calculate double precision result */
 
2268
    r = ((mp_word) t.dp[2*ix]) +
 
2269
        ((mp_word)a->dp[ix])*((mp_word)a->dp[ix]);
 
2270
 
 
2271
    /* store lower part in result */
 
2272
    t.dp[ix+ix] = (mp_digit) (r & ((mp_word) MP_MASK));
 
2273
 
 
2274
    /* get the carry */
 
2275
    u           = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
 
2276
 
 
2277
    /* left hand side of A[ix] * A[iy] */
 
2278
    tmpx        = a->dp[ix];
 
2279
 
 
2280
    /* alias for where to store the results */
 
2281
    tmpt        = t.dp + (2*ix + 1);
 
2282
    
 
2283
    for (iy = ix + 1; iy < pa; iy++) {
 
2284
      /* first calculate the product */
 
2285
      r       = ((mp_word)tmpx) * ((mp_word)a->dp[iy]);
 
2286
 
 
2287
      /* now calculate the double precision result, note we use
 
2288
       * addition instead of *2 since it's easier to optimize
 
2289
       */
 
2290
      r       = ((mp_word) *tmpt) + r + r + ((mp_word) u);
 
2291
 
 
2292
      /* store lower part */
 
2293
      *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
 
2294
 
 
2295
      /* get carry */
 
2296
      u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
 
2297
    }
 
2298
    /* propagate upwards */
 
2299
    while (u != ((mp_digit) 0)) {
 
2300
      r       = ((mp_word) *tmpt) + ((mp_word) u);
 
2301
      *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
 
2302
      u       = (mp_digit)(r >> ((mp_word) DIGIT_BIT));
 
2303
    }
 
2304
  }
 
2305
 
 
2306
  mp_clamp (&t);
 
2307
  mp_exch (&t, b);
 
2308
  mp_clear (&t);
 
2309
  return MP_OKAY;
 
2310
}
 
2311
 
 
2312
 
 
2313
/* multiplies |a| * |b| and does not compute the lower digs digits
 
2314
 * [meant to get the higher part of the product]
 
2315
 */
 
2316
static int s_mp_mul_high_digs (mp_int * a, mp_int * b, mp_int * c, int digs)
 
2317
{
 
2318
  mp_int  t;
 
2319
  int     res, pa, pb, ix, iy;
 
2320
  mp_digit u;
 
2321
  mp_word r;
 
2322
  mp_digit tmpx, *tmpt, *tmpy;
 
2323
 
 
2324
  /* can we use the fast multiplier? */
 
2325
#ifdef BN_FAST_S_MP_MUL_HIGH_DIGS_C
 
2326
  if (((a->used + b->used + 1) < MP_WARRAY)
 
2327
      && MIN (a->used, b->used) < (1 << ((CHAR_BIT * sizeof (mp_word)) - (2 * DIGIT_BIT)))) {
 
2328
    return fast_s_mp_mul_high_digs (a, b, c, digs);
 
2329
  }
 
2330
#endif
 
2331
 
 
2332
  if ((res = mp_init_size (&t, a->used + b->used + 1)) != MP_OKAY) {
 
2333
    return res;
 
2334
  }
 
2335
  t.used = a->used + b->used + 1;
 
2336
 
 
2337
  pa = a->used;
 
2338
  pb = b->used;
 
2339
  for (ix = 0; ix < pa; ix++) {
 
2340
    /* clear the carry */
 
2341
    u = 0;
 
2342
 
 
2343
    /* left hand side of A[ix] * B[iy] */
 
2344
    tmpx = a->dp[ix];
 
2345
 
 
2346
    /* alias to the address of where the digits will be stored */
 
2347
    tmpt = &(t.dp[digs]);
 
2348
 
 
2349
    /* alias for where to read the right hand side from */
 
2350
    tmpy = b->dp + (digs - ix);
 
2351
 
 
2352
    for (iy = digs - ix; iy < pb; iy++) {
 
2353
      /* calculate the double precision result */
 
2354
      r       = ((mp_word)*tmpt) +
 
2355
                ((mp_word)tmpx) * ((mp_word)*tmpy++) +
 
2356
                ((mp_word) u);
 
2357
 
 
2358
      /* get the lower part */
 
2359
      *tmpt++ = (mp_digit) (r & ((mp_word) MP_MASK));
 
2360
 
 
2361
      /* carry the carry */
 
2362
      u       = (mp_digit) (r >> ((mp_word) DIGIT_BIT));
 
2363
    }
 
2364
    *tmpt = u;
 
2365
  }
 
2366
  mp_clamp (&t);
 
2367
  mp_exch (&t, c);
 
2368
  mp_clear (&t);
 
2369
  return MP_OKAY;
 
2370
}