Coverage Report

Created: 2018-08-14 02:14

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/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
450M
#define CHECK(TEST)   assert(TEST)
73
123M
#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
51.8M
#define MP_VALUE_DIGITS(V) \
103
51.8M
((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
104
105
/* Round precision P to nearest word boundary */
106
21.9M
#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
107
108
/* Set array P of S digits to zero */
109
10.4M
#define ZERO(P, S) \
110
10.4M
do{ \
111
10.4M
  mp_size i__ = (S) * sizeof(mp_digit); \
112
10.4M
  mp_digit *p__ = (P); \
113
10.4M
  memset(p__, 0, i__); \
114
10.4M
} while(0)
115
116
/* Copy S digits from array P to array Q */
117
133M
#define COPY(P, Q, S) \
118
133M
do{ \
119
120M
  mp_size i__ = (S) * sizeof(mp_digit); \
120
120M
  mp_digit *p__ = (P), *q__ = (Q); \
121
120M
  memcpy(q__, p__, i__); \
122
120M
} 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_) { \
129
0
    T xch = *u_; \
130
0
    *u_++ = *v_; \
131
0
    *v_-- = xch; \
132
0
  } \
133
0
} while(0)
134
135
183M
#define CLAMP(Z) \
136
183M
do{ \
137
183M
  mp_int z_ = (Z); \
138
183M
  mp_size uz_ = MP_USED(z_); \
139
183M
  mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \
140
217M
  while (uz_ > 1 && 
(*dz_-- == 0)148M
) \
141
183M
    
--uz_33.8M
; \
142
183M
  MP_USED(z_) = uz_; \
143
183M
} while(0)
144
145
/* Select min/max.  Do not provide expressions for which multiple
146
   evaluation would be problematic, e.g. x++ */
147
2.78M
#define MIN(A, B) ((B)<(A)?
(B)550k
:
(A)2.23M
)
148
107M
#define MAX(A, B) ((B)>(A)?
(B)47.1M
:
(A)60.1M
)
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.72M
#define SWAP(T, A, B) \
153
7.72M
do{ \
154
7.72M
  T t_ = (A); \
155
7.72M
  A = (B); \
156
7.72M
  B = t_; \
157
7.72M
} while(0)
158
159
/* Used to set up and access simple temp stacks within functions. */
160
#define DECLARE_TEMP(N) \
161
4.22M
  mpz_t temp[(N)]; \
162
4.22M
  int last__ = 0
163
#define CLEANUP_TEMP() \
164
2.70M
 CLEANUP: \
165
4.42M
  while (--last__ >= 0) \
166
2.70M
    
mp_int_clear(1.72M
TEMP1.72M
(last__))
167
3.44M
#define TEMP(K) (temp + (K))
168
1.72M
#define LAST_TEMP() TEMP(last__)
169
1.72M
#define SETUP(E) \
170
1.72M
do{ \
171
1.72M
  if ((res = (E)) != MP_OK) \
172
1.72M
    
goto CLEANUP0
; \
173
1.72M
  ++(last__); \
174
1.72M
} while(0)
175
176
/* Compare value to zero. */
177
243M
#define CMPZ(Z) \
178
243M
(((Z)->used==1&&
(Z)->digits[0]==092.7M
)?
04.64M
:
((Z)->sign==MP_NEG)?238M
-1203M
:
134.9M
)
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_USED(X), ub_ = MP_USED(Y); \
185
0
  mp_size o_ = ua_ + ub_; \
186
0
  ZERO(MP_DIGITS(Z), o_); \
