Coverage Report

Created: 2017-06-23 12:40

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