Coverage Report

Created: 2019-04-21 11:35

/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
335M
#define CHECK(TEST)   assert(TEST)
73
67.7M
#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
29.4M
#define MP_VALUE_DIGITS(V) \
103
29.4M
((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
104
105
/* Round precision P to nearest word boundary */
106
17.3M
#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
107
108
/* Set array P of S digits to zero */
109
8.13M
#define ZERO(P, S) \
110
8.13M
do{ \
111
8.13M
  mp_size i__ = (S) * sizeof(mp_digit); \
112
8.13M
  mp_digit *p__ = (P); \
113
8.13M
  memset(p__, 0, i__); \
114
8.13M
} while(0)
115
116
/* Copy S digits from array P to array Q */
117
111M
#define COPY(P, Q, S) \
118
111M
do{ \
119
100M
  mp_size i__ = (S) * sizeof(mp_digit); \
120
100M
  mp_digit *p__ = (P), *q__ = (Q); \
121
100M
  memcpy(q__, p__, i__); \
122
100M
} 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
161M
#define CLAMP(Z) \
136
161M
do{ \
137
161M
  mp_int z_ = (Z); \
138
161M
  mp_size uz_ = MP_USED(z_); \
139
161M
  mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \
140
188M
  while (uz_ > 1 && 
(*dz_-- == 0)137M
) \
141
161M
    
--uz_27.2M
; \
142
161M
  MP_USED(z_) = uz_; \
143
161M
} while(0)
144
145
/* Select min/max.  Do not provide expressions for which multiple
146
   evaluation would be problematic, e.g. x++ */
147
2.31M
#define MIN(A, B) ((B)<(A)?
(B)486k
:
(A)1.82M
)
148
88.7M
#define MAX(A, B) ((B)>(A)?
(B)44.2M
:
(A)44.4M
)
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
5.85M
#define SWAP(T, A, B) \
153
5.85M
do{ \
154
5.85M
  T t_ = (A); \
155
5.85M
  A = (B); \
156
5.85M
  B = t_; \
157
5.85M
} while(0)
158
159
/* Used to set up and access simple temp stacks within functions. */
160
#define DECLARE_TEMP(N) \
161
3.48M
  mpz_t temp[(N)]; \
162
3.48M
  int last__ = 0
163
#define CLEANUP_TEMP() \
164
2.24M
 CLEANUP: \
165
3.85M
  while (--last__ >= 0) \
166
2.24M
    
