Coverage Report

Created: 2017-08-21 19:50

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/imath/imath.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
  Name:     imath.c
3
  Purpose:  Arbitrary precision integer arithmetic routines.
4
  Author:   M. J. Fromberger <http://spinning-yarns.org/michael/>
5
6
  Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
7
8
  Permission is hereby granted, free of charge, to any person obtaining a copy
9
  of this software and associated documentation files (the "Software"), to deal
10
  in the Software without restriction, including without limitation the rights
11
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12
  copies of the Software, and to permit persons to whom the Software is
13
  furnished to do so, subject to the following conditions:
14
15
  The above copyright notice and this permission notice shall be included in
16
  all copies or substantial portions of the Software.
17
18
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL THE
21
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24
  SOFTWARE.
25
 */
26
27
#include "imath.h"
28
29
#if DEBUG
30
#include <stdio.h>
31
#endif
32
33
#include <stdlib.h>
34
#include <string.h>
35
#include <ctype.h>
36
37
#include <assert.h>
38
39
#if DEBUG
40
#define STATIC /* public */
41
#else
42
#define STATIC static
43
#endif
44
45
const mp_result MP_OK     = 0;  /* no error, all is well  */
46
const mp_result MP_FALSE  = 0;  /* boolean false          */
47
const mp_result MP_TRUE   = -1; /* boolean true           */
48
const mp_result MP_MEMORY = -2; /* out of memory          */
49
const mp_result MP_RANGE  = -3; /* argument out of range  */
50
const mp_result MP_UNDEF  = -4; /* result undefined       */
51
const mp_result MP_TRUNC  = -5; /* output truncated       */
52
const mp_result MP_BADARG = -6; /* invalid null argument  */
53
const mp_result MP_MINERR = -6;
54
55
const mp_sign   MP_NEG  = 1;    /* value is strictly negative */
56
const mp_sign   MP_ZPOS = 0;    /* value is non-negative      */
57
58
STATIC const char *s_unknown_err = "unknown result code";
59
STATIC const char *s_error_msg[] = {
60
  "error code 0",
61
  "boolean true",
62
  "out of memory",
63
  "argument out of range",
64
  "result undefined",
65
  "output truncated",
66
  "invalid argument",
67
  NULL
68
};
69
70
/* Argument checking macros
71
   Use CHECK() where a return value is required; NRCHECK() elsewhere */
72
418M
#define CHECK(TEST)   assert(TEST)
73
113M
#define NRCHECK(TEST) assert(TEST)
74
75
/* The ith entry of this table gives the value of log_i(2).
76
77
   An integer value n requires ceil(log_i(n)) digits to be represented
78
   in base i.  Since it is easy to compute lg(n), by counting bits, we
79
   can compute log_i(n) = lg(n) * log_i(2).
80
81
   The use of this table eliminates a dependency upon linkage against
82
   the standard math libraries.
83
84
   If MP_MAX_RADIX is increased, this table should be expanded too.
85
 */
86
STATIC const double s_log2[] = {
87
   0.000000000, 0.000000000, 1.000000000, 0.630929754,  /* (D)(D) 2  3 */
88
   0.500000000, 0.430676558, 0.386852807, 0.356207187,  /*  4  5  6  7 */
89
   0.333333333, 0.315464877, 0.301029996, 0.289064826,  /*  8  9 10 11 */
90
   0.278942946, 0.270238154, 0.262649535, 0.255958025,  /* 12 13 14 15 */
91
   0.250000000, 0.244650542, 0.239812467, 0.235408913,  /* 16 17 18 19 */
92
   0.231378213, 0.227670249, 0.224243824, 0.221064729,  /* 20 21 22 23 */
93
   0.218104292, 0.215338279, 0.212746054, 0.210309918,  /* 24 25 26 27 */
94
   0.208014598, 0.205846832, 0.203795047, 0.201849087,  /* 28 29 30 31 */
95
   0.200000000, 0.198239863, 0.196561632, 0.194959022,  /* 32 33 34 35 */
96
   0.193426404,                                         /* 36          */
97
};
98
99
100
101
/* Return the number of digits needed to represent a static value */
102
47.8M
#define MP_VALUE_DIGITS(V) \
103
47.8M
((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
104
105
/* Round precision P to nearest word boundary */
106
20.7M
#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
107
108
/* Set array P of S digits to zero */
109
10.0M
#define ZERO(P, S) \
110
10.0M
do{ \
111
10.0M
  mp_size i__ = (S) * sizeof(mp_digit); \
112
10.0M
  mp_digit *p__ = (P); \
113
10.0M
  memset(p__, 0, i__); \
114
10.0M
} while(0)
115
116
/* Copy S digits from array P to array Q */
117
111M
#define COPY(P, Q, S) \
118
111M
do{ \
119
111M
  mp_size i__ = (S) * sizeof(mp_digit); \
120
111M
  mp_digit *p__ = (P), *q__ = (Q); \
121
111M
  memcpy(q__, p__, i__); \
122
111M
} while(0)
123
124
/* Reverse N elements of type T in array A */
125
0
#define REV(T, A, N) \
126
0
do{ \
127
0
  T *u_ = (A), *v_ = u_ + (N) - 1; \
128
0
  while (
u_ < v_0
)
{ \0
129
0
    T xch = *u_; \
130
0
    *u_++ = *v_; \
131
0
    *v_-- = xch; \
132
0
  } \
133
0
} while(0)
134
135
171M
#define CLAMP(Z) \
136
171M
do{ \
137
171M
  mp_int z_ = (Z); \
138
171M
  mp_size uz_ = MP_USED(z_); \
139
171M
  mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \
140
203M
  while (
uz_ > 1 && 203M
(*dz_-- == 0)145M
) \
141
32.4M
    --uz_; \
142
171M
  MP_USED(z_) = uz_; \
143
171M
} while(0)
144
145
/* Select min/max.  Do not provide expressions for which multiple
146
   evaluation would be problematic, e.g. x++ */
147
2.51M
#define MIN(A, B) 
((B)<(A)?2.51M
(B)510k
:
(A)2.00M
)
148
100M
#define MAX(A, B) 
((B)>(A)?100M
(B)46.4M
:
(A)53.8M
)
149
150
/* Exchange lvalues A and B of type T, e.g.
151
   SWAP(int, x, y) where x and y are variables of type int. */
152
7.45M
#define SWAP(T, A, B) \
153
7.45M
do{ \
154
7.45M
  T t_ = (A); \
155
7.45M
  A = (B); \
156
7.45M
  B = t_; \
157
7.45M
} while(0)
158
159
/* Used to set up and access simple temp stacks within functions. */
160
#define DECLARE_TEMP(N) \
161
3.91M
  mpz_t temp[(N)]; \
162
3.91M
  int last__ = 0
163
#define CLEANUP_TEMP() \
164
2.45M
 CLEANUP: \
165
4.09M
  while (--last__ >= 0) \
166
2.45M
    
mp_int_clear(1.64M
TEMP1.64M
(last__))
167
3.28M
#define TEMP(K) (temp + (K))
168
1.64M
#define LAST_TEMP() TEMP(last__)
169
1.64M
#define SETUP(E) \
170
1.64M
do{ \
171
1.64M
  if ((res = (E)) != MP_OK) \
172
0
    goto CLEANUP; \
173
1.64M
  ++(last__); \
174
1.64M
} while(0)
175
176
/* Compare value to zero. */
177
227M
#define CMPZ(Z) \
178
227M
(((Z)->used==1&&
(Z)->digits[0]==077.6M
)?
04.23M
:
((Z)->sign==MP_NEG)?223M
-1190M
:
132.3M
)
179
180
/* Multiply X by Y into Z, ignoring signs.  Requires that Z have
181
   enough storage preallocated to hold the result. */
182
0
#define UMUL(X, Y, Z) \
183
0
do{ \
184
0
  mp_size ua_ = 
MP_USED0
(X), ub_ =
MP_USED0
(Y); \
185
0
  mp_size o_ = ua_ + ub_; \
186
0
  ZERO(MP_DIGITS(Z), o_); \
187
0
  (void) s_kmul(
MP_DIGITS0
(X),
MP_DIGITS0
(Y),
MP_DIGITS0
(Z), ua_, ub_); \
188
0
  MP_USED(Z) = o_; \
189
0
  CLAMP(Z); \
190
0
} while(0)
191
192
/* Square X into Z.  Requires that Z have enough storage to hold the
193
   result. */
194
0
#define USQR(X, Z) \
195
0
do{ \
196
0
  mp_size ua_ = MP_USED(X), o_ = ua_ + ua_; \
197
0
  ZERO(MP_DIGITS(Z), o_); \
198
0
  (void) s_ksqr(
MP_DIGITS0
(X),
MP_DIGITS0
(Z), ua_); \
199
0
  MP_USED(Z) = o_; \
200
0
  CLAMP(Z); \
201
0
} while(0)
202
203
266M
#define UPPER_HALF(W)           
((mp_word)((W) >> 266M
MP_DIGIT_BIT266M
))
204
266M
#define LOWER_HALF(W)           ((mp_digit)(W))
205
0
#define HIGH_BIT_SET(W)         
((W) >> (0
MP_WORD_BIT0
- 1))
206
0
#define ADD_WILL_OVERFLOW(W, V) 
((0
MP_WORD_MAX0
- (V)) < (W))
207
208
209
210
/* Default number of digits allocated to a new mp_int */
211
#if IMATH_TEST
212
mp_size default_precision = MP_DEFAULT_PREC;
213
#else
214
STATIC const mp_size default_precision = MP_DEFAULT_PREC;
215
#endif
216
217
/* Minimum number of digits to invoke recursive multiply */
218
#if IMATH_TEST
219
mp_size multiply_threshold = MP_MULT_THRESH;
220
#else
221
STATIC const mp_size multiply_threshold = MP_MULT_THRESH;
222
#endif
223
224
/* Allocate a buffer of (at least) num digits, or return
225
   NULL if that couldn't be done.  */
226
STATIC mp_digit *s_alloc(mp_size num);
227
228
/* Release a buffer of digits allocated by s_alloc(). */
229
STATIC void s_free(void *ptr);
230
231
/* Insure that z has at least min digits allocated, resizing if
232
   necessary.  Returns true if successful, false if out of memory. */
233
STATIC int  s_pad(mp_int z, mp_size min);
234
235
/* Fill in a "fake" mp_int on the stack with a given value */
236
STATIC void      s_fake(mp_int z, mp_small value, mp_digit vbuf[]);
237
STATIC void      s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[]);
238
239
/* Compare two runs of digits of given length, returns <0, 0, >0 */
240
STATIC int       s_cdig(mp_digit *da, mp_digit *db, mp_size len);
241
242
/* Pack the unsigned digits of v into array t */
243
STATIC int       s_uvpack(mp_usmall v, mp_digit t[]);
244
245
/* Compare magnitudes of a and b, returns <0, 0, >0 */
246
STATIC int       s_ucmp(mp_int a, mp_int b);
247
248
/* Compare magnitudes of a and v, returns <0, 0, >0 */
249
STATIC int       s_vcmp(mp_int a, mp_small v);
250
STATIC int       s_uvcmp(mp_int a, mp_usmall uv);
251
252
/* Unsigned magnitude addition; assumes dc is big enough.
253
   Carry out is returned (no memory allocated). */
254
STATIC mp_digit  s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
255
            mp_size size_a, mp_size size_b);
256
257
/* Unsigned magnitude subtraction.  Assumes dc is big enough. */
258
STATIC void      s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
259
            mp_size size_a, mp_size size_b);
260
261
/* Unsigned recursive multiplication.  Assumes dc is big enough. */
262
STATIC int       s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
263
      mp_size size_a, mp_size size_b);
264
265
/* Unsigned magnitude multiplication.  Assumes dc is big enough. */
266
STATIC void      s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
267
      mp_size size_a, mp_size size_b);
268
269
/* Unsigned recursive squaring.  Assumes dc is big enough. */
270
STATIC int       s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
271
272
/* Unsigned magnitude squaring.  Assumes dc is big enough. */
273
STATIC void      s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
274
275
/* Single digit addition.  Assumes a is big enough. */
276
STATIC void      s_dadd(mp_int a, mp_digit b);
277
278
/* Single digit multiplication.  Assumes a is big enough. */
279
STATIC void      s_dmul(mp_int a, mp_digit b);
280
281
/* Single digit multiplication on buffers; assumes dc is big enough. */
282
STATIC void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
283
       mp_size size_a);
284
285
/* Single digit division.  Replaces a with the quotient,
286
   returns the remainder.  */
287
STATIC mp_digit  s_ddiv(mp_int a, mp_digit b);
288
289
/* Quick division by a power of 2, replaces z (no allocation) */
290
STATIC void      s_qdiv(mp_int z, mp_size p2);
291
292
/* Quick remainder by a power of 2, replaces z (no allocation) */
293
STATIC void      s_qmod(mp_int z, mp_size p2);
294
295
/* Quick multiplication by a power of 2, replaces z.
296
   Allocates if necessary; returns false in case this fails. */
297
STATIC int       s_qmul(mp_int z, mp_size p2);
298
299
/* Quick subtraction from a power of 2, replaces z.
300
   Allocates if necessary; returns false in case this fails. */
301
STATIC int       s_qsub(mp_int z, mp_size p2);
302
303
/* Return maximum k such that 2^k divides z. */
304
STATIC int       s_dp2k(mp_int z);
305
306
/* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
307
STATIC int       s_isp2(mp_int z);
308
309
/* Set z to 2^k.  May allocate; returns false in case this fails. */
310
STATIC int       s_2expt(mp_int z, mp_small k);
311
312
/* Normalize a and b for division, returns normalization constant */
313
STATIC int       s_norm(mp_int a, mp_int b);
314
315
/* Compute constant mu for Barrett reduction, given modulus m, result
316
   replaces z, m is untouched. */
317
STATIC mp_result s_brmu(mp_int z, mp_int m);
318
319
/* Reduce a modulo m, using Barrett's algorithm. */
320
STATIC int       s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
321
322
/* Modular exponentiation, using Barrett reduction */
323
STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
324
325
/* Unsigned magnitude division.  Assumes |a| > |b|.  Allocates temporaries;
326
   overwrites a with quotient, b with remainder. */
327
STATIC mp_result s_udiv_knuth(mp_int a, mp_int b);
328
329
/* Compute the number of digits in radix r required to represent the given
330
   value.  Does not account for sign flags, terminators, etc. */
331
STATIC int       s_outlen(mp_int z, mp_size r);
332
333
/* Guess how many digits of precision will be needed to represent a radix r
334
   value of the specified number of digits.  Returns a value guaranteed to be
335
   no smaller than the actual number required. */
336
STATIC mp_size   s_inlen(int len, mp_size r);
337
338
/* Convert a character to a digit value in radix r, or
339
   -1 if out of range */
340
STATIC int       s_ch2val(char c, int r);
341
342
/* Convert a digit value to a character */
343
STATIC char      s_val2ch(int v, int caps);
344
345
/* Take 2's complement of a buffer in place */
346
STATIC void      s_2comp(unsigned char *buf, int len);
347
348
/* Convert a value to binary, ignoring sign.  On input, *limpos is the bound on
349
   how many bytes should be written to buf; on output, *limpos is set to the
350
   number of bytes actually written. */