187
0
  (void) s_kmul(MP_DIGITS(X), MP_DIGITS(Y), MP_DIGITS(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_DIGITS(X), MP_DIGITS(Z), ua_); \
199
0
  MP_USED(Z) = o_; \
200
0
  CLAMP(Z); \
201
0
} while(0)
202
203
273M
#define UPPER_HALF(W)           ((mp_word)((W) >> MP_DIGIT_BIT))
204
273M
#define LOWER_HALF(W)           ((mp_digit)(W))
205
0
#define HIGH_BIT_SET(W)         ((W) >> (MP_WORD_BIT - 1))
206
0
#define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (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
40.6M
{
361
40.6M
  if (z == NULL)
362
40.6M
    
return MP_BADARG0
;
363
40.6M
364
40.6M
  z->single = 0;
365
40.6M
  z->digits = &(z->single);
366
40.6M
  z->alloc  = 1;
367
40.6M
  z->used   = 1;
368
40.6M
  z->sign   = MP_ZPOS;
369
40.6M
370
40.6M
  return MP_OK;
371
40.6M
}
372
373
mp_int    mp_int_alloc(void)
374
33.4M
{
375
33.4M
  mp_int out = malloc(sizeof(mpz_t));
376
33.4M
377
33.4M
  if (out != NULL)
378
33.4M
    mp_int_init(out);
379
33.4M
380
33.4M
  return out;
381
33.4M
}
382
383
mp_result mp_int_init_size(mp_int z, mp_size prec)
384
6.18M
{
385
6.18M
  CHECK(z != NULL);
386
6.18M
387
6.18M
  if (prec == 0)
388
0
    prec = default_precision;
389
6.18M
  else if (prec == 1)
390
338k
    return mp_int_init(z);
391
5.84M
  else
392
5.84M
    prec = (mp_size) ROUND_PREC(prec);
393
6.18M
394
6.18M
  
if (5.84M
(5.84M
MP_DIGITS5.84M
(z) = s_alloc(prec)) == NULL)
395
5.84M
    
return MP_MEMORY0
;
396
5.84M
397
5.84M
  z->digits[0] = 0;
398
5.84M
  MP_USED(z) = 1;
399
5.84M
  MP_ALLOC(z) = prec;
400
5.84M
  MP_SIGN(z) = MP_ZPOS;
401
5.84M
402
5.84M
  return MP_OK;
403
5.84M
}
404
405
mp_result mp_int_init_copy(mp_int z, mp_int old)
406
7.54M
{
407
7.54M
  mp_result res;
408
7.54M
  mp_size uold;
409
7.54M
410
7.54M
  CHECK(z != NULL && old != NULL);
411
7.54M
412
7.54M
  uold = MP_USED(old);
413
7.54M
  if (uold == 1) {
414
3.06M
    mp_int_init(z);
415
3.06M
  }
416
4.47M
  else {
417
4.47M
    mp_size target = MAX(uold, default_precision);
418
4.47M
419
4.47M
    if ((res = mp_int_init_size(z, target)) != MP_OK)
420
0
      return res;
421
7.54M
  }
422
7.54M
423
7.54M
  MP_USED(z) = uold;
424
7.54M
  MP_SIGN(z) = MP_SIGN(old);
425
7.54M
  COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
426
7.54M
427
7.54M
  return MP_OK;
428
7.54M
}
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
168k
{
441
168k
  mpz_t    vtmp;
442
168k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
443
168k
444
168k
  s_ufake(&vtmp, uvalue, vbuf);
445
168k
  return mp_int_init_copy(z, &vtmp);
446
168k
}
447
448
mp_result  mp_int_set_value(mp_int z, mp_small value)
449
972k
{
450
972k
  mpz_t    vtmp;
451
972k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
452
972k
453
972k
  s_fake(&vtmp, value, vbuf);
454
972k
  return mp_int_copy(&vtmp, z);
455
972k
}
456
457
mp_result  mp_int_set_uvalue(mp_int z, mp_usmall uvalue)
458
141k
{
459
141k
  mpz_t    vtmp;
460
141k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
461
141k
462
141k
  s_ufake(&vtmp, uvalue, vbuf);
463
141k
  return mp_int_copy(&vtmp, z);
464
141k
}
465
466
void      mp_int_clear(mp_int z)
467
46.4M
{
468
46.4M
  if (z == NULL)
469
46.4M
    
return0
;
470
46.4M
471
46.4M
  if (MP_DIGITS(z) != NULL) {
472
46.4M
    if (MP_DIGITS(z) != &(z->single))
473
19.8M
      s_free(MP_DIGITS(z));
474
46.4M
475
46.4M
    MP_DIGITS(z) = NULL;
476
46.4M
  }
477
46.4M
}
478
479
void      mp_int_free(mp_int z)
480
33.4M
{
481
33.4M
  NRCHECK(z != NULL);
482
33.4M
483
33.4M
  mp_int_clear(z);
484
33.4M
  free(z); /* note: NOT s_free() */
485
33.4M
}
486
487
mp_result mp_int_copy(mp_int a, mp_int c)
488
104M
{
489
104M
  CHECK(a != NULL && c != NULL);
490
104M
491
104M
  if (a != c) {
492
98.9M
    mp_size ua = MP_USED(a);
493
98.9M
    mp_digit *da, *dc;
494
98.9M
495
98.9M
    if (!s_pad(c, ua))
496
0
      return MP_MEMORY;
497
98.9M
498
98.9M
    da = MP_DIGITS(a); dc = MP_DIGITS(c);
499
98.9M
    COPY(da, dc, ua);
500
98.9M
501
98.9M
    MP_USED(c) = ua;
502
98.9M
    MP_SIGN(c) = MP_SIGN(a);
503
98.9M
  }
504
104M
505
104M
  return MP_OK;
506
104M
}
507
508
void      mp_int_swap(mp_int a, mp_int c)
509
0
{
510
0
  if (a != c) {
511
0
    mpz_t tmp = *a;
512
0
513
0
    *a = *c;
514
0
    *c = tmp;
515
0
516
0
    if (MP_DIGITS(a) == &(c->single))
517
0
      MP_DIGITS(a) = &(a->single);
518
0
    if (MP_DIGITS(c) == &(a->single))
519
0
      MP_DIGITS(c) = &(c->single);
520
0
  }
521
0
}
522
523
void      mp_int_zero(mp_int z)
524
21.2M
{
525
21.2M
  NRCHECK(z != NULL);
526
21.2M
527
21.2M
  z->digits[0] = 0;
528
21.2M
  MP_USED(z) = 1;
529
21.2M
  MP_SIGN(z) = MP_ZPOS;
530
21.2M
}
531
532
mp_result mp_int_abs(mp_int a, mp_int c)
533
3.52M
{
534
3.52M
  mp_result res;
535
3.52M
536
3.52M
  CHECK(a != NULL && c != NULL);
537
3.52M
538
3.52M
  if ((res = mp_int_copy(a, c)) != MP_OK)
539
0
    return res;
540
3.52M
541
3.52M
  MP_SIGN(c) = MP_ZPOS;
542
3.52M
  return MP_OK;
543
3.52M
}
544
545
mp_result mp_int_neg(mp_int a, mp_int c)
546
71.1M
{
547
71.1M
  mp_result res;
548
71.1M
549
71.1M
  CHECK(a != NULL && c != NULL);
550
71.1M
551
71.1M
  if ((res = mp_int_copy(a, c)) != MP_OK)
552
0
    return res;
553
71.1M
554
71.1M
  if (CMPZ(c) != 0)
555
71.1M
    
MP_SIGN71.0M
(c) = 1 - 71.0M
MP_SIGN71.0M
(a);
556
71.1M
557
71.1M
  return MP_OK;
558
71.1M
}
559
560
mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
561
14.0M
{
562
14.0M
  mp_size ua, ub, uc, max;
563
14.0M
564
14.0M
  CHECK(a != NULL && b != NULL && c != NULL);
565
14.0M
566
14.0M
  ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
567
14.0M
  max = MAX(ua, ub);
568
14.0M
569
14.0M
  if (MP_SIGN(a) == MP_SIGN(b)) {
570
8.12M
    /* Same sign -- add magnitudes, preserve sign of addends */
571
8.12M
    mp_digit carry;
572
8.12M
573
8.12M
    if (!s_pad(c, max))
574
0
      return MP_MEMORY;
575
8.12M
576
8.12M
    carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
577
8.12M
    uc = max;
578
8.12M
579
8.12M
    if (carry) {
580
49.8k
      if (!s_pad(c, max + 1))
581
0
  return MP_MEMORY;
582
49.8k
583
49.8k
      c->digits[max] = carry;
584
49.8k
      ++uc;
585
49.8k
    }
586
8.12M
587
8.12M
    MP_USED(c) = uc;
588
8.12M
    MP_SIGN(c) = MP_SIGN(a);
589
8.12M
590
8.12M
  }
591
5.89M
  else {
592
5.89M
    /* Different signs -- subtract magnitudes, preserve sign of greater */
593
5.89M
    mp_int  x, y;
594
5.89M
    int     cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
595
5.89M
596
5.89M
    /* Set x to max(a, b), y to min(a, b) to simplify later code.
597
5.89M
       A special case yields zero for equal magnitudes.
598
5.89M
    */
599
5.89M
    if (cmp == 0) {
600
819k
      mp_int_zero(c);
601
819k
      return MP_OK;
602
819k
    }
603
5.07M
    else if (cmp < 0) {
604
1.36M
      x = b; y = a;
605
1.36M
    }
606
3.70M
    else {
607
3.70M
      x = a; y = b;
608
3.70M
    }
609
5.89M
610
5.89M
    
if (5.07M
!s_pad(c, 5.07M
MP_USED5.07M
(x)))
611
0
      return MP_MEMORY;
612
5.07M
613
5.07M
    /* Subtract smaller from larger */
614
5.07M
    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
615
5.07M
    MP_USED(c) = MP_USED(x);
616
5.07M
    CLAMP(c);
617
5.07M
618
5.07M
    /* Give result the sign of the larger */
619
5.07M
    MP_SIGN(c) = MP_SIGN(x);
620
5.07M
  }
621
14.0M
622
14.0M
  
return MP_OK13.1M
;
623
14.0M
}
624
625
mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c)
626
3.38k
{
627
3.38k
  mpz_t    vtmp;
628
3.38k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
629
3.38k
630
3.38k
  s_fake(&vtmp, value, vbuf);
631
3.38k
632
3.38k
  return mp_int_add(a, &vtmp, c);
633
3.38k
}
634
635
mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
636
79.5M
{
637
79.5M
  mp_size ua, ub, uc, max;
638
79.5M
639
79.5M
  CHECK(a != NULL && b != NULL && c != NULL);
640
79.5M
641
79.5M
  ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
642
79.5M
  max = MAX(ua, ub);
643
79.5M
644
79.5M
  if (MP_SIGN(a) != MP_SIGN(b)) {
645
267k
    /* Different signs -- add magnitudes and keep sign of a */
646
267k
    mp_digit carry;
647
267k
648
267k
    if (!s_pad(c, max))
649
0
      return MP_MEMORY;
650
267k
651
267k
    carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
652
267k
    uc = max;
653
267k
654
267k
    if (carry) {
655
2.09k
      if (!s_pad(c, max + 1))
656
0
  return MP_MEMORY;
657
2.09k
658
2.09k
      c->digits[max] = carry;
659
2.09k
      ++uc;
660
2.09k
    }
661
267k
662
267k
    MP_USED(c) = uc;
663
267k
    MP_SIGN(c) = MP_SIGN(a);
664
267k
665
267k
  }
666
79.2M
  else {
667
79.2M
    /* Same signs -- subtract magnitudes */
668
79.2M
    mp_int  x, y;
669
79.2M
    mp_sign osign;
670
79.2M
    int     cmp = s_ucmp(a, b);
671
79.2M
672
79.2M
    if (!s_pad(c, max))
673
0
      return MP_MEMORY;
674
79.2M
675
79.2M
    if (cmp >= 0) {
676
13.4M
      x = a; y = b; osign = MP_ZPOS;
677
13.4M
    }
678
65.8M
    else {
679
65.8M
      x = b; y = a; osign = MP_NEG;
680
65.8M
    }
681
79.2M
682
79.2M
    if (MP_SIGN(a) == MP_NEG && 
cmp != 0110k
)
683
92.2k
      osign = 1 - osign;
684
79.2M
685
79.2M
    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
686
79.2M
    MP_USED(c) = MP_USED(x);
687
79.2M
    CLAMP(c);
688
79.2M
689
79.2M
    MP_SIGN(c) = osign;
690
79.2M
  }
691
79.5M
692
79.5M
  return MP_OK;
693
79.5M
}
694
695
mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c)
696
15.6k
{
697
15.6k
  mpz_t    vtmp;
698
15.6k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
699
15.6k
700
15.6k
  s_fake(&vtmp, value, vbuf);
701
15.6k
702
15.6k
  return mp_int_sub(a, &vtmp, c);
703
15.6k
}
704
705
mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
706
26.4M
{
707
26.4M
  mp_digit *out;
708
26.4M
  mp_size   osize, ua, ub, p = 0;
709
26.4M
  mp_sign   osign;
710
26.4M
711
26.4M
  CHECK(a != NULL && b != NULL && c != NULL);
712
26.4M
713
26.4M
  /* If either input is zero, we can shortcut multiplication */
714
26.4M
  if (mp_int_compare_zero(a) == 0 || 
mp_int_compare_zero(b) == 020.7M
) {
715
18.8M
    mp_int_zero(c);
716
18.8M
    return MP_OK;
717
18.8M
  }
718
7.57M
719
7.57M
  /* Output is positive if inputs have same sign, otherwise negative */
720
7.57M
  osign = (MP_SIGN(a) == MP_SIGN(b)) ? 
MP_ZPOS4.56M
:
MP_NEG3.01M
;
721
7.57M
722
7.57M
  /* If the output is not identical to any of the inputs, we'll write the
723
7.57M
     results directly; otherwise, allocate a temporary space. */
724
7.57M
  ua = MP_USED(a); ub = MP_USED(b);
725
7.57M
  osize = MAX(ua, ub);
726
7.57M
  osize = 4 * ((osize + 1) / 2);
727
7.57M
728
7.57M
  if (c == a || 
c == b5.86M
) {
729
1.70M
    p = ROUND_PREC(osize);
730
1.70M
    p = MAX(p, default_precision);
731
1.70M
732
1.70M
    if ((out = s_alloc(p)) == NULL)
733
1.70M
      
return MP_MEMORY0
;
734
5.86M
  }
735
5.86M
  else {
736
5.86M
    if (!s_pad(c, osize))
737
0
      return MP_MEMORY;
738
5.86M
739
5.86M
    out = MP_DIGITS(c);
740
5.86M
  }
741
7.57M
  ZERO(out, osize);
742
7.57M
743
7.57M
  if (!s_kmul(MP_DIGITS(a), MP_DIGITS(b), out, ua, ub))
744
0
    return MP_MEMORY;
745
7.57M
746
7.57M
  /* If we allocated a new buffer, get rid of whatever memory c was already
747
7.57M
     using, and fix up its fields to reflect that.
748
7.57M
   */
749
7.57M
  if (out != MP_DIGITS(c)) {
750
1.70M
    if ((void *) MP_DIGITS(c) != (void *) c)
751
1.45M
      s_free(MP_DIGITS(c));
752
1.70M
    MP_DIGITS(c) = out;
753
1.70M
    MP_ALLOC(c) = p;
754
1.70M
  }
755
7.57M
756
7.57M
  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
757
7.57M
  CLAMP(c);           /* ... right here */
758
7.57M
  MP_SIGN(c) = osign;
759
7.57M
760
7.57M
  return MP_OK;
761
7.57M
}
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.83k
{
775
3.83k
  mp_result res;
776
3.83k
  CHECK(a != NULL && c != NULL && p2 >= 0);
777
3.83k
778
3.83k
  if ((res = mp_int_copy(a, c)) != MP_OK)
779
0
    return res;
780
3.83k
781
3.83k
  if (s_qmul(c, (mp_size) p2))
782
3.83k
    return MP_OK;
783
0
  else
784
0
    return MP_MEMORY;
785
3.83k
}
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 == c) {
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 = MP_DIGITS(c);
808
0
  }
809
0
  ZERO(out, osize);