mp_int_clear(1.61M
TEMP1.61M
(last__))
167
3.23M
#define TEMP(K) (temp + (K))
168
1.61M
#define LAST_TEMP() TEMP(last__)
169
1.61M
#define SETUP(E) \
170
1.61M
do{ \
171
1.61M
  if ((res = (E)) != MP_OK) \
172
1.61M
    
goto CLEANUP0
; \
173
1.61M
  ++(last__); \
174
1.61M
} while(0)
175
176
/* Compare value to zero. */
177
217M
#define CMPZ(Z) \
178
217M
(((Z)->used==1&&
(Z)->digits[0]==070.1M
)?
03.93M
:
((Z)->sign==MP_NEG)?213M
-1182M
:
130.6M
)
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
250M
#define UPPER_HALF(W)           ((mp_word)((W) >> MP_DIGIT_BIT))
204
250M
#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
23.4M
{
361
23.4M
  if (z == NULL)
362
23.4M
    
return MP_BADARG0
;
363
23.4M
364
23.4M
  z->single = 0;
365
23.4M
  z->digits = &(z->single);
366
23.4M
  z->alloc  = 1;
367
23.4M
  z->used   = 1;
368
23.4M
  z->sign   = MP_ZPOS;
369
23.4M
370
23.4M
  return MP_OK;
371
23.4M
}
372
373
mp_int    mp_int_alloc(void)
374
18.0M
{
375
18.0M
  mp_int out = malloc(sizeof(mpz_t));
376
18.0M
377
18.0M
  if (out != NULL)
378
18.0M
    mp_int_init(out);
379
18.0M
380
18.0M
  return out;
381
18.0M
}
382
383
mp_result mp_int_init_size(mp_int z, mp_size prec)
384
5.71M
{
385
5.71M
  CHECK(z != NULL);
386
5.71M
387
5.71M
  if (prec == 0)
388
0
    prec = default_precision;
389
5.71M
  else if (prec == 1)
390
334k
    return mp_int_init(z);
391
5.37M
  else
392
5.37M
    prec = (mp_size) ROUND_PREC(prec);
393
5.71M
394
5.71M
  
if (5.37M
(5.37M
MP_DIGITS5.37M
(z) = s_alloc(prec)) == NULL)
395
5.37M
    
return MP_MEMORY0
;
396
5.37M
397
5.37M
  z->digits[0] = 0;
398
5.37M
  MP_USED(z) = 1;
399
5.37M
  MP_ALLOC(z) = prec;
400
5.37M
  MP_SIGN(z) = MP_ZPOS;
401
5.37M
402
5.37M
  return MP_OK;
403
5.37M
}
404
405
mp_result mp_int_init_copy(mp_int z, mp_int old)
406
6.30M
{
407
6.30M
  mp_result res;
408
6.30M
  mp_size uold;
409
6.30M
410
6.30M
  CHECK(z != NULL && old != NULL);
411
6.30M
412
6.30M
  uold = MP_USED(old);
413
6.30M
  if (uold == 1) {
414
2.27M
    mp_int_init(z);
415
2.27M
  }
416
4.02M
  else {
417
4.02M
    mp_size target = MAX(uold, default_precision);
418
4.02M
419
4.02M
    if ((res = mp_int_init_size(z, target)) != MP_OK)
420
0
      return res;
421
6.30M
  }
422
6.30M
423
6.30M
  MP_USED(z) = uold;
424
6.30M
  MP_SIGN(z) = MP_SIGN(old);
425
6.30M
  COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
426
6.30M
427
6.30M
  return MP_OK;
428
6.30M
}
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
49.2k
{
441
49.2k
  mpz_t    vtmp;
442
49.2k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
443
49.2k
444
49.2k
  s_ufake(&vtmp, uvalue, vbuf);
445
49.2k
  return mp_int_init_copy(z, &vtmp);
446
49.2k
}
447
448
mp_result  mp_int_set_value(mp_int z, mp_small value)
449
792k
{
450
792k
  mpz_t    vtmp;
451
792k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
452
792k
453
792k
  s_fake(&vtmp, value, vbuf);
454
792k
  return mp_int_copy(&vtmp, z);
455
792k
}
456
457
mp_result  mp_int_set_uvalue(mp_int z, mp_usmall uvalue)
458
9.76k
{
459
9.76k
  mpz_t    vtmp;
460
9.76k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
461
9.76k
462
9.76k
  s_ufake(&vtmp, uvalue, vbuf);
463
9.76k
  return mp_int_copy(&vtmp, z);
464
9.76k
}
465
466
void      mp_int_clear(mp_int z)
467
28.8M
{
468
28.8M
  if (z == NULL)
469
28.8M
    
return0
;
470
28.8M
471
28.8M
  if (MP_DIGITS(z) != NULL) {
472
28.8M
    if (MP_DIGITS(z) != &(z->single))
473
15.7M
      s_free(MP_DIGITS(z));
474
28.8M
475
28.8M
    MP_DIGITS(z) = NULL;
476
28.8M
  }
477
28.8M
}
478
479
void      mp_int_free(mp_int z)
480
18.0M
{
481
18.0M
  NRCHECK(z != NULL);
482
18.0M
483
18.0M
  mp_int_clear(z);
484
18.0M
  free(z); /* note: NOT s_free() */
485
18.0M
}
486
487
mp_result mp_int_copy(mp_int a, mp_int c)
488
88.9M
{
489
88.9M
  CHECK(a != NULL && c != NULL);
490
88.9M
491
88.9M
  if (a != c) {
492
84.3M
    mp_size ua = MP_USED(a);
493
84.3M
    mp_digit *da, *dc;
494
84.3M
495
84.3M
    if (!s_pad(c, ua))
496
0
      return MP_MEMORY;
497
84.3M
498
84.3M
    da = MP_DIGITS(a); dc = MP_DIGITS(c);
499
84.3M
    COPY(da, dc, ua);
500
84.3M
501
84.3M
    MP_USED(c) = ua;
502
84.3M
    MP_SIGN(c) = MP_SIGN(a);
503
84.3M
  }
504
88.9M
505
88.9M
  return MP_OK;
506
88.9M
}
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
9.75M
{
525
9.75M
  NRCHECK(z != NULL);
526
9.75M
527
9.75M
  z->digits[0] = 0;
528
9.75M
  MP_USED(z) = 1;
529
9.75M
  MP_SIGN(z) = MP_ZPOS;
530
9.75M
}
531
532
mp_result mp_int_abs(mp_int a, mp_int c)
533
3.00M
{
534
3.00M
  mp_result res;
535
3.00M
536
3.00M
  CHECK(a != NULL && c != NULL);
537
3.00M
538
3.00M
  if ((res = mp_int_copy(a, c)) != MP_OK)
539
0
    return res;
540
3.00M
541
3.00M
  MP_SIGN(c) = MP_ZPOS;
542
3.00M
  return MP_OK;
543
3.00M
}
544
545
mp_result mp_int_neg(mp_int a, mp_int c)
546
63.2M
{
547
63.2M
  mp_result res;
548
63.2M
549
63.2M
  CHECK(a != NULL && c != NULL);
550
63.2M
551
63.2M
  if ((res = mp_int_copy(a, c)) != MP_OK)
552
0
    return res;
553
63.2M
554
63.2M
  if (CMPZ(c) != 0)
555
63.2M
    
MP_SIGN63.2M
(c) = 1 - 63.2M
MP_SIGN63.2M
(a);
556
63.2M
557
63.2M
  return MP_OK;
558
63.2M
}
559
560
mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
561
6.81M
{
562
6.81M
  mp_size ua, ub, uc, max;
563
6.81M
564
6.81M
  CHECK(a != NULL && b != NULL && c != NULL);
565
6.81M
566
6.81M
  ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
567
6.81M
  max = MAX(ua, ub);
568
6.81M
569
6.81M
  if (MP_SIGN(a) == MP_SIGN(b)) {
570
3.65M
    /* Same sign -- add magnitudes, preserve sign of addends */
571
3.65M
    mp_digit carry;
572
3.65M
573
3.65M
    if (!s_pad(c, max))
574
0
      return MP_MEMORY;
575
3.65M
576
3.65M
    carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
577
3.65M
    uc = max;
578
3.65M
579
3.65M
    if (carry) {
580
16.0k
      if (!s_pad(c, max + 1))
581
0
  return MP_MEMORY;
582
16.0k
583
16.0k
      c->digits[max] = carry;
584
16.0k
      ++uc;
585
16.0k
    }
586
3.65M
587
3.65M
    MP_USED(c) = uc;
588
3.65M
    MP_SIGN(c) = MP_SIGN(a);
589
3.65M
590
3.65M
  }
591
3.15M
  else {
592
3.15M
    /* Different signs -- subtract magnitudes, preserve sign of greater */
593
3.15M
    mp_int  x, y;
594
3.15M
    int     cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
595
3.15M
596
3.15M
    /* Set x to max(a, b), y to min(a, b) to simplify later code.
597
3.15M
       A special case yields zero for equal magnitudes.
598
3.15M
    */
599
3.15M
    if (cmp == 0) {
600
578k
      mp_int_zero(c);
601
578k
      return MP_OK;
602
578k
    }
603
2.57M
    else if (cmp < 0) {
604
966k
      x = b; y = a;
605
966k
    }
606
1.61M
    else {
607
1.61M
      x = a; y = b;
608
1.61M
    }
609
3.15M
610
3.15M
    
if (2.57M
!s_pad(c, 2.57M
MP_USED2.57M
(x)))
611
0
      return MP_MEMORY;
612
2.57M
613
2.57M
    /* Subtract smaller from larger */
614
2.57M
    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
615
2.57M
    MP_USED(c) = MP_USED(x);
616
2.57M
    CLAMP(c);
617
2.57M
618
2.57M
    /* Give result the sign of the larger */
619
2.57M
    MP_SIGN(c) = MP_SIGN(x);
620
2.57M
  }
621
6.81M
622
6.81M
  
return MP_OK6.23M
;
623
6.81M
}
624
625
mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c)
626
436
{
627
436
  mpz_t    vtmp;
628
436
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
629
436
630
436
  s_fake(&vtmp, value, vbuf);
631
436
632
436
  return mp_int_add(a, &vtmp, c);
633
436
}
634
635
mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
636
71.3M
{
637
71.3M
  mp_size ua, ub, uc, max;
638
71.3M
639
71.3M
  CHECK(a != NULL && b != NULL && c != NULL);
640
71.3M
641
71.3M
  ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
642
71.3M
  max = MAX(ua, ub);
643
71.3M
644
71.3M
  if (MP_SIGN(a) != MP_SIGN(b)) {
645
208k
    /* Different signs -- add magnitudes and keep sign of a */
646
208k
    mp_digit carry;
647
208k
648
208k
    if (!s_pad(c, max))
649
0
      return MP_MEMORY;
650
208k
651
208k
    carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
652
208k
    uc = max;
653
208k
654
208k
    if (carry) {
655
1.80k
      if (!s_pad(c, max + 1))
656
0
  return MP_MEMORY;
657
1.80k
658
1.80k
      c->digits[max] = carry;
659
1.80k
      ++uc;
660
1.80k
    }
661
208k
662
208k
    MP_USED(c) = uc;
663
208k
    MP_SIGN(c) = MP_SIGN(a);
664
208k
665
208k
  }
666
71.1M
  else {
667
71.1M
    /* Same signs -- subtract magnitudes */
668
71.1M
    mp_int  x, y;
669
71.1M
    mp_sign osign;
670
71.1M
    int     cmp = s_ucmp(a, b);
671
71.1M
672
71.1M
    if (!s_pad(c, max))
673
0
      return MP_MEMORY;
674
71.1M
675
71.1M
    if (cmp >= 0) {
676
12.0M
      x = a; y = b; osign = MP_ZPOS;
677
12.0M
    }
678
59.0M
    else {
679
59.0M
      x = b; y = a; osign = MP_NEG;
680
59.0M
    }
681
71.1M
682
71.1M
    if (MP_SIGN(a) == MP_NEG && 
cmp != 064.9k
)
683
52.3k
      osign = 1 - osign;
684
71.1M
685
71.1M
    s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
686
71.1M
    MP_USED(c) = MP_USED(x);
687
71.1M
    CLAMP(c);
688
71.1M
689
71.1M
    MP_SIGN(c) = osign;
690
71.1M
  }
691
71.3M
692
71.3M
  return MP_OK;
693
71.3M
}
694
695
mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c)
696
12.2k
{
697
12.2k
  mpz_t    vtmp;
698
12.2k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
699
12.2k
700
12.2k
  s_fake(&vtmp, value, vbuf);
701
12.2k
702
12.2k
  return mp_int_sub(a, &vtmp, c);
703
12.2k
}
704
705
mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
706
13.2M
{
707
13.2M
  mp_digit *out;
708
13.2M
  mp_size   osize, ua, ub, p = 0;
709
13.2M
  mp_sign   osign;
710
13.2M
711
13.2M
  CHECK(a != NULL && b != NULL && c != NULL);
712
13.2M
713
13.2M
  /* If either input is zero, we can shortcut multiplication */
714
13.2M
  if (mp_int_compare_zero(a) == 0 || 
mp_int_compare_zero(b) == 011.0M
) {
715
7.93M
    mp_int_zero(c);
716
7.93M
    return MP_OK;
717
7.93M
  }
718
5.35M
719
5.35M
  /* Output is positive if inputs have same sign, otherwise negative */
720
5.35M
  osign = (MP_SIGN(a) == MP_SIGN(b)) ? 
MP_ZPOS3.13M
:
MP_NEG2.22M
;
721
5.35M
722
5.35M
  /* If the output is not identical to any of the inputs, we'll write the
723
5.35M
     results directly; otherwise, allocate a temporary space. */
724
5.35M
  ua = MP_USED(a); ub = MP_USED(b);
725
5.35M
  osize = MAX(ua, ub);
726
5.35M
  osize = 4 * ((osize + 1) / 2);
727
5.35M
728
5.35M
  if (c == a || 
c == b4.18M
) {
729
1.16M
    p = ROUND_PREC(osize);
730
1.16M
    p = MAX(p, default_precision);
731
1.16M
732
1.16M
    if ((out = s_alloc(p)) == NULL)
733
1.16M
      
return MP_MEMORY0
;
734
4.18M
  }
735
4.18M
  else {
736
4.18M
    if (!s_pad(c, osize))
737
0
      return MP_MEMORY;
738
4.18M
739
4.18M
    out = MP_DIGITS(c);
740
4.18M
  }
741
5.35M
  ZERO(out, osize);
742
5.35M
743
5.35M
  if (!s_kmul(MP_DIGITS(a), MP_DIGITS(b), out, ua, ub))
744
0
    return MP_MEMORY;
745
5.35M
746
5.35M
  /* If we allocated a new buffer, get rid of whatever memory c was already
747
5.35M
     using, and fix up its fields to reflect that.
748
5.35M
   */
749
5.35M
  if (out != MP_DIGITS(c)) {
750
1.16M
    if ((void *) MP_DIGITS(c) != (void *) c)
751
1.02M
      s_free(MP_DIGITS(c));
752
1.16M
    MP_DIGITS(c) = out;
753
1.16M
    MP_ALLOC(c) = p;
754
1.16M
  }
755
5.35M
756
5.35M
  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
757
5.35M
  CLAMP(c);           /* ... right here */
758
5.35M
  MP_SIGN(c) = osign;
759
5.35M
760
5.35M
  return MP_OK;
761
5.35M
}
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
1.19k
{
775
1.19k
  mp_result res;
776
1.19k
  CHECK(a != NULL && c != NULL && p2 >= 0);
777
1.19k
778
1.19k
  if ((res = mp_int_copy(a, c)) != MP_OK)
779
0
    return res;
780
1.19k
781
1.19k
  if (s_qmul(c, (mp_size) p2))
782
1.19k
    return MP_OK;
783
0
  else
784
0
    return MP_MEMORY;
785
1.19k
}
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
3.48M
{
832
3.48M
  int cmp, lg;
833
3.48M
  mp_result res = MP_OK;
834
3.48M
  mp_int qout, rout;
835
3.48M
  mp_sign sa = MP_SIGN(a), sb = MP_SIGN(b);
836
3.48M
  DECLARE_TEMP(2);
837
3.48M
838
3.48M
  CHECK(a != NULL && b != NULL && q != r);
839
3.48M
840
3.48M
  if (CMPZ(b) == 0)
841
0
    return MP_UNDEF;
842
3.48M
  else if ((cmp = s_ucmp(a, b)) < 0) {
843
640k
    /* If |a| < |b|, no division is required:
844
640k
       q = 0, r = a
845
640k
     */
846
640k
    if (r && 
(res = mp_int_copy(a, r)) != MP_OK19.7k
)
847
0
      return res;
848
640k
849
640k
    if (q)
850
634k
      mp_int_zero(q);
851
640k
852
640k
    return MP_OK;
853
640k
  }
854
2.84M
  else if (cmp == 0) {
855
601k
    /* If |a| = |b|, no division is required:
856
601k
       q = 1 or -1, r = 0
857
601k
     */
858
601k
    if (r)
859
2.08k
      mp_int_zero(r);
860
601k
861
601k
    if (q) {
862
600k
      mp_int_zero(q);
863
600k
      q->digits[0] = 1;
864
600k
865
600k
      if (sa != sb)
866
600k
  
MP_SIGN60.7k
(q) = MP_NEG60.7k
;
867
600k
    }
868
601k
869
601k
    return MP_OK;
870
601k
  }
871
2.24M
872
2.24M
  /* When |a| > |b|, real division is required.  We need someplace to store
873
2.24M
     quotient and remainder, but q and r are allowed to be NULL or to overlap
874
2.24M
     with the inputs.
875
2.24M
   */
876
2.24M
  if ((lg = s_isp2(b)) < 0) {
877
1.61M
    if (q && 
b != q1.60M
) {
878
1.50M
      if ((res = mp_int_copy(a, q)) != MP_OK)
879
0
  goto CLEANUP;
880
1.50M
      else
881
1.50M
  qout = q;
882
1.50M
    }
883
115k
    else {
884
115k
      qout = LAST_TEMP();
885
115k
      SETUP(mp_int_init_copy(LAST_TEMP(), a));
886
115k
    }
887
1.61M
888
1.61M
    if (r && 
a != r116k
) {
889
116k
      if ((res = mp_int_copy(b, r)) != MP_OK)
890
0
  goto CLEANUP;
891
116k
      else
892
116k
  rout = r;
893
116k
    }
894
1.50M
    else {
895
1.50M
      rout = LAST_TEMP();
896
1.50M
      SETUP(mp_int_init_copy(LAST_TEMP(), b));
897
1.50M
    }
898
1.61M
899
1.61M
    if ((res = s_udiv_knuth(qout, rout)) != MP_OK) 
goto CLEANUP0
;
900
622k
  }
901
622k
  else {
902
622k
    if (q && 
(res = mp_int_copy(a, q)) != MP_OK607k
)
goto CLEANUP0
;
903
622k
    if (r && 
(res = mp_int_copy(a, r)) != MP_OK59.5k
)
goto CLEANUP0
;
904
622k
905
622k
    if (q) 
s_qdiv(q, (mp_size) lg)607k
; qout = q;
906
622k
    if (r) 
s_qmod(r, (mp_size) lg)59.5k
; rout = r;
907
622k
  }
908
2.24M
909
2.24M
  /* Recompute signs for output */
910
2.24M
  if (rout) {
911
1.67M
    MP_SIGN(rout) = sa;
912
1.67M
    if (CMPZ(rout) == 0)
913
1.67M
      
MP_SIGN1.60M
(rout) = MP_ZPOS1.60M
;
914
1.67M
  }
915
2.24M
  if (qout) {
916
2.22M
    MP_SIGN(qout) = (sa == sb) ? 
MP_ZPOS1.55M
:
MP_NEG676k
;
917
2.22M
    if (CMPZ(qout) == 0)
918
2.22M
      
MP_SIGN0
(qout) = MP_ZPOS0
;
919
2.22M
  }
920
2.24M
921
2.24M
  if (q && 
(res = mp_int_copy(qout, q)) != MP_OK2.21M
)
goto CLEANUP0
;
922
2.24M
  if (r && 
(res = mp_int_copy(rout, r)) != MP_OK175k
)
goto CLEANUP0
;
923
2.24M
924
2.24M
  CLEANUP_TEMP();
925
2.24M
  return res;
926
2.24M
}
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
24.4k
{
959
24.4k
  mpz_t     vtmp, rtmp;
960
24.4k
  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
961
24.4k
  mp_result res;
962
24.4k
963
24.4k
  mp_int_init(&rtmp);
964
24.4k
  s_fake(&vtmp, value, vbuf);
965
24.4k
966
24.4k
  if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
967
0
    goto CLEANUP;
968
24.4k
969
24.4k
  if (r)
970
21.0k
    (void) mp_int_to_int(&rtmp, r); /* can't fail */
971
24.4k
972
24.4k
 CLEANUP:
973
24.4k
  mp_int_clear(&rtmp);
974
24.4k
  return res;
975
24.4k
}
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
159k
{
1094
159k
  mp_sign sa;
1095
159k
1096
159k
  CHECK(a != NULL && b != NULL);
1097
159k
1098
159k
  sa = MP_SIGN(a);
1099
159k
  if (sa == MP_SIGN(b)) {
1100
141k
    int cmp = s_ucmp(a, b);
1101
141k
1102
141k
    /* If they're both zero or positive, the normal comparison applies; if both
1103
141k
       negative, the sense is reversed. */
1104
141k
    if (sa == MP_ZPOS)
1105
98.1k
      return cmp;
1106
43.4k
    else
1107
43.4k
      return -cmp;
1108
17.8k
1109
17.8k
  }
1110
17.8k
  else {
1111
17.8k
    if (sa == MP_ZPOS)
1112
8.48k
      return 1;
1113
9.39k
    else
1114
9.39k
      return -1;
1115
17.8k
  }
1116
159k
}
1117
1118
int       mp_int_compare_unsigned(mp_int a, mp_int b)
1119
2.24M
{
1120
2.24M
  NRCHECK(a != NULL && b != NULL);
1121
2.24M
1122
2.24M
  return s_ucmp(a, b);
1123
2.24M
}
1124
1125
int       mp_int_compare_zero(mp_int z)
1126
37.6M
{
1127
37.6M
  NRCHECK(z != NULL);
1128
37.6M
1129
37.6M
  if (MP_USED(z) == 1 && 
z->digits[0] == 017.2M
)
1130
8.06M
    return 0;
1131
29.6M
  else if (MP_SIGN(z) == MP_ZPOS)
1132
21.1M
    return 1;
1133
8.45M
  else
1134
8.45M
    return -1;
1135
37.6M
}
1136
1137
int       mp_int_compare_value(mp_int z, mp_small value)
1138
45.0M
{
1139
45.0M
  mp_sign vsign = (value < 0) ? 
MP_NEG21.6M
:
MP_ZPOS23.3M
;
1140
45.0M
  int cmp;
1141
45.0M
1142
45.0M
  CHECK(z != NULL);
1143
45.0M
1144
45.0M
  if (vsign == MP_SIGN(z)) {
1145
28.5M
    cmp = s_vcmp(z, value);
1146
28.5M
1147
28.5M
    return (vsign == MP_ZPOS) ? 
cmp23.3M
:
-cmp5.25M
;
1148
28.5M
  }
1149
16.4M
  else {
1150
16.4M
    return (value < 0) ? 
116.4M
:
-118.1k
;
1151
16.4M
  }
1152
45.0M
}
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.32M
{
1314
2.32M
  int ca, cb, k = 0;
1315
2.32M
  mpz_t u, v, t;
1316
2.32M
  mp_result res;
1317
2.32M
1318
2.32M
  CHECK(a != NULL && b != NULL && c != NULL);
1319
2.32M
1320
2.32M
  ca = CMPZ(a);
1321
2.32M
  cb = CMPZ(b);
1322
2.32M
  if (ca == 0 && 
cb == 03
)
1323
0
    return MP_UNDEF;
1324
2.32M
  else if (ca == 0)
1325
3
    return mp_int_abs(b, c);
1326
2.32M
  else if (cb == 0)
1327
9.43k
    return mp_int_abs(a, c);
1328
2.31M
1329
2.31M
  mp_int_init(&t);
1330
2.31M
  if ((res = mp_int_init_copy(&u, a)) != MP_OK)
1331
0
    goto U;
1332
2.31M
  if ((res = mp_int_init_copy(&v, b)) != MP_OK)
1333
0
    goto V;
1334
2.31M
1335
2.31M
  MP_SIGN(&u) = MP_ZPOS; MP_SIGN(&v) = MP_ZPOS;
1336
2.31M
1337
2.31M
  { /* Divide out common factors of 2 from u and v */
1338
2.31M
    int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
1339
2.31M
1340
2.31M
    k = MIN(div2_u, div2_v);
1341
2.31M
    s_qdiv(&u, (mp_size) k);
1342
2.31M
    s_qdiv(&v, (mp_size) k);
1343
2.31M
  }
1344
2.31M
1345
2.31M
  if (mp_int_is_odd(&u)) {
1346
1.82M
    if ((res = mp_int_neg(&v, &t)) != MP_OK)
1347
0
      goto CLEANUP;
1348
486k
  }
1349
486k
  else {
1350
486k
    if ((res = mp_int_copy(&u, &t)) != MP_OK)
1351
0
      goto CLEANUP;
1352
2.31M
  }
1353
2.31M
1354
70.9M
  
for (;;)2.31M
{
1355
70.9M
    s_qdiv(&t, s_dp2k(&t));
1356
70.9M
1357
70.9M
    if (CMPZ(&t) > 0) {
1358
10.0M
      if ((res = mp_int_copy(&t, &u)) != MP_OK)
1359
0
  goto CLEANUP;
1360
60.8M
    }
1361
60.8M
    else {
1362
60.8M
      if ((res = mp_int_neg(&t, &v)) != MP_OK)
1363
0
  goto CLEANUP;
1364
70.9M
    }
1365
70.9M
1366
70.9M
    if ((res = mp_int_sub(&u, &v, &t)) != MP_OK)
1367
0
      goto CLEANUP;
1368
70.9M
1369
70.9M
    if (CMPZ(&t) == 0)
1370
2.31M
      break;
1371
70.9M
  }
1372
2.31M
1373
2.31M
  if ((res = mp_int_abs(&u, c)) != MP_OK)
1374
0
    goto CLEANUP;
1375
2.31M
  if (!s_qmul(c, (mp_size) k))
1376
0
    res = MP_MEMORY;
1377
2.31M
1378
2.31M
 CLEANUP:
1379
2.31M
  mp_int_clear(&v);
1380
2.31M
 V: mp_int_clear(&u);
1381
2.31M
 U: mp_int_clear(&t);
1382
2.31M
1383
2.31M
  return res;
1384
2.31M
}
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
250k
{
1500
250k
  mpz_t lcm;
1501
250k
  mp_result res;
1502
250k
1503
250k
  CHECK(a != NULL && b != NULL && c != NULL);
1504
250k
1505
250k
  /* Since a * b = gcd(a, b) * lcm(a, b), we can compute
1506
250k
     lcm(a, b) = (a / gcd(a, b)) * b.
1507
250k
1508
250k
     This formulation insures everything works even if the input
1509
250k
     variables share space.
1510
250k
   */
1511
250k
  if ((res = mp_int_init(&lcm)) != MP_OK)
1512
0
    return res;
1513
250k
  if ((res = mp_int_gcd(a, b, &lcm)) != MP_OK)
1514
0
    goto CLEANUP;
1515
250k
  if ((res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK)
1516
0
    goto CLEANUP;
1517
250k
  if ((res = mp_int_mul(&lcm, b, &lcm)) != MP_OK)
1518
0
    goto CLEANUP;
1519
250k
1520
250k
  res = mp_int_copy(&lcm, c);
1521
250k
1522
250k
  CLEANUP:
1523
250k
    mp_int_clear(&lcm);
1524
250k
1525
250k
  return res;
1526
250k
}
1527
1528
int       mp_int_divisible_value(mp_int a, mp_small v)
1529
21.0k
{
1530
21.0k
  mp_small rem = 0;
1531
21.0k
1532
21.0k
  if (mp_int_div_value(a, v, NULL, &rem) != MP_OK)
1533
0
    return 0;
1534
21.0k
1535
21.0k
  return rem == 0;
1536
21.0k
}
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
25.7M
{
1615
25.7M
  mp_usmall uv = 0;
1616
25.7M
  mp_size   uz;
1617
25.7M
  mp_digit *dz;
1618
25.7M
  mp_sign   sz;
1619
25.7M
1620
25.7M
  CHECK(z != NULL);
1621
25.7M
1622
25.7M
  /* Make sure the value is representable as a small integer */
1623
25.7M
  sz = MP_SIGN(z);
1624
25.7M
  if ((sz == MP_ZPOS && 
mp_int_compare_value(z, 20.5M
MP_SMALL_MAX20.5M
) > 0) ||
1625
25.7M
      
mp_int_compare_value(z, 21.6M
MP_SMALL_MIN21.6M
) < 0)
1626
6.23M
    return MP_RANGE;
1627
19.5M
1628
19.5M
  uz = MP_USED(z);
1629
19.5M
  dz = MP_DIGITS(z) + uz - 1;
1630
19.5M
1631
44.5M
  while (uz > 0) {
1632
24.9M
    uv <<= MP_DIGIT_BIT/2;
1633
24.9M
    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1634
24.9M
    --uz;
1635
24.9M
  }
1636
19.5M
1637
19.5M
  if (out)
1638
19.5M
    *out = (mp_small)((sz == MP_NEG) ? 
-uv3.14M
:
uv16.4M
);
1639
19.5M
1640
19.5M
  return MP_OK;
1641
19.5M
}
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
571
{
1675
571
  mp_result res;
1676
571
  int       cmp = 0;
1677
571
1678
571
  CHECK(z != NULL && str != NULL && limit >= 2);
1679
571
1680
571
  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1681
571
    
return MP_RANGE0
;
1682
571
1683
571
  if (CMPZ(z) == 0) {
1684
0
    *str++ = s_val2ch(0, 1);
1685
0
  }
1686
571
  else {
1687
571
    mpz_t tmp;
1688
571
    char  *h, *t;
1689
571
1690
571
    if ((res = mp_int_init_copy(&tmp, z)) != MP_OK)
1691
0
      return res;
1692
571
1693
571
    if (MP_SIGN(z) == MP_NEG) {
1694
303
      *str++ = '-';
1695
303
      --limit;
1696
303
    }
1697
571
    h = str;
1698
571
1699
571
    /* Generate digits in reverse order until finished or limit reached */
1700
11.6k
    for (/* */; limit > 0; 
--limit11.0k
) {
1701
11.6k
      mp_digit d;
1702
11.6k
1703
11.6k
      if ((cmp = CMPZ(&tmp)) == 0)
1704
571
  break;
1705
11.0k
1706
11.0k
      d = s_ddiv(&tmp, (mp_digit)radix);
1707
11.0k
      *str++ = s_val2ch(d, 1);
1708
11.0k
    }
1709
571
    t = str - 1;
1710
571
1711
571
    /* Put digits back in correct output order */
1712
5.87k
    while (h < t) {
1713
5.30k
      char tc = *h;
1714
5.30k
      *h++ = *t;
1715
5.30k
      *t-- = tc;
1716
5.30k
    }
1717
571
1718
571
    mp_int_clear(&tmp);
1719
571
  }
1720
571
1721
571
  *str = '\0';
1722
571
  if (cmp == 0)
1723
571
    return MP_OK;
1724
0
  else
1725
0
    return MP_TRUNC;
1726
571
}
1727
1728
mp_result mp_int_string_len(mp_int z, mp_size radix)
1729
593
{
1730
593
  int  len;
1731
593
1732
593
  CHECK(z != NULL);
1733
593
1734
593
  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1735
593
    
return MP_RANGE0
;
1736
593
1737
593
  len = s_outlen(z, radix) + 1; /* for terminator */
1738
593
1739
593
  /* Allow for sign marker on negatives */
1740
593
  if (MP_SIGN(z) == MP_NEG)
1741
310
    len += 1;
1742
593
1743
593
  return len;
1744
593
}
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
72
{
1749
72
  return mp_int_read_cstring(z, radix, str, NULL);
1750
72
}
1751
1752
mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
1753
72
{
1754
72
  int ch;
1755
72
1756
72
  CHECK(z != NULL && str != NULL);
1757
72
1758
72
  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1759
72
    
return MP_RANGE0
;
1760
72
1761
72
  /* Skip leading whitespace */
1762
72
  while (isspace((int)*str))
1763
0
    ++str;
1764
72
1765
72
  /* Handle leading sign tag (+/-, positive default) */
1766
72
  switch (*str) {
1767
72
  case '-':
1768
16
    MP_SIGN(z) = MP_NEG;
1769
16
    ++str;
1770
16
    break;
1771
72
  case '+':
1772
0
    ++str; /* fallthrough */
1773
56
  default:
1774
56
    MP_SIGN(z) = MP_ZPOS;
1775
56
    break;
1776
72
  }
1777
72
1778
72
  /* Skip leading zeroes */
1779
72
  while ((ch = s_ch2val(*str, radix)) == 0)
1780
0
    ++str;
1781
72
1782
72
  /* Make sure there is enough space for the value */
1783
72
  if (!s_pad(z, s_inlen(strlen(str), radix)))
1784
0
    return MP_MEMORY;
1785
72
1786
72
  MP_USED(z) = 1; z->digits[0] = 0;
1787
72
1788
828
  while (*str != '\0' && 
((ch = s_ch2val(*str, radix)) >= 0)756
) {
1789
756
    s_dmul(z, (mp_digit)radix);
1790
756
    s_dadd(z, (mp_digit)ch);
1791
756
    ++str;
1792
756
  }
1793
72
1794
72
  CLAMP(z);
1795
72
1796
72
  /* Override sign for zero, even if negative specified. */
1797
72
  if (CMPZ(z) == 0)
1798
72
    
MP_SIGN0
(z) = MP_ZPOS0
;
1799
72
1800
72
  if (end != NULL)
1801
72
    
*end = (char *)str0
;
1802
72
1803
72
  /* Return a truncation error if the string has unprocessed characters
1804
72
     remaining, so the caller can tell if the whole string was done */
1805
72
  if (*str != '\0')
1806
0
    return MP_TRUNC;
1807
72
  else
1808
72
    return MP_OK;
1809
72
}
1810
1811
mp_result mp_int_count_bits(mp_int z)
1812
1.19k
{
1813
1.19k
  mp_size  nbits = 0, uz;
1814
1.19k
  mp_digit d;
1815
1.19k
1816
1.19k
  CHECK(z != NULL);
1817
1.19k
1818
1.19k
  uz = MP_USED(z);
1819
1.19k
  if (uz == 1 && 
z->digits[0] == 0708
)
1820
0
    return 1;
1821
1.19k
1822
1.19k
  --uz;
1823
1.19k
  nbits = uz * MP_DIGIT_BIT;
1824
1.19k
  d = z->digits[uz];
1825
1.19k
1826
21.9k
  while (d != 0) {
1827
20.7k
    d >>= 1;
1828
20.7k
    ++nbits;
1829
20.7k
  }
1830
1.19k
1831
1.19k
  return nbits;
1832
1.19k
}
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
604
{
1938
604
  mp_result  res = mp_int_count_bits(z);
1939
604
  int        bytes;
1940
604
1941
604
  if (res <= 0)
1942
0
    return res;
1943
604
1944
604
  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
1945
604
1946
604
  return bytes;
1947
604
}
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
16.7M
{
1970
16.7M
  mp_digit *out = malloc(num * sizeof(mp_digit));
1971
16.7M
1972
16.7M
  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
16.7M
  return out;
1984
16.7M
}
1985
1986
STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
1987
605k
{
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
605k
1999
605k
  assert(new != NULL); /* for debugging */
2000
605k
#endif
2001
605k
  return new;
2002
605k
}
2003
2004
STATIC void s_free(void *ptr)
2005
16.7M
{
2006
16.7M
  free(ptr);
2007
16.7M
}
2008
2009
STATIC int      s_pad(mp_int z, mp_size min)
2010
169M
{
2011
169M
  if (MP_ALLOC(z) < min) {
2012
10.7M
    mp_size nsize = ROUND_PREC(min);
2013
10.7M
    mp_digit *tmp;
2014
10.7M
2015
10.7M
    if ((void *)z->digits == (void *)z) {
2016
10.1M
      if ((tmp = s_alloc(nsize)) == NULL)
2017
10.1M
        
return 00
;
2018
10.1M
2019
10.1M
      COPY(MP_DIGITS(z), tmp, MP_USED(z));
2020
10.1M
    }
2021
605k
    else if ((tmp = s_realloc(MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL)
2022
605k
      
return 00
;
2023
10.7M
2024
10.7M
    MP_DIGITS(z) = tmp;
2025
10.7M
    MP_ALLOC(z) = nsize;
2026
10.7M
  }
2027
169M
2028
169M
  return 1;
2029
169M
}
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
829k
{
2034
829k
  mp_usmall uv = (mp_usmall) (value < 0) ? 
-value1.67k
:
value827k
;
2035
829k
  s_ufake(z, uv, vbuf);
2036
829k
  if (value < 0)
2037
1.67k
    z->sign = MP_NEG;
2038
829k
}
2039
2040
STATIC void      s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[])
2041
29.4M
{
2042
29.4M
  mp_size ndig = (mp_size) s_uvpack(value, vbuf);
2043
29.4M
2044
29.4M
  z->used = ndig;
2045
29.4M
  z->alloc = MP_VALUE_DIGITS(value);
2046
29.4M
  z->sign = MP_ZPOS;
2047
29.4M
  z->digits = vbuf;
2048
29.4M
}
2049
2050
STATIC int      s_cdig(mp_digit *da, mp_digit *db, mp_size len)
2051
43.7M
{
2052
43.7M
  mp_digit *dat = da + len - 1, *dbt = db + len - 1;
2053
43.7M
2054
58.5M
  for (/* */; len != 0; 
--len, --dat, --dbt14.7M
) {
2055
52.7M
    if (*dat > *dbt)
2056
9.40M
      return 1;
2057
43.3M
    else if (*dat < *dbt)
2058
28.5M
      return -1;
2059
52.7M
  }
2060
43.7M
2061
43.7M
  
return 05.77M
;
2062
43.7M
}
2063
2064
STATIC int       s_uvpack(mp_usmall uv, mp_digit t[])
2065
29.4M
{
2066
29.4M
  int ndig = 0;
2067
29.4M
2068
29.4M
  if (uv == 0)
2069
765k
    t[ndig++] = 0;
2070
28.6M
  else {
2071
83.1M
    while (uv != 0) {
2072
54.4M
      t[ndig++] = (mp_digit) uv;
2073
54.4M
      uv >>= MP_DIGIT_BIT/2;
2074
54.4M
      uv >>= MP_DIGIT_BIT/2;
2075
54.4M
    }
2076
28.6M
  }
2077
29.4M
2078
29.4M
  return ndig;
2079
29.4M
}
2080
2081
STATIC int      s_ucmp(mp_int a, mp_int b)
2082
110M
{
2083
110M
  mp_size  ua = MP_USED(a), ub = MP_USED(b);
2084
110M
2085
110M
  if (ua > ub)
2086
14.8M
    return 1;
2087
95.9M
  else if (ub > ua)
2088
52.2M
    return -1;
2089
43.7M
  else
2090
43.7M
    return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua);
2091
110M
}
2092
2093
STATIC int      s_vcmp(mp_int a, mp_small v)
2094
28.5M
{
2095
28.5M
  mp_usmall uv = (v < 0) ? 
-(mp_usmall) v5.25M
:
(mp_usmall) v23.3M
;
2096
28.5M
  return s_uvcmp(a, uv);
2097
28.5M
}
2098
2099
STATIC int      s_uvcmp(mp_int a, mp_usmall uv)
2100
28.5M
{
2101
28.5M
  mpz_t vtmp;
2102
28.5M
  mp_digit vdig[MP_VALUE_DIGITS(uv)];
2103
28.5M
2104
28.5M
  s_ufake(&vtmp, uv, vdig);
2105
28.5M
  return s_ucmp(a, &vtmp);
2106
28.5M
}
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
3.86M
{
2111
3.86M
  mp_size pos;
2112
3.86M
  mp_word w = 0;
2113
3.86M
2114
3.86M
  /* Insure that da is the longer of the two to simplify later code */
2115
3.86M
  if (size_b > size_a) {
2116
600k
    SWAP(mp_digit *, da, db);
2117
600k
    SWAP(mp_size, size_a, size_b);
2118
600k
  }
2119
3.86M
2120
3.86M
  /* Add corresponding digits until the shorter number runs out */
2121
8.70M
  for (pos = 0; pos < size_b; 
++pos, ++da, ++db, ++dc4.83M
) {
2122
4.83M
    w = w + (mp_word) *da + (mp_word) *db;
2123
4.83M
    *dc = LOWER_HALF(w);
2124
4.83M
    w = UPPER_HALF(w);
2125
4.83M
  }
2126
3.86M
2127
3.86M
  /* Propagate carries as far as necessary */
2128
7.31M
  for (/* */; pos < size_a; 
++pos, ++da, ++dc3.45M
) {
2129
3.45M
    w = w + *da;
2130
3.45M
2131
3.45M
    *dc = LOWER_HALF(w);
2132
3.45M
    w = UPPER_HALF(w);
2133
3.45M
  }
2134
3.86M
2135
3.86M
  /* Return carry out */
2136
3.86M
  return (mp_digit)w;
2137
3.86M
}
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
75.6M
{
2142
75.6M
  mp_size pos;
2143
75.6M
  mp_word w = 0;
2144
75.6M
2145
75.6M
  /* We assume that |a| >= |b| so this should definitely hold */
2146
75.6M
  assert(size_a >= size_b);
2147
75.6M
2148
75.6M
  /* Subtract corresponding digits and propagate borrow */
2149
205M
  for (pos = 0; pos < size_b; 
++pos, ++da, ++db, ++dc130M
) {
2150
130M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2151
130M
   (mp_word)*da) - w - (mp_word)*db;
2152
130M
2153
130M
    *dc = LOWER_HALF(w);
2154
130M
    w = (UPPER_HALF(w) == 0);
2155
130M
  }
2156
75.6M
2157
75.6M
  /* Finish the subtraction for remaining upper digits of da */
2158
159M
  for (/* */; pos < size_a; 
++pos, ++da, ++dc83.4M
) {
2159
83.4M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2160
83.4M
   (mp_word)*da) - w;
2161
83.4M
2162
83.4M
    *dc = LOWER_HALF(w);
2163
83.4M
    w = (UPPER_HALF(w) == 0);
2164
83.4M
  }