351
STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad);
352
353
#if DEBUG
354
/* Dump a representation of the mp_int to standard output */
355
void      s_print(char *tag, mp_int z);
356
void      s_print_buf(char *tag, mp_digit *buf, mp_size num);
357
#endif
358
359
mp_result mp_int_init(mp_int z)
360
36.8M
{
361
36.8M
  if (z == NULL)
362
0
    return MP_BADARG;
363
36.8M
364
36.8M
  z->single = 0;
365
36.8M
  z->digits = &(z->single);
366
36.8M
  z->alloc  = 1;
367
36.8M
  z->used   = 1;
368
36.8M
  z->sign   = MP_ZPOS;
369
36.8M
370
36.8M
  return MP_OK;
371
36.8M
}
372
373
mp_int    mp_int_alloc(void)
374
30.6M
{
375
30.6M
  mp_int out = malloc(sizeof(mpz_t));
376
30.6M
377
30.6M
  if (out != NULL)
378
30.6M
    mp_int_init(out);
379
30.6M
380
30.6M
  return out;
381
30.6M
}
382
383
mp_result mp_int_init_size(mp_int z, mp_size prec)
384
5.95M
{
385
5.95M
  CHECK(z != NULL);
386
5.95M
387
5.95M
  if (prec == 0)
388
0
    prec = default_precision;
389
5.95M
  else 
if (5.95M
prec == 15.95M
)
390
332k
    return mp_int_init(z);
391
5.95M
  else
392
5.62M
    
prec = (mp_size) 5.62M
ROUND_PREC5.62M
(prec);
393
5.95M
394
5.62M
  
if (5.62M
(5.62M
MP_DIGITS5.62M
(z) = s_alloc(prec)) == NULL)
395
0
    return MP_MEMORY;
396
5.62M
397
5.62M
  z->digits[0] = 0;
398
5.62M
  MP_USED(z) = 1;
399
5.62M
  MP_ALLOC(z) = prec;
400
5.62M
  MP_SIGN(z) = MP_ZPOS;
401
5.62M
402
5.62M
  return MP_OK;
403
5.95M
}
404
405
mp_result mp_int_init_copy(mp_int z, mp_int old)
406
6.86M
{
407
6.86M
  mp_result res;
408
6.86M
  mp_size uold;
409
6.86M
410
6.86M
  CHECK(z != NULL && old != NULL);
411
6.86M
412
6.86M
  uold = MP_USED(old);
413
6.86M
  if (
uold == 16.86M
)
{2.60M
414
2.60M
    mp_int_init(z);
415
6.86M
  }
416
6.86M
  else {
417
4.25M
    mp_size target = MAX(uold, default_precision);
418
4.25M
419
4.25M
    if ((res = mp_int_init_size(z, target)) != MP_OK)
420
0
      return res;
421
6.86M
  }
422
6.86M
423
6.86M
  
MP_USED6.86M
(z) = uold;6.86M
424
6.86M
  
MP_SIGN6.86M
(z) =
MP_SIGN6.86M
(old);
425
6.86M
  COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
426
6.86M
427
6.86M
  return MP_OK;
428
6.86M
}
429
430
mp_result mp_int_init_value(mp_int z, mp_small value)
431
0
{
432
0
  mpz_t    vtmp;
433
0
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
434
0
435
0
  s_fake(&vtmp, value, vbuf);
436
0
  return mp_int_init_copy(z, &vtmp);
437
0
}
438
439
mp_result mp_int_init_uvalue(mp_int z, mp_usmall uvalue)
440
165k
{
441
165k
  mpz_t    vtmp;
442
165k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
443
165k
444
165k
  s_ufake(&vtmp, uvalue, vbuf);
445
165k
  return mp_int_init_copy(z, &vtmp);
446
165k
}
447
448
mp_result  mp_int_set_value(mp_int z, mp_small value)
449
852k
{
450
852k
  mpz_t    vtmp;
451
852k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
452
852k
453
852k
  s_fake(&vtmp, value, vbuf);
454
852k
  return mp_int_copy(&vtmp, z);
455
852k
}
456
457
mp_result  mp_int_set_uvalue(mp_int z, mp_usmall uvalue)
458
48.0k
{
459
48.0k
  mpz_t    vtmp;
460
48.0k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
461
48.0k
462
48.0k
  s_ufake(&vtmp, uvalue, vbuf);
463
48.0k
  return mp_int_copy(&vtmp, z);
464
48.0k
}
465
466
void      mp_int_clear(mp_int z)
467
42.4M
{
468
42.4M
  if (z == NULL)
469
0
    return;
470
42.4M
471
42.4M
  
if (42.4M
MP_DIGITS42.4M
(z) != NULL42.4M
)
{42.4M
472
42.4M
    if (
MP_DIGITS42.4M
(z) != &(z->single)42.4M
)
473
18.8M
      
s_free(18.8M
MP_DIGITS18.8M
(z));
474
42.4M
475
42.4M
    MP_DIGITS(z) = NULL;
476
42.4M
  }
477
42.4M
}
478
479
void      mp_int_free(mp_int z)
480
30.6M
{
481
30.6M
  NRCHECK(z != NULL);
482
30.6M
483
30.6M
  mp_int_clear(z);
484
30.6M
  free(z); /* note: NOT s_free() */
485
30.6M
}
486
487
mp_result mp_int_copy(mp_int a, mp_int c)
488
96.9M
{
489
96.9M
  CHECK(a != NULL && c != NULL);
490
96.9M
491
96.9M
  if (
a != c96.9M
)
{91.9M
492
91.9M
    mp_size ua = MP_USED(a);
493
91.9M
    mp_digit *da, *dc;
494
91.9M
495
91.9M
    if (!s_pad(c, ua))
496
0
      return MP_MEMORY;
497
91.9M
498
91.9M
    
da = 91.9M
MP_DIGITS91.9M
(a); dc =
MP_DIGITS91.9M
(c);
499
91.9M
    COPY(da, dc, ua);
500
91.9M
501
91.9M
    MP_USED(c) = ua;
502
91.9M
    
MP_SIGN91.9M
(c) =
MP_SIGN91.9M
(a);
503
96.9M
  }
504
96.9M
505
96.9M
  return MP_OK;
506
96.9M
}
507
508
void      mp_int_swap(mp_int a, mp_int c)
509
0
{
510
0
  if (
a != c0
)
{0
511
0
    mpz_t tmp = *a;
512
0
513
0
    *a = *c;
514
0
    *c = tmp;
515
0
516
0
    if (
MP_DIGITS0
(a) == &(c->single)0
)
517
0
      
MP_DIGITS0
(a) = &(a->single)0
;
518
0
    if (
MP_DIGITS0
(c) == &(a->single)0
)
519
0
      
MP_DIGITS0
(c) = &(c->single)0
;
520
0
  }
521
0
}
522
523
void      mp_int_zero(mp_int z)
524
19.3M
{
525
19.3M
  NRCHECK(z != NULL);
526
19.3M
527
19.3M
  z->digits[0] = 0;
528
19.3M
  MP_USED(z) = 1;
529
19.3M
  MP_SIGN(z) = MP_ZPOS;
530
19.3M
}
531
532
mp_result mp_int_abs(mp_int a, mp_int c)
533
3.24M
{
534
3.24M
  mp_result res;
535
3.24M
536
3.24M
  CHECK(a != NULL && c != NULL);
537
3.24M
538
3.24M
  if ((res = mp_int_copy(a, c)) != MP_OK)
539
0
    return res;
540
3.24M
541
3.24M
  
MP_SIGN3.24M
(c) = MP_ZPOS;3.24M
542
3.24M
  return MP_OK;
543
3.24M
}
544
545
mp_result mp_int_neg(mp_int a, mp_int c)
546
66.4M
{
547
66.4M
  mp_result res;
548
66.4M
549
66.4M
  CHECK(a != NULL && c != NULL);
550
66.4M
551
66.4M
  if ((res = mp_int_copy(a, c)) != MP_OK)
552
0
    return res;
553
66.4M
554
66.4M
  
if (66.4M
CMPZ66.4M
(c) != 066.4M
)
555
66.4M
    
MP_SIGN66.4M
(c) = 1 - 66.4M
MP_SIGN66.4M
(a);
556
66.4M
557
66.4M
  return MP_OK;
558
66.4M
}
559
560
mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
561
12.8M
{
562
12.8M
  mp_size ua, ub, uc, max;
563
12.8M
564
12.8M
  CHECK(a != NULL && b != NULL && c != NULL);
565
12.8M
566
12.8M
  ua = 
MP_USED12.8M
(a); ub =
MP_USED12.8M
(b); uc =
MP_USED12.8M
(c);
567
12.8M
  max = MAX(ua, ub);
568
12.8M
569
12.8M
  if (
MP_SIGN12.8M
(a) == 12.8M
MP_SIGN12.8M
(b))
{7.58M
570
7.58M
    /* Same sign -- add magnitudes, preserve sign of addends */
571
7.58M
    mp_digit carry;
572
7.58M
573
7.58M
    if (!s_pad(c, max))
574
0
      return MP_MEMORY;
575
7.58M
576
7.58M
    
carry = s_uadd(7.58M
MP_DIGITS7.58M
(a),
MP_DIGITS7.58M
(b),
MP_DIGITS7.58M
(c), ua, ub);
577
7.58M
    uc = max;
578
7.58M
579
7.58M
    if (
carry7.58M
)
{46.0k
580
46.0k
      if (!s_pad(c, max + 1))
581
0
  return MP_MEMORY;
582
46.0k
583
46.0k
      c->digits[max] = carry;
584
46.0k
      ++uc;
585
7.58M
    }
586
7.58M
587
7.58M
    
MP_USED7.58M
(c) = uc;7.58M
588
7.58M
    
MP_SIGN7.58M
(c) =
MP_SIGN7.58M
(a);
589
7.58M
590
12.8M
  }
591
12.8M
  else {
592
5.23M
    /* Different signs -- subtract magnitudes, preserve sign of greater */
593
5.23M
    mp_int  x, y;
594
5.23M
    int     cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
595
5.23M
596
5.23M
    /* Set x to max(a, b), y to min(a, b) to simplify later code.
597
5.23M
       A special case yields zero for equal magnitudes.
598
5.23M
    */
599
5.23M
    if (
cmp == 05.23M
)
{786k
600
786k
      mp_int_zero(c);
601
786k
      return MP_OK;
602
5.23M
    }
603
4.44M
    else 
if (4.44M
cmp < 04.44M
)
{1.28M
604
1.28M
      x = b; y = a;
605
4.44M
    }
606
4.44M
    else {
607
3.15M
      x = a; y = b;
608
5.23M
    }
609
5.23M
610
4.44M
    
if (4.44M
!s_pad(c, 4.44M
MP_USED4.44M
(x)))
611
0
      return MP_MEMORY;
612
4.44M
613
4.44M
    /* Subtract smaller from larger */
614
4.44M
    
s_usub(4.44M
MP_DIGITS4.44M
(x),
MP_DIGITS4.44M
(y),
MP_DIGITS4.44M
(c),
MP_USED4.44M
(x),
MP_USED4.44M
(y));
615
4.44M
    
MP_USED4.44M
(c) =
MP_USED4.44M
(x);
616
4.44M
    CLAMP(c);
617
4.44M
618
4.44M
    /* Give result the sign of the larger */
619
4.44M
    
MP_SIGN4.44M
(c) =
MP_SIGN4.44M
(x);
620
12.8M
  }
621
12.8M
622
12.0M
  return MP_OK;
623
12.8M
}
624
625
mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c)
626
1.26k
{
627
1.26k
  mpz_t    vtmp;
628
1.26k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
629
1.26k
630
1.26k
  s_fake(&vtmp, value, vbuf);
631
1.26k
632
1.26k
  return mp_int_add(a, &vtmp, c);
633
1.26k
}
634
635
mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
636
74.5M
{
637
74.5M
  mp_size ua, ub, uc, max;
638
74.5M
639
74.5M
  CHECK(a != NULL && b != NULL && c != NULL);
640
74.5M
641
74.5M
  ua = 
MP_USED74.5M
(a); ub =
MP_USED74.5M
(b); uc =
MP_USED74.5M
(c);
642
74.5M
  max = MAX(ua, ub);
643
74.5M
644
74.5M
  if (
MP_SIGN74.5M
(a) != 74.5M
MP_SIGN74.5M
(b))
{264k
645
264k
    /* Different signs -- add magnitudes and keep sign of a */
646
264k
    mp_digit carry;
647
264k
648
264k
    if (!s_pad(c, max))
649
0
      return MP_MEMORY;
650
264k
651
264k
    
carry = s_uadd(264k
MP_DIGITS264k
(a),
MP_DIGITS264k
(b),
MP_DIGITS264k
(c), ua, ub);
652
264k
    uc = max;
653
264k
654
264k
    if (
carry264k
)
{2.13k
655
2.13k
      if (!s_pad(c, max + 1))
656
0
  return MP_MEMORY;
657
2.13k
658
2.13k
      c->digits[max] = carry;
659
2.13k
      ++uc;
660
264k
    }
661
264k
662
264k
    
MP_USED264k
(c) = uc;264k
663
264k
    
MP_SIGN264k
(c) =
MP_SIGN264k
(a);
664
264k
665
74.5M
  }
666
74.5M
  else {
667
74.2M
    /* Same signs -- subtract magnitudes */
668
74.2M
    mp_int  x, y;
669
74.2M
    mp_sign osign;
670
74.2M
    int     cmp = s_ucmp(a, b);
671
74.2M
672
74.2M
    if (!s_pad(c, max))
673
0
      return MP_MEMORY;
674
74.2M
675
74.2M
    
if (74.2M
cmp >= 074.2M
)
{12.5M
676
12.5M
      x = a; y = b; osign = MP_ZPOS;
677
74.2M
    }
678
74.2M
    else {
679
61.7M
      x = b; y = a; osign = MP_NEG;
680
74.2M
    }
681
74.2M
682
74.2M
    if (
MP_SIGN74.2M
(a) == MP_NEG && 74.2M
cmp != 090.9k
)
683
74.1k
      osign = 1 - osign;
684
74.2M
685
74.2M
    s_usub(
MP_DIGITS74.2M
(x),
MP_DIGITS74.2M
(y),
MP_DIGITS74.2M
(c),
MP_USED74.2M
(x),
MP_USED74.2M
(y));
686
74.2M
    
MP_USED74.2M
(c) =
MP_USED74.2M
(x);
687
74.2M
    CLAMP(c);
688
74.2M
689
74.2M
    MP_SIGN(c) = osign;
690
74.5M
  }
691
74.5M
692
74.5M
  return MP_OK;
693
74.5M
}
694
695
mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c)
696
14.1k
{
697
14.1k
  mpz_t    vtmp;
698
14.1k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
699
14.1k
700
14.1k
  s_fake(&vtmp, value, vbuf);
701
14.1k
702
14.1k
  return mp_int_sub(a, &vtmp, c);
703
14.1k
}
704
705
mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
706
24.3M
{
707
24.3M
  mp_digit *out;
708
24.3M
  mp_size   osize, ua, ub, p = 0;
709
24.3M
  mp_sign   osign;
710
24.3M
711
24.3M
  CHECK(a != NULL && b != NULL && c != NULL);
712
24.3M
713
24.3M
  /* If either input is zero, we can shortcut multiplication */
714
24.3M
  if (
mp_int_compare_zero(a) == 0 || 24.3M
mp_int_compare_zero(b) == 019.4M
)
{17.1M
715
17.1M
    mp_int_zero(c);
716
17.1M
    return MP_OK;
717
24.3M
  }
718
24.3M
719
24.3M
  /* Output is positive if inputs have same sign, otherwise negative */
720
7.18M
  
osign = (7.18M
MP_SIGN7.18M
(a) ==
MP_SIGN7.18M
(b)) ?
MP_ZPOS4.28M
:
MP_NEG2.89M
;
721
7.18M
722
7.18M
  /* If the output is not identical to any of the inputs, we'll write the
723
7.18M
     results directly; otherwise, allocate a temporary space. */
724
7.18M
  ua = 
MP_USED7.18M
(a); ub =
MP_USED7.18M
(b);
725
7.18M
  osize = MAX(ua, ub);
726
7.18M
  osize = 4 * ((osize + 1) / 2);
727
7.18M
728
7.18M
  if (
c == a || 7.18M
c == b5.70M
)
{1.47M
729
1.47M
    p = ROUND_PREC(osize);
730
1.47M
    p = MAX(p, default_precision);
731
1.47M
732
1.47M
    if ((out = s_alloc(p)) == NULL)
733
0
      return MP_MEMORY;
734
7.18M
  }
735
7.18M
  else {
736
5.70M
    if (!s_pad(c, osize))
737
0
      return MP_MEMORY;
738
5.70M
739
5.70M
    
out = 5.70M
MP_DIGITS5.70M
(c);
740
7.18M
  }
741
7.18M
  
ZERO7.18M
(out, osize);7.18M
742
7.18M
743
7.18M
  if (
!s_kmul(7.18M
MP_DIGITS7.18M
(a),
MP_DIGITS7.18M
(b), out, ua, ub))
744
0
    return MP_MEMORY;
745
7.18M
746
7.18M
  /* If we allocated a new buffer, get rid of whatever memory c was already
747
7.18M
     using, and fix up its fields to reflect that.
748
7.18M
   */
749
7.18M
  
if (7.18M
out != 7.18M
MP_DIGITS7.18M
(c))
{1.47M
750
1.47M
    if (
(void *) 1.47M
MP_DIGITS1.47M
(c) != (void *) c)
751
1.27M
      
s_free(1.27M
MP_DIGITS1.27M
(c));
752
1.47M
    MP_DIGITS(c) = out;
753
1.47M
    MP_ALLOC(c) = p;
754
7.18M
  }
755
7.18M
756
7.18M
  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
757
7.18M
  CLAMP(c);           /* ... right here */
758
7.18M
  MP_SIGN(c) = osign;
759
7.18M
760
7.18M
  return MP_OK;
761
24.3M
}
762
763
mp_result mp_int_mul_value(mp_int a, mp_small value, mp_int c)
764
0
{
765
0
  mpz_t    vtmp;
766
0
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
767
0
768
0
  s_fake(&vtmp, value, vbuf);
769
0
770
0
  return mp_int_mul(a, &vtmp, c);
771
0
}
772
773
mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c)
774
3.94k
{
775
3.94k
  mp_result res;
776
3.94k
  CHECK(a != NULL && c != NULL && p2 >= 0);
777
3.94k
778
3.94k
  if ((res = mp_int_copy(a, c)) != MP_OK)
779
0
    return res;
780
3.94k
781
3.94k
  
if (3.94k
s_qmul(c, (mp_size) p2)3.94k
)
782
3.94k
    return MP_OK;
783
3.94k
  else
784
0
    return MP_MEMORY;
785
3.94k
}
786
787
mp_result mp_int_sqr(mp_int a, mp_int c)
788
0
{
789
0
  mp_digit *out;
790
0
  mp_size   osize, p = 0;
791
0
792
0
  CHECK(a != NULL && c != NULL);
793
0
794
0
  /* Get a temporary buffer big enough to hold the result */
795
0
  osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2);
796
0
  if (
a == c0
)
{0
797
0
    p = ROUND_PREC(osize);
798
0
    p = MAX(p, default_precision);
799
0
800
0
    if ((out = s_alloc(p)) == NULL)
801
0
      return MP_MEMORY;
802
0
  }
803
0
  else {
804
0
    if (!s_pad(c, osize))
805
0
      return MP_MEMORY;
806
0
807
0
    
out = 0
MP_DIGITS0
(c);
808
0
  }
809
0
  
ZERO0
(out, osize);0
810
0
811
0
  s_ksqr(
MP_DIGITS0
(a), out,
MP_USED0
(a));
812
0
813
0
  /* Get rid of whatever memory c was already using, and fix up its fields to
814
0
     reflect the new digit array it's using
815
0
   */