810
0
811
0
  s_ksqr(MP_DIGITS(a), out, MP_USED(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 != MP_DIGITS(c)) {
817
0
    if ((void *) MP_DIGITS(c) != (void *) c)
818
0
      s_free(MP_DIGITS(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
4.22M
{
832
4.22M
  int cmp, lg;
833
4.22M
  mp_result res = MP_OK;
834
4.22M
  mp_int qout, rout;
835
4.22M
  mp_sign sa = MP_SIGN(a), sb = MP_SIGN(b);
836
4.22M
  DECLARE_TEMP(2);
837
4.22M
838
4.22M
  CHECK(a != NULL && b != NULL && q != r);
839
4.22M
840
4.22M
  if (CMPZ(b) == 0)
841
0
    return MP_UNDEF;
842
4.22M
  else if ((cmp = s_ucmp(a, b)) < 0) {
843
795k
    /* If |a| < |b|, no division is required:
844
795k
       q = 0, r = a
845
795k
     */
846
795k
    if (r && 
(res = mp_int_copy(a, r)) != MP_OK85.4k
)
847
0
      return res;
848
795k
849
795k
    if (q)
850
786k
      mp_int_zero(q);
851
795k
852
795k
    return MP_OK;
853
795k
  }
854
3.42M
  else if (cmp == 0) {
855
723k
    /* If |a| = |b|, no division is required:
856
723k
       q = 1 or -1, r = 0
857
723k
     */
858
723k
    if (r)
859
5.41k
      mp_int_zero(r);
860
723k
861
723k
    if (q) {
862
720k
      mp_int_zero(q);
863
720k
      q->digits[0] = 1;
864
720k
865
720k
      if (sa != sb)
866
720k
  
MP_SIGN86.3k
(q) = MP_NEG86.3k
;
867
720k
    }
868
723k
869
723k
    return MP_OK;
870
723k
  }
871
2.70M
872
2.70M
  /* When |a| > |b|, real division is required.  We need someplace to store
873
2.70M
     quotient and remainder, but q and r are allowed to be NULL or to overlap
874
2.70M
     with the inputs.
875
2.70M
   */
876
2.70M
  if ((lg = s_isp2(b)) < 0) {
877
1.73M
    if (q && 
b != q1.71M
) {
878
1.61M
      if ((res = mp_int_copy(a, q)) != MP_OK)
879
0
  goto CLEANUP;
880
1.61M
      else
881
1.61M
  qout = q;
882
1.61M
    }
883
117k
    else {
884
117k
      qout = LAST_TEMP();
885
117k
      SETUP(mp_int_init_copy(LAST_TEMP(), a));
886
117k
    }
887
1.73M
888
1.73M
    if (r && 
a != r125k
) {
889
125k
      if ((res = mp_int_copy(b, r)) != MP_OK)
890
0
  goto CLEANUP;
891
125k
      else
892
125k
  rout = r;
893
125k
    }
894
1.60M
    else {
895
1.60M
      rout = LAST_TEMP();
896
1.60M
      SETUP(mp_int_init_copy(LAST_TEMP(), b));
897
1.60M
    }
898
1.73M
899
1.73M
    if ((res = s_udiv_knuth(qout, rout)) != MP_OK) 
goto CLEANUP0
;
900
970k
  }
901
970k
  else {
902
970k
    if (q && 
(res = mp_int_copy(a, q)) != MP_OK923k
)
goto CLEANUP0
;
903
970k
    if (r && 
(res = mp_int_copy(a, r)) != MP_OK140k
)
goto CLEANUP0
;
904
970k
905
970k
    if (q) 
s_qdiv(q, (mp_size) lg)923k
; qout = q;
906
970k
    if (r) 
s_qmod(r, (mp_size) lg)140k
; rout = r;
907
970k
  }
908
2.70M
909
2.70M
  /* Recompute signs for output */
910
2.70M
  if (rout) {
911
1.87M
    MP_SIGN(rout) = sa;
912
1.87M
    if (CMPZ(rout) == 0)
913
1.87M
      
MP_SIGN1.78M
(rout) = MP_ZPOS1.78M
;
914
1.87M
  }
915
2.70M
  if (qout) {
916
2.65M
    MP_SIGN(qout) = (sa == sb) ? 
MP_ZPOS1.87M
:
MP_NEG779k
;
917
2.65M
    if (CMPZ(qout) == 0)
918
2.65M
      
MP_SIGN0
(qout) = MP_ZPOS0
;
919
2.65M
  }
920
2.70M
921
2.70M
  if (q && 
(res = mp_int_copy(qout, q)) != MP_OK2.64M
)
goto CLEANUP0
;
922
2.70M
  if (r && 
(res = mp_int_copy(rout, r)) != MP_OK265k
)
goto CLEANUP0
;
923
2.70M
924
2.70M
  CLEANUP_TEMP();
925
2.70M
  return res;
926
2.70M
}
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 == c) {
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 (CMPZ(out) < 0)
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
65.1k
{
959
65.1k
  mpz_t     vtmp, rtmp;
960
65.1k
  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
961
65.1k
  mp_result res;
962
65.1k
963
65.1k
  mp_int_init(&rtmp);
964
65.1k
  s_fake(&vtmp, value, vbuf);
965
65.1k
966
65.1k
  if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
967
0
    goto CLEANUP;
968
65.1k
969
65.1k
  if (r)
970
54.1k
    (void) mp_int_to_int(&rtmp, r); /* can't fail */
971
65.1k
972
65.1k
 CLEANUP:
973
65.1k
  mp_int_clear(&rtmp);
974
65.1k
  return res;
975
65.1k
}
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 && (res = mp_int_copy(a, q)) == MP_OK)
984
0
    s_qdiv(q, (mp_size) p2);
985
0
986
0
  if (res == MP_OK && r != NULL && (res = mp_int_copy(a, r)) == MP_OK)
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 ((res = mp_int_init_copy(&t, a)) != MP_OK)
1003
0
    return res;
1004
0
1005
0
  (void) mp_int_set_value(c, 1);
1006
0
  while (v != 0) {
1007
0
    if (v & 1) {
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 == 0) break;
1014
0
1015
0
    if ((res = mp_int_sqr(&t, &t)) != MP_OK)
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 ((res = mp_int_init_value(&t, a)) != MP_OK)
1035
0
    return res;
1036
0
1037
0
  (void) mp_int_set_value(c, 1);
1038
0
  while (v != 0) {
1039
0
    if (v & 1) {
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 == 0) break;
1046
0
1047
0
    if ((res = mp_int_sqr(&t, &t)) != MP_OK)
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_SIGN(b) == MP_NEG)
1064
0
    return MP_RANGE;
1065
0
1066
0
  if ((res = mp_int_init_copy(&t, a)) != MP_OK)
1067
0
    return res;
1068
0
1069
0
  (void) mp_int_set_value(c, 1);
1070
0
  for (ix = 0; ix < MP_USED(b); ++ix) {
1071
0
    mp_digit d = b->digits[ix];
1072
0
1073
0
    for (jx = 0; jx < MP_DIGIT_BIT; ++jx) {
1074
0
      if (d & 1) {
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 && ix + 1 == MP_USED(b))
1081
0
  break;
1082
0
      if ((res = mp_int_sqr(&t, &t)) != MP_OK)
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
269k
{
1094
269k
  mp_sign sa;
1095
269k
1096
269k
  CHECK(a != NULL && b != NULL);
1097
269k
1098
269k
  sa = MP_SIGN(a);
1099
269k
  if (sa == MP_SIGN(b)) {
1100
235k
    int cmp = s_ucmp(a, b);
1101
235k
1102
235k
    /* If they're both zero or positive, the normal comparison applies; if both
1103
235k
       negative, the sense is reversed. */
1104
235k
    if (sa == MP_ZPOS)
1105
161k
      return cmp;
1106
74.0k
    else
1107
74.0k
      return -cmp;
1108
34.3k
1109
34.3k
  }
1110
34.3k
  else {
1111
34.3k
    if (sa == MP_ZPOS)
1112
15.6k
      return 1;
1113
18.6k
    else
1114
18.6k
      return -1;
1115
34.3k
  }
1116
269k
}
1117
1118
int       mp_int_compare_unsigned(mp_int a, mp_int b)
1119
3.14M
{
1120
3.14M
  NRCHECK(a != NULL && b != NULL);
1121
3.14M
1122
3.14M
  return s_ucmp(a, b);
1123
3.14M
}
1124
1125
int       mp_int_compare_zero(mp_int z)
1126
65.7M
{
1127
65.7M
  NRCHECK(z != NULL);
1128
65.7M
1129
65.7M
  if (MP_USED(z) == 1 && 
z->digits[0] == 038.3M
)
1130
19.1M
    return 0;
1131
46.5M
  else if (MP_SIGN(z) == MP_ZPOS)
1132
33.0M
    return 1;
1133
13.4M
  else
1134
13.4M
    return -1;
1135
65.7M
}
1136
1137
int       mp_int_compare_value(mp_int z, mp_small value)
1138
83.2M
{
1139
83.2M
  mp_sign vsign = (value < 0) ? 
MP_NEG41.2M
:
MP_ZPOS41.9M
;
1140
83.2M
  int cmp;
1141
83.2M
1142
83.2M
  CHECK(z != NULL);
1143
83.2M
1144
83.2M
  if (vsign == MP_SIGN(z)) {
1145
50.4M
    cmp = s_vcmp(z, value);
1146
50.4M
1147
50.4M
    return (vsign == MP_ZPOS) ? 
cmp41.9M
:
-cmp8.56M
;
1148
50.4M
  }
1149
32.7M
  else {
1150
32.7M
    return (value < 0) ? 
132.7M
:
-132.5k
;
1151
32.7M
  }
1152
83.2M
}
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_SIGN(z) == MP_NEG)
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 (CMPZ(m) == 0)
1175
0
    return MP_UNDEF;
1176
0
  if (CMPZ(b) < 0)
1177
0
    return MP_RANGE;
1178
0
1179
0
  um = MP_USED(m);
1180
0
  SETUP(mp_int_init_size(TEMP(0), 2 * um));
1181
0
  SETUP(mp_int_init_size(TEMP(1), 2 * um));
1182
0
1183
0
  if (c == b || c == m) {
1184
0
    SETUP(mp_int_init_size(TEMP(2), 2 * um));
1185
0
    s = TEMP(2);
1186
0
  }
1187
0
  else {
1188
0
    s = c;
1189
0
  }
1190
0
1191
0
  if ((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;
1192
0
1193
0
  if ((res = s_brmu(TEMP(1), m)) != MP_OK) goto CLEANUP;
1194
0
1195
0
  if ((res = s_embar(TEMP(0), b, m, TEMP(1), s)) != MP_OK)
1196
0
    goto CLEANUP;
1197
0
1198
0
  res = mp_int_copy(s, c);
1199
0
1200
0
  CLEANUP_TEMP();
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 (CMPZ(m) == 0)
1236
0
    return MP_UNDEF;
1237
0
  if (CMPZ(b) < 0)
1238
0
    return MP_RANGE;
1239
0
1240
0
  um = MP_USED(m);
1241
0
  SETUP(mp_int_init_size(TEMP(0), 2 * um));
1242
0
1243
0
  if (c == b || c == m) {
1244
0
    SETUP(mp_int_init_size(TEMP(1), 2 * um));
1245
0
    s = TEMP(1);
1246
0
  }
1247
0
  else {
1248
0
    s = c;
1249
0
  }
1250
0
1251
0
  if ((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;
1252
0
1253
0
  if ((res = s_embar(TEMP(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_TEMP();
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 (CMPZ(a) == 0 || CMPZ(m) <= 0)
1278
0
    return MP_RANGE;
1279
0
1280
0
  sa = MP_SIGN(a); /* need this for the result later */
1281
0
1282
0
  for (last__ = 0; last__ < 2; ++last__)
1283
0
    mp_int_init(LAST_TEMP());
1284
0
1285
0
  if ((res = mp_int_egcd(a, m, TEMP(0), TEMP(1), NULL)) != MP_OK)
1286
0
    goto CLEANUP;
1287
0
1288
0
  if (mp_int_compare_value(TEMP(0), 1) != 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 ((res = mp_int_mod(TEMP(1), m, TEMP(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 (sa == MP_NEG)
1303
0
    res = mp_int_sub(m, TEMP(1), c);
1304
0
  else
1305
0
    res = mp_int_copy(TEMP(1), c);
1306
0
1307
0
  CLEANUP_TEMP();
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.79M
{
1314
2.79M
  int ca, cb, k = 0;
1315
2.79M
  mpz_t u, v, t;
1316
2.79M
  mp_result res;
1317
2.79M
1318
2.79M
  CHECK(a != NULL && b != NULL && c != NULL);
1319
2.79M
1320
2.79M
  ca = CMPZ(a);
1321
2.79M
  cb = CMPZ(b);
1322
2.79M
  if (ca == 0 && 
cb == 018
)
1323
0
    return MP_UNDEF;
1324
2.79M
  else if (ca == 0)
1325
18
    return mp_int_abs(b, c);
1326
2.79M
  else if (cb == 0)
1327
9.43k
    return mp_int_abs(a, c);
1328
2.78M
1329
2.78M
  mp_int_init(&t);
1330
2.78M
  if ((res = mp_int_init_copy(&u, a)) != MP_OK)
1331
0
    goto U;
1332
2.78M
  if ((res = mp_int_init_copy(&v, b)) != MP_OK)
1333
0
    goto V;
1334
2.78M
1335
2.78M
  MP_SIGN(&u) = MP_ZPOS; MP_SIGN(&v) = MP_ZPOS;
1336
2.78M
1337
2.78M
  { /* Divide out common factors of 2 from u and v */
1338
2.78M
    int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
1339
2.78M
1340
2.78M
    k = MIN(div2_u, div2_v);
1341
2.78M
    s_qdiv(&u, (mp_size) k);
1342
2.78M
    s_qdiv(&v, (mp_size) k);
1343
2.78M
  }
1344
2.78M
1345
2.78M
  if (mp_int_is_odd(&u)) {
1346
2.23M
    if ((res = mp_int_neg(&v, &t)) != MP_OK)
1347
0
      goto CLEANUP;
1348
550k
  }
1349
550k
  else {
1350
550k
    if ((res = mp_int_copy(&u, &t)) != MP_OK)
1351
0
      goto CLEANUP;
1352
2.78M
  }
1353
2.78M
1354
78.9M
  
for (;;)2.78M
{
1355
78.9M
    s_qdiv(&t, s_dp2k(&t));
1356
78.9M
1357
78.9M
    if (CMPZ(&t) > 0) {
1358
10.9M
      if ((res = mp_int_copy(&t, &u)) != MP_OK)
1359
0
  goto CLEANUP;
1360
67.9M
    }
1361
67.9M
    else {
1362
67.9M
      if ((res = mp_int_neg(&t, &v)) != MP_OK)
1363
0
  goto CLEANUP;
1364
78.9M
    }
1365
78.9M
1366
78.9M
    if ((res = mp_int_sub(&u, &v, &t)) != MP_OK)
1367
0
      goto CLEANUP;
1368
78.9M
1369
78.9M
    if (CMPZ(&t) == 0)
1370
2.78M
      break;
1371
78.9M
  }
1372
2.78M
1373
2.78M
  if ((res = mp_int_abs(&u, c)) != MP_OK)
1374
0
    goto CLEANUP;
1375
2.78M
  if (!s_qmul(c, (mp_size) k))
1376
0
    res = MP_MEMORY;
1377
2.78M
1378
2.78M
 CLEANUP:
1379
2.78M
  mp_int_clear(&v);
1380
2.78M
 V: mp_int_clear(&u);
1381
2.78M
 U: mp_int_clear(&t);
1382
2.78M
1383
2.78M
  return res;
1384
2.78M
}
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 && cb == 0)
1403
0
    return MP_UNDEF;
1404
0
  else if (ca == 0) {
1405
0
    if ((res = mp_int_abs(b, c)) != MP_OK) return res;
1406
0
    mp_int_zero(x); (void) mp_int_set_value(y, 1); return MP_OK;
1407
0
  }
1408
0
  else if (cb == 0) {
1409
0
    if ((res = mp_int_abs(a, c)) != MP_OK) return res;
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; last__ < 4; ++last__)
1416
0
    mp_int_init(LAST_TEMP());
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
  SETUP(mp_int_init_copy(TEMP(5), b));
1422
0
1423
0
  /* We will work with absolute values here */
1424
0
  MP_SIGN(TEMP(4)) = MP_ZPOS;
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(TEMP(4)), div2_v = s_dp2k(TEMP(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
  SETUP(mp_int_init_copy(TEMP(7), TEMP(5)));
1437
0
1438
0
  for (;;) {
1439
0
    while (mp_int_is_even(TEMP(4))) {
1440
0
      s_qdiv(TEMP(4), 1);
1441
0
1442
0
      if (mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) {
1443
0
  if ((res = mp_int_add(TEMP(0), TEMP(7), TEMP(0))) != MP_OK)
1444
0
    goto CLEANUP;
1445
0
  if ((res = mp_int_sub(TEMP(1), TEMP(6), TEMP(1))) != MP_OK)
1446
0
    goto CLEANUP;
1447
0
      }
1448
0
1449
0
      s_qdiv(TEMP(0), 1);
1450
0
      s_qdiv(TEMP(1), 1);
1451
0
    }
1452
0
1453
0
    while (mp_int_is_even(TEMP(5))) {
1454
0
      s_qdiv(TEMP(5), 1);
1455
0
1456
0
      if (mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) {
1457
0
  if ((res = mp_int_add(TEMP(2), TEMP(7), TEMP(2))) != MP_OK)
1458
0
    goto CLEANUP;
1459
0
  if ((res = mp_int_sub(TEMP(3), TEMP(6), TEMP(3))) != MP_OK)
1460
0
    goto CLEANUP;
1461
0
      }
1462
0
1463
0
      s_qdiv(TEMP(2), 1);
1464
0
      s_qdiv(TEMP(3), 1);
1465
0
    }
1466
0
1467
0
    if (mp_int_compare(TEMP(4), TEMP(5)) >= 0) {
1468
0
      if ((res = mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != MP_OK) goto CLEANUP;
1469
0
      if ((res = mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != MP_OK) goto CLEANUP;
1470
0
      if ((res = mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != MP_OK) goto CLEANUP;
1471
0
    }
1472
0
    else {
1473
0
      if ((res = mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != MP_OK) goto CLEANUP;
1474
0
      if ((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK) goto CLEANUP;
1475
0
      if ((res = mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != MP_OK) goto CLEANUP;
1476
0
    }
1477
0
1478
0
    if (CMPZ(TEMP(4)) == 0) {
1479
0
      if (x && (res = mp_int_copy(TEMP(2), x)) != MP_OK) goto CLEANUP;
1480
0
      if (y && (res = mp_int_copy(TEMP(3), y)) != MP_OK) goto CLEANUP;
1481
0
      if (c) {
1482
0
  if (!s_qmul(TEMP(5), k)) {
1483
0
    res = MP_MEMORY;
1484
0
    goto CLEANUP;
1485
0
  }
1486
0
1487
0
  res = mp_int_copy(TEMP(5), c);
1488
0
      }
1489
0
1490
0
      break;
1491
0
    }
1492
0
  }
1493
0
1494
0
  CLEANUP_TEMP();
1495
0
  return res;
1496
0
}
1497
1498
mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c)
1499
269k
{
1500
269k
  mpz_t lcm;
1501
269k
  mp_result res;
1502
269k
1503
269k
  CHECK(a != NULL && b != NULL && c != NULL);
1504
269k
1505
269k
  /* Since a * b = gcd(a, b) * lcm(a, b), we can compute
1506
269k
     lcm(a, b) = (a / gcd(a, b)) * b.
1507
269k
1508
269k
     This formulation insures everything works even if the input
1509
269k
     variables share space.
1510
269k
   */
1511
269k
  if ((res = mp_int_init(&lcm)) != MP_OK)
1512
0
    return res;
1513
269k
  if ((res = mp_int_gcd(a, b, &lcm)) != MP_OK)
1514
0
    goto CLEANUP;
1515
269k
  if ((res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK)
1516
0
    goto CLEANUP;
1517
269k
  if ((res = mp_int_mul(&lcm, b, &lcm)) != MP_OK)
1518
0
    goto CLEANUP;
1519
269k
1520
269k
  res = mp_int_copy(&lcm, c);
1521
269k
1522
269k
  CLEANUP:
1523
269k
    mp_int_clear(&lcm);
1524
269k
1525
269k
  return res;
1526
269k
}
1527
1528
int       mp_int_divisible_value(mp_int a, mp_small v)
1529
54.1k
{
1530
54.1k
  mp_small rem = 0;
1531
54.1k
1532
54.1k
  if (mp_int_div_value(a, v, NULL, &rem) != MP_OK)
1533
0
    return 0;
1534
54.1k
1535
54.1k
  return rem == 0;
1536
54.1k
}
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 == 1) {
1558
0
    return mp_int_copy(a, c);
1559
0
  }
1560
0
  if (MP_SIGN(a) == MP_NEG) {
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
  SETUP(mp_int_init_copy(LAST_TEMP(), a));
1568
0
  SETUP(mp_int_init_copy(LAST_TEMP(), a));
1569
0
  SETUP(mp_int_init(LAST_TEMP()));
1570
0
  SETUP(mp_int_init(LAST_TEMP()));
1571
0
  SETUP(mp_int_init(LAST_TEMP()));
1572
0
1573
0
  (void) mp_int_abs(TEMP(0), TEMP(0));
1574
0
  (void) mp_int_abs(TEMP(1), TEMP(1));
1575
0
1576
0
  for (;;) {
1577
0
    if ((res = mp_int_expt(TEMP(1), b, TEMP(2))) != MP_OK)
1578
0
      goto CLEANUP;
1579
0
1580
0
    if (mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0)
1581
0
      break;
1582
0
1583
0
    if ((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK)
1584
0
      goto CLEANUP;
1585
0
    if ((res = mp_int_expt(TEMP(1), b - 1, TEMP(3))) != MP_OK)
1586
0
      goto CLEANUP;
1587
0
    if ((res = mp_int_mul_value(TEMP(3), b, TEMP(3))) != MP_OK)
1588
0
      goto CLEANUP;
1589
0
    if ((res = mp_int_div(TEMP(2), TEMP(3), TEMP(4), NULL)) != MP_OK)
1590
0
      goto CLEANUP;
1591
0
    if ((res = mp_int_sub(TEMP(1), TEMP(4), TEMP(4))) != MP_OK)
1592
0
      goto CLEANUP;
1593
0
1594
0
    if (mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) {
1595
0
      if ((res = mp_int_sub_value(TEMP(4), 1, TEMP(4))) != MP_OK)
1596
0
  goto CLEANUP;
1597
0
    }
1598
0
    if ((res = mp_int_copy(TEMP(4), TEMP(1))) != MP_OK)
1599
0
      goto CLEANUP;
1600
0
  }
1601
0
1602
0
  if ((res = mp_int_copy(TEMP(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 (flips)
1607
0
    (void) mp_int_neg(c, c); /* cannot fail */
1608
0
1609
0
  CLEANUP_TEMP();
1610
0
  return res;
1611
0
}
1612
1613
mp_result mp_int_to_int(mp_int z, mp_small *out)
1614
47.1M
{
1615
47.1M
  mp_usmall uv = 0;
1616
47.1M
  mp_size   uz;
1617
47.1M
  mp_digit *dz;
1618
47.1M
  mp_sign   sz;
1619
47.1M
1620
47.1M
  CHECK(z != NULL);
1621
47.1M
1622
47.1M
  /* Make sure the value is representable as a small integer */
1623
47.1M
  sz = MP_SIGN(z);
1624
47.1M
  if ((sz == MP_ZPOS && 
mp_int_compare_value(z, 38.6M
MP_SMALL_MAX38.6M
) > 0) ||
1625
47.1M
      
mp_int_compare_value(z, 41.2M
MP_SMALL_MIN41.2M
) < 0)
1626
8.33M
    return MP_RANGE;
1627
38.8M
1628
38.8M
  uz = MP_USED(z);
1629
38.8M
  dz = MP_DIGITS(z) + uz - 1;
1630
38.8M
1631
85.7M
  while (uz > 0) {
1632
46.9M
    uv <<= MP_DIGIT_BIT/2;
1633
46.9M
    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1634
46.9M
    --uz;
1635
46.9M
  }
1636
38.8M
1637
38.8M
  if (out)
1638
38.8M
    *out = (mp_small)((sz == MP_NEG) ? 
-uv6.11M
:
uv32.6M
);
1639
38.8M
1640
38.8M
  return MP_OK;
1641
38.8M
}
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 || mp_int_compare_uvalue(z, MP_USMALL_MAX) > 0)
1655
0
    return MP_RANGE;
1656
0
     
1657
0
  uz = MP_USED(z);
1658
0
  dz = MP_DIGITS(z) + uz - 1;
1659
0
  
1660
0
  while (uz > 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.45k
{
1675
1.45k
  mp_result res;
1676
1.45k
  int       cmp = 0;
1677
1.45k
1678
1.45k
  CHECK(z != NULL && str != NULL && limit >= 2);
1679
1.45k
1680
1.45k
  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1681
1.45k
    
return MP_RANGE0
;
1682
1.45k
1683
1.45k
  if (CMPZ(z) == 0) {
1684
0
    *str++ = s_val2ch(0, 1);
1685
0
  }
1686
1.45k
  else {
1687
1.45k
    mpz_t tmp;
1688
1.45k
    char  *h, *t;
1689
1.45k
1690
1.45k
    if ((res = mp_int_init_copy(&tmp, z)) != MP_OK)
1691
0
      return res;
1692
1.45k
1693
1.45k
    if (MP_SIGN(z) == MP_NEG) {
1694
747
      *str++ = '-';
1695
747
      --limit;
1696
747
    }
1697
1.45k
    h = str;
1698
1.45k
1699
1.45k
    /* Generate digits in reverse order until finished or limit reached */
1700
27.8k
    for (/* */; limit > 0; 
--limit26.3k
) {
1701
27.8k
      mp_digit d;
1702
27.8k
1703
27.8k
      if ((cmp = CMPZ(&tmp)) == 0)
1704
1.45k
  break;
1705
26.3k
1706
26.3k
      d = s_ddiv(&tmp, (mp_digit)radix);
1707
26.3k
      *str++ = s_val2ch(d, 1);
1708
26.3k
    }
1709
1.45k
    t = str - 1;
1710
1.45k
1711
1.45k
    /* Put digits back in correct output order */
1712
14.0k
    while (h < t) {
1713
12.5k
      char tc = *h;
1714
12.5k
      *h++ = *t;
1715
12.5k
      *t-- = tc;
1716
12.5k
    }
1717
1.45k
1718
1.45k
    mp_int_clear(&tmp);
1719
1.45k
  }
1720
1.45k
1721
1.45k
  *str = '\0';
1722
1.45k
  if (cmp == 0)
1723
1.45k
    return MP_OK;
1724
0
  else
1725
0
    return MP_TRUNC;
1726
1.45k
}
1727
1728
mp_result mp_int_string_len(mp_int z, mp_size radix)
1729
1.56k
{
1730
1.56k
  int  len;
1731
1.56k
1732
1.56k
  CHECK(z != NULL);
1733
1.56k
1734
1.56k
  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1735
1.56k
    
return MP_RANGE0
;
1736
1.56k
1737
1.56k
  len = s_outlen(z, radix) + 1; /* for terminator */
1738
1.56k
1739
1.56k
  /* Allow for sign marker on negatives */
1740
1.56k
  if (MP_SIGN(z) == MP_NEG)
1741
776
    len += 1;
1742
1.56k
1743
1.56k
  return len;
1744
1.56k
}
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
162
{
1749
162
  return mp_int_read_cstring(z, radix, str, NULL);
1750
162
}
1751
1752
mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
1753
162
{
1754
162
  int ch;
1755
162
1756
162
  CHECK(z != NULL && str != NULL);
1757
162
1758
162
  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1759
162
    
return MP_RANGE0
;
1760
162
1761
162
  /* Skip leading whitespace */
1762
162
  while (isspace((int)*str))
1763
0
    ++str;
1764
162
1765
162
  /* Handle leading sign tag (+/-, positive default) */
1766
162
  switch (*str) {
1767
162
  case '-':
1768
61
    MP_SIGN(z) = MP_NEG;
1769
61
    ++str;
1770
61
    break;
1771
162
  case '+':
1772
0
    ++str; /* fallthrough */
1773
101
  default:
1774
101
    MP_SIGN(z) = MP_ZPOS;
1775
101
    break;
1776
162
  }
1777
162
1778
162
  /* Skip leading zeroes */
1779
162
  while ((ch = s_ch2val(*str, radix)) == 0)
1780
0
    ++str;
1781
162
1782
162
  /* Make sure there is enough space for the value */
1783
162
  if (!s_pad(z, s_inlen(strlen(str), radix)))
1784
0
    return MP_MEMORY;
1785
162
1786
162
  MP_USED(z) = 1; z->digits[0] = 0;
1787
162
1788
1.96k
  while (*str != '\0' && 
((ch = s_ch2val(*str, radix)) >= 0)1.80k
) {
1789
1.80k
    s_dmul(z, (mp_digit)radix);
1790
1.80k
    s_dadd(z, (mp_digit)ch);
1791
1.80k
    ++str;
1792
1.80k
  }
1793
162
1794
162
  CLAMP(z);
1795
162
1796
162
  /* Override sign for zero, even if negative specified. */
1797
162
  if (CMPZ(z) == 0)
1798
162
    
MP_SIGN0
(z) = MP_ZPOS0
;
1799
162
1800
162
  if (end != NULL)
1801
162
    
*end = (char *)str0
;
1802
162
1803
162
  /* Return a truncation error if the string has unprocessed characters
1804
162
     remaining, so the caller can tell if the whole string was done */
1805
162
  if (*str != '\0')
1806
0
    return MP_TRUNC;
1807
162
  else
1808
162
    return MP_OK;
1809
162
}
1810
1811
mp_result mp_int_count_bits(mp_int z)
1812
4.00k
{
1813
4.00k
  mp_size  nbits = 0, uz;
1814
4.00k
  mp_digit d;
1815
4.00k
1816
4.00k
  CHECK(z != NULL);
1817
4.00k
1818
4.00k
  uz = MP_USED(z);
1819
4.00k
  if (uz == 1 && 
z->digits[0] == 02.65k
)
1820
0
    return 1;
1821
4.00k
1822
4.00k
  --uz;
1823
4.00k
  nbits = uz * MP_DIGIT_BIT;
1824
4.00k
  d = z->digits[uz];
1825
4.00k
1826
63.1k
  while (d != 0) {
1827
59.1k
    d >>= 1;
1828
59.1k
    ++nbits;
1829
59.1k
  }
1830
4.00k
1831
4.00k
  return nbits;
1832
4.00k
}
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_SIGN(z) == MP_NEG)
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_BIT - 1)) / MP_DIGIT_BIT;
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)) {
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 > 0; --i, ++tmp) {
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_SIGN(z) == MP_NEG)
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_BIT - 1)) / MP_DIGIT_BIT;
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 > 0; --i, ++tmp) {
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
2.44k
{
1938
2.44k
  mp_result  res = mp_int_count_bits(z);
1939
2.44k
  int        bytes;
1940
2.44k
1941
2.44k
  if (res <= 0)
1942
0
    return res;
1943
2.44k
1944
2.44k
  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
1945
2.44k
1946
2.44k
  return bytes;
1947
2.44k
}
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 && s_error_msg[ix] != NULL; ++ix)
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
21.2M
{
1970
21.2M
  mp_digit *out = malloc(num * sizeof(mp_digit));
1971
21.2M
1972
21.2M
  assert(out != NULL); /* for debugging */
1973
#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
1983
21.2M
  return out;
1984
21.2M
}
1985
1986
STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
1987
650k
{
1988
#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
  mp_digit *new = realloc(old, nsize * sizeof(mp_digit));
1998
650k
1999
650k
  assert(new != NULL); /* for debugging */
2000
650k
#endif
2001
650k
  return new;
2002
650k
}
2003
2004
STATIC void s_free(void *ptr)
2005
21.2M
{
2006
21.2M
  free(ptr);
2007
21.2M
}
2008
2009
STATIC int      s_pad(mp_int z, mp_size min)
2010
201M
{
2011
201M
  if (MP_ALLOC(z) < min) {
2012
14.3M
    mp_size nsize = ROUND_PREC(min);
2013
14.3M
    mp_digit *tmp;
2014
14.3M
2015
14.3M
    if ((void *)z->digits == (void *)z) {
2016
13.7M
      if ((tmp = s_alloc(nsize)) == NULL)
2017
13.7M
        
return 00
;
2018
13.7M
2019
13.7M
      COPY(MP_DIGITS(z), tmp, MP_USED(z));
2020
13.7M
    }
2021
650k
    else if ((tmp = s_realloc(MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL)
2022
650k
      
return 00
;
2023
14.3M
2024
14.3M
    MP_DIGITS(z) = tmp;
2025
14.3M
    MP_ALLOC(z) = nsize;
2026
14.3M
  }
2027
201M
2028
201M
  return 1;
2029
201M
}
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
1.05M
{
2034
1.05M
  mp_usmall uv = (mp_usmall) (value < 0) ? 
-value6.58k
:
value1.05M
;
2035
1.05M
  s_ufake(z, uv, vbuf);
2036
1.05M
  if (value < 0)
2037
6.58k
    z->sign = MP_NEG;
2038
1.05M
}
2039
2040
STATIC void      s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[])
2041
51.8M
{
2042
51.8M
  mp_size ndig = (mp_size) s_uvpack(value, vbuf);
2043
51.8M
2044
51.8M
  z->used = ndig;
2045
51.8M
  z->alloc = MP_VALUE_DIGITS(value);
2046
51.8M
  z->sign = MP_ZPOS;
2047
51.8M
  z->digits = vbuf;
2048
51.8M
}
2049
2050
STATIC int      s_cdig(mp_digit *da, mp_digit *db, mp_size len)
2051
58.2M
{
2052
58.2M
  mp_digit *dat = da + len - 1, *dbt = db + len - 1;
2053
58.2M
2054
77.5M
  for (/* */; len != 0; 
--len, --dat, --dbt19.3M
) {
2055
69.4M
    if (*dat > *dbt)
2056
13.8M
      return 1;
2057
55.5M
    else if (*dat < *dbt)
2058
36.1M
      return -1;
2059
69.4M
  }
2060
58.2M
2061
58.2M
  
return 08.17M
;
2062
58.2M
}
2063
2064
STATIC int       s_uvpack(mp_usmall uv, mp_digit t[])
2065
51.8M
{
2066
51.8M
  int ndig = 0;
2067
51.8M
2068
51.8M
  if (uv == 0)
2069
868k
    t[ndig++] = 0;
2070
50.9M
  else {
2071
149M
    while (uv != 0) {
2072
98.1M
      t[ndig++] = (mp_digit) uv;
2073
98.1M
      uv >>= MP_DIGIT_BIT/2;
2074
98.1M
      uv >>= MP_DIGIT_BIT/2;
2075
98.1M
    }
2076
50.9M
  }
2077
51.8M
2078
51.8M
  return ndig;
2079
51.8M
}
2080
2081
STATIC int      s_ucmp(mp_int a, mp_int b)
2082
145M
{
2083
145M
  mp_size  ua = MP_USED(a), ub = MP_USED(b);
2084
145M
2085
145M
  if (ua > ub)
2086
17.1M
    return 1;
2087
128M
  else if (ub > ua)
2088
69.9M
    return -1;
2089
58.2M
  else
2090
58.2M
    return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua);
2091
145M
}
2092
2093
STATIC int      s_vcmp(mp_int a, mp_small v)
2094
50.4M
{
2095
50.4M
  mp_usmall uv = (v < 0) ? 
-(mp_usmall) v8.56M
:
(mp_usmall) v41.9M
;
2096
50.4M
  return s_uvcmp(a, uv);
2097
50.4M
}
2098
2099
STATIC int      s_uvcmp(mp_int a, mp_usmall uv)
2100
50.4M
{
2101
50.4M
  mpz_t vtmp;
2102
50.4M
  mp_digit vdig[MP_VALUE_DIGITS(uv)];
2103
50.4M
2104
50.4M
  s_ufake(&vtmp, uv, vdig);
2105
50.4M
  return s_ucmp(a, &vtmp);
2106
50.4M
}
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
8.39M
{
2111
8.39M
  mp_size pos;
2112
8.39M
  mp_word w = 0;
2113
8.39M
2114
8.39M
  /* Insure that da is the longer of the two to simplify later code */
2115
8.39M
  if (size_b > size_a) {
2116
828k
    SWAP(mp_digit *, da, db);
2117
828k
    SWAP(mp_size, size_a, size_b);
2118
828k
  }
2119
8.39M
2120
8.39M
  /* Add corresponding digits until the shorter number runs out */
2121
18.1M
  for (pos = 0; pos < size_b; 
++pos, ++da, ++db, ++dc9.74M
) {
2122
9.74M
    w = w + (mp_word) *da + (mp_word) *db;
2123
9.74M
    *dc = LOWER_HALF(w);
2124
9.74M
    w = UPPER_HALF(w);
2125
9.74M
  }
2126
8.39M
2127
8.39M
  /* Propagate carries as far as necessary */
2128
13.8M
  for (/* */; pos < size_a; 
++pos, ++da, ++dc5.49M
) {
2129
5.49M
    w = w + *da;
2130
5.49M
2131
5.49M
    *dc = LOWER_HALF(w);
2132
5.49M
    w = UPPER_HALF(w);
2133
5.49M
  }
2134
8.39M
2135
8.39M
  /* Return carry out */
2136
8.39M
  return (mp_digit)w;
2137
8.39M
}
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
86.2M
{
2142
86.2M
  mp_size pos;
2143
86.2M
  mp_word w = 0;
2144
86.2M
2145
86.2M
  /* We assume that |a| >= |b| so this should definitely hold */
2146
86.2M
  assert(size_a >= size_b);
2147
86.2M
2148
86.2M
  /* Subtract corresponding digits and propagate borrow */
2149
227M
  for (pos = 0; pos < size_b; 
++pos, ++da, ++db, ++dc141M
) {
2150
141M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2151
141M
   (mp_word)*da) - w - (mp_word)*db;
2152
141M
2153
141M
    *dc = LOWER_HALF(w);
2154
141M
    w = (UPPER_HALF(w) == 0);
2155
141M
  }
2156
86.2M
2157
86.2M
  /* Finish the subtraction for remaining upper digits of da */
2158
171M
  for (/* */; pos < size_a; 
++pos, ++da, ++dc85.1M
) {
2159
85.1M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2160
85.1M
   (mp_word)*da) - w;
2161
85.1M
2162
85.1M
    *dc = LOWER_HALF(w);
2163
85.1M
    w = (UPPER_HALF(w) == 0);
2164
85.1M
  }
2165
86.2M
2166
86.2M
  /* If there is a borrow out at the end, it violates the precondition */
2167
86.2M
  assert(w == 0);
2168
86.2M
}
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.57M
{
2173
7.57M
  mp_size  bot_size;
2174
7.57M
2175
7.57M
  /* Make sure b is the smaller of the two input values */
2176
7.57M
  if (size_b > size_a) {
2177
3.03M
    SWAP(mp_digit *, da, db);
2178
3.03M
    SWAP(mp_size, size_a, size_b);
2179
3.03M
  }
2180
7.57M
2181
7.57M
  /* Insure that the bottom is the larger half in an odd-length split; the code
2182
7.57M
     below relies on this being true.
2183
7.57M
   */
2184
7.57M
  bot_size = (size_a + 1) / 2;
2185
7.57M
2186
7.57M
  /* If the values are big enough to bother with recursion, use the Karatsuba
2187
7.57M
     algorithm to compute the product; otherwise use the normal multiplication
2188
7.57M
     algorithm
2189
7.57M
   */
2190
7.57M
  if (multiply_threshold &&
2191
7.57M
      size_a >= multiply_threshold &&
2192
7.57M
      
size_b > bot_size2.41k
) {
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)) == NULL) 
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
36
  }
2248
7.57M
  else {
2249
7.57M
    s_umul(da, db, dc, size_a, size_b);
2250
7.57M
  }
2251
7.57M
2252
7.57M
  return 1;
2253
7.57M
}
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.57M
{
2258
7.57M
  mp_size a, b;
2259
7.57M
  mp_word w;
2260
7.57M
2261
25.9M
  for (a = 0; a < size_a; 
++a, ++dc, ++da18.3M
) {
2262
18.3M
    mp_digit *dct = dc;
2263
18.3M
    mp_digit *dbt = db;
2264
18.3M
2265
18.3M
    if (*da == 0)
2266
1.29M
      continue;
2267
17.0M
2268
17.0M
    w = 0;
2269
40.5M
    for (b = 0; b < size_b; 
++b, ++dbt, ++dct23.4M
) {
2270
23.4M
      w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
2271
23.4M
2272
23.4M
      *dct = LOWER_HALF(w);
2273
23.4M
      w = UPPER_HALF(w);
2274
23.4M
    }
2275
17.0M
2276
17.0M
    *dct = (mp_digit)w;
2277
17.0M
  }
2278
7.57M
}
2279
2280
STATIC int       s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2281
0
{
2282
0
  if (multiply_threshold && size_a > multiply_threshold) {
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)) == NULL) return 0;
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 < top; ++i) {
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_a; ++i, dc += 2, ++da) {
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_a; ++j, ++dat, ++dct) {
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 (ov) {
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 = UPPER_HALF(w)) != 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.80k
{
2387
1.80k
  mp_word w = 0;
2388
1.80k
  mp_digit *da = MP_DIGITS(a);
2389
1.80k
  mp_size ua = MP_USED(a);
2390
1.80k
2391
1.80k
  w = (mp_word)*da + b;
2392
1.80k
  *da++ = LOWER_HALF(w);
2393
1.80k
  w = UPPER_HALF(w);
2394
1.80k
2395
2.00k
  for (ua -= 1; ua > 0; 
--ua, ++da200
) {
2396
200
    w = (mp_word)*da + w;
2397
200
2398
200
    *da = LOWER_HALF(w);
2399
200
    w = UPPER_HALF(w);
2400
200
  }
2401
1.80k
2402
1.80k
  if (w) {
2403
2
    *da = (mp_digit)w;
2404
2
    MP_USED(a) += 1;
2405
2
  }
2406
1.80k
}
2407
2408
STATIC void      s_dmul(mp_int a, mp_digit b)
2409
1.80k
{
2410
1.80k
  mp_word w = 0;
2411
1.80k
  mp_digit *da = MP_DIGITS(a);
2412
1.80k
  mp_size ua = MP_USED(a);
2413
1.80k
2414
3.78k
  while (ua > 0) {
2415
1.98k
    w = (mp_word)*da * b + w;
2416
1.98k
    *da++ = LOWER_HALF(w);
2417
1.98k
    w = UPPER_HALF(w);
2418
1.98k
    --ua;
2419
1.98k
  }
2420
1.80k
2421
1.80k
  if (w) {
2422
20
    *da = (mp_digit)w;
2423
20
    MP_USED(a) += 1;
2424
20
  }
2425
1.80k
}
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 > 0) {
2432
8.55M
    w = (mp_word)*da++ * (mp_word)b + w;
2433
8.55M
2434
8.55M
    *dc++ = LOWER_HALF(w);
2435
8.55M
    w = UPPER_HALF(w);
2436
8.55M
    --size_a;
2437
8.55M
  }
2438
2.08M
2439
2.08M
  if (w)
2440
0
    *dc = LOWER_HALF(w);
2441
2.08M
}
2442
2443
STATIC mp_digit  s_ddiv(mp_int a, mp_digit b)
2444
913k
{
2445
913k
  mp_word w = 0, qdigit;
2446
913k
  mp_size ua = MP_USED(a);
2447
913k
  mp_digit *da = MP_DIGITS(a) + ua - 1;
2448
913k
2449
3.09M
  for (/* */; ua > 0; 
--ua, --da2.17M
) {
2450
2.17M
    w = (w << MP_DIGIT_BIT) | *da;
2451
2.17M
2452
2.17M
    if (w >= b) {
2453
1.63M
      qdigit = w / b;
2454
1.63M
      w = w % b;
2455
1.63M
    }
2456
541k
    else {
2457
541k
      qdigit = 0;
2458
541k
    }
2459
2.17M
2460
2.17M
    *da = (mp_digit)qdigit;
2461
2.17M
  }
2462
913k
2463
913k
  CLAMP(a);
2464
913k
  return (mp_digit)w;
2465
913k
}
2466
2467
STATIC void     s_qdiv(mp_int z, mp_size p2)
2468
86.2M
{
2469
86.2M
  mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
2470
86.2M
  mp_size uz = MP_USED(z);
2471
86.2M
2472
86.2M
  if (ndig) {
2473
253k
    mp_size  mark;
2474
253k
    mp_digit *to, *from;
2475
253k
2476
253k
    if (ndig >= uz) {
2477
0
      mp_int_zero(z);
2478
0
      return;
2479
0
    }
2480
253k
2481
253k
    to = MP_DIGITS(z); from = to + ndig;
2482
253k
2483
586k
    for (mark = ndig; mark < uz; 
++mark332k
)
2484
332k
      *to++ = *from++;
2485
253k
2486
253k
    MP_USED(z) = uz - ndig;
2487
253k
  }
2488
86.2M
2489
86.2M
  if (nbits) {
2490
81.5M
    mp_digit d = 0, *dz, save;
2491
81.5M
    mp_size  up = MP_DIGIT_BIT - nbits;
2492
81.5M
2493
81.5M
    uz = MP_USED(z);
2494
81.5M
    dz = MP_DIGITS(z) + uz - 1;
2495
81.5M
2496
296M
    for (/* */; uz > 0; 
--uz, --dz214M
) {
2497
214M
      save = *dz;
2498
214M
2499
214M
      *dz = (*dz >> nbits) | (d << up);
2500
214M
      d = save;
2501
214M
    }
2502
81.5M
2503
81.5M
    CLAMP(z);
2504
81.5M
  }
2505
86.2M
2506
86.2M
  if (MP_USED(z) == 1 && 
z->digits[0] == 033.5M
)
2507
86.2M
    
MP_SIGN782k
(z) = MP_ZPOS782k
;
2508
86.2M
}
2509
2510
STATIC void     s_qmod(mp_int z, mp_size p2)
2511
140k
{
2512
140k
  mp_size start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
2513
140k
  mp_size uz = MP_USED(z);
2514
140k
  mp_digit mask = (1u << rest) - 1;
2515
140k
2516
140k
  if (start <= uz) {
2517
140k
    MP_USED(z) = start;
2518
140k
    z->digits[start - 1] &= mask;
2519
140k
    CLAMP(z);
2520
140k
  }
2521
140k
}
2522
2523
STATIC int      s_qmul(mp_int z, mp_size p2)
2524
4.45M
{
2525
4.45M
  mp_size   uz, need, rest, extra, i;
2526
4.45M
  mp_digit *from, *to, d;
2527
4.45M
2528
4.45M
  if (p2 == 0)
2529
1.47M
    return 1;
2530
2.97M
2531
2.97M
  uz = MP_USED(z); 
2532
2.97M
  need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT;
2533
2.97M
2534
2.97M
  /* Figure out if we need an extra digit at the top end; this occurs if the
2535
2.97M
     topmost `rest' bits of the high-order digit of z are not zero, meaning
2536
2.97M
     they will be shifted off the end if not preserved */
2537
2.97M
  extra = 0;
2538
2.97M
  if (rest != 0) {
2539
2.91M
    mp_digit *dz = MP_DIGITS(z) + uz - 1;
2540
2.91M
2541
2.91M
    if ((*dz >> (MP_DIGIT_BIT - rest)) != 0)
2542
586k
      extra = 1;
2543
2.91M
  }
2544
2.97M
2545
2.97M
  if (!s_pad(z, uz + need + extra))
2546
0
    return 0;
2547
2.97M
2548
2.97M
  /* If we need to shift by whole digits, do that in one pass, then
2549
2.97M
     to back and shift by partial digits.
2550
2.97M
   */
2551
2.97M
  if (need > 0) {
2552
83.9k
    from = MP_DIGITS(z) + uz - 1;
2553
83.9k
    to = from + need;
2554
83.9k
2555
168k
    for (i = 0; i < uz; 
++i84.0k
)
2556
84.0k
      *to-- = *from--;
2557
83.9k
2558
83.9k
    ZERO(MP_DIGITS(z), need);
2559
83.9k
    uz += need;
2560
83.9k
  }
2561
2.97M
2562
2.97M
  if (rest) {
2563
2.91M
    d = 0;
2564
10.4M
    for (i = need, from = 
MP_DIGITS2.91M
(z) + need; i < uz;
++i, ++from7.58M
) {
2565
7.58M
      mp_digit save = *from;
2566
7.58M
      
2567
7.58M
      *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
2568
7.58M
      d = save;
2569
7.58M
    }
2570
2.91M
2571
2.91M
    d >>= (MP_DIGIT_BIT - rest);
2572
2.91M
    if (d != 0) {
2573
586k
      *from = d;
2574
586k
      uz += extra;
2575
586k
    }
2576
2.91M
  }
2577
2.97M
2578
2.97M
  MP_USED(z) = uz;
2579
2.97M
  CLAMP(z);
2580
2.97M
2581
2.97M
  return 1;
2582
2.97M
}
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 = MP_DIGITS(z); pos < tdig; ++pos, ++zp) {
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_HALF(w) ? 0 : 1;
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
84.5M
{
2616
84.5M
  int       k = 0;
2617
84.5M
  mp_digit *dp = MP_DIGITS(z), d;
2618
84.5M
2619
84.5M
  if (MP_USED(z) == 1 && 
*dp == 030.4M
)
2620
0
    return 1;
2621
84.5M
2622
84.9M
  
while (84.5M
*dp == 0) {
2623
430k
    k += MP_DIGIT_BIT;
2624
430k
    ++dp;
2625
430k
  }
2626
84.5M
2627
84.5M
  d = *dp;
2628
245M
  while ((d & 1) == 0) {
2629
161M
    d >>= 1;
2630
161M
    ++k;
2631
161M
  }
2632
84.5M
2633
84.5M
  return k;
2634
84.5M
}
2635
2636
STATIC int       s_isp2(mp_int z)
2637
2.70M
{
2638
2.70M
  mp_size uz = MP_USED(z), k = 0;
2639
2.70M
  mp_digit *dz = MP_DIGITS(z), d;
2640
2.70M
2641
2.75M
  while (uz > 1) {
2642
900k
    if (*dz++ != 0)
2643
844k
      return -1;
2644
55.9k
    k += MP_DIGIT_BIT;
2645
55.9k
    --uz;
2646
55.9k
  }
2647
2.70M
2648
2.70M
  d = *dz;
2649
5.19M
  while (d > 1) {
2650
4.22M
    if (d & 1)
2651
888k
      return -1;
2652
3.33M
    ++k; d >>= 1;
2653
3.33M
  }
2654
1.85M
2655
1.85M
  
return (int) k970k
;
2656
1.85M
}
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_BIT) / MP_DIGIT_BIT;
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 = MP_DIGITS(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
846k
{
2679
846k
  mp_digit d = b->digits[MP_USED(b) - 1];
2680
846k
  int k = 0;
2681
846k
2682
15.1M
  while (d < (1u << (mp_digit)(MP_DIGIT_BIT - 1))) { /* d < (MP_RADIX / 2) */
2683
14.3M
    d <<= 1;
2684
14.3M
    ++k;
2685
14.3M
  }
2686
846k
2687
846k
  /* These multiplications can't fail */
2688
846k
  if (k != 0) {
2689
830k
    (void) s_qmul(a, (mp_size) k);
2690
830k
    (void) s_qmul(b, (mp_size) k);
2691
830k
  }
2692
846k
2693
846k
  return k;
2694
846k
}
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, MP_DIGIT_BIT * 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 ((CMPZ(x) < 0) && !s_qsub(x, umb_p1))
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 (mp_int_compare(x, m) >= 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_USED(mu); db = MP_DIGITS(b); dbt = db + MP_USED(b) - 1;
2759
0
2760
0
  while (last__ < 3) {
2761
0
    SETUP(mp_int_init_size(LAST_TEMP(), 4 * umu));
2762
0
    ZERO(MP_DIGITS(TEMP(last__ - 1)), MP_ALLOC(TEMP(last__ - 1)));
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 < dbt) {
2769
0
    int      i;
2770
0
2771
0
    for (d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) {
2772
0
      if (d & 1) {
2773
0
  /* The use of a second temporary avoids allocation */
2774
0
  UMUL(c, a, TEMP(0));
2775
0
  if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
2776
0
    res = MP_MEMORY; goto CLEANUP;
2777
0
  }
2778
0
  mp_int_copy(TEMP(0), c);
2779
0
      }
2780
0
2781
0
2782
0
      USQR(a, TEMP(0));
2783
0
      assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
2784
0
      if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
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 & 1) {
2798
0
      UMUL(c, a, TEMP(0));
2799
0
      if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
2800
0
  res = MP_MEMORY; goto CLEANUP;
2801
0
      }
2802
0
      mp_int_copy(TEMP(0), c);
2803
0
    }
2804
0
2805
0
    d >>= 1;
2806
0
    if (!d) break;
2807
0
2808
0
    USQR(a, TEMP(0));
2809
0
    if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
2810
0
      res = MP_MEMORY; goto CLEANUP;
2811
0
    }
2812
0
    (void) mp_int_copy(TEMP(0), a);
2813
0
  }
2814
0
2815
0
  CLEANUP_TEMP();
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.73M
STATIC mp_result s_udiv_knuth(mp_int u, mp_int v) {
2840
1.73M
  mpz_t q, r, t;
2841
1.73M
  mp_result
2842
1.73M
  res = MP_OK;
2843
1.73M
  int k,j;
2844
1.73M
  mp_size m,n;
2845
1.73M
2846
1.73M
  /* Force signs to positive */
2847
1.73M
  MP_SIGN(u) = MP_ZPOS;
2848
1.73M
  MP_SIGN(v) = MP_ZPOS;
2849
1.73M
2850
1.73M
  /* Use simple division algorithm when v is only one digit long */
2851
1.73M
  if (MP_USED(v) == 1) {
2852
886k
    mp_digit d, rem;
2853
886k
    d   = v->digits[0];
2854
886k
    rem = s_ddiv(u, d);
2855
886k
    mp_int_set_value(v, rem);
2856
886k
    return MP_OK;
2857
886k
  }
2858
846k
2859
846k
  /* Algorithm D
2860
846k
2861
846k
     The n and m variables are defined as used by Knuth.
2862
846k
     u is an n digit number with digits u_{n-1}..u_0.
2863
846k
     v is an n+m digit number with digits from v_{m+n-1}..v_0.
2864
846k
     We require that n > 1 and m >= 0
2865
846k
   */
2866
846k
  n = MP_USED(v);
2867
846k
  m = MP_USED(u) - n;
2868
846k
  assert(n > 1);
2869
846k
  assert(m >= 0);
2870
846k
2871
846k
  /* D1: Normalize.
2872
846k
     The normalization step provides the necessary condition for Theorem B,
2873
846k
     which states that the quotient estimate for q_j, call it qhat
2874
846k
2875
846k
       qhat = u_{j+n}u_{j+n-1} / v_{n-1}
2876
846k
2877
846k
     is bounded by
2878
846k
2879
846k
      qhat - 2 <= q_j <= qhat.
2880
846k
2881
846k
     That is, qhat is always greater than the actual quotient digit q,
2882
846k
     and it is never more than two larger than the actual quotient digit.
2883
846k
   */
2884
846k
  k = s_norm(u, v);
2885
846k
2886
846k
  /* Extend size of u by one if needed.
2887
846k
2888
846k
     The algorithm begins with a value of u that has one more digit of input.
2889
846k
     The normalization step sets u_{m+n}..u_0 = 2^k * u_{m+n-1}..u_0. If the
2890
846k
     multiplication did not increase the number of digits of u, we need to add
2891
846k
     a leading zero here.
2892
846k
   */
2893
846k
  if (k == 0 || 
MP_USED830k
(u) != m + n + 1830k
) {
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
318k
  }
2899
846k
2900
846k
  /* Add a leading 0 to v.
2901
846k
2902
846k
     The multiplication in step D4 multiplies qhat * 0v_{n-1}..v_0.  We need to
2903
846k
     add the leading zero to v here to ensure that the multiplication will
2904
846k
     produce the full n+1 digit result.
2905
846k
   */
2906
846k
  if (!s_pad(v, n+1)) 
return MP_MEMORY0
; v->digits[n] = 0;
2907
846k
2908
846k
  /* Initialize temporary variables q and t.
2909
846k
     q allocates space for m+1 digits to store the quotient digits
2910
846k
     t allocates space for n+1 digits to hold the result of q_j*v
2911
846k
   */
2912
846k
  if ((res = mp_int_init_size(&q, m + 1)) != MP_OK) 
return res0
;
2913
846k
  if ((res = mp_int_init_size(&t, n + 1)) != MP_OK) 
goto CLEANUP0
;
2914
846k
2915
846k
  /* D2: Initialize j */
2916
846k
  j = m;
2917
846k
  r.digits = MP_DIGITS(u) + j;  /* The contents of r are shared with u */
2918
846k
  r.used   = n + 1;
2919
846k
  r.sign   = MP_ZPOS;
2920
846k
  r.alloc  = MP_ALLOC(u);
2921
846k
  ZERO(t.digits, t.alloc);
2922
846k
2923
846k
  /* Calculate the m+1 digits of the quotient result */
2924
2.77M
  for (; j >= 0; 
j--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 > MP_DIGIT_MAX)
2939
1.92M
      
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) > 0) {   /* would the remainder be negative? */
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) > 0) { /* would the remainder be negative? */
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
4.39k
      }
2969
157k
      assert(s_ucmp(&t, &r) <=  0 && "The mathematics failed us.");
2970
157k
    }
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
846k
2999
846k
  /* Get rid of leading zeros in q */
3000
846k
  q.used = m + 1;
3001
846k
  CLAMP(&q);
3002
846k
3003
846k
  /* Denormalize the remainder */
3004
846k
  CLAMP(u); /* use u here because the r.digits pointer is off-by-one */
3005
846k
  if (k != 0)
3006
830k
    s_qdiv(u, k);
3007
846k
3008
846k
  mp_int_copy(u, v);  /* ok:  0 <= r < v */
3009
846k
  mp_int_copy(&q, u); /* ok:  q <= u     */
3010
846k
3011
846k
  mp_int_clear(&t);
3012
846k
 CLEANUP:
3013
846k
  mp_int_clear(&q);
3014
846k
  return res;
3015
846k
}
3016
3017
STATIC int       s_outlen(mp_int z, mp_size r)
3018
1.56k
{
3019
1.56k
  mp_result bits;
3020
1.56k
  double raw;
3021
1.56k
3022
1.56k
  assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
3023
1.56k
3024
1.56k
  bits = mp_int_count_bits(z);
3025
1.56k
  raw = (double)bits * s_log2[r];
3026
1.56k
3027
1.56k
  return (int)(raw + 0.999999);
3028
1.56k
}
3029
3030
STATIC mp_size   s_inlen(int len, mp_size r)
3031
162
{
3032
162
  double  raw = (double)len / s_log2[r];
3033
162
  mp_size bits = (mp_size)(raw + 0.5);
3034
162
3035
162
  return (mp_size)((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT) + 1;
3036
162
}
3037
3038
STATIC int       s_ch2val(char c, int r)
3039
1.96k
{
3040
1.96k
  int out;
3041
1.96k
3042
1.96k
  if (isdigit((unsigned char) c))
3043
1.96k
    out = c - '0';
3044
0
  else if (r > 10 && isalpha((unsigned char) c))
3045
0
    out = toupper(c) - 'A' + 10;
3046
0
  else
3047
0
    return -1;
3048
1.96k
3049
1.96k
  return (out >= r) ? 
-10
: out;
3050
1.96k
}
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
0
  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
0
  }
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 >= 0; --i) {
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_USED(z); dz = MP_DIGITS(z);
3093
0
  while (uz > 0 && pos < limit) {
3094
0
    mp_digit d = *dz++;
3095
0
    int i;
3096
0
3097
0
    for (i = sizeof(mp_digit); i > 0 && pos < limit; --i) {
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 && uz == 1)
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 > 0) break;
3108
0
3109
0
    --uz;
3110
0
  }
3111
0
3112
0
  if (pad != 0 && (buf[pos - 1] >> (CHAR_BIT - 1))) {
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_OK : MP_TRUNC;
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 */