2165
75.6M
2166
75.6M
  /* If there is a borrow out at the end, it violates the precondition */
2167
75.6M
  assert(w == 0);
2168
75.6M
}
2169
2170
STATIC int       s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
2171
      mp_size size_a, mp_size size_b)
2172
5.35M
{
2173
5.35M
  mp_size  bot_size;
2174
5.35M
2175
5.35M
  /* Make sure b is the smaller of the two input values */
2176
5.35M
  if (size_b > size_a) {
2177
2.32M
    SWAP(mp_digit *, da, db);
2178
2.32M
    SWAP(mp_size, size_a, size_b);
2179
2.32M
  }
2180
5.35M
2181
5.35M
  /* Insure that the bottom is the larger half in an odd-length split; the code
2182
5.35M
     below relies on this being true.
2183
5.35M
   */
2184
5.35M
  bot_size = (size_a + 1) / 2;
2185
5.35M
2186
5.35M
  /* If the values are big enough to bother with recursion, use the Karatsuba
2187
5.35M
     algorithm to compute the product; otherwise use the normal multiplication
2188
5.35M
     algorithm
2189
5.35M
   */
2190
5.35M
  if (multiply_threshold &&
2191
5.35M
      size_a >= multiply_threshold &&
2192
5.35M
      
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
5.35M
  else {
2249
5.35M
    s_umul(da, db, dc, size_a, size_b);
2250
5.35M
  }
2251
5.35M
2252
5.35M
  return 1;
2253
5.35M
}
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
5.35M
{
2258
5.35M
  mp_size a, b;
2259
5.35M
  mp_word w;
2260
5.35M
2261
19.7M
  for (a = 0; a < size_a; 
++a, ++dc, ++da14.4M
) {
2262
14.4M
    mp_digit *dct = dc;
2263
14.4M
    mp_digit *dbt = db;
2264
14.4M
2265
14.4M
    if (*da == 0)
2266
439k
      continue;
2267
13.9M
2268
13.9M
    w = 0;
2269
34.3M
    for (b = 0; b < size_b; 
++b, ++dbt, ++dct20.3M
) {
2270
20.3M
      w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
2271
20.3M
2272
20.3M
      *dct = LOWER_HALF(w);
2273
20.3M
      w = UPPER_HALF(w);
2274
20.3M
    }
2275
13.9M
2276
13.9M
    *dct = (mp_digit)w;
2277
13.9M
  }
2278
5.35M
}
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
756
{
2387
756
  mp_word w = 0;
2388
756
  mp_digit *da = MP_DIGITS(a);
2389
756
  mp_size ua = MP_USED(a);
2390
756
2391
756
  w = (mp_word)*da + b;
2392
756
  *da++ = LOWER_HALF(w);
2393
756
  w = UPPER_HALF(w);
2394
756
2395
796
  for (ua -= 1; ua > 0; 
--ua, ++da40
) {
2396
40
    w = (mp_word)*da + w;
2397
40
2398
40
    *da = LOWER_HALF(w);
2399
40
    w = UPPER_HALF(w);
2400
40
  }
2401
756
2402
756
  if (w) {
2403
2
    *da = (mp_digit)w;
2404
2
    MP_USED(a) += 1;
2405
2
  }
2406
756
}
2407
2408
STATIC void      s_dmul(mp_int a, mp_digit b)
2409
756
{
2410
756
  mp_word w = 0;
2411
756
  mp_digit *da = MP_DIGITS(a);
2412
756
  mp_size ua = MP_USED(a);
2413
756
2414
1.54k
  while (ua > 0) {
2415
792
    w = (mp_word)*da * b + w;
2416
792
    *da++ = LOWER_HALF(w);
2417
792
    w = UPPER_HALF(w);
2418
792
    --ua;
2419
792
  }
2420
756
2421
756
  if (w) {
2422
4
    *da = (mp_digit)w;
2423
4
    MP_USED(a) += 1;
2424
4
  }
2425
756
}
2426
2427
STATIC void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
2428
2.07M
{
2429
2.07M
  mp_word  w = 0;
2430
2.07M
2431
10.5M
  while (size_a > 0) {
2432
8.52M
    w = (mp_word)*da++ * (mp_word)b + w;
2433
8.52M
2434
8.52M
    *dc++ = LOWER_HALF(w);
2435
8.52M
    w = UPPER_HALF(w);
2436
8.52M
    --size_a;
2437
8.52M
  }
2438
2.07M
2439
2.07M
  if (w)
2440
0
    *dc = LOWER_HALF(w);
2441
2.07M
}
2442
2443
STATIC mp_digit  s_ddiv(mp_int a, mp_digit b)
2444
791k
{
2445
791k
  mp_word w = 0, qdigit;
2446
791k
  mp_size ua = MP_USED(a);
2447
791k
  mp_digit *da = MP_DIGITS(a) + ua - 1;
2448
791k
2449
2.74M
  for (/* */; ua > 0; 
--ua, --da1.94M
) {
2450
1.94M
    w = (w << MP_DIGIT_BIT) | *da;
2451
1.94M
2452
1.94M
    if (w >= b) {
2453
1.49M
      qdigit = w / b;
2454
1.49M
      w = w % b;
2455
1.49M
    }
2456
454k
    else {
2457
454k
      qdigit = 0;
2458
454k
    }
2459
1.94M
2460
1.94M
    *da = (mp_digit)qdigit;
2461
1.94M
  }
2462
791k
2463
791k
  CLAMP(a);
2464
791k
  return (mp_digit)w;
2465
791k
}
2466
2467
STATIC void     s_qdiv(mp_int z, mp_size p2)
2468
76.9M
{
2469
76.9M
  mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
2470
76.9M
  mp_size uz = MP_USED(z);
2471
76.9M
2472
76.9M
  if (ndig) {
2473
85.5k
    mp_size  mark;
2474
85.5k
    mp_digit *to, *from;
2475
85.5k
2476
85.5k
    if (ndig >= uz) {
2477
0
      mp_int_zero(z);
2478
0
      return;
2479
0
    }
2480
85.5k
2481
85.5k
    to = MP_DIGITS(z); from = to + ndig;
2482
85.5k
2483
203k
    for (mark = ndig; mark < uz; 
++mark117k
)
2484
117k
      *to++ = *from++;
2485
85.5k
2486
85.5k
    MP_USED(z) = uz - ndig;
2487
85.5k
  }
2488
76.9M
2489
76.9M
  if (nbits) {
2490
72.9M
    mp_digit d = 0, *dz, save;
2491
72.9M
    mp_size  up = MP_DIGIT_BIT - nbits;
2492
72.9M
2493
72.9M
    uz = MP_USED(z);
2494
72.9M
    dz = MP_DIGITS(z) + uz - 1;
2495
72.9M
2496
277M
    for (/* */; uz > 0; 
--uz, --dz204M
) {
2497
204M
      save = *dz;
2498
204M
2499
204M
      *dz = (*dz >> nbits) | (d << up);
2500
204M
      d = save;
2501
204M
    }
2502
72.9M
2503
72.9M
    CLAMP(z);
2504
72.9M
  }
2505
76.9M
2506
76.9M
  if (MP_USED(z) == 1 && 
z->digits[0] == 025.3M
)
2507
76.9M
    
MP_SIGN775k
(z) = MP_ZPOS775k
;
2508
76.9M
}
2509
2510
STATIC void     s_qmod(mp_int z, mp_size p2)
2511
59.5k
{
2512
59.5k
  mp_size start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
2513
59.5k
  mp_size uz = MP_USED(z);
2514
59.5k
  mp_digit mask = (1u << rest) - 1;
2515
59.5k
2516
59.5k
  if (start <= uz) {
2517
59.5k
    MP_USED(z) = start;
2518
59.5k
    z->digits[start - 1] &= mask;
2519
59.5k
    CLAMP(z);
2520
59.5k
  }
2521
59.5k
}
2522
2523
STATIC int      s_qmul(mp_int z, mp_size p2)
2524
3.96M
{
2525
3.96M
  mp_size   uz, need, rest, extra, i;
2526
3.96M
  mp_digit *from, *to, d;
2527
3.96M
2528
3.96M
  if (p2 == 0)
2529
1.31M
    return 1;
2530
2.64M
2531
2.64M
  uz = MP_USED(z); 
2532
2.64M
  need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT;
2533
2.64M
2534
2.64M
  /* Figure out if we need an extra digit at the top end; this occurs if the
2535
2.64M
     topmost `rest' bits of the high-order digit of z are not zero, meaning
2536
2.64M
     they will be shifted off the end if not preserved */
2537
2.64M
  extra = 0;
2538
2.64M
  if (rest != 0) {
2539
2.62M
    mp_digit *dz = MP_DIGITS(z) + uz - 1;
2540
2.62M
2541
2.62M
    if ((*dz >> (MP_DIGIT_BIT - rest)) != 0)
2542
577k
      extra = 1;
2543
2.62M
  }
2544
2.64M
2545
2.64M
  if (!s_pad(z, uz + need + extra))
2546
0
    return 0;
2547
2.64M
2548
2.64M
  /* If we need to shift by whole digits, do that in one pass, then
2549
2.64M
     to back and shift by partial digits.
2550
2.64M
   */
2551
2.64M
  if (need > 0) {
2552
27.9k
    from = MP_DIGITS(z) + uz - 1;
2553
27.9k
    to = from + need;
2554
27.9k
2555
55.9k
    for (i = 0; i < uz; 
++i28.0k
)
2556
28.0k
      *to-- = *from--;
2557
27.9k
2558
27.9k
    ZERO(MP_DIGITS(z), need);
2559
27.9k
    uz += need;
2560
27.9k
  }
2561
2.64M
2562
2.64M
  if (rest) {
2563
2.62M
    d = 0;
2564
9.89M
    for (i = need, from = 
MP_DIGITS2.62M
(z) + need; i < uz;
++i, ++from7.27M
) {
2565
7.27M
      mp_digit save = *from;
2566
7.27M
      
2567
7.27M
      *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
2568
7.27M
      d = save;
2569
7.27M
    }
2570
2.62M
2571
2.62M
    d >>= (MP_DIGIT_BIT - rest);
2572
2.62M
    if (d != 0) {
2573
577k
      *from = d;
2574
577k
      uz += extra;
2575
577k
    }
2576
2.62M
  }
2577
2.64M
2578
2.64M
  MP_USED(z) = uz;
2579
2.64M
  CLAMP(z);
2580
2.64M
2581
2.64M
  return 1;
2582
2.64M
}
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
75.5M
{
2616
75.5M
  int       k = 0;
2617
75.5M
  mp_digit *dp = MP_DIGITS(z), d;
2618
75.5M
2619
75.5M
  if (MP_USED(z) == 1 && 
*dp == 022.8M
)
2620
0
    return 1;
2621
75.5M
2622
75.6M
  
while (75.5M
*dp == 0) {
2623
146k
    k += MP_DIGIT_BIT;
2624
146k
    ++dp;
2625
146k
  }
2626
75.5M
2627
75.5M
  d = *dp;
2628
221M
  while ((d & 1) == 0) {
2629
145M
    d >>= 1;
2630
145M
    ++k;
2631
145M
  }
2632
75.5M
2633
75.5M
  return k;
2634
75.5M
}
2635
2636
STATIC int       s_isp2(mp_int z)
2637
2.24M
{
2638
2.24M
  mp_size uz = MP_USED(z), k = 0;
2639
2.24M
  mp_digit *dz = MP_DIGITS(z), d;
2640
2.24M
2641
2.26M
  while (uz > 1) {
2642
857k
    if (*dz++ != 0)
2643
838k
      return -1;
2644
19.5k
    k += MP_DIGIT_BIT;
2645
19.5k
    --uz;
2646
19.5k
  }
2647
2.24M
2648
2.24M
  d = *dz;
2649
2.99M
  while (d > 1) {
2650
2.37M
    if (d & 1)
2651
780k
      return -1;
2652
1.58M
    ++k; d >>= 1;
2653
1.58M
  }
2654
1.40M
2655
1.40M
  
return (int) k622k
;
2656
1.40M
}
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
839k
{
2679
839k
  mp_digit d = b->digits[MP_USED(b) - 1];
2680
839k
  int k = 0;
2681
839k
2682
14.9M
  while (d < (1u << (mp_digit)(MP_DIGIT_BIT - 1))) { /* d < (MP_RADIX / 2) */
2683
14.1M
    d <<= 1;
2684
14.1M
    ++k;
2685
14.1M
  }
2686
839k
2687
839k
  /* These multiplications can't fail */
2688
839k
  if (k != 0) {
2689
824k
    (void) s_qmul(a, (mp_size) k);
2690
824k
    (void) s_qmul(b, (mp_size) k);
2691
824k
  }
2692
839k
2693
839k
  return k;
2694
839k
}
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.61M
STATIC mp_result s_udiv_knuth(mp_int u, mp_int v) {
2840
1.61M
  mpz_t q, r, t;
2841
1.61M
  mp_result
2842
1.61M
  res = MP_OK;
2843
1.61M
  int k,j;
2844
1.61M
  mp_size m,n;
2845
1.61M
2846
1.61M
  /* Force signs to positive */
2847
1.61M
  MP_SIGN(u) = MP_ZPOS;
2848
1.61M
  MP_SIGN(v) = MP_ZPOS;
2849
1.61M
2850
1.61M
  /* Use simple division algorithm when v is only one digit long */
2851
1.61M
  if (MP_USED(v) == 1) {
2852
779k
    mp_digit d, rem;
2853
779k
    d   = v->digits[0];
2854
779k
    rem = s_ddiv(u, d);
2855
779k
    mp_int_set_value(v, rem);
2856
779k
    return MP_OK;
2857
779k
  }
2858
839k
2859
839k
  /* Algorithm D
2860
839k
2861
839k
     The n and m variables are defined as used by Knuth.
2862
839k
     u is an n digit number with digits u_{n-1}..u_0.
2863
839k
     v is an n+m digit number with digits from v_{m+n-1}..v_0.
2864
839k
     We require that n > 1 and m >= 0
2865
839k
   */
2866
839k
  n = MP_USED(v);
2867
839k
  m = MP_USED(u) - n;
2868
839k
  assert(n > 1);
2869
839k
  assert(m >= 0);
2870
839k
2871
839k
  /* D1: Normalize.
2872
839k
     The normalization step provides the necessary condition for Theorem B,
2873
839k
     which states that the quotient estimate for q_j, call it qhat
2874
839k
2875
839k
       qhat = u_{j+n}u_{j+n-1} / v_{n-1}
2876
839k
2877
839k
     is bounded by
2878
839k
2879
839k
      qhat - 2 <= q_j <= qhat.
2880
839k
2881
839k
     That is, qhat is always greater than the actual quotient digit q,
2882
839k
     and it is never more than two larger than the actual quotient digit.
2883
839k
   */
2884
839k
  k = s_norm(u, v);
2885
839k
2886
839k
  /* Extend size of u by one if needed.
2887
839k
2888
839k
     The algorithm begins with a value of u that has one more digit of input.
2889
839k
     The normalization step sets u_{m+n}..u_0 = 2^k * u_{m+n-1}..u_0. If the
2890
839k
     multiplication did not increase the number of digits of u, we need to add
2891
839k
     a leading zero here.
2892
839k
   */
2893
839k
  if (k == 0 || 
MP_USED824k
(u) != m + n + 1824k
) {
2894
316k
    if (!s_pad(u, m+n+1))
2895
0
      return MP_MEMORY;
2896
316k
    u->digits[m+n] = 0;
2897
316k
    u->used = m+n+1;
2898
316k
  }
2899
839k
2900
839k
  /* Add a leading 0 to v.
2901
839k
2902
839k
     The multiplication in step D4 multiplies qhat * 0v_{n-1}..v_0.  We need to
2903
839k
     add the leading zero to v here to ensure that the multiplication will
2904
839k
     produce the full n+1 digit result.
2905
839k
   */
2906
839k
  if (!s_pad(v, n+1)) 
return MP_MEMORY0
; v->digits[n] = 0;
2907
839k
2908
839k
  /* Initialize temporary variables q and t.
2909
839k
     q allocates space for m+1 digits to store the quotient digits
2910
839k
     t allocates space for n+1 digits to hold the result of q_j*v
2911
839k
   */
2912
839k
  if ((res = mp_int_init_size(&q, m + 1)) != MP_OK) 
return res0
;
2913
839k
  if ((res = mp_int_init_size(&t, n + 1)) != MP_OK) 
goto CLEANUP0
;
2914
839k
2915
839k
  /* D2: Initialize j */
2916
839k
  j = m;
2917
839k
  r.digits = MP_DIGITS(u) + j;  /* The contents of r are shared with u */
2918
839k
  r.used   = n + 1;
2919
839k
  r.sign   = MP_ZPOS;
2920
839k
  r.alloc  = MP_ALLOC(u);
2921
839k
  ZERO(t.digits, t.alloc);
2922
839k
2923
839k
  /* Calculate the m+1 digits of the quotient result */
2924
2.75M
  for (; j >= 0; 
j--1.91M
) {
2925
1.91M
    /* D3: Calculate q' */
2926
1.91M
    /* r->digits is aligned to position j of the number u */
2927
1.91M
    mp_word pfx, qhat;
2928
1.91M
    pfx   = r.digits[n];
2929
1.91M
    pfx <<= MP_DIGIT_BIT / 2;
2930
1.91M
    pfx <<= MP_DIGIT_BIT / 2;
2931
1.91M
    pfx |= r.digits[n-1]; /* pfx = u_{j+n}{j+n-1} */
2932
1.91M
2933
1.91M
    qhat = pfx / v->digits[n-1];
2934
1.91M
    /* Check to see if qhat > b, and decrease qhat if so.
2935
1.91M
       Theorem B guarantess that qhat is at most 2 larger than the
2936
1.91M
       actual value, so it is possible that qhat is greater than
2937
1.91M
       the maximum value that will fit in a digit */
2938
1.91M
    if (qhat > MP_DIGIT_MAX)
2939
1.91M
      
qhat = 55
MP_DIGIT_MAX55
;
2940
1.91M
2941
1.91M
    /* D4,D5,D6: Multiply qhat * v and test for a correct value of q
2942
1.91M
2943
1.91M
       We proceed a bit different than the way described by Knuth. This way is
2944
1.91M
       simpler but less efficent. Instead of doing the multiply and subtract
2945
1.91M
       then checking for underflow, we first do the multiply of qhat * v and
2946
1.91M
       see if it is larger than the current remainder r. If it is larger, we
2947
1.91M
       decrease qhat by one and try again. We may need to decrease qhat one
2948
1.91M
       more time before we get a value that is smaller than r.
2949
1.91M
2950
1.91M
       This way is less efficent than Knuth becuase we do more multiplies, but
2951
1.91M
       we do not need to worry about underflow this way.
2952
1.91M
     */
2953
1.91M
    /* t = qhat * v */
2954
1.91M
    s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1); t.used = n + 1;
2955
1.91M
    CLAMP(&t);
2956
1.91M
2957
1.91M
    /* Clamp r for the comparison. Comparisons do not like leading zeros. */
2958
1.91M
    CLAMP(&r);
2959
1.91M
    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.91M
    /* Unclamp r. The D algorithm expects r = u_{j+n}..u_j to always be n+1
2972
1.91M
       digits long. */
2973
1.91M
    r.used = n + 1;
2974
1.91M
2975
1.91M
    /* D4: Multiply and subtract
2976
1.91M
2977
1.91M
       Note: The multiply was completed above so we only need to subtract here.
2978
1.91M
     */
2979
1.91M
    s_usub(r.digits, t.digits, r.digits, r.used, t.used);
2980
1.91M
2981
1.91M
    /* D5: Test remainder
2982
1.91M
2983
1.91M
       Note: Not needed because we always check that qhat is the correct value
2984
1.91M
             before performing the subtract.  Value cast to mp_digit to prevent
2985
1.91M
             warning, qhat has been clamped to MP_DIGIT_MAX
2986
1.91M
     */
2987
1.91M
    q.digits[j] = (mp_digit)qhat;
2988
1.91M
2989
1.91M
    /* D6: Add back
2990
1.91M
       Note: Not needed because we always check that qhat is the correct value
2991
1.91M
             before performing the subtract.
2992
1.91M
     */
2993
1.91M
2994
1.91M
    /* D7: Loop on j */
2995
1.91M
    r.digits--;
2996
1.91M
    ZERO(t.digits, t.alloc);
2997
1.91M
  }
2998
839k
2999
839k
  /* Get rid of leading zeros in q */
3000
839k
  q.used = m + 1;
3001
839k
  CLAMP(&q);
3002
839k
3003
839k
  /* Denormalize the remainder */
3004
839k
  CLAMP(u); /* use u here because the r.digits pointer is off-by-one */
3005
839k
  if (k != 0)
3006
824k
    s_qdiv(u, k);
3007
839k
3008
839k
  mp_int_copy(u, v);  /* ok:  0 <= r < v */
3009
839k
  mp_int_copy(&q, u); /* ok:  q <= u     */
3010
839k
3011
839k
  mp_int_clear(&t);
3012
839k
 CLEANUP:
3013
839k
  mp_int_clear(&q);
3014
839k
  return res;
3015
839k
}
3016
3017
STATIC int       s_outlen(mp_int z, mp_size r)
3018
593
{
3019
593
  mp_result bits;
3020
593
  double raw;
3021
593
3022
593
  assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
3023
593
3024
593
  bits = mp_int_count_bits(z);
3025
593
  raw = (double)bits * s_log2[r];
3026
593
3027
593
  return (int)(raw + 0.999999);
3028
593
}
3029
3030
STATIC mp_size   s_inlen(int len, mp_size r)
3031
72
{
3032
72
  double  raw = (double)len / s_log2[r];
3033
72
  mp_size bits = (mp_size)(raw + 0.5);
3034
72
3035
72
  return (mp_size)((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT) + 1;
3036
72
}
3037
3038
STATIC int       s_ch2val(char c, int r)
3039
828
{
3040
828
  int out;
3041
828
3042
828
  if (isdigit((unsigned char) c))
3043
828
    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
828
3049
828
  return (out >= r) ? 
-10
: out;
3050
828
}
3051
3052
STATIC char      s_val2ch(int v, int caps)
3053
11.0k
{
3054
11.0k
  assert(v >= 0);
3055
11.0k
3056
11.0k
  if (v < 10)
3057
11.0k
    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
11.0k
}
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 */