816
0
  if (
out != 0
MP_DIGITS0
(c))
{0
817
0
    if (
(void *) 0
MP_DIGITS0
(c) != (void *) c)
818
0
      
s_free(0
MP_DIGITS0
(c));
819
0
    MP_DIGITS(c) = out;
820
0
    MP_ALLOC(c) = p;
821
0
  }
822
0
823
0
  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
824
0
  CLAMP(c);           /* ... right here */
825
0
  MP_SIGN(c) = MP_ZPOS;
826
0
827
0
  return MP_OK;
828
0
}
829
830
mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r)
831
3.91M
{
832
3.91M
  int cmp, lg;
833
3.91M
  mp_result res = MP_OK;
834
3.91M
  mp_int qout, rout;
835
3.91M
  mp_sign sa = 
MP_SIGN3.91M
(a), sb =
MP_SIGN3.91M
(b);
836
3.91M
  DECLARE_TEMP(2);
837
3.91M
838
3.91M
  CHECK(a != NULL && b != NULL && q != r);
839
3.91M
840
3.91M
  if (
CMPZ3.91M
(b) == 03.91M
)
841
0
    return MP_UNDEF;
842
3.91M
  else 
if (3.91M
(cmp = s_ucmp(a, b)) < 03.91M
)
{735k
843
735k
    /* If |a| < |b|, no division is required:
844
735k
       q = 0, r = a
845
735k
     */
846
735k
    if (
r && 735k
(res = mp_int_copy(a, r)) != MP_OK49.8k
)
847
0
      return res;
848
735k
849
735k
    
if (735k
q735k
)
850
726k
      mp_int_zero(q);
851
735k
852
735k
    return MP_OK;
853
3.91M
  }
854
3.17M
  else 
if (3.17M
cmp == 03.17M
)
{718k
855
718k
    /* If |a| = |b|, no division is required:
856
718k
       q = 1 or -1, r = 0
857
718k
     */
858
718k
    if (r)
859
4.77k
      mp_int_zero(r);
860
718k
861
718k
    if (
q718k
)
{715k
862
715k
      mp_int_zero(q);
863
715k
      q->digits[0] = 1;
864
715k
865
715k
      if (sa != sb)
866
84.1k
  
MP_SIGN84.1k
(q) = MP_NEG84.1k
;
867
718k
    }
868
718k
869
718k
    return MP_OK;
870
3.91M
  }
871
3.91M
872
3.91M
  /* When |a| > |b|, real division is required.  We need someplace to store
873
3.91M
     quotient and remainder, but q and r are allowed to be NULL or to overlap
874
3.91M
     with the inputs.
875
3.91M
   */
876
2.45M
  
if (2.45M
(lg = s_isp2(b)) < 02.45M
)
{1.64M
877
1.64M
    if (
q && 1.64M
b != q1.63M
)
{1.53M
878
1.53M
      if ((res = mp_int_copy(a, q)) != MP_OK)
879
0
  goto CLEANUP;
880
1.53M
      else
881
1.53M
  qout = q;
882
1.64M
    }
883
1.64M
    else {
884
116k
      qout = LAST_TEMP();
885
116k
      SETUP(mp_int_init_copy(LAST_TEMP(), a));
886
1.64M
    }
887
1.64M
888
1.64M
    
if (1.64M
r && 1.64M
a != r120k
)
{120k
889
120k
      if ((res = mp_int_copy(b, r)) != MP_OK)
890
0
  goto CLEANUP;
891
120k
      else
892
120k
  rout = r;
893
1.64M
    }
894
1.64M
    else {
895
1.52M
      rout = LAST_TEMP();
896
1.52M
      SETUP(mp_int_init_copy(LAST_TEMP(), b));
897
1.64M
    }
898
1.64M
899
1.64M
    
if (1.64M
(res = s_udiv_knuth(qout, rout)) != MP_OK1.64M
)
goto CLEANUP0
;
900
2.45M
  }
901
2.45M
  else {
902
808k
    if (
q && 808k
(res = mp_int_copy(a, q)) != MP_OK764k
)
goto CLEANUP0
;
903
808k
    
if (808k
r && 808k
(res = mp_int_copy(a, r)) != MP_OK126k
)
goto CLEANUP0
;
904
808k
905
808k
    
if (808k
q808k
)
s_qdiv(q, (mp_size) lg)764k
; qout = q;
906
808k
    if (
r808k
)
s_qmod(r, (mp_size) lg)126k
; rout = r;
907
2.45M
  }
908
2.45M
909
2.45M
  /* Recompute signs for output */
910
2.45M
  
if (2.45M
rout2.45M
)
{1.77M
911
1.77M
    MP_SIGN(rout) = sa;
912
1.77M
    if (
CMPZ1.77M
(rout) == 01.77M
)
913
1.69M
      
MP_SIGN1.69M
(rout) = MP_ZPOS1.69M
;
914
2.45M
  }
915
2.45M
  if (
qout2.45M
)
{2.41M
916
2.41M
    
MP_SIGN2.41M
(qout) = (sa == sb) ?
MP_ZPOS1.69M
:
MP_NEG719k
;
917
2.41M
    if (
CMPZ2.41M
(qout) == 02.41M
)
918
0
      
MP_SIGN0
(qout) = MP_ZPOS0
;
919
2.45M
  }
920
2.45M
921
2.45M
  if (
q && 2.45M
(res = mp_int_copy(qout, q)) != MP_OK2.39M
)
goto CLEANUP0
;
922
2.45M
  
if (2.45M
r && 2.45M
(res = mp_int_copy(rout, r)) != MP_OK247k
)
goto CLEANUP0
;
923
2.45M
924
2.45M
  
CLEANUP_TEMP2.45M
();2.45M
925
2.45M
  return res;
926
3.91M
}
927
928
mp_result mp_int_mod(mp_int a, mp_int m, mp_int c)
929
0
{
930
0
  mp_result res;
931
0
  mpz_t     tmp;
932
0
  mp_int    out;
933
0
934
0
  if (
m == c0
)
{0
935
0
    mp_int_init(&tmp);
936
0
    out = &tmp;
937
0
  }
938
0
  else {
939
0
    out = c;
940
0
  }
941
0
942
0
  if ((res = mp_int_div(a, m, NULL, out)) != MP_OK)
943
0
    goto CLEANUP;
944
0
945
0
  
if (0
CMPZ0
(out) < 00
)
946
0
    res = mp_int_add(out, m, c);
947
0
  else
948
0
    res = mp_int_copy(out, c);
949
0
950
0
 CLEANUP:
951
0
  if (out != c)
952
0
    mp_int_clear(&tmp);
953
0
954
0
  return res;
955
0
}
956
957
mp_result mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r)
958
61.6k
{
959
61.6k
  mpz_t     vtmp, rtmp;
960
61.6k
  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
961
61.6k
  mp_result res;
962
61.6k
963
61.6k
  mp_int_init(&rtmp);
964
61.6k
  s_fake(&vtmp, value, vbuf);
965
61.6k
966
61.6k
  if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
967
0
    goto CLEANUP;
968
61.6k
969
61.6k
  
if (61.6k
r61.6k
)
970
50.7k
    (void) mp_int_to_int(&rtmp, r); /* can't fail */
971
61.6k
972
61.6k
 CLEANUP:
973
61.6k
  mp_int_clear(&rtmp);
974
61.6k
  return res;
975
61.6k
}
976
977
mp_result mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r)
978
0
{
979
0
  mp_result res = MP_OK;
980
0
981
0
  CHECK(a != NULL && p2 >= 0 && q != r);
982
0
983
0
  if (
q != NULL && 0
(res = mp_int_copy(a, q)) == MP_OK0
)
984
0
    s_qdiv(q, (mp_size) p2);
985
0
986
0
  if (
res == MP_OK && 0
r != NULL0
&&
(res = mp_int_copy(a, r)) == MP_OK0
)
987
0
    s_qmod(r, (mp_size) p2);
988
0
989
0
  return res;
990
0
}
991
992
mp_result mp_int_expt(mp_int a, mp_small b, mp_int c)
993
0
{
994
0
  mpz_t t;
995
0
  mp_result res;
996
0
  unsigned int v = labs(b);
997
0
998
0
  CHECK(c != NULL);
999
0
  if (b < 0)
1000
0
    return MP_RANGE;
1001
0
1002
0
  
if (0
(res = mp_int_init_copy(&t, a)) != MP_OK0
)
1003
0
    return res;
1004
0
1005
0
  (void) mp_int_set_value(c, 1);
1006
0
  while (
v != 00
)
{0
1007
0
    if (
v & 10
)
{0
1008
0
      if ((res = mp_int_mul(c, &t, c)) != MP_OK)
1009
0
  goto CLEANUP;
1010
0
    }
1011
0
1012
0
    v >>= 1;
1013
0
    if (
v == 00
)
break0
;
1014
0
1015
0
    
if (0
(res = mp_int_sqr(&t, &t)) != MP_OK0
)
1016
0
      goto CLEANUP;
1017
0
  }
1018
0
1019
0
 CLEANUP:
1020
0
  mp_int_clear(&t);
1021
0
  return res;
1022
0
}
1023
1024
mp_result mp_int_expt_value(mp_small a, mp_small b, mp_int c)
1025
0
{
1026
0
  mpz_t     t;
1027
0
  mp_result res;
1028
0
  unsigned int v = labs(b);
1029
0
1030
0
  CHECK(c != NULL);
1031
0
  if (b < 0)
1032
0
    return MP_RANGE;
1033
0
1034
0
  
if (0
(res = mp_int_init_value(&t, a)) != MP_OK0
)
1035
0
    return res;
1036
0
1037
0
  (void) mp_int_set_value(c, 1);
1038
0
  while (
v != 00
)
{0
1039
0
    if (
v & 10
)
{0
1040
0
      if ((res = mp_int_mul(c, &t, c)) != MP_OK)
1041
0
  goto CLEANUP;
1042
0
    }
1043
0
1044
0
    v >>= 1;
1045
0
    if (
v == 00
)
break0
;
1046
0
1047
0
    
if (0
(res = mp_int_sqr(&t, &t)) != MP_OK0
)
1048
0
      goto CLEANUP;
1049
0
  }
1050
0
1051
0
 CLEANUP:
1052
0
  mp_int_clear(&t);
1053
0
  return res;
1054
0
}
1055
1056
mp_result mp_int_expt_full(mp_int a, mp_int b, mp_int c)
1057
0
{
1058
0
  mpz_t t;
1059
0
  mp_result res;
1060
0
  unsigned ix, jx;
1061
0
1062
0
  CHECK(a != NULL && b != NULL && c != NULL);
1063
0
  if (
MP_SIGN0
(b) == MP_NEG0
)
1064
0
    return MP_RANGE;
1065
0
1066
0
  
if (0
(res = mp_int_init_copy(&t, a)) != MP_OK0
)
1067
0
    return res;
1068
0
1069
0
  (void) mp_int_set_value(c, 1);
1070
0
  for (ix = 0; 
ix < 0
MP_USED0
(b);
++ix0
)
{0
1071
0
    mp_digit d = b->digits[ix];
1072
0
1073
0
    for (jx = 0; 
jx < 0
MP_DIGIT_BIT0
;
++jx0
)
{0
1074
0
      if (
d & 10
)
{0
1075
0
  if ((res = mp_int_mul(c, &t, c)) != MP_OK)
1076
0
    goto CLEANUP;
1077
0
      }
1078
0
1079
0
      d >>= 1;
1080
0
      if (
d == 0 && 0
ix + 1 == 0
MP_USED0
(b))
1081
0
  break;
1082
0
      
if (0
(res = mp_int_sqr(&t, &t)) != MP_OK0
)
1083
0
  goto CLEANUP;
1084
0
    }
1085
0
  }
1086
0
1087
0
 CLEANUP:
1088
0
  mp_int_clear(&t);
1089
0
  return res;
1090
0
}
1091
1092
int       mp_int_compare(mp_int a, mp_int b)
1093
261k
{
1094
261k
  mp_sign sa;
1095
261k
1096
261k
  CHECK(a != NULL && b != NULL);
1097
261k
1098
261k
  sa = MP_SIGN(a);
1099
261k
  if (
sa == 261k
MP_SIGN261k
(b))
{228k
1100
228k
    int cmp = s_ucmp(a, b);
1101
228k
1102
228k
    /* If they're both zero or positive, the normal comparison applies; if both
1103
228k
       negative, the sense is reversed. */
1104
228k
    if (sa == MP_ZPOS)
1105
162k
      return cmp;
1106
228k
    else
1107
66.0k
      return -cmp;
1108
228k
1109
261k
  }
1110
261k
  else {
1111
33.4k
    if (sa == MP_ZPOS)
1112
15.1k
      return 1;
1113
33.4k
    else
1114
18.2k
      return -1;
1115
261k
  }
1116
261k
}
1117
1118
int       mp_int_compare_unsigned(mp_int a, mp_int b)
1119
2.83M
{
1120
2.83M
  NRCHECK(a != NULL && b != NULL);
1121
2.83M
1122
2.83M
  return s_ucmp(a, b);
1123
2.83M
}
1124
1125
int       mp_int_compare_zero(mp_int z)
1126
60.9M
{
1127
60.9M
  NRCHECK(z != NULL);
1128
60.9M
1129
60.9M
  if (
MP_USED60.9M
(z) == 1 && 60.9M
z->digits[0] == 034.3M
)
1130
17.3M
    return 0;
1131
43.5M
  else 
if (43.5M
MP_SIGN43.5M
(z) == MP_ZPOS43.5M
)
1132
30.7M
    return 1;
1133
43.5M
  else
1134
12.8M
    return -1;
1135
60.9M
}
1136
1137
int       mp_int_compare_value(mp_int z, mp_small value)
1138
76.5M
{
1139
76.5M
  mp_sign vsign = (value < 0) ? 
MP_NEG37.6M
:
MP_ZPOS38.9M
;
1140
76.5M
  int cmp;
1141
76.5M
1142
76.5M
  CHECK(z != NULL);
1143
76.5M
1144
76.5M
  if (
vsign == 76.5M
MP_SIGN76.5M
(z))
{46.6M
1145
46.6M
    cmp = s_vcmp(z, value);
1146
46.6M
1147
46.6M
    return (vsign == MP_ZPOS) ? 
cmp38.8M
:
-cmp7.77M
;
1148
76.5M
  }
1149
76.5M
  else {
1150
29.8M
    return (value < 0) ? 
129.8M
:
-132.6k
;
1151
76.5M
  }
1152
76.5M
}
1153
1154
int       mp_int_compare_uvalue(mp_int z, mp_usmall uv)
1155
0
{
1156
0
  CHECK(z != NULL);
1157
0
1158
0
  if (
MP_SIGN0
(z) == MP_NEG0
)
1159
0
    return -1;
1160
0
  else
1161
0
    return s_uvcmp(z, uv);
1162
0
}
1163
1164
mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c)
1165
0
{
1166
0
  mp_result res;
1167
0
  mp_size um;
1168
0
  mp_int s;
1169
0
  DECLARE_TEMP(3);
1170
0
1171
0
  CHECK(a != NULL && b != NULL && c != NULL && m != NULL);
1172
0
1173
0
  /* Zero moduli and negative exponents are not considered. */
1174
0
  if (
CMPZ0
(m) == 00
)
1175
0
    return MP_UNDEF;
1176
0
  
if (0
CMPZ0
(b) < 00
)
1177
0
    return MP_RANGE;
1178
0
1179
0
  
um = 0
MP_USED0
(m);
1180
0
  SETUP(mp_int_init_size(TEMP(0), 2 * um));
1181
0
  
SETUP0
(mp_int_init_size(TEMP(1), 2 * um));0
1182
0
1183
0
  
if (0
c == b || 0
c == m0
)
{0
1184
0
    SETUP(mp_int_init_size(TEMP(2), 2 * um));
1185
0
    
s = 0
TEMP0
(2);
1186
0
  }
1187
0
  else {
1188
0
    s = c;
1189
0
  }
1190
0
1191
0
  
if (0
(res = mp_int_mod(a, m, 0
TEMP0
(0))) != MP_OK)
goto CLEANUP0
;
1192
0
1193
0
  
if (0
(res = s_brmu(0
TEMP0
(1), m)) != MP_OK)
goto CLEANUP0
;
1194
0
1195
0
  
if (0
(res = s_embar(0
TEMP0
(0), b, m,
TEMP0
(1), s)) != MP_OK)
1196
0
    goto CLEANUP;
1197
0
1198
0
  res = mp_int_copy(s, c);
1199
0
1200
0
  
CLEANUP_TEMP0
();0
1201
0
  return res;
1202
0
}
1203
1204
mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c)
1205
0
{
1206
0
  mpz_t vtmp;
1207
0
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1208
0
1209
0
  s_fake(&vtmp, value, vbuf);
1210
0
1211
0
  return mp_int_exptmod(a, &vtmp, m, c);
1212
0
}
1213
1214
mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b,
1215
        mp_int m, mp_int c)
1216
0
{
1217
0
  mpz_t vtmp;
1218
0
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1219
0
1220
0
  s_fake(&vtmp, value, vbuf);
1221
0
1222
0
  return mp_int_exptmod(&vtmp, b, m, c);
1223
0
}
1224
1225
mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
1226
0
{
1227
0
  mp_result res;
1228
0
  mp_size um;
1229
0
  mp_int s;
1230
0
  DECLARE_TEMP(2);
1231
0
1232
0
  CHECK(a && b && m && c);
1233
0
1234
0
  /* Zero moduli and negative exponents are not considered. */
1235
0
  if (
CMPZ0
(m) == 00
)
1236
0
    return MP_UNDEF;
1237
0
  
if (0
CMPZ0
(b) < 00
)
1238
0
    return MP_RANGE;
1239
0
1240
0
  
um = 0
MP_USED0
(m);
1241
0
  SETUP(mp_int_init_size(TEMP(0), 2 * um));
1242
0
1243
0
  
if (0
c == b || 0
c == m0
)
{0
1244
0
    SETUP(mp_int_init_size(TEMP(1), 2 * um));
1245
0
    
s = 0
TEMP0
(1);
1246
0
  }
1247
0
  else {
1248
0
    s = c;
1249
0
  }
1250
0
1251
0
  
if (0
(res = mp_int_mod(a, m, 0
TEMP0
(0))) != MP_OK)
goto CLEANUP0
;
1252
0
1253
0
  
if (0
(res = s_embar(0
TEMP0
(0), b, m, mu, s)) != MP_OK)
1254
0
    goto CLEANUP;
1255
0
1256
0
  res = mp_int_copy(s, c);
1257
0
1258
0
  
CLEANUP_TEMP0
();0
1259
0
  return res;
1260
0
}
1261
1262
mp_result mp_int_redux_const(mp_int m, mp_int c)
1263
0
{
1264
0
  CHECK(m != NULL && c != NULL && m != c);
1265
0
1266
0
  return s_brmu(c, m);
1267
0
}
1268
1269
mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c)
1270
0
{
1271
0
  mp_result res;
1272
0
  mp_sign sa;
1273
0
  DECLARE_TEMP(2);
1274
0
1275
0
  CHECK(a != NULL && m != NULL && c != NULL);
1276
0
1277
0
  if (
CMPZ0
(a) == 0 || 0
CMPZ0
(m) <= 00
)
1278
0
    return MP_RANGE;
1279
0
1280
0
  
sa = 0
MP_SIGN0
(a); /* need this for the result later */
1281
0
1282
0
  for (last__ = 0; 
last__ < 20
;
++last__0
)
1283
0
    
mp_int_init(0
LAST_TEMP0
());
1284
0
1285
0
  if (
(res = mp_int_egcd(a, m, 0
TEMP0
(0),
TEMP0
(1), NULL)) != MP_OK)
1286
0
    goto CLEANUP;
1287
0
1288
0
  
if (0
mp_int_compare_value(0
TEMP0
(0), 1) != 0)
{0
1289
0
    res = MP_UNDEF;
1290
0
    goto CLEANUP;
1291
0
  }
1292
0
1293
0
  /* It is first necessary to constrain the value to the proper range */
1294
0
  
if (0
(res = mp_int_mod(0
TEMP0
(1), m,
TEMP0
(1))) != MP_OK)
1295
0
    goto CLEANUP;
1296
0
1297
0
  /* Now, if 'a' was originally negative, the value we have is actually the
1298
0
     magnitude of the negative representative; to get the positive value we
1299
0
     have to subtract from the modulus.  Otherwise, the value is okay as it
1300
0
     stands.
1301
0
   */
1302
0
  
if (0
sa == MP_NEG0
)
1303
0
    
res = mp_int_sub(m, 0
TEMP0
(1), c);
1304
0
  else
1305
0
    
res = mp_int_copy(0
TEMP0
(1), c);
1306
0
1307
0
  
CLEANUP_TEMP0
();0
1308
0
  return res;
1309
0
}
1310
1311
/* Binary GCD algorithm due to Josef Stein, 1961 */
1312
mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c)
1313
2.51M
{
1314
2.51M
  int ca, cb, k = 0;
1315
2.51M
  mpz_t u, v, t;
1316
2.51M
  mp_result res;
1317
2.51M
1318
2.51M
  CHECK(a != NULL && b != NULL && c != NULL);
1319
2.51M
1320
2.51M
  ca = CMPZ(a);
1321
2.51M
  cb = CMPZ(b);
1322
2.51M
  if (
ca == 0 && 2.51M
cb == 021
)
1323
0
    return MP_UNDEF;
1324
2.51M
  else 
if (2.51M
ca == 02.51M
)
1325
21
    return mp_int_abs(b, c);
1326
2.51M
  else 
if (2.51M
cb == 02.51M
)
1327
9.43k
    return mp_int_abs(a, c);
1328
2.51M
1329
2.51M
  mp_int_init(&t);
1330
2.51M
  if ((res = mp_int_init_copy(&u, a)) != MP_OK)
1331
0
    goto U;
1332
2.51M
  
if (2.51M
(res = mp_int_init_copy(&v, b)) != MP_OK2.51M
)
1333
0
    goto V;
1334
2.51M
1335
2.51M
  
MP_SIGN2.51M
(&u) = MP_ZPOS; 2.51M
MP_SIGN2.51M
(&v) = MP_ZPOS;
1336
2.51M
1337
2.51M
  { /* Divide out common factors of 2 from u and v */
1338
2.51M
    int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
1339
2.51M
1340
2.51M
    k = MIN(div2_u, div2_v);
1341
2.51M
    s_qdiv(&u, (mp_size) k);
1342
2.51M
    s_qdiv(&v, (mp_size) k);
1343
2.51M
  }
1344
2.51M
1345
2.51M
  if (
mp_int_is_odd2.51M
(&u))
{2.00M
1346
2.00M
    if ((res = mp_int_neg(&v, &t)) != MP_OK)
1347
0
      goto CLEANUP;
1348
2.51M
  }
1349
2.51M
  else {
1350
510k
    if ((res = mp_int_copy(&u, &t)) != MP_OK)
1351
0
      goto CLEANUP;
1352
2.51M
  }
1353
2.51M
1354
2.51M
  
for (;;) 2.51M
{73.9M
1355
73.9M
    s_qdiv(&t, s_dp2k(&t));
1356
73.9M
1357
73.9M
    if (
CMPZ73.9M
(&t) > 073.9M
)
{10.3M
1358
10.3M
      if ((res = mp_int_copy(&t, &u)) != MP_OK)
1359
0
  goto CLEANUP;
1360
73.9M
    }
1361
73.9M
    else {
1362
63.6M
      if ((res = mp_int_neg(&t, &v)) != MP_OK)
1363
0
  goto CLEANUP;
1364
73.9M
    }
1365
73.9M
1366
73.9M
    
if (73.9M
(res = mp_int_sub(&u, &v, &t)) != MP_OK73.9M
)
1367
0
      goto CLEANUP;
1368
73.9M
1369
73.9M
    
if (73.9M
CMPZ73.9M
(&t) == 073.9M
)
1370
2.51M
      break;
1371
73.9M
  }
1372
2.51M
1373
2.51M
  
if (2.51M
(res = mp_int_abs(&u, c)) != MP_OK2.51M
)
1374
0
    goto CLEANUP;
1375
2.51M
  
if (2.51M
!s_qmul(c, (mp_size) k)2.51M
)
1376
0
    res = MP_MEMORY;
1377
2.51M
1378
2.51M
 CLEANUP:
1379
2.51M
  mp_int_clear(&v);
1380
2.51M
 V: mp_int_clear(&u);
1381
2.51M
 U: mp_int_clear(&t);
1382
2.51M
1383
2.51M
  return res;
1384
2.51M
}
1385
1386
/* This is the binary GCD algorithm again, but this time we keep track of the
1387
   elementary matrix operations as we go, so we can get values x and y
1388
   satisfying c = ax + by.
1389
 */
1390
mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c,
1391
          mp_int x, mp_int y)
1392
0
{
1393
0
  int k, ca, cb;
1394
0
  mp_result res;
1395
0
  DECLARE_TEMP(8);
1396
0
1397
0
  CHECK(a != NULL && b != NULL && c != NULL &&
1398
0
  (x != NULL || y != NULL));
1399
0
1400
0
  ca = CMPZ(a);
1401
0
  cb = CMPZ(b);
1402
0
  if (
ca == 0 && 0
cb == 00
)
1403
0
    return MP_UNDEF;
1404
0
  else 
if (0
ca == 00
)
{0
1405
0
    if (
(res = mp_int_abs(b, c)) != MP_OK0
)
return res0
;
1406
0
    mp_int_zero(x); (void) mp_int_set_value(y, 1); return MP_OK;
1407
0
  }
1408
0
  else 
if (0
cb == 00
)
{0
1409
0
    if (
(res = mp_int_abs(a, c)) != MP_OK0
)
return res0
;
1410
0
    (void) mp_int_set_value(x, 1); mp_int_zero(y); return MP_OK;
1411
0
  }
1412
0
1413
0
  /* Initialize temporaries:
1414
0
     A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */
1415
0
  
for (last__ = 0; 0
last__ < 40
;
++last__0
)
1416
0
    
mp_int_init(0
LAST_TEMP0
());
1417
0
  TEMP(0)->digits[0] = 1;
1418
0
  TEMP(3)->digits[0] = 1;
1419
0
1420
0
  SETUP(mp_int_init_copy(TEMP(4), a));
1421
0
  
SETUP0
(mp_int_init_copy(TEMP(5), b));0
1422
0
1423
0
  /* We will work with absolute values here */
1424
0
  
MP_SIGN0
(TEMP(4)) = MP_ZPOS;0
1425
0
  MP_SIGN(TEMP(5)) = MP_ZPOS;
1426
0
1427
0
  { /* Divide out common factors of 2 from u and v */
1428
0
    int  div2_u = s_dp2k(
TEMP0
(4)), div2_v = s_dp2k(
TEMP0
(5));
1429
0
1430
0
    k = MIN(div2_u, div2_v);
1431
0
    s_qdiv(TEMP(4), k);
1432
0
    s_qdiv(TEMP(5), k);
1433
0
  }
1434
0
1435
0
  SETUP(mp_int_init_copy(TEMP(6), TEMP(4)));
1436
0
  
SETUP0
(mp_int_init_copy(TEMP(7), TEMP(5)));0
1437
0
1438
0
  
for (;;) 0
{0
1439
0
    while (
mp_int_is_even0
(TEMP(4)))
{0
1440
0
      s_qdiv(TEMP(4), 1);
1441
0
1442
0
      if (
mp_int_is_odd0
(TEMP(0)) || 0
mp_int_is_odd0
(TEMP(1)))
{0
1443
0
  if (
(res = mp_int_add(0
TEMP0
(0),
TEMP0
(7),
TEMP0
(0))) != MP_OK)
1444
0
    goto CLEANUP;
1445
0
  
if (0
(res = mp_int_sub(0
TEMP0
(1),
TEMP0
(6),
TEMP0
(1))) != MP_OK)
1446
0
    goto CLEANUP;
1447
0
      }
1448
0
1449
0
      
s_qdiv(0
TEMP0
(0), 1);
1450
0
      s_qdiv(TEMP(1), 1);
1451
0
    }
1452
0
1453
0
    
while (0
mp_int_is_even0
(TEMP(5)))
{0
1454
0
      s_qdiv(TEMP(5), 1);
1455
0
1456
0
      if (
mp_int_is_odd0
(TEMP(2)) || 0
mp_int_is_odd0
(TEMP(3)))
{0
1457
0
  if (
(res = mp_int_add(0
TEMP0
(2),
TEMP0
(7),
TEMP0
(2))) != MP_OK)
1458
0
    goto CLEANUP;
1459
0
  
if (0
(res = mp_int_sub(0
TEMP0
(3),
TEMP0
(6),
TEMP0
(3))) != MP_OK)
1460
0
    goto CLEANUP;
1461
0
      }
1462
0
1463
0
      
s_qdiv(0
TEMP0
(2), 1);
1464
0
      s_qdiv(TEMP(3), 1);
1465
0
    }
1466
0
1467
0
    
if (0
mp_int_compare(0
TEMP0
(4),
TEMP0
(5)) >= 0)
{0
1468
0
      if (
(res = mp_int_sub(0
TEMP0
(4),
TEMP0
(5),
TEMP0
(4))) != MP_OK)
goto CLEANUP0
;
1469
0
      
if (0
(res = mp_int_sub(0
TEMP0
(0),
TEMP0
(2),
TEMP0
(0))) != MP_OK)
goto CLEANUP0
;
1470
0
      
if (0
(res = mp_int_sub(0
TEMP0
(1),
TEMP0
(3),
TEMP0
(1))) != MP_OK)
goto CLEANUP0
;
1471
0
    }
1472
0
    else {
1473
0
      if (
(res = mp_int_sub(0
TEMP0
(5),
TEMP0
(4),
TEMP0
(5))) != MP_OK)
goto CLEANUP0
;
1474
0
      
if (0
(res = mp_int_sub(0
TEMP0
(2),
TEMP0
(0),
TEMP0
(2))) != MP_OK)
goto CLEANUP0
;
1475
0
      
if (0
(res = mp_int_sub(0
TEMP0
(3),
TEMP0
(1),
TEMP0
(3))) != MP_OK)
goto CLEANUP0
;
1476
0
    }
1477
0
1478
0
    
if (0
CMPZ0
(TEMP(4)) == 00
)
{0
1479
0
      if (
x && 0
(res = mp_int_copy(0
TEMP0
(2), x)) != MP_OK)
goto CLEANUP0
;
1480
0
      
if (0
y && 0
(res = mp_int_copy(0
TEMP0
(3), y)) != MP_OK)
goto CLEANUP0
;
1481
0
      
if (0
c0
)
{0
1482
0
  if (
!s_qmul(0
TEMP0
(5), k))
{0
1483
0
    res = MP_MEMORY;
1484
0
    goto CLEANUP;
1485
0
  }
1486
0
1487
0
  
res = mp_int_copy(0
TEMP0
(5), c);
1488
0
      }
1489
0
1490
0
      break;
1491
0
    }
1492
0
  }
1493
0
1494
0
  
CLEANUP_TEMP0
();0
1495
0
  return res;
1496
0
}
1497
1498
mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c)
1499
270k
{
1500
270k
  mpz_t lcm;
1501
270k
  mp_result res;
1502
270k
1503
270k
  CHECK(a != NULL && b != NULL && c != NULL);
1504
270k
1505
270k
  /* Since a * b = gcd(a, b) * lcm(a, b), we can compute
1506
270k
     lcm(a, b) = (a / gcd(a, b)) * b.
1507
270k
1508
270k
     This formulation insures everything works even if the input
1509
270k
     variables share space.
1510
270k
   */
1511
270k
  if ((res = mp_int_init(&lcm)) != MP_OK)
1512
0
    return res;
1513
270k
  
if (270k
(res = mp_int_gcd(a, b, &lcm)) != MP_OK270k
)
1514
0
    goto CLEANUP;
1515
270k
  
if (270k
(res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK270k
)
1516
0
    goto CLEANUP;
1517
270k
  
if (270k
(res = mp_int_mul(&lcm, b, &lcm)) != MP_OK270k
)
1518
0
    goto CLEANUP;
1519
270k
1520
270k
  res = mp_int_copy(&lcm, c);
1521
270k
1522
270k
  CLEANUP:
1523
270k
    mp_int_clear(&lcm);
1524
270k
1525
270k
  return res;
1526
270k
}
1527
1528
int       mp_int_divisible_value(mp_int a, mp_small v)
1529
50.7k
{
1530
50.7k
  mp_small rem = 0;
1531
50.7k
1532
50.7k
  if (mp_int_div_value(a, v, NULL, &rem) != MP_OK)
1533
0
    return 0;
1534
50.7k
1535
50.7k
  return rem == 0;
1536
50.7k
}
1537
1538
int       mp_int_is_pow2(mp_int z)
1539
0
{
1540
0
  CHECK(z != NULL);
1541
0
1542
0
  return s_isp2(z);
1543
0
}
1544
1545
/* Implementation of Newton's root finding method, based loosely on a patch
1546
   contributed by Hal Finkel <half@halssoftware.com>
1547
   modified by M. J. Fromberger.
1548
 */
1549
mp_result mp_int_root(mp_int a, mp_small b, mp_int c)
1550
0
{
1551
0
  mp_result res = MP_OK;
1552
0
  int flips = 0;
1553
0
  DECLARE_TEMP(5);
1554
0
1555
0
  CHECK(a != NULL && c != NULL && b > 0);
1556
0
1557
0
  if (
b == 10
)
{0
1558
0
    return mp_int_copy(a, c);
1559
0
  }
1560
0
  
if (0
MP_SIGN0
(a) == MP_NEG0
)
{0
1561
0
    if (b % 2 == 0)
1562
0
      return MP_UNDEF; /* root does not exist for negative a with even b */
1563
0
    else
1564
0
      flips = 1;
1565
0
  }
1566
0
1567
0
  
SETUP0
(mp_int_init_copy(LAST_TEMP(), a));0
1568
0
  
SETUP0
(mp_int_init_copy(LAST_TEMP(), a));0
1569
0
  
SETUP0
(mp_int_init(LAST_TEMP()));0
1570
0
  
SETUP0
(mp_int_init(LAST_TEMP()));0
1571
0
  
SETUP0
(mp_int_init(LAST_TEMP()));0
1572
0
1573
0
  
(void) mp_int_abs(0
TEMP0
(0),
TEMP0
(0));
1574
0
  (void) mp_int_abs(
TEMP0
(1),
TEMP0
(1));
1575
0
1576
0
  for (;;) {
1577
0
    if (
(res = mp_int_expt(0
TEMP0
(1), b,
TEMP0
(2))) != MP_OK)
1578
0
      goto CLEANUP;
1579
0
1580
0
    
if (0
mp_int_compare_unsigned(0
TEMP0
(2),
TEMP0
(0)) <= 0)
1581
0
      break;
1582
0
1583
0
    
if (0
(res = mp_int_sub(0
TEMP0
(2),
TEMP0
(0),
TEMP0
(2))) != MP_OK)
1584
0
      goto CLEANUP;
1585
0
    
if (0
(res = mp_int_expt(0
TEMP0
(1), b - 1,
TEMP0
(3))) != MP_OK)
1586
0
      goto CLEANUP;
1587
0
    
if (0
(res = mp_int_mul_value(0
TEMP0
(3), b,
TEMP0
(3))) != MP_OK)
1588
0
      goto CLEANUP;
1589
0
    
if (0
(res = mp_int_div(0
TEMP0
(2),
TEMP0
(3),
TEMP0
(4), NULL)) != MP_OK)
1590
0
      goto CLEANUP;
1591
0
    
if (0
(res = mp_int_sub(0
TEMP0
(1),
TEMP0
(4),
TEMP0
(4))) != MP_OK)
1592
0
      goto CLEANUP;
1593
0
1594
0
    
if (0
mp_int_compare_unsigned(0
TEMP0
(1),
TEMP0
(4)) == 0)
{0
1595
0
      if (
(res = mp_int_sub_value(0
TEMP0
(4), 1,
TEMP0
(4))) != MP_OK)
1596
0
  goto CLEANUP;
1597
0
    }
1598
0
    
if (0
(res = mp_int_copy(0
TEMP0
(4),
TEMP0
(1))) != MP_OK)
1599
0
      goto CLEANUP;
1600
0
  }
1601
0
1602
0
  
if (0
(res = mp_int_copy(0
TEMP0
(1), c)) != MP_OK)
1603
0
    goto CLEANUP;
1604
0
1605
0
  /* If the original value of a was negative, flip the output sign. */
1606
0
  
if (0
flips0
)
1607
0
    (void) mp_int_neg(c, c); /* cannot fail */
1608
0
1609
0
  
CLEANUP_TEMP0
();0
1610
0
  return res;
1611
0
}
1612
1613
mp_result mp_int_to_int(mp_int z, mp_small *out)
1614
43.4M
{
1615
43.4M
  mp_usmall uv = 0;
1616
43.4M
  mp_size   uz;
1617
43.4M
  mp_digit *dz;
1618
43.4M
  mp_sign   sz;
1619
43.4M
1620
43.4M
  CHECK(z != NULL);
1621
43.4M
1622
43.4M
  /* Make sure the value is representable as a small integer */
1623
43.4M
  sz = MP_SIGN(z);
1624
43.4M
  if (
(sz == MP_ZPOS && 43.4M
mp_int_compare_value(z, 35.7M
MP_SMALL_MAX35.7M
) > 0) ||
1625
37.5M
      
mp_int_compare_value(z, 37.5M
MP_SMALL_MIN37.5M
) < 0)
1626
8.33M
    return MP_RANGE;
1627
43.4M
1628
35.1M
  
uz = 35.1M
MP_USED35.1M
(z);
1629
35.1M
  dz = MP_DIGITS(z) + uz - 1;
1630
35.1M
1631
78.0M
  while (
uz > 078.0M
)
{42.8M
1632
42.8M
    uv <<= MP_DIGIT_BIT/2;
1633
42.8M
    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1634
42.8M
    --uz;
1635
42.8M
  }
1636
35.1M
1637
35.1M
  if (out)
1638
35.1M
    
*out = (mp_small)((sz == MP_NEG) ? 35.1M
-uv5.33M
:
uv29.8M
);
1639
35.1M
1640
43.4M
  return MP_OK;
1641
43.4M
}
1642
1643
mp_result mp_int_to_uint(mp_int z, mp_usmall *out)
1644
0
{
1645
0
  mp_usmall uv = 0;
1646
0
  mp_size   uz;
1647
0
  mp_digit *dz;
1648
0
  mp_sign   sz;
1649
0
  
1650
0
  CHECK(z != NULL);
1651
0
1652
0
  /* Make sure the value is representable as an unsigned small integer */
1653
0
  sz = MP_SIGN(z);
1654
0
  if (
sz == MP_NEG || 0
mp_int_compare_uvalue(z, 0
MP_USMALL_MAX0
) > 0)
1655
0
    return MP_RANGE;
1656
0
     
1657
0
  
uz = 0
MP_USED0
(z);
1658
0
  dz = MP_DIGITS(z) + uz - 1;
1659
0
  
1660
0
  while (
uz > 00
)
{0
1661
0
    uv <<= MP_DIGIT_BIT/2;
1662
0
    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1663
0
    --uz;
1664
0
  }
1665
0
  
1666
0
  if (out)
1667
0
    *out = uv;
1668
0
  
1669
0
  return MP_OK;
1670
0
}
1671
1672
mp_result mp_int_to_string(mp_int z, mp_size radix,
1673
         char *str, int limit)
1674
1.44k
{
1675
1.44k
  mp_result res;
1676
1.44k
  int       cmp = 0;
1677
1.44k
1678
1.44k
  CHECK(z != NULL && str != NULL && limit >= 2);
1679
1.44k
1680
1.44k
  if (
radix < 1.44k
MP_MIN_RADIX1.44k
||
radix > 1.44k
MP_MAX_RADIX1.44k
)
1681
0
    return MP_RANGE;
1682
1.44k
1683
1.44k
  
if (1.44k
CMPZ1.44k
(z) == 01.44k
)
{0
1684
0
    *str++ = s_val2ch(0, 1);
1685
1.44k
  }
1686
1.44k
  else {
1687
1.44k
    mpz_t tmp;
1688
1.44k
    char  *h, *t;
1689
1.44k
1690
1.44k
    if ((res = mp_int_init_copy(&tmp, z)) != MP_OK)
1691
0
      return res;
1692
1.44k
1693
1.44k
    
if (1.44k
MP_SIGN1.44k
(z) == MP_NEG1.44k
)
{740
1694
740
      *str++ = '-';
1695
740
      --limit;
1696
1.44k
    }
1697
1.44k
    h = str;
1698
1.44k
1699
1.44k
    /* Generate digits in reverse order until finished or limit reached */
1700
27.7k
    for (/* */; 
limit > 027.7k
;
--limit26.3k
)
{27.7k
1701
27.7k
      mp_digit d;
1702
27.7k
1703
27.7k
      if (
(cmp = 27.7k
CMPZ27.7k
(&tmp)) == 0)
1704
1.44k
  break;
1705
27.7k
1706
27.7k
      d = s_ddiv(&tmp, (mp_digit)radix);
1707
26.3k
      *str++ = s_val2ch(d, 1);
1708
26.3k
    }
1709
1.44k
    t = str - 1;
1710
1.44k
1711
1.44k
    /* Put digits back in correct output order */
1712
14.0k
    while (
h < t14.0k
)
{12.5k
1713
12.5k
      char tc = *h;
1714
12.5k
      *h++ = *t;
1715
12.5k
      *t-- = tc;
1716
12.5k
    }
1717
1.44k
1718
1.44k
    mp_int_clear(&tmp);
1719
1.44k
  }
1720
1.44k
1721
1.44k
  *str = '\0';
1722
1.44k
  if (cmp == 0)
1723
1.44k
    return MP_OK;
1724
1.44k
  else
1725
0
    return MP_TRUNC;
1726
1.44k
}
1727
1728
mp_result mp_int_string_len(mp_int z, mp_size radix)
1729
1.52k
{
1730
1.52k
  int  len;
1731
1.52k
1732
1.52k
  CHECK(z != NULL);
1733
1.52k
1734
1.52k
  if (
radix < 1.52k
MP_MIN_RADIX1.52k
||
radix > 1.52k
MP_MAX_RADIX1.52k
)
1735
0
    return MP_RANGE;
1736
1.52k
1737
1.52k
  len = s_outlen(z, radix) + 1; /* for terminator */
1738
1.52k
1739
1.52k
  /* Allow for sign marker on negatives */
1740
1.52k
  if (
MP_SIGN1.52k
(z) == MP_NEG1.52k
)
1741
766
    len += 1;
1742
1.52k
1743
1.52k
  return len;
1744
1.52k
}
1745
1746
/* Read zero-terminated string into z */
1747
mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str)
1748
148
{
1749
148
  return mp_int_read_cstring(z, radix, str, NULL);
1750
148
}
1751
1752
mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
1753
148
{
1754
148
  int ch;
1755
148
1756
148
  CHECK(z != NULL && str != NULL);
1757
148
1758
148
  if (
radix < 148
MP_MIN_RADIX148
||
radix > 148
MP_MAX_RADIX148
)
1759
0
    return MP_RANGE;
1760
148
1761
148
  /* Skip leading whitespace */
1762
148
  
while (148
isspace((int)*str)148
)
1763
0
    ++str;
1764
148
1765
148
  /* Handle leading sign tag (+/-, positive default) */
1766
148
  switch (*str) {
1767
148
  case '-':
1768
54
    MP_SIGN(z) = MP_NEG;
1769
54
    ++str;
1770
148
    break;
1771
148
  case '+':
1772
0
    ++str; /* fallthrough */
1773
94
  default:
1774
94
    MP_SIGN(z) = MP_ZPOS;
1775
148
    break;
1776
148
  }
1777
148
1778
148
  /* Skip leading zeroes */
1779
148
  
while (148
(ch = s_ch2val(*str, radix)) == 0148
)
1780
0
    ++str;
1781
148
1782
148
  /* Make sure there is enough space for the value */
1783
148
  if (!s_pad(z, s_inlen(strlen(str), radix)))
1784
0
    return MP_MEMORY;
1785
148
1786
148
  
MP_USED148
(z) = 1; z->digits[0] = 0;148
1787
148
1788
1.80k
  while (
*str != '\0' && 1.80k
((ch = s_ch2val(*str, radix)) >= 0)1.66k
)
{1.66k
1789
1.66k
    s_dmul(z, (mp_digit)radix);
1790
1.66k
    s_dadd(z, (mp_digit)ch);
1791
1.66k
    ++str;
1792
1.66k
  }
1793
148
1794
148
  CLAMP(z);
1795
148
1796
148
  /* Override sign for zero, even if negative specified. */
1797
148
  if (
CMPZ148
(z) == 0148
)
1798
0
    
MP_SIGN0
(z) = MP_ZPOS0
;
1799
148
1800
148
  if (end != NULL)
1801
0
    *end = (char *)str;
1802
148
1803
148
  /* Return a truncation error if the string has unprocessed characters
1804
148
     remaining, so the caller can tell if the whole string was done */
1805
148
  if (*str != '\0')
1806
0
    return MP_TRUNC;
1807
148
  else
1808
148
    return MP_OK;
1809
148
}
1810
1811
mp_result mp_int_count_bits(mp_int z)
1812
3.36k
{
1813
3.36k
  mp_size  nbits = 0, uz;
1814
3.36k
  mp_digit d;
1815
3.36k
1816
3.36k
  CHECK(z != NULL);
1817
3.36k
1818
3.36k
  uz = MP_USED(z);
1819
3.36k
  if (
uz == 1 && 3.36k
z->digits[0] == 02.04k
)
1820
0
    return 1;
1821
3.36k
1822
3.36k
  --uz;
1823
3.36k
  nbits = uz * MP_DIGIT_BIT;
1824
3.36k
  d = z->digits[uz];
1825
3.36k
1826
59.0k
  while (
d != 059.0k
)
{55.6k
1827
55.6k
    d >>= 1;
1828
55.6k
    ++nbits;
1829
55.6k
  }
1830
3.36k
1831
3.36k
  return nbits;
1832
3.36k
}
1833
1834
mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit)
1835
0
{
1836
0
  static const int PAD_FOR_2C = 1;
1837
0
1838
0
  mp_result res;
1839
0
  int limpos = limit;
1840
0
1841
0
  CHECK(z != NULL && buf != NULL);
1842
0
1843
0
  res = s_tobin(z, buf, &limpos, PAD_FOR_2C);
1844
0
1845
0
  if (
MP_SIGN0
(z) == MP_NEG0
)
1846
0
    s_2comp(buf, limpos);
1847
0
1848
0
  return res;
1849
0
}
1850
1851
mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len)
1852
0
{
1853
0
  mp_size need, i;
1854
0
  unsigned char *tmp;
1855
0
  mp_digit *dz;
1856
0
1857
0
  CHECK(z != NULL && buf != NULL && len > 0);
1858
0
1859
0
  /* Figure out how many digits are needed to represent this value */
1860
0
  need = ((len * CHAR_BIT) + (
MP_DIGIT_BIT0
- 1)) /
MP_DIGIT_BIT0
;
1861
0
  if (!s_pad(z, need))
1862
0
    return MP_MEMORY;
1863
0
1864
0
  mp_int_zero(z);
1865
0
1866
0
  /* If the high-order bit is set, take the 2's complement before reading the
1867
0
     value (it will be restored afterward) */
1868
0
  if (
buf[0] >> (CHAR_BIT - 1)0
)
{0
1869
0
    MP_SIGN(z) = MP_NEG;
1870
0
    s_2comp(buf, len);
1871
0
  }
1872
0
1873
0
  dz = MP_DIGITS(z);
1874
0
  for (tmp = buf, i = len; 
i > 00
;
--i, ++tmp0
)
{0
1875
0
    s_qmul(z, (mp_size) CHAR_BIT);
1876
0
    *dz |= *tmp;
1877
0
  }
1878
0
1879
0
  /* Restore 2's complement if we took it before */
1880
0
  if (
MP_SIGN0
(z) == MP_NEG0
)
1881
0
    s_2comp(buf, len);
1882
0
1883
0
  return MP_OK;
1884
0
}
1885
1886
mp_result mp_int_binary_len(mp_int z)
1887
0
{
1888
0
  mp_result  res = mp_int_count_bits(z);
1889
0
  int        bytes = mp_int_unsigned_len(z);
1890
0
1891
0
  if (res <= 0)
1892
0
    return res;
1893
0
1894
0
  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
1895
0
1896
0
  /* If the highest-order bit falls exactly on a byte boundary, we need to pad
1897
0
     with an extra byte so that the sign will be read correctly when reading it
1898
0
     back in. */
1899
0
  if (bytes * CHAR_BIT == res)
1900
0
    ++bytes;
1901
0
1902
0
  return bytes;
1903
0
}
1904
1905
mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit)
1906
0
{
1907
0
  static const int NO_PADDING = 0;
1908
0
1909
0
  CHECK(z != NULL && buf != NULL);
1910
0
1911
0
  return s_tobin(z, buf, &limit, NO_PADDING);
1912
0
}
1913
1914
mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len)
1915
0
{
1916
0
  mp_size need, i;
1917
0
  unsigned char *tmp;
1918
0
1919
0
  CHECK(z != NULL && buf != NULL && len > 0);
1920
0
1921
0
  /* Figure out how many digits are needed to represent this value */
1922
0
  need = ((len * CHAR_BIT) + (
MP_DIGIT_BIT0
- 1)) /
MP_DIGIT_BIT0
;
1923
0
  if (!s_pad(z, need))
1924
0
    return MP_MEMORY;
1925
0
1926
0
  mp_int_zero(z);
1927
0
1928
0
  for (tmp = buf, i = len; 
i > 00
;
--i, ++tmp0
)
{0
1929
0
    (void) s_qmul(z, CHAR_BIT);
1930
0
    *MP_DIGITS(z) |= *tmp;
1931
0
  }
1932
0
1933
0
  return MP_OK;
1934
0
}
1935
1936
mp_result mp_int_unsigned_len(mp_int z)
1937
1.84k
{
1938
1.84k
  mp_result  res = mp_int_count_bits(z);
1939
1.84k
  int        bytes;
1940
1.84k
1941
1.84k
  if (res <= 0)
1942
0
    return res;
1943
1.84k
1944
1.84k
  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
1945
1.84k
1946
1.84k
  return bytes;
1947
1.84k
}
1948
1949
const char *mp_error_string(mp_result res)
1950
0
{
1951
0
  int ix;
1952
0
  if (res > 0)
1953
0
    return s_unknown_err;
1954
0
1955
0
  res = -res;
1956
0
  for (ix = 0; 
ix < res && 0
s_error_msg[ix] != NULL0
;
++ix0
)
1957
0
    ;
1958
0
1959
0
  if (s_error_msg[ix] != NULL)
1960
0
    return s_error_msg[ix];
1961
0
  else
1962
0
    return s_unknown_err;
1963
0
}
1964
1965
/*------------------------------------------------------------------------*/
1966
/* Private functions for internal use.  These make assumptions.           */
1967
1968
STATIC mp_digit *s_alloc(mp_size num)
1969
20.1M
{
1970
20.1M
  mp_digit *out = malloc(num * sizeof(mp_digit));
1971
20.1M
1972
20.1M
  assert(out != NULL); /* for debugging */
1973
20.1M
#if DEBUG > 1
1974
  {
1975
    mp_digit v = (mp_digit) 0xdeadbeef;
1976
    int      ix;
1977
1978
    for (ix = 0; ix < num; ++ix)
1979
      out[ix] = v;
1980
  }
1981
#endif
1982
20.1M
1983
20.1M
  return out;
1984
20.1M
}
1985
1986
STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
1987
644k
{
1988
644k
#if DEBUG > 1
1989
  mp_digit *new = s_alloc(nsize);
1990
  int       ix;
1991
1992
  for (ix = 0; ix < nsize; ++ix)
1993
    new[ix] = (mp_digit) 0xdeadbeef;
1994
1995
  memcpy(new, old, osize * sizeof(mp_digit));
1996
#else
1997
644k
  mp_digit *new = realloc(old, nsize * sizeof(mp_digit));
1998
644k
1999
644k
  assert(new != NULL); /* for debugging */
2000
644k
#endif
2001
644k
  return new;
2002
644k
}
2003
2004
STATIC void s_free(void *ptr)
2005
20.1M
{
2006
20.1M
  free(ptr);
2007
20.1M
}
2008
2009
STATIC int      s_pad(mp_int z, mp_size min)
2010
188M
{
2011
188M
  if (
MP_ALLOC188M
(z) < min188M
)
{13.6M
2012
13.6M
    mp_size nsize = ROUND_PREC(min);
2013
13.6M
    mp_digit *tmp;
2014
13.6M
2015
13.6M
    if (
(void *)z->digits == (void *)z13.6M
)
{13.0M
2016
13.0M
      if ((tmp = s_alloc(nsize)) == NULL)
2017
0
        return 0;
2018
13.0M
2019
13.0M
      
COPY13.0M
(MP_DIGITS(z), tmp, MP_USED(z));13.0M
2020
13.6M
    }
2021
644k
    else 
if (644k
(tmp = s_realloc(644k
MP_DIGITS644k
(z),
MP_ALLOC644k
(z), nsize)) == NULL)
2022
0
      return 0;
2023
13.6M
2024
13.6M
    
MP_DIGITS13.6M
(z) = tmp;13.6M
2025
13.6M
    MP_ALLOC(z) = nsize;
2026
188M
  }
2027
188M
2028
188M
  return 1;
2029
188M
}
2030
2031
/* Note: This will not work correctly when value == MP_SMALL_MIN */
2032
STATIC void      s_fake(mp_int z, mp_small value, mp_digit vbuf[])
2033
929k
{
2034
929k
  mp_usmall uv = (mp_usmall) (value < 0) ? 
-value5.53k
:
value924k
;
2035
929k
  s_ufake(z, uv, vbuf);
2036
929k
  if (value < 0)
2037
5.53k
    z->sign = MP_NEG;
2038
929k
}
2039
2040
STATIC void      s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[])
2041
47.8M
{
2042
47.8M
  mp_size ndig = (mp_size) s_uvpack(value, vbuf);
2043
47.8M
2044
47.8M
  z->used = ndig;
2045
47.8M
  z->alloc = MP_VALUE_DIGITS(value);
2046
47.8M
  z->sign = MP_ZPOS;
2047
47.8M
  z->digits = vbuf;
2048
47.8M
}
2049
2050
STATIC int      s_cdig(mp_digit *da, mp_digit *db, mp_size len)
2051
52.0M
{
2052
52.0M
  mp_digit *dat = da + len - 1, *dbt = db + len - 1;
2053
52.0M
2054
70.9M
  for (/* */; 
len != 070.9M
;
--len, --dat, --dbt18.9M
)
{63.1M
2055
63.1M
    if (*dat > *dbt)
2056
12.5M
      return 1;
2057
50.6M
    else 
if (50.6M
*dat < *dbt50.6M
)
2058
31.6M
      return -1;
2059
63.1M
  }
2060
52.0M
2061
7.79M
  return 0;
2062
52.0M
}
2063
2064
STATIC int       s_uvpack(mp_usmall uv, mp_digit t[])
2065
47.8M
{
2066
47.8M
  int ndig = 0;
2067
47.8M
2068
47.8M
  if (uv == 0)
2069
791k
    t[ndig++] = 0;
2070
47.8M
  else {
2071
137M
    while (
uv != 0137M
)
{90.4M
2072
90.4M
      t[ndig++] = (mp_digit) uv;
2073
90.4M
      uv >>= MP_DIGIT_BIT/2;
2074
90.4M
      uv >>= MP_DIGIT_BIT/2;
2075
90.4M
    }
2076
47.8M
  }
2077
47.8M
2078
47.8M
  return ndig;
2079
47.8M
}
2080
2081
STATIC int      s_ucmp(mp_int a, mp_int b)
2082
135M
{
2083
135M
  mp_size  ua = 
MP_USED135M
(a), ub =
MP_USED135M
(b);
2084
135M
2085
135M
  if (ua > ub)
2086
16.6M
    return 1;
2087
118M
  else 
if (118M
ub > ua118M
)
2088
66.5M
    return -1;
2089
118M
  else
2090
52.0M
    
return s_cdig(52.0M
MP_DIGITS52.0M
(a),
MP_DIGITS52.0M
(b), ua);
2091
135M
}
2092
2093
STATIC int      s_vcmp(mp_int a, mp_small v)
2094
46.6M
{
2095
46.6M
  mp_usmall uv = (v < 0) ? 
-(mp_usmall) v7.77M
:
(mp_usmall) v38.8M
;
2096
46.6M
  return s_uvcmp(a, uv);
2097
46.6M
}
2098
2099
STATIC int      s_uvcmp(mp_int a, mp_usmall uv)
2100
46.6M
{
2101
46.6M
  mpz_t vtmp;
2102
46.6M
  mp_digit vdig[MP_VALUE_DIGITS(uv)];
2103
46.6M
2104
46.6M
  s_ufake(&vtmp, uv, vdig);
2105
46.6M
  return s_ucmp(a, &vtmp);
2106
46.6M
}
2107
2108
STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
2109
           mp_size size_a, mp_size size_b)
2110
7.84M
{
2111
7.84M
  mp_size pos;
2112
7.84M
  mp_word w = 0;
2113
7.84M
2114
7.84M
  /* Insure that da is the longer of the two to simplify later code */
2115
7.84M
  if (
size_b > size_a7.84M
)
{749k
2116
749k
    SWAP(mp_digit *, da, db);
2117
749k
    SWAP(mp_size, size_a, size_b);
2118
7.84M
  }
2119
7.84M
2120
7.84M
  /* Add corresponding digits until the shorter number runs out */
2121
17.0M
  for (pos = 0; 
pos < size_b17.0M
;
++pos, ++da, ++db, ++dc9.18M
)
{9.18M
2122
9.18M
    w = w + (mp_word) *da + (mp_word) *db;
2123
9.18M
    *dc = LOWER_HALF(w);
2124
9.18M
    w = UPPER_HALF(w);
2125
9.18M
  }
2126
7.84M
2127
7.84M
  /* Propagate carries as far as necessary */
2128
13.2M
  for (/* */; 
pos < size_a13.2M
;
++pos, ++da, ++dc5.36M
)
{5.36M
2129
5.36M
    w = w + *da;
2130
5.36M
2131
5.36M
    *dc = LOWER_HALF(w);
2132
5.36M
    w = UPPER_HALF(w);
2133
7.84M
  }
2134
7.84M
2135
7.84M
  /* Return carry out */
2136
7.84M
  return (mp_digit)w;
2137
7.84M
}
2138
2139
STATIC void     s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
2140
           mp_size size_a, mp_size size_b)
2141
80.6M
{
2142
80.6M
  mp_size pos;
2143
80.6M
  mp_word w = 0;
2144
80.6M
2145
80.6M
  /* We assume that |a| >= |b| so this should definitely hold */
2146
80.6M
  assert(size_a >= size_b);
2147
80.6M
2148
80.6M
  /* Subtract corresponding digits and propagate borrow */
2149
215M
  for (pos = 0; 
pos < size_b215M
;
++pos, ++da, ++db, ++dc135M
)
{135M
2150
135M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2151
135M
   (mp_word)*da) - w - (mp_word)*db;
2152
135M
2153
135M
    *dc = LOWER_HALF(w);
2154
135M
    w = (UPPER_HALF(w) == 0);
2155
135M
  }
2156
80.6M
2157
80.6M
  /* Finish the subtraction for remaining upper digits of da */
2158
165M
  for (/* */; 
pos < size_a165M
;
++pos, ++da, ++dc84.9M
)
{84.9M
2159
84.9M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2160
84.9M
   (mp_word)*da) - w;
2161
84.9M
2162
84.9M
    *dc = LOWER_HALF(w);
2163
84.9M
    w = (UPPER_HALF(w) == 0);
2164
84.9M
  }
2165
80.6M
2166
80.6M
  /* If there is a borrow out at the end, it violates the precondition */
2167
80.6M
  assert(w == 0);
2168
80.6M
}
2169
2170
STATIC int       s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
2171
      mp_size size_a, mp_size size_b)
2172
7.18M
{
2173
7.18M
  mp_size  bot_size;
2174
7.18M
2175
7.18M
  /* Make sure b is the smaller of the two input values */
2176
7.18M
  if (
size_b > size_a7.18M
)
{2.97M
2177
2.97M
    SWAP(mp_digit *, da, db);
2178
2.97M
    SWAP(mp_size, size_a, size_b);
2179
7.18M
  }
2180
7.18M
2181
7.18M
  /* Insure that the bottom is the larger half in an odd-length split; the code
2182
7.18M
     below relies on this being true.
2183
7.18M
   */
2184
7.18M
  bot_size = (size_a + 1) / 2;
2185
7.18M
2186
7.18M
  /* If the values are big enough to bother with recursion, use the Karatsuba
2187
7.18M
     algorithm to compute the product; otherwise use the normal multiplication
2188
7.18M
     algorithm
2189
7.18M
   */
2190
7.18M
  if (multiply_threshold &&
2191
7.18M
      size_a >= multiply_threshold &&
2192
7.18M
      
size_b > bot_size2.41k
)
{36
2193
36
2194
36
    mp_digit *t1, *t2, *t3, carry;
2195
36
2196
36
    mp_digit *a_top = da + bot_size;
2197
36
    mp_digit *b_top = db + bot_size;
2198
36
2199
36
    mp_size  at_size = size_a - bot_size;
2200
36
    mp_size  bt_size = size_b - bot_size;
2201
36
    mp_size  buf_size = 2 * bot_size;
2202
36
2203
36
    /* Do a single allocation for all three temporary buffers needed; each
2204
36
       buffer must be big enough to hold the product of two bottom halves, and
2205
36
       one buffer needs space for the completed product; twice the space is
2206
36
       plenty.
2207
36
     */
2208
36
    if (
(t1 = s_alloc(4 * buf_size)) == NULL36
)
return 00
;
2209
36
    t2 = t1 + buf_size;
2210
36
    t3 = t2 + buf_size;
2211
36
    ZERO(t1, 4 * buf_size);
2212
36
2213
36
    /* t1 and t2 are initially used as temporaries to compute the inner product
2214
36
       (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
2215
36
     */
2216
36
    carry = s_uadd(da, a_top, t1, bot_size, at_size);      /* t1 = a1 + a0 */
2217
36
    t1[bot_size] = carry;
2218
36
2219
36
    carry = s_uadd(db, b_top, t2, bot_size, bt_size);      /* t2 = b1 + b0 */
2220
36
    t2[bot_size] = carry;
2221
36
2222
36
    (void) s_kmul(t1, t2, t3, bot_size + 1, bot_size + 1); /* t3 = t1 * t2 */
2223
36
2224
36
    /* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that
2225
36
       we're left with only the pieces we want:  t3 = a1b0 + a0b1
2226
36
     */
2227
36
    ZERO(t1, buf_size);
2228
36
    ZERO(t2, buf_size);
2229
36
    (void) s_kmul(da, db, t1, bot_size, bot_size);     /* t1 = a0 * b0 */
2230
36
    (void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */
2231
36
2232
36
    /* Subtract out t1 and t2 to get the inner product */
2233
36
    s_usub(t3, t1, t3, buf_size + 2, buf_size);
2234
36
    s_usub(t3, t2, t3, buf_size + 2, buf_size);
2235
36
2236
36
    /* Assemble the output value */
2237
36
    COPY(t1, dc, buf_size);
2238
36
    carry = s_uadd(t3, dc + bot_size, dc + bot_size,
2239
36
       buf_size + 1, buf_size);
2240
36
    assert(carry == 0);
2241
36
2242
36
    carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
2243
36
       buf_size, buf_size);
2244
36
    assert(carry == 0);
2245
36
2246
36
    s_free(t1); /* note t2 and t3 are just internal pointers to t1 */
2247
7.18M
  }
2248
7.18M
  else {
2249
7.18M
    s_umul(da, db, dc, size_a, size_b);
2250
7.18M
  }
2251
7.18M
2252
7.18M
  return 1;
2253
7.18M
}
2254
2255
STATIC void     s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
2256
           mp_size size_a, mp_size size_b)
2257
7.18M
{
2258
7.18M
  mp_size a, b;
2259
7.18M
  mp_word w;
2260
7.18M
2261
25.0M
  for (a = 0; 
a < size_a25.0M
;
++a, ++dc, ++da17.8M
)
{17.8M
2262
17.8M
    mp_digit *dct = dc;
2263
17.8M
    mp_digit *dbt = db;
2264
17.8M
2265
17.8M
    if (*da == 0)
2266
1.30M
      continue;
2267
17.8M
2268
17.8M
    w = 0;
2269
39.4M
    for (b = 0; 
b < size_b39.4M
;
++b, ++dbt, ++dct22.9M
)
{22.9M
2270
22.9M
      w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
2271
22.9M
2272
22.9M
      *dct = LOWER_HALF(w);
2273
22.9M
      w = UPPER_HALF(w);
2274
22.9M
    }
2275
16.5M
2276
16.5M
    *dct = (mp_digit)w;
2277
16.5M
  }
2278
7.18M
}
2279
2280
STATIC int       s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2281
0
{
2282
0
  if (
multiply_threshold && 0
size_a > multiply_threshold0
)
{0
2283
0
    mp_size  bot_size = (size_a + 1) / 2;
2284
0
    mp_digit *a_top = da + bot_size;
2285
0
    mp_digit *t1, *t2, *t3, carry;
2286
0
    mp_size  at_size = size_a - bot_size;
2287
0
    mp_size  buf_size = 2 * bot_size;
2288
0
2289
0
    if (
(t1 = s_alloc(4 * buf_size)) == NULL0
)
return 00
;
2290
0
    t2 = t1 + buf_size;
2291
0
    t3 = t2 + buf_size;
2292
0
    ZERO(t1, 4 * buf_size);
2293
0
2294
0
    (void) s_ksqr(da, t1, bot_size);    /* t1 = a0 ^ 2 */
2295
0
    (void) s_ksqr(a_top, t2, at_size);  /* t2 = a1 ^ 2 */
2296
0
2297
0
    (void) s_kmul(da, a_top, t3, bot_size, at_size);  /* t3 = a0 * a1 */
2298
0
2299
0
    /* Quick multiply t3 by 2, shifting left (can't overflow) */
2300
0
    {
2301
0
      int i, top = bot_size + at_size;
2302
0
      mp_word w, save = 0;
2303
0
2304
0
      for (i = 0; 
i < top0
;
++i0
)
{0
2305
0
  w = t3[i];
2306
0
  w = (w << 1) | save;
2307
0
  t3[i] = LOWER_HALF(w);
2308
0
  save = UPPER_HALF(w);
2309
0
      }
2310
0
      t3[i] = LOWER_HALF(save);
2311
0
    }
2312
0
2313
0
    /* Assemble the output value */
2314
0
    COPY(t1, dc, 2 * bot_size);
2315
0
    carry = s_uadd(t3, dc + bot_size, dc + bot_size,
2316
0
       buf_size + 1, buf_size);
2317
0
    assert(carry == 0);
2318
0
2319
0
    carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
2320
0
       buf_size, buf_size);
2321
0
    assert(carry == 0);
2322
0
2323
0
    s_free(t1); /* note that t2 and t2 are internal pointers only */
2324
0
2325
0
  } 
2326
0
  else {
2327
0
    s_usqr(da, dc, size_a);
2328
0
  }
2329
0
2330
0
  return 1;
2331
0
}
2332
2333
STATIC void      s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2334
0
{
2335
0
  mp_size i, j;
2336
0
  mp_word w;
2337
0
2338
0
  for (i = 0; 
i < size_a0
;
++i, dc += 2, ++da0
)
{0
2339
0
    mp_digit  *dct = dc, *dat = da;
2340
0
2341
0
    if (*da == 0)
2342
0
      continue;
2343
0
2344
0
    /* Take care of the first digit, no rollover */
2345
0
    w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct;
2346
0
    *dct = LOWER_HALF(w);
2347
0
    w = UPPER_HALF(w);
2348
0
    ++dat; ++dct;
2349
0
2350
0
    for (j = i + 1; 
j < size_a0
;
++j, ++dat, ++dct0
)
{0
2351
0
      mp_word  t = (mp_word)*da * (mp_word)*dat;
2352
0
      mp_word  u = w + (mp_word)*dct, ov = 0;
2353
0
2354
0
      /* Check if doubling t will overflow a word */
2355
0
      if (HIGH_BIT_SET(t))
2356
0
  ov = 1;
2357
0
2358
0
      w = t + t;
2359
0
2360
0
      /* Check if adding u to w will overflow a word */
2361
0
      if (ADD_WILL_OVERFLOW(w, u))
2362
0
  ov = 1;
2363
0
2364
0
      w += u;
2365
0
2366
0
      *dct = LOWER_HALF(w);
2367
0
      w = UPPER_HALF(w);
2368
0
      if (
ov0
)
{0
2369
0
  w += MP_DIGIT_MAX; /* MP_RADIX */
2370
0
  ++w;
2371
0
      }
2372
0
    }
2373
0
2374
0
    w = w + *dct;
2375
0
    *dct = (mp_digit)w;
2376
0
    while (
(w = 0
UPPER_HALF0
(w)) != 0)
{0
2377
0
      ++dct; w = w + *dct;
2378
0
      *dct = LOWER_HALF(w);
2379
0
    }
2380
0
2381
0
    assert(w == 0);
2382
0
  }
2383
0
}
2384
2385
STATIC void      s_dadd(mp_int a, mp_digit b)
2386
1.66k
{
2387
1.66k
  mp_word w = 0;
2388
1.66k
  mp_digit *da = MP_DIGITS(a);
2389
1.66k
  mp_size ua = MP_USED(a);
2390
1.66k
2391
1.66k
  w = (mp_word)*da + b;
2392
1.66k
  *da++ = LOWER_HALF(w);
2393
1.66k
  w = UPPER_HALF(w);
2394
1.66k
2395
1.86k
  for (ua -= 1; 
ua > 01.86k
;
--ua, ++da200
)
{200
2396
200
    w = (mp_word)*da + w;
2397
200
2398
200
    *da = LOWER_HALF(w);
2399
200
    w = UPPER_HALF(w);
2400
1.66k
  }
2401
1.66k
2402
1.66k
  if (
w1.66k
)
{2
2403
2
    *da = (mp_digit)w;
2404
2
    MP_USED(a) += 1;
2405
1.66k
  }
2406
1.66k
}
2407
2408
STATIC void      s_dmul(mp_int a, mp_digit b)
2409
1.66k
{
2410
1.66k
  mp_word w = 0;
2411
1.66k
  mp_digit *da = MP_DIGITS(a);
2412
1.66k
  mp_size ua = MP_USED(a);
2413
1.66k
2414
3.50k
  while (
ua > 03.50k
)
{1.84k
2415
1.84k
    w = (mp_word)*da * b + w;
2416
1.84k
    *da++ = LOWER_HALF(w);
2417
1.84k
    w = UPPER_HALF(w);
2418
1.84k
    --ua;
2419
1.84k
  }
2420
1.66k
2421
1.66k
  if (
w1.66k
)
{20
2422
20
    *da = (mp_digit)w;
2423
20
    MP_USED(a) += 1;
2424
1.66k
  }
2425
1.66k
}
2426
2427
STATIC void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
2428
2.08M
{
2429
2.08M
  mp_word  w = 0;
2430
2.08M
2431
10.6M
  while (
size_a > 010.6M
)
{8.54M
2432
8.54M
    w = (mp_word)*da++ * (mp_word)b + w;
2433
8.54M
2434
8.54M
    *dc++ = LOWER_HALF(w);
2435
8.54M
    w = UPPER_HALF(w);
2436
8.54M
    --size_a;
2437
8.54M
  }
2438
2.08M
2439
2.08M
  if (w)
2440
0
    
*dc = 0
LOWER_HALF0
(w);
2441
2.08M
}
2442
2443
STATIC mp_digit  s_ddiv(mp_int a, mp_digit b)
2444
832k
{
2445
832k
  mp_word w = 0, qdigit;
2446
832k
  mp_size ua = MP_USED(a);
2447
832k
  mp_digit *da = MP_DIGITS(a) + ua - 1;
2448
832k
2449
2.86M
  for (/* */; 
ua > 02.86M
;
--ua, --da2.02M
)
{2.02M
2450
2.02M
    w = (w << MP_DIGIT_BIT) | *da;
2451
2.02M
2452
2.02M
    if (
w >= b2.02M
)
{1.55M
2453
1.55M
      qdigit = w / b;
2454
1.55M
      w = w % b;
2455
2.02M
    }
2456
2.02M
    else {
2457
476k
      qdigit = 0;
2458
2.02M
    }
2459
2.02M
2460
2.02M
    *da = (mp_digit)qdigit;
2461
2.02M
  }
2462
832k
2463
832k
  CLAMP(a);
2464
832k
  return (mp_digit)w;
2465
832k
}
2466
2467
STATIC void     s_qdiv(mp_int z, mp_size p2)
2468
80.5M
{
2469
80.5M
  mp_size ndig = p2 / 
MP_DIGIT_BIT80.5M
, nbits = p2 %
MP_DIGIT_BIT80.5M
;
2470
80.5M
  mp_size uz = MP_USED(z);
2471
80.5M
2472
80.5M
  if (
ndig80.5M
)
{260k
2473
260k
    mp_size  mark;
2474
260k
    mp_digit *to, *from;
2475
260k
2476
260k
    if (
ndig >= uz260k
)
{0
2477
0
      mp_int_zero(z);
2478
0
      return;
2479
260k
    }
2480
260k
2481
260k
    
to = 260k
MP_DIGITS260k
(z); from = to + ndig;
2482
260k
2483
603k
    for (mark = ndig; 
mark < uz603k
;
++mark342k
)
2484
342k
      *to++ = *from++;
2485
260k
2486
260k
    MP_USED(z) = uz - ndig;
2487
80.5M
  }
2488
80.5M
2489
80.5M
  
if (80.5M
nbits80.5M
)
{76.1M
2490
76.1M
    mp_digit d = 0, *dz, save;
2491
76.1M
    mp_size  up = MP_DIGIT_BIT - nbits;
2492
76.1M
2493
76.1M
    uz = MP_USED(z);
2494
76.1M
    dz = MP_DIGITS(z) + uz - 1;
2495
76.1M
2496
285M
    for (/* */; 
uz > 0285M
;
--uz, --dz208M
)
{208M
2497
208M
      save = *dz;
2498
208M
2499
208M
      *dz = (*dz >> nbits) | (d << up);
2500
208M
      d = save;
2501
208M
    }
2502
76.1M
2503
76.1M
    CLAMP(z);
2504
80.5M
  }
2505
80.5M
2506
80.5M
  if (
MP_USED80.5M
(z) == 1 && 80.5M
z->digits[0] == 028.0M
)
2507
777k
    
MP_SIGN777k
(z) = MP_ZPOS777k
;
2508
80.5M
}
2509
2510
STATIC void     s_qmod(mp_int z, mp_size p2)
2511
126k
{
2512
126k
  mp_size start = p2 / 
MP_DIGIT_BIT126k
+ 1, rest = p2 %
MP_DIGIT_BIT126k
;
2513
126k
  mp_size uz = MP_USED(z);
2514
126k
  mp_digit mask = (1u << rest) - 1;
2515
126k
2516
126k
  if (
start <= uz126k
)
{126k
2517
126k
    MP_USED(z) = start;
2518
126k
    z->digits[start - 1] &= mask;
2519
126k
    CLAMP(z);
2520
126k
  }
2521
126k
}
2522
2523
STATIC int      s_qmul(mp_int z, mp_size p2)
2524
4.16M
{
2525
4.16M
  mp_size   uz, need, rest, extra, i;
2526
4.16M
  mp_digit *from, *to, d;
2527
4.16M
2528
4.16M
  if (p2 == 0)
2529
1.38M
    return 1;
2530
4.16M
2531
2.77M
  
uz = 2.77M
MP_USED2.77M
(z);
2532
2.77M
  need = p2 / 
MP_DIGIT_BIT2.77M
; rest = p2 %
MP_DIGIT_BIT2.77M
;
2533
2.77M
2534
2.77M
  /* Figure out if we need an extra digit at the top end; this occurs if the
2535
2.77M
     topmost `rest' bits of the high-order digit of z are not zero, meaning
2536
2.77M
     they will be shifted off the end if not preserved */
2537
2.77M
  extra = 0;
2538
2.77M
  if (
rest != 02.77M
)
{2.71M
2539
2.71M
    mp_digit *dz = MP_DIGITS(z) + uz - 1;
2540
2.71M
2541
2.71M
    if (
(*dz >> (2.71M
MP_DIGIT_BIT2.71M
- rest)) != 0)
2542
578k
      extra = 1;
2543
2.77M
  }
2544
2.77M
2545
2.77M
  if (!s_pad(z, uz + need + extra))
2546
0
    return 0;
2547
2.77M
2548
2.77M
  /* If we need to shift by whole digits, do that in one pass, then
2549
2.77M
     to back and shift by partial digits.
2550
2.77M
   */
2551
2.77M
  
if (2.77M
need > 02.77M
)
{86.9k
2552
86.9k
    from = MP_DIGITS(z) + uz - 1;
2553
86.9k
    to = from + need;
2554
86.9k
2555
173k
    for (i = 0; 
i < uz173k
;
++i87.0k
)
2556
87.0k
      *to-- = *from--;
2557
86.9k
2558
86.9k
    ZERO(MP_DIGITS(z), need);
2559
86.9k
    uz += need;
2560
2.77M
  }
2561
2.77M
2562
2.77M
  if (
rest2.77M
)
{2.71M
2563
2.71M
    d = 0;
2564
10.0M
    for (i = need, from = 
MP_DIGITS2.71M
(z) + need;
i < uz10.0M
;
++i, ++from7.37M
)
{7.37M
2565
7.37M
      mp_digit save = *from;
2566
7.37M
      
2567
7.37M
      *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
2568
7.37M
      d = save;
2569
7.37M
    }
2570
2.71M
2571
2.71M
    d >>= (MP_DIGIT_BIT - rest);
2572
2.71M
    if (
d != 02.71M
)
{578k
2573
578k
      *from = d;
2574
578k
      uz += extra;
2575
2.71M
    }
2576
2.77M
  }
2577
2.77M
2578
2.77M
  MP_USED(z) = uz;
2579
2.77M
  CLAMP(z);
2580
2.77M
2581
2.77M
  return 1;
2582
4.16M
}
2583
2584
/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
2585
   The sign of the result is always zero/positive.
2586
 */
2587
STATIC int       s_qsub(mp_int z, mp_size p2)
2588
0
{
2589
0
  mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
2590
0
  mp_size  tdig = (p2 / MP_DIGIT_BIT), pos;
2591
0
  mp_word  w = 0;
2592
0
2593
0
  if (!s_pad(z, tdig + 1))
2594
0
    return 0;
2595
0
2596
0
  
for (pos = 0, zp = 0
MP_DIGITS0
(z);
pos < tdig0
;
++pos, ++zp0
)
{0
2597
0
    w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp;
2598
0
2599
0
    *zp = LOWER_HALF(w);
2600
0
    w = 
UPPER_HALF0
(w) ?
00
:
10
;
2601
0
  }
2602
0
2603
0
  w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp;
2604
0
  *zp = LOWER_HALF(w);
2605
0
2606
0
  assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
2607
0
2608
0
  MP_SIGN(z) = MP_ZPOS;
2609
0
  CLAMP(z);
2610
0
2611
0
  return 1;
2612
0
}
2613
2614
STATIC int      s_dp2k(mp_int z)
2615
78.9M
{
2616
78.9M
  int       k = 0;
2617
78.9M
  mp_digit *dp = MP_DIGITS(z), d;
2618
78.9M
2619
78.9M
  if (
MP_USED78.9M
(z) == 1 && 78.9M
*dp == 025.2M
)
2620
0
    return 1;
2621
78.9M
2622
79.4M
  
while (78.9M
*dp == 079.4M
)
{425k
2623
425k
    k += MP_DIGIT_BIT;
2624
425k
    ++dp;
2625
78.9M
  }
2626
78.9M
2627
78.9M
  d = *dp;
2628
231M
  while (
(d & 1) == 0231M
)
{152M
2629
152M
    d >>= 1;
2630
152M
    ++k;
2631
152M
  }
2632
78.9M
2633
78.9M
  return k;
2634
78.9M
}
2635
2636
STATIC int       s_isp2(mp_int z)
2637
2.45M
{
2638
2.45M
  mp_size uz = MP_USED(z), k = 0;
2639
2.45M
  mp_digit *dz = MP_DIGITS(z), d;
2640
2.45M
2641
2.51M
  while (
uz > 12.51M
)
{895k
2642
895k
    if (*dz++ != 0)
2643
838k
      return -1;
2644
56.9k
    
k += 56.9k
MP_DIGIT_BIT56.9k
;
2645
56.9k
    --uz;
2646
2.45M
  }
2647
2.45M
2648
2.45M
  d = *dz;
2649
3.99M
  while (
d > 13.99M
)
{3.18M
2650
3.18M
    if (d & 1)
2651
808k
      return -1;
2652
3.18M
    ++k; d >>= 1;
2653
2.37M
  }
2654
1.61M
2655
808k
  return (int) k;
2656
2.45M
}
2657
2658
STATIC int       s_2expt(mp_int z, mp_small k)
2659
0
{
2660
0
  mp_size  ndig, rest;
2661
0
  mp_digit *dz;
2662
0
2663
0
  ndig = (k + 
MP_DIGIT_BIT0
) /
MP_DIGIT_BIT0
;
2664
0
  rest = k % MP_DIGIT_BIT;
2665
0
2666
0
  if (!s_pad(z, ndig))
2667
0
    return 0;
2668
0
2669
0
  
dz = 0
MP_DIGITS0
(z);
2670
0
  ZERO(dz, ndig);
2671
0
  *(dz + ndig - 1) = (1 << rest);
2672
0
  MP_USED(z) = ndig;
2673
0
2674
0
  return 1;
2675
0
}
2676
2677
STATIC int      s_norm(mp_int a, mp_int b)
2678
840k
{
2679
840k
  mp_digit d = b->digits[MP_USED(b) - 1];
2680
840k
  int k = 0;
2681
840k
2682
15.0M
  while (
d < (1u << (mp_digit)(15.0M
MP_DIGIT_BIT15.0M
- 1)))
{ /* d < (MP_RADIX / 2) */14.1M
2683
14.1M
    d <<= 1;
2684
14.1M
    ++k;
2685
14.1M
  }
2686
840k
2687
840k
  /* These multiplications can't fail */
2688
840k
  if (
k != 0840k
)
{825k
2689
825k
    (void) s_qmul(a, (mp_size) k);
2690
825k
    (void) s_qmul(b, (mp_size) k);
2691
840k
  }
2692
840k
2693
840k
  return k;
2694
840k
}
2695
2696
STATIC mp_result s_brmu(mp_int z, mp_int m)
2697
0
{
2698
0
  mp_size um = MP_USED(m) * 2;
2699
0
2700
0
  if (!s_pad(z, um))
2701
0
    return MP_MEMORY;
2702
0
2703
0
  
s_2expt(z, 0
MP_DIGIT_BIT0
* um);
2704
0
  return mp_int_div(z, m, z, NULL);
2705
0
}
2706
2707
STATIC int       s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
2708
0
{
2709
0
  mp_size   um = MP_USED(m), umb_p1, umb_m1;
2710
0
2711
0
  umb_p1 = (um + 1) * MP_DIGIT_BIT;
2712
0
  umb_m1 = (um - 1) * MP_DIGIT_BIT;
2713
0
2714
0
  if (mp_int_copy(x, q1) != MP_OK)
2715
0
    return 0;
2716
0
2717
0
  /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
2718
0
  s_qdiv(q1, umb_m1);
2719
0
  UMUL(q1, mu, q2);
2720
0
  s_qdiv(q2, umb_p1);
2721
0
2722
0
  /* Set x = x mod b^(k+1) */
2723
0
  s_qmod(x, umb_p1);
2724
0
2725
0
  /* Now, q is a guess for the quotient a / m.
2726
0
     Compute x - q * m mod b^(k+1), replacing x.  This may be off
2727
0
     by a factor of 2m, but no more than that.
2728
0
   */
2729
0
  UMUL(q2, m, q1);
2730
0
  s_qmod(q1, umb_p1);
2731
0
  (void) mp_int_sub(x, q1, x); /* can't fail */
2732
0
2733
0
  /* The result may be < 0; if it is, add b^(k+1) to pin it in the proper
2734
0
     range. */
2735
0
  if (
(0
CMPZ0
(x) < 0) &&
!s_qsub(x, umb_p1)0
)
2736
0
    return 0;
2737
0
2738
0
  /* If x > m, we need to back it off until it is in range.  This will be
2739
0
     required at most twice.  */
2740
0
  
if (0
mp_int_compare(x, m) >= 00
)
{0
2741
0
    (void) mp_int_sub(x, m, x);
2742
0
    if (mp_int_compare(x, m) >= 0)
2743
0
      (void) mp_int_sub(x, m, x);
2744
0
  }
2745
0
2746
0
  /* At this point, x has been properly reduced. */
2747
0
  return 1;
2748
0
}
2749
2750
/* Perform modular exponentiation using Barrett's method, where mu is the
2751
   reduction constant for m.  Assumes a < m, b > 0. */
2752
STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
2753
0
{
2754
0
  mp_digit  *db, *dbt, umu, d;
2755
0
  mp_result res;
2756
0
  DECLARE_TEMP(3);
2757
0
2758
0
  umu = 
MP_USED0
(mu); db =
MP_DIGITS0
(b); dbt = db +
MP_USED0
(b) - 1;
2759
0
2760
0
  while (
last__ < 30
)
{0
2761
0
    SETUP(mp_int_init_size(LAST_TEMP(), 4 * umu));
2762
0
    
ZERO0
(MP_DIGITS(TEMP(last__ - 1)), MP_ALLOC(TEMP(last__ - 1)));0
2763
0
  }
2764
0
2765
0
  (void) mp_int_set_value(c, 1);
2766
0
2767
0
  /* Take care of low-order digits */
2768
0
  while (
db < dbt0
)
{0
2769
0
    int      i;
2770
0
2771
0
    for (d = *db, i = 
MP_DIGIT_BIT0
;
i > 00
;
--i, d >>= 10
)
{0
2772
0
      if (
d & 10
)
{0
2773
0
  /* The use of a second temporary avoids allocation */
2774
0
  UMUL(c, a, TEMP(0));
2775
0
  if (
!s_reduce(0
TEMP0
(0), m, mu,
TEMP0
(1),
TEMP0
(2)))
{0
2776
0
    res = MP_MEMORY; goto CLEANUP;
2777
0
  }
2778
0
  
mp_int_copy(0
TEMP0
(0), c);
2779
0
      }
2780
0
2781
0
2782
0
      
USQR0
(a, TEMP(0));0
2783
0
      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
2784
0
      if (
!s_reduce(0
TEMP0
(0), m, mu,
TEMP0
(1),
TEMP0
(2)))
{0
2785
0
  res = MP_MEMORY; goto CLEANUP;
2786
0
      }
2787
0
      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
2788
0
      mp_int_copy(TEMP(0), a);
2789
0
    }
2790
0
2791
0
    ++db;
2792
0
  }
2793
0
2794
0
  /* Take care of highest-order digit */
2795
0
  d = *dbt;
2796
0
  for (;;) {
2797
0
    if (
d & 10
)
{0
2798
0
      UMUL(c, a, TEMP(0));
2799
0
      if (
!s_reduce(0
TEMP0
(0), m, mu,
TEMP0
(1),
TEMP0
(2)))
{0
2800
0
  res = MP_MEMORY; goto CLEANUP;
2801
0
      }
2802
0
      
mp_int_copy(0
TEMP0
(0), c);
2803
0
    }
2804
0
2805
0
    d >>= 1;
2806
0
    if (
!d0
)
break0
;
2807
0
2808
0
    
USQR0
(a, TEMP(0));0
2809
0
    if (
!s_reduce(0
TEMP0
(0), m, mu,
TEMP0
(1),
TEMP0
(2)))
{0
2810
0
      res = MP_MEMORY; goto CLEANUP;
2811
0
    }
2812
0
    
(void) mp_int_copy(0
TEMP0
(0), a);
2813
0
  }
2814
0
2815
0
  
CLEANUP_TEMP0
();0
2816
0
  return res;
2817
0
}
2818
2819
/* Division of nonnegative integers
2820
2821
   This function implements division algorithm for unsigned multi-precision
2822
   integers. The algorithm is based on Algorithm D from Knuth's "The Art of
2823
   Computer Programming", 3rd ed. 1998, pg 272-273.
2824
2825
   We diverge from Knuth's algorithm in that we do not perform the subtraction
2826
   from the remainder until we have determined that we have the correct
2827
   quotient digit. This makes our algorithm less efficient that Knuth because
2828
   we might have to perform multiple multiplication and comparison steps before
2829
   the subtraction. The advantage is that it is easy to implement and ensure
2830
   correctness without worrying about underflow from the subtraction.
2831
2832
   inputs: u   a n+m digit integer in base b (b is 2^MP_DIGIT_BIT)
2833
           v   a n   digit integer in base b (b is 2^MP_DIGIT_BIT)
2834
           n >= 1
2835
           m >= 0
2836
  outputs: u / v stored in u
2837
           u % v stored in v
2838
 */
2839
1.64M
STATIC mp_result s_udiv_knuth(mp_int u, mp_int v) {
2840
1.64M
  mpz_t q, r, t;
2841
1.64M
  mp_result
2842
1.64M
  res = MP_OK;
2843
1.64M
  int k,j;
2844
1.64M
  mp_size m,n;
2845
1.64M
2846
1.64M
  /* Force signs to positive */
2847
1.64M
  MP_SIGN(u) = MP_ZPOS;
2848
1.64M
  MP_SIGN(v) = MP_ZPOS;
2849
1.64M
2850
1.64M
  /* Use simple division algorithm when v is only one digit long */
2851
1.64M
  if (
MP_USED1.64M
(v) == 11.64M
)
{806k
2852
806k
    mp_digit d, rem;
2853
806k
    d   = v->digits[0];
2854
806k
    rem = s_ddiv(u, d);
2855
806k
    mp_int_set_value(v, rem);
2856
806k
    return MP_OK;
2857
1.64M
  }
2858
1.64M
2859
1.64M
  /* Algorithm D
2860
1.64M
2861
1.64M
     The n and m variables are defined as used by Knuth.
2862
1.64M
     u is an n digit number with digits u_{n-1}..u_0.
2863
1.64M
     v is an n+m digit number with digits from v_{m+n-1}..v_0.
2864
1.64M
     We require that n > 1 and m >= 0
2865
1.64M
   */
2866
840k
  
n = 840k
MP_USED840k
(v);
2867
840k
  m = MP_USED(u) - n;
2868
840k
  assert(n > 1);
2869
840k
  assert(m >= 0);
2870
840k
2871
840k
  /* D1: Normalize.
2872
840k
     The normalization step provides the necessary condition for Theorem B,
2873
840k
     which states that the quotient estimate for q_j, call it qhat
2874
840k
2875
840k
       qhat = u_{j+n}u_{j+n-1} / v_{n-1}
2876
840k
2877
840k
     is bounded by
2878
840k
2879
840k
      qhat - 2 <= q_j <= qhat.
2880
840k
2881
840k
     That is, qhat is always greater than the actual quotient digit q,
2882
840k
     and it is never more than two larger than the actual quotient digit.
2883
840k
   */
2884
840k
  k = s_norm(u, v);
2885
840k
2886
840k
  /* Extend size of u by one if needed.
2887
840k
2888
840k
     The algorithm begins with a value of u that has one more digit of input.
2889
840k
     The normalization step sets u_{m+n}..u_0 = 2^k * u_{m+n-1}..u_0. If the
2890
840k
     multiplication did not increase the number of digits of u, we need to add
2891
840k
     a leading zero here.
2892
840k
   */
2893
840k
  if (
k == 0 || 840k
MP_USED825k
(u) != m + n + 1825k
)
{318k
2894
318k
    if (!s_pad(u, m+n+1))
2895
0
      return MP_MEMORY;
2896
318k
    u->digits[m+n] = 0;
2897
318k
    u->used = m+n+1;
2898
840k
  }
2899
840k
2900
840k
  /* Add a leading 0 to v.
2901
840k
2902
840k
     The multiplication in step D4 multiplies qhat * 0v_{n-1}..v_0.  We need to
2903
840k
     add the leading zero to v here to ensure that the multiplication will
2904
840k
     produce the full n+1 digit result.
2905
840k
   */
2906
840k
  
if (840k
!s_pad(v, n+1)840k
)
return MP_MEMORY0
;
v->digits[n] = 0;840k
2907
840k
2908
840k
  /* Initialize temporary variables q and t.
2909
840k
     q allocates space for m+1 digits to store the quotient digits
2910
840k
     t allocates space for n+1 digits to hold the result of q_j*v
2911
840k
   */
2912
840k
  if (
(res = mp_int_init_size(&q, m + 1)) != MP_OK840k
)
return res0
;
2913
840k
  
if (840k
(res = mp_int_init_size(&t, n + 1)) != MP_OK840k
)
goto CLEANUP0
;
2914
840k
2915
840k
  /* D2: Initialize j */
2916
840k
  j = m;
2917
840k
  r.digits = MP_DIGITS(u) + j;  /* The contents of r are shared with u */
2918
840k
  r.used   = n + 1;
2919
840k
  r.sign   = MP_ZPOS;
2920
840k
  r.alloc  = MP_ALLOC(u);
2921
840k
  ZERO(t.digits, t.alloc);
2922
840k
2923
840k
  /* Calculate the m+1 digits of the quotient result */
2924
2.76M
  for (; 
j >= 02.76M
;
j--1.92M
)
{1.92M
2925
1.92M
    /* D3: Calculate q' */
2926
1.92M
    /* r->digits is aligned to position j of the number u */
2927
1.92M
    mp_word pfx, qhat;
2928
1.92M
    pfx   = r.digits[n];
2929
1.92M
    pfx <<= MP_DIGIT_BIT / 2;
2930
1.92M
    pfx <<= MP_DIGIT_BIT / 2;
2931
1.92M
    pfx |= r.digits[n-1]; /* pfx = u_{j+n}{j+n-1} */
2932
1.92M
2933
1.92M
    qhat = pfx / v->digits[n-1];
2934
1.92M
    /* Check to see if qhat > b, and decrease qhat if so.
2935
1.92M
       Theorem B guarantess that qhat is at most 2 larger than the
2936
1.92M
       actual value, so it is possible that qhat is greater than
2937
1.92M
       the maximum value that will fit in a digit */
2938
1.92M
    if (
qhat > 1.92M
MP_DIGIT_MAX1.92M
)
2939
86
      
qhat = 86
MP_DIGIT_MAX86
;
2940
1.92M
2941
1.92M
    /* D4,D5,D6: Multiply qhat * v and test for a correct value of q
2942
1.92M
2943
1.92M
       We proceed a bit different than the way described by Knuth. This way is
2944
1.92M
       simpler but less efficent. Instead of doing the multiply and subtract
2945
1.92M
       then checking for underflow, we first do the multiply of qhat * v and
2946
1.92M
       see if it is larger than the current remainder r. If it is larger, we
2947
1.92M
       decrease qhat by one and try again. We may need to decrease qhat one
2948
1.92M
       more time before we get a value that is smaller than r.
2949
1.92M
2950
1.92M
       This way is less efficent than Knuth becuase we do more multiplies, but
2951
1.92M
       we do not need to worry about underflow this way.
2952
1.92M
     */
2953
1.92M
    /* t = qhat * v */
2954
1.92M
    s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1); t.used = n + 1;
2955
1.92M
    CLAMP(&t);
2956
1.92M
2957
1.92M
    /* Clamp r for the comparison. Comparisons do not like leading zeros. */
2958
1.92M
    CLAMP(&r);
2959
1.92M
    if (
s_ucmp(&t, &r) > 01.92M
)
{ /* would the remainder be negative? */157k
2960
157k
      qhat -= 1;   /* try a smaller q */
2961
157k
      s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1);
2962
157k
      t.used = n + 1; CLAMP(&t);
2963
157k
      if (
s_ucmp(&t, &r) > 0157k
)
{ /* would the remainder be negative? */4.39k
2964
4.39k
        assert(qhat > 0);
2965
4.39k
        qhat -= 1; /* try a smaller q */
2966
4.39k
        s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1);
2967
4.39k
        t.used = n + 1; CLAMP(&t);
2968
157k
      }
2969
157k
      assert(s_ucmp(&t, &r) <=  0 && "The mathematics failed us.");
2970
1.92M
    }
2971
1.92M
    /* Unclamp r. The D algorithm expects r = u_{j+n}..u_j to always be n+1
2972
1.92M
       digits long. */
2973
1.92M
    r.used = n + 1;
2974
1.92M
2975
1.92M
    /* D4: Multiply and subtract
2976
1.92M
2977
1.92M
       Note: The multiply was completed above so we only need to subtract here.
2978
1.92M
     */
2979
1.92M
    s_usub(r.digits, t.digits, r.digits, r.used, t.used);
2980
1.92M
2981
1.92M
    /* D5: Test remainder
2982
1.92M
2983
1.92M
       Note: Not needed because we always check that qhat is the correct value
2984
1.92M
             before performing the subtract.  Value cast to mp_digit to prevent
2985
1.92M
             warning, qhat has been clamped to MP_DIGIT_MAX
2986
1.92M
     */
2987
1.92M
    q.digits[j] = (mp_digit)qhat;
2988
1.92M
2989
1.92M
    /* D6: Add back
2990
1.92M
       Note: Not needed because we always check that qhat is the correct value
2991
1.92M
             before performing the subtract.
2992
1.92M
     */
2993
1.92M
2994
1.92M
    /* D7: Loop on j */
2995
1.92M
    r.digits--;
2996
1.92M
    ZERO(t.digits, t.alloc);
2997
1.92M
  }
2998
840k
2999
840k
  /* Get rid of leading zeros in q */
3000
840k
  q.used = m + 1;
3001
840k
  CLAMP(&q);
3002
840k
3003
840k
  /* Denormalize the remainder */
3004
840k
  CLAMP(u); /* use u here because the r.digits pointer is off-by-one */
3005
840k
  if (k != 0)
3006
825k
    s_qdiv(u, k);
3007
840k
3008
840k
  mp_int_copy(u, v);  /* ok:  0 <= r < v */
3009
840k
  mp_int_copy(&q, u); /* ok:  q <= u     */
3010
840k
3011
840k
  mp_int_clear(&t);
3012
840k
 CLEANUP:
3013
840k
  mp_int_clear(&q);
3014
840k
  return res;
3015
1.64M
}
3016
3017
STATIC int       s_outlen(mp_int z, mp_size r)
3018
1.52k
{
3019
1.52k
  mp_result bits;
3020
1.52k
  double raw;
3021
1.52k
3022
1.52k
  assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
3023
1.52k
3024
1.52k
  bits = mp_int_count_bits(z);
3025
1.52k
  raw = (double)bits * s_log2[r];
3026
1.52k
3027
1.52k
  return (int)(raw + 0.999999);
3028
1.52k
}
3029
3030
STATIC mp_size   s_inlen(int len, mp_size r)
3031
148
{
3032
148
  double  raw = (double)len / s_log2[r];
3033
148
  mp_size bits = (mp_size)(raw + 0.5);
3034
148
3035
148
  return (mp_size)((bits + (
MP_DIGIT_BIT148
- 1)) /
MP_DIGIT_BIT148
) + 1;
3036
148
}
3037
3038
STATIC int       s_ch2val(char c, int r)
3039
1.80k
{
3040
1.80k
  int out;
3041
1.80k
3042
1.80k
  if (isdigit((unsigned char) c))
3043
1.80k
    out = c - '0';
3044
0
  else 
if (0
r > 10 && 0
isalpha((unsigned char) c)0
)
3045
0
    out = toupper(c) - 'A' + 10;
3046
0
  else
3047
0
    return -1;
3048
1.80k
3049
1.80k
  
return (out >= r) ? 1.80k
-10
:
out1.80k
;
3050
1.80k
}
3051
3052
STATIC char      s_val2ch(int v, int caps)
3053
26.3k
{
3054
26.3k
  assert(v >= 0);
3055
26.3k
3056
26.3k
  if (v < 10)
3057
26.3k
    return v + '0';
3058
26.3k
  else {
3059
0
    char out = (v - 10) + 'a';
3060
0
3061
0
    if (caps)
3062
0
      return toupper(out);
3063
0
    else
3064
0
      return out;
3065
26.3k
  }
3066
26.3k
}
3067
3068
STATIC void      s_2comp(unsigned char *buf, int len)
3069
0
{
3070
0
  int i;
3071
0
  unsigned short s = 1;
3072
0
3073
0
  for (i = len - 1; 
i >= 00
;
--i0
)
{0
3074
0
    unsigned char c = ~buf[i];
3075
0
3076
0
    s = c + s;
3077
0
    c = s & UCHAR_MAX;
3078
0
    s >>= CHAR_BIT;
3079
0
3080
0
    buf[i] = c;
3081
0
  }
3082
0
3083
0
  /* last carry out is ignored */
3084
0
}
3085
3086
STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
3087
0
{
3088
0
  mp_size uz;
3089
0
  mp_digit *dz;
3090
0
  int pos = 0, limit = *limpos;
3091
0
3092
0
  uz = 
MP_USED0
(z); dz =
MP_DIGITS0
(z);
3093
0
  while (
uz > 0 && 0
pos < limit0
)
{0
3094
0
    mp_digit d = *dz++;
3095
0
    int i;
3096
0
3097
0
    for (i = sizeof(mp_digit); 
i > 0 && 0
pos < limit0
;
--i0
)
{0
3098
0
      buf[pos++] = (unsigned char)d;
3099
0
      d >>= CHAR_BIT;
3100
0
3101
0
      /* Don't write leading zeroes */
3102
0
      if (
d == 0 && 0
uz == 10
)
3103
0
  i = 0; /* exit loop without signaling truncation */
3104
0
    }
3105
0
3106
0
    /* Detect truncation (loop exited with pos >= limit) */
3107
0
    if (
i > 00
)
break0
;
3108
0
3109
0
    --uz;
3110
0
  }
3111
0
3112
0
  if (
pad != 0 && 0
(buf[pos - 1] >> (CHAR_BIT - 1))0
)
{0
3113
0
    if (pos < limit)
3114
0
      buf[pos++] = 0;
3115
0
    else
3116
0
      uz = 1;
3117
0
  }
3118
0
3119
0
  /* Digits are in reverse order, fix that */
3120
0
  REV(unsigned char, buf, pos);
3121
0
3122
0
  /* Return the number of bytes actually written */
3123
0
  *limpos = pos;
3124
0
3125
0
  return (uz == 0) ? 
MP_OK0
:
MP_TRUNC0
;
3126
0
}
3127
3128
#if DEBUG
3129
void      s_print(char *tag, mp_int z)
3130
{
3131
  int  i;
3132
3133
  fprintf(stderr, "%s: %c ", tag,
3134
    (MP_SIGN(z) == MP_NEG) ? '-' : '+');
3135
3136
  for (i = MP_USED(z) - 1; i >= 0; --i)
3137
    fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), z->digits[i]);
3138
3139
  fputc('\n', stderr);
3140
3141
}
3142
3143
void      s_print_buf(char *tag, mp_digit *buf, mp_size num)
3144
{
3145
  int i;
3146
3147
  fprintf(stderr, "%s: ", tag);
3148
3149
  for (i = num - 1; i >= 0; --i)
3150
    fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), buf[i]);
3151
3152
  fputc('\n', stderr);
3153
}
3154
#endif
3155
3156
/* Here there be dragons */