Coverage Report

Created: 2017-03-28 09:59

/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
370M
#define CHECK(TEST)   assert(TEST)
73
82.2M
#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
36.0M
#define MP_VALUE_DIGITS(V) \
103
36.0M
((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
104
105
/* Round precision P to nearest word boundary */
106
20.2M
#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
107
108
/* Set array P of S digits to zero */
109
9.62M
#define ZERO(P, S) \
110
9.62M
do{ \
111
9.62M
  mp_size i__ = (S) * sizeof(mp_digit); \
112
9.62M
  mp_digit *p__ = (P); \
113
9.62M
  memset(p__, 0, i__); \
114
9.62M
} while(0)
115
116
/* Copy S digits from array P to array Q */
117
108M
#define COPY(P, Q, S) \
118
108M
do{ \
119
108M
  mp_size i__ = (S) * sizeof(mp_digit); \
120
108M
  mp_digit *p__ = (P), *q__ = (Q); \
121
108M
  memcpy(q__, p__, i__); \
122
108M
} 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
169M
#define CLAMP(Z) \
136
169M
do{ \
137
169M
  mp_int z_ = (Z); \
138
169M
  mp_size uz_ = MP_USED(z_); \
139
169M
  mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \
140
201M
  while (
uz_ > 1 && 201M
(*dz_-- == 0)144M
) \
141
31.2M
    --uz_; \
142
169M
  MP_USED(z_) = uz_; \
143
169M
} while(0)
144
145
/* Select min/max.  Do not provide expressions for which multiple
146
   evaluation would be problematic, e.g. x++ */
147
2.49M
#define MIN(A, B) 
((B)<(A)?2.49M
(B)504k
:
(A)1.98M
)
148
96.3M
#define MAX(A, B) 
((B)>(A)?96.3M
(B)46.3M
:
(A)50.0M
)
149
150
/* Exchange lvalues A and B of type T, e.g.
151
   SWAP(int, x, y) where x and y are variables of type int. */
152
7.38M
#define SWAP(T, A, B) \
153
7.38M
do{ \
154
7.38M
  T t_ = (A); \
155
7.38M
  A = (B); \
156
7.38M
  B = t_; \
157
7.38M
} while(0)
158
159
/* Used to set up and access simple temp stacks within functions. */
160
#define DECLARE_TEMP(N) \
161
3.87M
  mpz_t temp[(N)]; \
162
3.87M
  int last__ = 0
163
#define CLEANUP_TEMP() \
164
2.44M
 CLEANUP: \
165
4.08M
  while (--last__ >= 0) \
166
1.63M
    
mp_int_clear(1.63M
TEMP1.63M
(last__))
167
3.27M
#define TEMP(K) (temp + (K))
168
1.63M
#define LAST_TEMP() TEMP(last__)
169
1.63M
#define SETUP(E) \
170
1.63M
do{ \
171
1.63M
  if ((res = (E)) != MP_OK) \
172
0
    goto CLEANUP; \
173
1.63M
  ++(last__); \
174
1.63M
} while(0)
175
176
/* Compare value to zero. */
177
227M
#define CMPZ(Z) \
178
223M
(((Z)->used==1&&
(Z)->digits[0]==077.6M
)?
04.20M
:
((Z)->sign==MP_NEG)?223M
-1190M
:
132.2M
)
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
262M
#define UPPER_HALF(W)           
((mp_word)((W) >> 262M
MP_DIGIT_BIT262M
))
204
262M
#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
28.0M
{
361
28.0M
  if (z == NULL)
362
0
    return MP_BADARG;
363
28.0M
364
28.0M
  z->single = 0;
365
28.0M
  z->digits = &(z->single);
366
28.0M
  z->alloc  = 1;
367
28.0M
  z->used   = 1;
368
28.0M
  z->sign   = MP_ZPOS;
369
28.0M
370
28.0M
  return MP_OK;
371
28.0M
}
372
373
mp_int    mp_int_alloc(void)
374
21.9M
{
375
21.9M
  mp_int out = malloc(sizeof(mpz_t));
376
21.9M
377
21.9M
  if (out != NULL)
378
21.9M
    mp_int_init(out);
379
21.9M
380
21.9M
  return out;
381
21.9M
}
382
383
mp_result mp_int_init_size(mp_int z, mp_size prec)
384
5.93M
{
385
5.93M
  CHECK(z != NULL);
386
5.93M
387
5.93M
  if (prec == 0)
388
0
    prec = default_precision;
389
5.93M
  else 
if (5.93M
prec == 15.93M
)
390
332k
    return mp_int_init(z);
391
5.93M
  else
392
5.60M
    
prec = (mp_size) 5.60M
ROUND_PREC5.60M
(prec);
393
5.93M
394
5.60M
  
if (5.60M
(5.60M
MP_DIGITS5.60M
(z) = s_alloc(prec)) == NULL)
395
0
    return MP_MEMORY;
396
5.60M
397
5.60M
  z->digits[0] = 0;
398
5.60M
  MP_USED(z) = 1;
399
5.60M
  MP_ALLOC(z) = prec;
400
5.60M
  MP_SIGN(z) = MP_ZPOS;
401
5.60M
402
5.60M
  return MP_OK;
403
5.60M
}
404
405
mp_result mp_int_init_copy(mp_int z, mp_int old)
406
6.81M
{
407
6.81M
  mp_result res;
408
6.81M
  mp_size uold;
409
6.81M
410
6.81M
  CHECK(z != NULL && old != NULL);
411
6.81M
412
6.81M
  uold = MP_USED(old);
413
6.81M
  if (
uold == 16.81M
)
{2.56M
414
2.56M
    mp_int_init(z);
415
2.56M
  }
416
4.24M
  else {
417
4.24M
    mp_size target = MAX(uold, default_precision);
418
4.24M
419
4.24M
    if ((res = mp_int_init_size(z, target)) != MP_OK)
420
0
      return res;
421
4.24M
  }
422
6.81M
423
6.81M
  
MP_USED6.81M
(z) = uold;6.81M
424
6.81M
  
MP_SIGN6.81M
(z) =
MP_SIGN6.81M
(old);
425
6.81M
  COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
426
6.81M
427
6.81M
  return MP_OK;
428
6.81M
}
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
155k
{
441
155k
  mpz_t    vtmp;
442
155k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
443
155k
444
155k
  s_ufake(&vtmp, uvalue, vbuf);
445
155k
  return mp_int_init_copy(z, &vtmp);
446
155k
}
447
448
mp_result  mp_int_set_value(mp_int z, mp_small value)
449
840k
{
450
840k
  mpz_t    vtmp;
451
840k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
452
840k
453
840k
  s_fake(&vtmp, value, vbuf);
454
840k
  return mp_int_copy(&vtmp, z);
455
840k
}
456
457
mp_result  mp_int_set_uvalue(mp_int z, mp_usmall uvalue)
458
38.4k
{
459
38.4k
  mpz_t    vtmp;
460
38.4k
  mp_digit vbuf[MP_VALUE_DIGITS(uvalue)];
461
38.4k
462
38.4k
  s_ufake(&vtmp, uvalue, vbuf);
463
38.4k
  return mp_int_copy(&vtmp, z);
464
38.4k
}
465
466
void      mp_int_clear(mp_int z)
467
33.6M
{
468
33.6M
  if (z == NULL)
469
0
    return;
470
33.6M
471
33.6M
  
if (33.6M
MP_DIGITS33.6M
(z) != NULL33.6M
)
{33.6M
472
33.6M
    if (
MP_DIGITS33.6M
(z) != &(z->single)33.6M
)
473
18.3M
      
s_free(18.3M
MP_DIGITS18.3M
(z));
474
33.6M
475
33.6M
    MP_DIGITS(z) = NULL;
476
33.6M
  }
477
33.6M
}
478
479
void      mp_int_free(mp_int z)
480
21.9M
{
481
21.9M
  NRCHECK(z != NULL);
482
21.9M
483
21.9M
  mp_int_clear(z);
484
21.9M
  free(z); /* note: NOT s_free() */
485
21.9M
}
486
487
mp_result mp_int_copy(mp_int a, mp_int c)
488
94.3M
{
489
94.3M
  CHECK(a != NULL && c != NULL);
490
94.3M
491
94.3M
  if (
a != c94.3M
)
{89.3M
492
89.3M
    mp_size ua = MP_USED(a);
493
89.3M
    mp_digit *da, *dc;
494
89.3M
495
89.3M
    if (!s_pad(c, ua))
496
0
      return MP_MEMORY;
497
89.3M
498
89.3M
    
da = 89.3M
MP_DIGITS89.3M
(a); dc =
MP_DIGITS89.3M
(c);
499
89.3M
    COPY(da, dc, ua);
500
89.3M
501
89.3M
    MP_USED(c) = ua;
502
89.3M
    
MP_SIGN89.3M
(c) =
MP_SIGN89.3M
(a);
503
89.3M
  }
504
94.3M
505
94.3M
  return MP_OK;
506
94.3M
}
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
11.5M
{
525
11.5M
  NRCHECK(z != NULL);
526
11.5M
527
11.5M
  z->digits[0] = 0;
528
11.5M
  MP_USED(z) = 1;
529
11.5M
  MP_SIGN(z) = MP_ZPOS;
530
11.5M
}
531
532
mp_result mp_int_abs(mp_int a, mp_int c)
533
3.22M
{
534
3.22M
  mp_result res;
535
3.22M
536
3.22M
  CHECK(a != NULL && c != NULL);
537
3.22M
538
3.22M
  if ((res = mp_int_copy(a, c)) != MP_OK)
539
0
    return res;
540
3.22M
541
3.22M
  
MP_SIGN3.22M
(c) = MP_ZPOS;3.22M
542
3.22M
  return MP_OK;
543
3.22M
}
544
545
mp_result mp_int_neg(mp_int a, mp_int c)
546
66.4M
{
547
66.4M
  mp_result res;
548
66.4M
549
66.4M
  CHECK(a != NULL && c != NULL);
550
66.4M
551
66.4M
  if ((res = mp_int_copy(a, c)) != MP_OK)
552
0
    return res;
553
66.4M
554
66.4M
  
if (66.4M
CMPZ66.4M
(c) != 066.4M
)
555
66.4M
    
MP_SIGN66.4M
(c) = 1 - 66.4M
MP_SIGN66.4M
(a);
556
66.4M
557
66.4M
  return MP_OK;
558
66.4M
}
559
560
mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
561
9.37M
{
562
9.37M
  mp_size ua, ub, uc, max;
563
9.37M
564
9.37M
  CHECK(a != NULL && b != NULL && c != NULL);
565
9.37M
566
9.37M
  ua = 
MP_USED9.37M
(a); ub =
MP_USED9.37M
(b); uc =
MP_USED9.37M
(c);
567
9.37M
  max = MAX(ua, ub);
568
9.37M
569
9.37M
  if (
MP_SIGN9.37M
(a) == 9.37M
MP_SIGN9.37M
(b))
{5.46M
570
5.46M
    /* Same sign -- add magnitudes, preserve sign of addends */
571
5.46M
    mp_digit carry;
572
5.46M
573
5.46M
    if (!s_pad(c, max))
574
0
      return MP_MEMORY;
575
5.46M
576
5.46M
    
carry = s_uadd(5.46M
MP_DIGITS5.46M
(a),
MP_DIGITS5.46M
(b),
MP_DIGITS5.46M
(c), ua, ub);
577
5.46M
    uc = max;
578
5.46M
579
5.46M
    if (
carry5.46M
)
{44.5k
580
44.5k
      if (!s_pad(c, max + 1))
581
0
  return MP_MEMORY;
582
44.5k
583
44.5k
      c->digits[max] = carry;
584
44.5k
      ++uc;
585
44.5k
    }
586
5.46M
587
5.46M
    
MP_USED5.46M
(c) = uc;5.46M
588
5.46M
    
MP_SIGN5.46M
(c) =
MP_SIGN5.46M
(a);
589
5.46M
590
5.46M
  }
591
3.91M
  else {
592
3.91M
    /* Different signs -- subtract magnitudes, preserve sign of greater */
593
3.91M
    mp_int  x, y;
594
3.91M
    int     cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
595
3.91M
596
3.91M
    /* Set x to max(a, b), y to min(a, b) to simplify later code.
597
3.91M
       A special case yields zero for equal magnitudes.
598
3.91M
    */
599
3.91M
    if (
cmp == 03.91M
)
{727k
600
727k
      mp_int_zero(c);
601
727k
      return MP_OK;
602
727k
    }
603
3.18M
    else 
if (3.18M
cmp < 03.18M
)
{1.21M
604
1.21M
      x = b; y = a;
605
1.21M
    }
606
1.97M
    else {
607
1.97M
      x = a; y = b;
608
1.97M
    }
609
3.91M
610
3.18M
    
if (3.18M
!s_pad(c, 3.18M
MP_USED3.18M
(x)))
611
0
      return MP_MEMORY;
612
3.18M
613
3.18M
    /* Subtract smaller from larger */
614
3.18M
    
s_usub(3.18M
MP_DIGITS3.18M
(x),
MP_DIGITS3.18M
(y),
MP_DIGITS3.18M
(c),
MP_USED3.18M
(x),
MP_USED3.18M
(y));
615
3.18M
    
MP_USED3.18M
(c) =
MP_USED3.18M
(x);
616
3.18M
    CLAMP(c);
617
3.18M
618
3.18M
    /* Give result the sign of the larger */
619
3.18M
    
MP_SIGN3.18M
(c) =
MP_SIGN3.18M
(x);
620
3.18M
  }
621
9.37M
622
8.65M
  return MP_OK;
623
9.37M
}
624
625
mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c)
626
1.32k
{
627
1.32k
  mpz_t    vtmp;
628
1.32k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
629
1.32k
630
1.32k
  s_fake(&vtmp, value, vbuf);
631
1.32k
632
1.32k
  return mp_int_add(a, &vtmp, c);
633
1.32k
}
634
635
mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
636
74.5M
{
637
74.5M
  mp_size ua, ub, uc, max;
638
74.5M
639
74.5M
  CHECK(a != NULL && b != NULL && c != NULL);
640
74.5M
641
74.5M
  ua = 
MP_USED74.5M
(a); ub =
MP_USED74.5M
(b); uc =
MP_USED74.5M
(c);
642
74.5M
  max = MAX(ua, ub);
643
74.5M
644
74.5M
  if (
MP_SIGN74.5M
(a) != 74.5M
MP_SIGN74.5M
(b))
{260k
645
260k
    /* Different signs -- add magnitudes and keep sign of a */
646
260k
    mp_digit carry;
647
260k
648
260k
    if (!s_pad(c, max))
649
0
      return MP_MEMORY;
650
260k
651
260k
    
carry = s_uadd(260k
MP_DIGITS260k
(a),
MP_DIGITS260k
(b),
MP_DIGITS260k
(c), ua, ub);
652
260k
    uc = max;
653
260k
654
260k
    if (
carry260k
)
{2.08k
655
2.08k
      if (!s_pad(c, max + 1))
656
0
  return MP_MEMORY;
657
2.08k
658
2.08k
      c->digits[max] = carry;
659
2.08k
      ++uc;
660
2.08k
    }
661
260k
662
260k
    
MP_USED260k
(c) = uc;260k
663
260k
    
MP_SIGN260k
(c) =
MP_SIGN260k
(a);
664
260k
665
260k
  }
666
74.2M
  else {
667
74.2M
    /* Same signs -- subtract magnitudes */
668
74.2M
    mp_int  x, y;
669
74.2M
    mp_sign osign;
670
74.2M
    int     cmp = s_ucmp(a, b);
671
74.2M
672
74.2M
    if (!s_pad(c, max))
673
0
      return MP_MEMORY;
674
74.2M
675
74.2M
    
if (74.2M
cmp >= 074.2M
)
{12.5M
676
12.5M
      x = a; y = b; osign = MP_ZPOS;
677
12.5M
    }
678
61.7M
    else {
679
61.7M
      x = b; y = a; osign = MP_NEG;
680
61.7M
    }
681
74.2M
682
74.2M
    if (
MP_SIGN74.2M
(a) == MP_NEG && 74.2M
cmp != 091.5k
)
683
75.2k
      osign = 1 - osign;
684
74.2M
685
74.2M
    s_usub(
MP_DIGITS74.2M
(x),
MP_DIGITS74.2M
(y),
MP_DIGITS74.2M
(c),
MP_USED74.2M
(x),
MP_USED74.2M
(y));
686
74.2M
    
MP_USED74.2M
(c) =
MP_USED74.2M
(x);
687
74.2M
    CLAMP(c);
688
74.2M
689
74.2M
    MP_SIGN(c) = osign;
690
74.2M
  }
691
74.5M
692
74.5M
  return MP_OK;
693
74.5M
}
694
695
mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c)
696
13.6k
{
697
13.6k
  mpz_t    vtmp;
698
13.6k
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
699
13.6k
700
13.6k
  s_fake(&vtmp, value, vbuf);
701
13.6k
702
13.6k
  return mp_int_sub(a, &vtmp, c);
703
13.6k
}
704
705
mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
706
16.1M
{
707
16.1M
  mp_digit *out;
708
16.1M
  mp_size   osize, ua, ub, p = 0;
709
16.1M
  mp_sign   osign;
710
16.1M
711
16.1M
  CHECK(a != NULL && b != NULL && c != NULL);
712
16.1M
713
16.1M
  /* If either input is zero, we can shortcut multiplication */
714
16.1M
  if (
mp_int_compare_zero(a) == 0 || 16.1M
mp_int_compare_zero(b) == 013.6M
)
{9.36M
715
9.36M
    mp_int_zero(c);
716
9.36M
    return MP_OK;
717
9.36M
  }
718
16.1M
719
16.1M
  /* Output is positive if inputs have same sign, otherwise negative */
720
6.78M
  
osign = (6.78M
MP_SIGN6.78M
(a) ==
MP_SIGN6.78M
(b)) ?
MP_ZPOS4.02M
:
MP_NEG2.76M
;
721
6.78M
722
6.78M
  /* If the output is not identical to any of the inputs, we'll write the
723
6.78M
     results directly; otherwise, allocate a temporary space. */
724
6.78M
  ua = 
MP_USED6.78M
(a); ub =
MP_USED6.78M
(b);
725
6.78M
  osize = MAX(ua, ub);
726
6.78M
  osize = 4 * ((osize + 1) / 2);
727
6.78M
728
6.78M
  if (
c == a || 6.78M
c == b5.33M
)
{1.44M
729
1.44M
    p = ROUND_PREC(osize);
730
1.44M
    p = MAX(p, default_precision);
731
1.44M
732
1.44M
    if ((out = s_alloc(p)) == NULL)
733
0
      return MP_MEMORY;
734
1.44M
  }
735
5.33M
  else {
736
5.33M
    if (!s_pad(c, osize))
737
0
      return MP_MEMORY;
738
5.33M
739
5.33M
    
out = 5.33M
MP_DIGITS5.33M
(c);
740
5.33M
  }
741
6.78M
  
ZERO6.78M
(out, osize);6.78M
742
6.78M
743
6.78M
  if (
!s_kmul(6.78M
MP_DIGITS6.78M
(a),
MP_DIGITS6.78M
(b), out, ua, ub))
744
0
    return MP_MEMORY;
745
6.78M
746
6.78M
  /* If we allocated a new buffer, get rid of whatever memory c was already
747
6.78M
     using, and fix up its fields to reflect that.
748
6.78M
   */
749
6.78M
  
if (6.78M
out != 6.78M
MP_DIGITS6.78M
(c))
{1.44M
750
1.44M
    if (
(void *) 1.44M
MP_DIGITS1.44M
(c) != (void *) c)
751
1.26M
      
s_free(1.26M
MP_DIGITS1.26M
(c));
752
1.44M
    MP_DIGITS(c) = out;
753
1.44M
    MP_ALLOC(c) = p;
754
1.44M
  }
755
6.78M
756
6.78M
  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
757
6.78M
  CLAMP(c);           /* ... right here */
758
6.78M
  MP_SIGN(c) = osign;
759
6.78M
760
6.78M
  return MP_OK;
761
6.78M
}
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.85k
{
775
3.85k
  mp_result res;
776
3.85k
  CHECK(a != NULL && c != NULL && p2 >= 0);
777
3.85k
778
3.85k
  if ((res = mp_int_copy(a, c)) != MP_OK)
779
0
    return res;
780
3.85k
781
3.85k
  
if (3.85k
s_qmul(c, (mp_size) p2)3.85k
)
782
3.85k
    return MP_OK;
783
3.85k
  else
784
0
    return MP_MEMORY;
785
3.85k
}
786
787
mp_result mp_int_sqr(mp_int a, mp_int c)
788
0
{
789
0
  mp_digit *out;
790
0
  mp_size   osize, p = 0;
791
0
792
0
  CHECK(a != NULL && c != NULL);
793
0
794
0
  /* Get a temporary buffer big enough to hold the result */
795
0
  osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2);
796
0
  if (
a == c0
)
{0
797
0
    p = ROUND_PREC(osize);
798
0
    p = MAX(p, default_precision);
799
0
800
0
    if ((out = s_alloc(p)) == NULL)
801
0
      return MP_MEMORY;
802
0
  }
803
0
  else {
804
0
    if (!s_pad(c, osize))
805
0
      return MP_MEMORY;
806
0
807
0
    
out = 0
MP_DIGITS0
(c);
808
0
  }
809
0
  
ZERO0
(out, osize);0
810
0
811
0
  s_ksqr(
MP_DIGITS0
(a), out,
MP_USED0
(a));
812
0
813
0
  /* Get rid of whatever memory c was already using, and fix up its fields to
814
0
     reflect the new digit array it's using
815
0
   */
816
0
  if (
out != 0
MP_DIGITS0
(c))
{0
817
0
    if (
(void *) 0
MP_DIGITS0
(c) != (void *) c)
818
0
      
s_free(0
MP_DIGITS0
(c));
819
0
    MP_DIGITS(c) = out;
820
0
    MP_ALLOC(c) = p;
821
0
  }
822
0
823
0
  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
824
0
  CLAMP(c);           /* ... right here */
825
0
  MP_SIGN(c) = MP_ZPOS;
826
0
827
0
  return MP_OK;
828
0
}
829
830
mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r)
831
3.87M
{
832
3.87M
  int cmp, lg;
833
3.87M
  mp_result res = MP_OK;
834
3.87M
  mp_int qout, rout;
835
3.87M
  mp_sign sa = 
MP_SIGN3.87M
(a), sb =
MP_SIGN3.87M
(b);
836
3.87M
  DECLARE_TEMP(2);
837
3.87M
838
3.87M
  CHECK(a != NULL && b != NULL && q != r);
839
3.87M
840
3.87M
  if (
CMPZ3.87M
(b) == 03.87M
)
841
0
    return MP_UNDEF;
842
3.87M
  else 
if (3.87M
(cmp = s_ucmp(a, b)) < 03.87M
)
{724k
843
724k
    /* If |a| < |b|, no division is required:
844
724k
       q = 0, r = a
845
724k
     */
846
724k
    if (
r && 724k
(res = mp_int_copy(a, r)) != MP_OK46.3k
)
847
0
      return res;
848
724k
849
724k
    
if (724k
q724k
)
850
715k
      mp_int_zero(q);
851
724k
852
724k
    return MP_OK;
853
724k
  }
854
3.14M
  else 
if (3.14M
cmp == 03.14M
)
{705k
855
705k
    /* If |a| = |b|, no division is required:
856
705k
       q = 1 or -1, r = 0
857
705k
     */
858
705k
    if (r)
859
4.45k
      mp_int_zero(r);
860
705k
861
705k
    if (
q705k
)
{702k
862
702k
      mp_int_zero(q);
863
702k
      q->digits[0] = 1;
864
702k
865
702k
      if (sa != sb)
866
81.0k
  
MP_SIGN81.0k
(q) = MP_NEG81.0k
;
867
702k
    }
868
705k
869
705k
    return MP_OK;
870
705k
  }
871
3.87M
872
3.87M
  /* When |a| > |b|, real division is required.  We need someplace to store
873
3.87M
     quotient and remainder, but q and r are allowed to be NULL or to overlap
874
3.87M
     with the inputs.
875
3.87M
   */
876
2.44M
  
if (2.44M
(lg = s_isp2(b)) < 02.44M
)
{1.64M
877
1.64M
    if (
q && 1.64M
b != q1.62M
)
{1.52M
878
1.52M
      if ((res = mp_int_copy(a, q)) != MP_OK)
879
0
  goto CLEANUP;
880
1.52M
      else
881
1.52M
  qout = q;
882
1.52M
    }
883
116k
    else {
884
116k
      qout = LAST_TEMP();
885
116k
      SETUP(mp_int_init_copy(LAST_TEMP(), a));
886
116k
    }
887
1.64M
888
1.64M
    
if (1.64M
r && 1.64M
a != r118k
)
{118k
889
118k
      if ((res = mp_int_copy(b, r)) != MP_OK)
890
0
  goto CLEANUP;
891
118k
      else
892
118k
  rout = r;
893
118k
    }
894
1.52M
    else {
895
1.52M
      rout = LAST_TEMP();
896
1.52M
      SETUP(mp_int_init_copy(LAST_TEMP(), b));
897
1.52M
    }
898
1.64M
899
1.64M
    
if (1.64M
(res = s_udiv_knuth(qout, rout)) != MP_OK1.64M
)
goto CLEANUP0
;
900
1.64M
  }
901
799k
  else {
902
799k
    if (
q && 799k
(res = mp_int_copy(a, q)) != MP_OK756k
)
goto CLEANUP0
;
903
799k
    
if (799k
r && 799k
(res = mp_int_copy(a, r)) != MP_OK122k
)
goto CLEANUP0
;
904
799k
905
799k
    
if (799k
q799k
)
s_qdiv(q, (mp_size) lg)756k
; qout = q;
906
799k
    if (
r799k
)
s_qmod(r, (mp_size) lg)122k
; rout = r;
907
799k
  }
908
2.44M
909
2.44M
  /* Recompute signs for output */
910
2.44M
  
if (2.44M
rout2.44M
)
{1.76M
911
1.76M
    MP_SIGN(rout) = sa;
912
1.76M
    if (
CMPZ1.76M
(rout) == 01.76M
)
913
1.68M
      
MP_SIGN1.68M
(rout) = MP_ZPOS1.68M
;
914
1.76M
  }
915
2.44M
  if (
qout2.44M
)
{2.39M
916
2.39M
    
MP_SIGN2.39M
(qout) = (sa == sb) ?
MP_ZPOS1.68M
:
MP_NEG717k
;
917
2.39M
    if (
CMPZ2.39M
(qout) == 02.39M
)
918
0
      
MP_SIGN0
(qout) = MP_ZPOS0
;
919
2.39M
  }
920
2.44M
921
2.44M
  if (
q && 2.44M
(res = mp_int_copy(qout, q)) != MP_OK2.38M
)
goto CLEANUP0
;
922
2.44M
  
if (2.44M
r && 2.44M
(res = mp_int_copy(rout, r)) != MP_OK241k
)
goto CLEANUP0
;
923
2.44M
924
2.44M
  
CLEANUP_TEMP2.44M
();2.44M
925
2.44M
  return res;
926
2.44M
}
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
59.3k
{
959
59.3k
  mpz_t     vtmp, rtmp;
960
59.3k
  mp_digit  vbuf[MP_VALUE_DIGITS(value)];
961
59.3k
  mp_result res;
962
59.3k
963
59.3k
  mp_int_init(&rtmp);
964
59.3k
  s_fake(&vtmp, value, vbuf);
965
59.3k
966
59.3k
  if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
967
0
    goto CLEANUP;
968
59.3k
969
59.3k
  
if (59.3k
r59.3k
)
970
49.4k
    (void) mp_int_to_int(&rtmp, r); /* can't fail */
971
59.3k
972
59.3k
 CLEANUP:
973
59.3k
  mp_int_clear(&rtmp);
974
59.3k
  return res;
975
59.3k
}
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
246k
{
1094
246k
  mp_sign sa;
1095
246k
1096
246k
  CHECK(a != NULL && b != NULL);
1097
246k
1098
246k
  sa = MP_SIGN(a);
1099
246k
  if (
sa == 246k
MP_SIGN246k
(b))
{214k
1100
214k
    int cmp = s_ucmp(a, b);
1101
214k
1102
214k
    /* If they're both zero or positive, the normal comparison applies; if both
1103
214k
       negative, the sense is reversed. */
1104
214k
    if (sa == MP_ZPOS)
1105
150k
      return cmp;
1106
214k
    else
1107
63.4k
      return -cmp;
1108
214k
1109
214k
  }
1110
31.9k
  else {
1111
31.9k
    if (sa == MP_ZPOS)
1112
14.5k
      return 1;
1113
31.9k
    else
1114
17.3k
      return -1;
1115
31.9k
  }
1116
246k
}
1117
1118
int       mp_int_compare_unsigned(mp_int a, mp_int b)
1119
2.76M
{
1120
2.76M
  NRCHECK(a != NULL && b != NULL);
1121
2.76M
1122
2.76M
  return s_ucmp(a, b);
1123
2.76M
}
1124
1125
int       mp_int_compare_zero(mp_int z)
1126
45.9M
{
1127
45.9M
  NRCHECK(z != NULL);
1128
45.9M
1129
45.9M
  if (
MP_USED45.9M
(z) == 1 && 45.9M
z->digits[0] == 019.5M
)
1130
9.56M
    return 0;
1131
36.4M
  else 
if (36.4M
MP_SIGN36.4M
(z) == MP_ZPOS36.4M
)
1132
26.3M
    return 1;
1133
36.4M
  else
1134
10.0M
    return -1;
1135
45.9M
}
1136
1137
int       mp_int_compare_value(mp_int z, mp_small value)
1138
54.5M
{
1139
28.5M
  mp_sign vsign = (value < 0) ? 
MP_NEG25.9M
:
MP_ZPOS28.5M
;
1140
54.5M
  int cmp;
1141
54.5M
1142
54.5M
  CHECK(z != NULL);
1143
54.5M
1144
54.5M
  if (
vsign == 54.5M
MP_SIGN54.5M
(z))
{34.9M
1145
34.9M
    cmp = s_vcmp(z, value);
1146
34.9M
1147
28.5M
    return (vsign == MP_ZPOS) ? 
cmp28.5M
:
-cmp6.38M
;
1148
34.9M
  }
1149
19.6M
  else {
1150
19.5M
    return (value < 0) ? 
119.5M
:
-131.1k
;
1151
19.6M
  }
1152
54.5M
}
1153
1154
int       mp_int_compare_uvalue(mp_int z, mp_usmall uv)
1155
0
{
1156
0
  CHECK(z != NULL);
1157
0
1158
0
  if (
MP_SIGN0
(z) == MP_NEG0
)
1159
0
    return -1;
1160
0
  else
1161
0
    return s_uvcmp(z, uv);
1162
0
}
1163
1164
mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c)
1165
0
{
1166
0
  mp_result res;
1167
0
  mp_size um;
1168
0
  mp_int s;
1169
0
  DECLARE_TEMP(3);
1170
0
1171
0
  CHECK(a != NULL && b != NULL && c != NULL && m != NULL);
1172
0
1173
0
  /* Zero moduli and negative exponents are not considered. */
1174
0
  if (
CMPZ0
(m) == 00
)
1175
0
    return MP_UNDEF;
1176
0
  
if (0
CMPZ0
(b) < 00
)
1177
0
    return MP_RANGE;
1178
0
1179
0
  
um = 0
MP_USED0
(m);
1180
0
  SETUP(mp_int_init_size(TEMP(0), 2 * um));
1181
0
  
SETUP0
(mp_int_init_size(TEMP(1), 2 * um));0
1182
0
1183
0
  
if (0
c == b || 0
c == m0
)
{0
1184
0
    SETUP(mp_int_init_size(TEMP(2), 2 * um));
1185
0
    
s = 0
TEMP0
(2);
1186
0
  }
1187
0
  else {
1188
0
    s = c;
1189
0
  }
1190
0
1191
0
  
if (0
(res = mp_int_mod(a, m, 0
TEMP0
(0))) != MP_OK)
goto CLEANUP0
;
1192
0
1193
0
  
if (0
(res = s_brmu(0
TEMP0
(1), m)) != MP_OK)
goto CLEANUP0
;
1194
0
1195
0
  
if (0
(res = s_embar(0
TEMP0
(0), b, m,
TEMP0
(1), s)) != MP_OK)
1196
0
    goto CLEANUP;
1197
0
1198
0
  res = mp_int_copy(s, c);
1199
0
1200
0
  
CLEANUP_TEMP0
();0
1201
0
  return res;
1202
0
}
1203
1204
mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c)
1205
0
{
1206
0
  mpz_t vtmp;
1207
0
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1208
0
1209
0
  s_fake(&vtmp, value, vbuf);
1210
0
1211
0
  return mp_int_exptmod(a, &vtmp, m, c);
1212
0
}
1213
1214
mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b,
1215
        mp_int m, mp_int c)
1216
0
{
1217
0
  mpz_t vtmp;
1218
0
  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1219
0
1220
0
  s_fake(&vtmp, value, vbuf);
1221
0
1222
0
  return mp_int_exptmod(&vtmp, b, m, c);
1223
0
}
1224
1225
mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
1226
0
{
1227
0
  mp_result res;
1228
0
  mp_size um;
1229
0
  mp_int s;
1230
0
  DECLARE_TEMP(2);
1231
0
1232
0
  CHECK(a && b && m && c);
1233
0
1234
0
  /* Zero moduli and negative exponents are not considered. */
1235
0
  if (
CMPZ0
(m) == 00
)
1236
0
    return MP_UNDEF;
1237
0
  
if (0
CMPZ0
(b) < 00
)
1238
0
    return MP_RANGE;
1239
0
1240
0
  
um = 0
MP_USED0
(m);
1241
0
  SETUP(mp_int_init_size(TEMP(0), 2 * um));
1242
0
1243
0
  
if (0
c == b || 0
c == m0
)
{0
1244
0
    SETUP(mp_int_init_size(TEMP(1), 2 * um));
1245
0
    
s = 0
TEMP0
(1);
1246
0
  }
1247
0
  else {
1248
0
    s = c;
1249
0
  }
1250
0
1251
0
  
if (0
(res = mp_int_mod(a, m, 0
TEMP0
(0))) != MP_OK)
goto CLEANUP0
;
1252
0
1253
0
  
if (0
(res = s_embar(0
TEMP0
(0), b, m, mu, s)) != MP_OK)
1254
0
    goto CLEANUP;
1255
0
1256
0
  res = mp_int_copy(s, c);
1257
0
1258
0
  
CLEANUP_TEMP0
();0
1259
0
  return res;
1260
0
}
1261
1262
mp_result mp_int_redux_const(mp_int m, mp_int c)
1263
0
{
1264
0
  CHECK(m != NULL && c != NULL && m != c);
1265
0
1266
0
  return s_brmu(c, m);
1267
0
}
1268
1269
mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c)
1270
0
{
1271
0
  mp_result res;
1272
0
  mp_sign sa;
1273
0
  DECLARE_TEMP(2);
1274
0
1275
0
  CHECK(a != NULL && m != NULL && c != NULL);
1276
0
1277
0
  if (
CMPZ0
(a) == 0 || 0
CMPZ0
(m) <= 00
)
1278
0
    return MP_RANGE;
1279
0
1280
0
  
sa = 0
MP_SIGN0
(a); /* need this for the result later */
1281
0
1282
0
  for (last__ = 0; 
last__ < 20
;
++last__0
)
1283
0
    
mp_int_init(0
LAST_TEMP0
());
1284
0
1285
0
  if (
(res = mp_int_egcd(a, m, 0
TEMP0
(0),
TEMP0
(1), NULL)) != MP_OK)
1286
0
    goto CLEANUP;
1287
0
1288
0
  
if (0
mp_int_compare_value(0
TEMP0
(0), 1) != 0)
{0
1289
0
    res = MP_UNDEF;
1290
0
    goto CLEANUP;
1291
0
  }
1292
0
1293
0
  /* It is first necessary to constrain the value to the proper range */
1294
0
  
if (0
(res = mp_int_mod(0
TEMP0
(1), m,
TEMP0
(1))) != MP_OK)
1295
0
    goto CLEANUP;
1296
0
1297
0
  /* Now, if 'a' was originally negative, the value we have is actually the
1298
0
     magnitude of the negative representative; to get the positive value we
1299
0
     have to subtract from the modulus.  Otherwise, the value is okay as it
1300
0
     stands.
1301
0
   */
1302
0
  
if (0
sa == MP_NEG0
)
1303
0
    
res = mp_int_sub(m, 0
TEMP0
(1), c);
1304
0
  else
1305
0
    
res = mp_int_copy(0
TEMP0
(1), c);
1306
0
1307
0
  
CLEANUP_TEMP0
();0
1308
0
  return res;
1309
0
}
1310
1311
/* Binary GCD algorithm due to Josef Stein, 1961 */
1312
mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c)
1313
2.50M
{
1314
2.50M
  int ca, cb, k = 0;
1315
2.50M
  mpz_t u, v, t;
1316
2.50M
  mp_result res;
1317
2.50M
1318
2.50M
  CHECK(a != NULL && b != NULL && c != NULL);
1319
2.50M
1320
2.50M
  ca = CMPZ(a);
1321
2.50M
  cb = CMPZ(b);
1322
2.50M
  if (
ca == 0 && 2.50M
cb == 021
)
1323
0
    return MP_UNDEF;
1324
2.50M
  else 
if (2.50M
ca == 02.50M
)
1325
21
    return mp_int_abs(b, c);
1326
2.50M
  else 
if (2.50M
cb == 02.50M
)
1327
9.43k
    return mp_int_abs(a, c);
1328
2.50M
1329
2.49M
  mp_int_init(&t);
1330
2.49M
  if ((res = mp_int_init_copy(&u, a)) != MP_OK)
1331
0
    goto U;
1332
2.49M
  
if (2.49M
(res = mp_int_init_copy(&v, b)) != MP_OK2.49M
)
1333
0
    goto V;
1334
2.49M
1335
2.49M
  
MP_SIGN2.49M
(&u) = MP_ZPOS; 2.49M
MP_SIGN2.49M
(&v) = MP_ZPOS;
1336
2.49M
1337
2.49M
  { /* Divide out common factors of 2 from u and v */
1338
2.49M
    int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
1339
2.49M
1340
2.49M
    k = MIN(div2_u, div2_v);
1341
2.49M
    s_qdiv(&u, (mp_size) k);
1342
2.49M
    s_qdiv(&v, (mp_size) k);
1343
2.49M
  }
1344
2.49M
1345
2.49M
  if (
mp_int_is_odd2.49M
(&u))
{1.98M
1346
1.98M
    if ((res = mp_int_neg(&v, &t)) != MP_OK)
1347
0
      goto CLEANUP;
1348
1.98M
  }
1349
504k
  else {
1350
504k
    if ((res = mp_int_copy(&u, &t)) != MP_OK)
1351
0
      goto CLEANUP;
1352
504k
  }
1353
2.49M
1354
73.9M
  
for (;;) 2.49M
{73.9M
1355
73.9M
    s_qdiv(&t, s_dp2k(&t));
1356
73.9M
1357
73.9M
    if (
CMPZ73.9M
(&t) > 073.9M
)
{10.3M
1358
10.3M
      if ((res = mp_int_copy(&t, &u)) != MP_OK)
1359
0
  goto CLEANUP;
1360
10.3M
    }
1361
63.6M
    else {
1362
63.6M
      if ((res = mp_int_neg(&t, &v)) != MP_OK)
1363
0
  goto CLEANUP;
1364
63.6M
    }
1365
73.9M
1366
73.9M
    
if (73.9M
(res = mp_int_sub(&u, &v, &t)) != MP_OK73.9M
)
1367
0
      goto CLEANUP;
1368
73.9M
1369
73.9M
    
if (73.9M
CMPZ73.9M
(&t) == 073.9M
)
1370
2.49M
      break;
1371
73.9M
  }
1372
2.49M
1373
2.49M
  
if (2.49M
(res = mp_int_abs(&u, c)) != MP_OK2.49M
)
1374
0
    goto CLEANUP;
1375
2.49M
  
if (2.49M
!s_qmul(c, (mp_size) k)2.49M
)
1376
0
    res = MP_MEMORY;
1377
2.49M
1378
2.49M
 CLEANUP:
1379
2.49M
  mp_int_clear(&v);
1380
2.49M
 V: mp_int_clear(&u);
1381
2.49M
 U: mp_int_clear(&t);
1382
2.49M
1383
2.49M
  return res;
1384
2.49M
}
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
267k
{
1500
267k
  mpz_t lcm;
1501
267k
  mp_result res;
1502
267k
1503
267k
  CHECK(a != NULL && b != NULL && c != NULL);
1504
267k
1505
267k
  /* Since a * b = gcd(a, b) * lcm(a, b), we can compute
1506
267k
     lcm(a, b) = (a / gcd(a, b)) * b.
1507
267k
1508
267k
     This formulation insures everything works even if the input
1509
267k
     variables share space.
1510
267k
   */
1511
267k
  if ((res = mp_int_init(&lcm)) != MP_OK)
1512
0
    return res;
1513
267k
  
if (267k
(res = mp_int_gcd(a, b, &lcm)) != MP_OK267k
)
1514
0
    goto CLEANUP;
1515
267k
  
if (267k
(res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK267k
)
1516
0
    goto CLEANUP;
1517
267k
  
if (267k
(res = mp_int_mul(&lcm, b, &lcm)) != MP_OK267k
)
1518
0
    goto CLEANUP;
1519
267k
1520
267k
  res = mp_int_copy(&lcm, c);
1521
267k
1522
267k
  CLEANUP:
1523
267k
    mp_int_clear(&lcm);
1524
267k
1525
267k
  return res;
1526
267k
}
1527
1528
int       mp_int_divisible_value(mp_int a, mp_small v)
1529
49.4k
{
1530
49.4k
  mp_small rem = 0;
1531
49.4k
1532
49.4k
  if (mp_int_div_value(a, v, NULL, &rem) != MP_OK)
1533
0
    return 0;
1534
49.4k
1535
49.4k
  return rem == 0;
1536
49.4k
}
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
31.8M
{
1615
31.8M
  mp_usmall uv = 0;
1616
31.8M
  mp_size   uz;
1617
31.8M
  mp_digit *dz;
1618
31.8M
  mp_sign   sz;
1619
31.8M
1620
31.8M
  CHECK(z != NULL);
1621
31.8M
1622
31.8M
  /* Make sure the value is representable as a small integer */
1623
31.8M
  sz = MP_SIGN(z);
1624
31.8M
  if (
(sz == MP_ZPOS && 31.8M
mp_int_compare_value(z, 25.4M
MP_SMALL_MAX25.4M
) > 0) ||
1625
25.9M
      
mp_int_compare_value(z, 25.9M
MP_SMALL_MIN25.9M
) < 0)
1626
8.28M
    return MP_RANGE;
1627
31.8M
1628
23.5M
  
uz = 23.5M
MP_USED23.5M
(z);
1629
23.5M
  dz = MP_DIGITS(z) + uz - 1;
1630
23.5M
1631
54.6M
  while (
uz > 054.6M
)
{31.1M
1632
31.1M
    uv <<= MP_DIGIT_BIT/2;
1633
31.1M
    uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1634
31.1M
    --uz;
1635
31.1M
  }
1636
23.5M
1637
23.5M
  if (out)
1638
23.5M
    
*out = (sz == MP_NEG) ? 23.5M
-(mp_small)uv3.95M
:
(mp_small)uv19.5M
;
1639
23.5M
1640
23.5M
  return MP_OK;
1641
31.8M
}
1642
1643
mp_result mp_int_to_uint(mp_int z, mp_usmall *out)
1644
0
{
1645
0
  mp_usmall uv = 0;
1646
0
  mp_size   uz;
1647
0
  mp_digit *dz;
1648
0
  mp_sign   sz;
1649
0
  
1650
0
  CHECK(z != NULL);
1651
0
1652
0
  /* Make sure the value is representable as an unsigned small integer */
1653
0
  sz = MP_SIGN(z);
1654
0
  if (
sz == MP_NEG || 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.35k
{
1675
1.35k
  mp_result res;
1676
1.35k
  int       cmp = 0;
1677
1.35k
1678
1.35k
  CHECK(z != NULL && str != NULL && limit >= 2);
1679
1.35k
1680
1.35k
  if (
radix < 1.35k
MP_MIN_RADIX1.35k
||
radix > 1.35k
MP_MAX_RADIX1.35k
)
1681
0
    return MP_RANGE;
1682
1.35k
1683
1.35k
  
if (1.35k
CMPZ1.35k
(z) == 01.35k
)
{0
1684
0
    *str++ = s_val2ch(0, 1);
1685
0
  }
1686
1.35k
  else {
1687
1.35k
    mpz_t tmp;
1688
1.35k
    char  *h, *t;
1689
1.35k
1690
1.35k
    if ((res = mp_int_init_copy(&tmp, z)) != MP_OK)
1691
0
      return res;
1692
1.35k
1693
1.35k
    
if (1.35k
MP_SIGN1.35k
(z) == MP_NEG1.35k
)
{677
1694
677
      *str++ = '-';
1695
677
      --limit;
1696
677
    }
1697
1.35k
    h = str;
1698
1.35k
1699
1.35k
    /* Generate digits in reverse order until finished or limit reached */
1700
26.3k
    for (/* */; 
limit > 026.3k
;
--limit24.9k
)
{26.3k
1701
26.3k
      mp_digit d;
1702
26.3k
1703
26.3k
      if (
(cmp = 26.3k
CMPZ26.3k
(&tmp)) == 0)
1704
1.35k
  break;
1705
26.3k
1706
24.9k
      d = s_ddiv(&tmp, (mp_digit)radix);
1707
24.9k
      *str++ = s_val2ch(d, 1);
1708
24.9k
    }
1709
1.35k
    t = str - 1;
1710
1.35k
1711
1.35k
    /* Put digits back in correct output order */
1712
13.2k
    while (
h < t13.2k
)
{11.9k
1713
11.9k
      char tc = *h;
1714
11.9k
      *h++ = *t;
1715
11.9k
      *t-- = tc;
1716
11.9k
    }
1717
1.35k
1718
1.35k
    mp_int_clear(&tmp);
1719
1.35k
  }
1720
1.35k
1721
1.35k
  *str = '\0';
1722
1.35k
  if (cmp == 0)
1723
1.35k
    return MP_OK;
1724
1.35k
  else
1725
0
    return MP_TRUNC;
1726
1.35k
}
1727
1728
mp_result mp_int_string_len(mp_int z, mp_size radix)
1729
1.43k
{
1730
1.43k
  int  len;
1731
1.43k
1732
1.43k
  CHECK(z != NULL);
1733
1.43k
1734
1.43k
  if (
radix < 1.43k
MP_MIN_RADIX1.43k
||
radix > 1.43k
MP_MAX_RADIX1.43k
)
1735
0
    return MP_RANGE;
1736
1.43k
1737
1.43k
  len = s_outlen(z, radix) + 1; /* for terminator */
1738
1.43k
1739
1.43k
  /* Allow for sign marker on negatives */
1740
1.43k
  if (
MP_SIGN1.43k
(z) == MP_NEG1.43k
)
1741
707
    len += 1;
1742
1.43k
1743
1.43k
  return len;
1744
1.43k
}
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
78
{
1749
78
  return mp_int_read_cstring(z, radix, str, NULL);
1750
78
}
1751
1752
mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
1753
78
{
1754
78
  int ch;
1755
78
1756
78
  CHECK(z != NULL && str != NULL);
1757
78
1758
78
  if (
radix < 78
MP_MIN_RADIX78
||
radix > 78
MP_MAX_RADIX78
)
1759
0
    return MP_RANGE;
1760
78
1761
78
  /* Skip leading whitespace */
1762
78
  
while (78
isspace((int)*str)78
)
1763
0
    ++str;
1764
78
1765
78
  /* Handle leading sign tag (+/-, positive default) */
1766
78
  switch (*str) {
1767
19
  case '-':
1768
19
    MP_SIGN(z) = MP_NEG;
1769
19
    ++str;
1770
19
    break;
1771
0
  case '+':
1772
0
    ++str; /* fallthrough */
1773
59
  default:
1774
59
    MP_SIGN(z) = MP_ZPOS;
1775
59
    break;
1776
78
  }
1777
78
1778
78
  /* Skip leading zeroes */
1779
78
  
while (78
(ch = s_ch2val(*str, radix)) == 078
)
1780
0
    ++str;
1781
78
1782
78
  /* Make sure there is enough space for the value */
1783
78
  if (!s_pad(z, s_inlen(strlen(str), radix)))
1784
0
    return MP_MEMORY;
1785
78
1786
78
  
MP_USED78
(z) = 1; z->digits[0] = 0;78
1787
78
1788
1.03k
  while (
*str != '\0' && 1.03k
((ch = s_ch2val(*str, radix)) >= 0)960
)
{960
1789
960
    s_dmul(z, (mp_digit)radix);
1790
960
    s_dadd(z, (mp_digit)ch);
1791
960
    ++str;
1792
960
  }
1793
78
1794
78
  CLAMP(z);
1795
78
1796
78
  /* Override sign for zero, even if negative specified. */
1797
78
  if (
CMPZ78
(z) == 078
)
1798
0
    
MP_SIGN0
(z) = MP_ZPOS0
;
1799
78
1800
78
  if (end != NULL)
1801
0
    *end = (char *)str;
1802
78
1803
78
  /* Return a truncation error if the string has unprocessed characters
1804
78
     remaining, so the caller can tell if the whole string was done */
1805
78
  if (*str != '\0')
1806
0
    return MP_TRUNC;
1807
78
  else
1808
78
    return MP_OK;
1809
78
}
1810
1811
mp_result mp_int_count_bits(mp_int z)
1812
2.73k
{
1813
2.73k
  mp_size  nbits = 0, uz;
1814
2.73k
  mp_digit d;
1815
2.73k
1816
2.73k
  CHECK(z != NULL);
1817
2.73k
1818
2.73k
  uz = MP_USED(z);
1819
2.73k
  if (
uz == 1 && 2.73k
z->digits[0] == 01.45k
)
1820
0
    return 1;
1821
2.73k
1822
2.73k
  --uz;
1823
2.73k
  nbits = uz * MP_DIGIT_BIT;
1824
2.73k
  d = z->digits[uz];
1825
2.73k
1826
53.4k
  while (
d != 053.4k
)
{50.7k
1827
50.7k
    d >>= 1;
1828
50.7k
    ++nbits;
1829
50.7k
  }
1830
2.73k
1831
2.73k
  return nbits;
1832
2.73k
}
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.29k
{
1938
1.29k
  mp_result  res = mp_int_count_bits(z);
1939
1.29k
  int        bytes;
1940
1.29k
1941
1.29k
  if (res <= 0)
1942
0
    return res;
1943
1.29k
1944
1.29k
  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
1945
1.29k
1946
1.29k
  return bytes;
1947
1.29k
}
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
19.6M
{
1970
19.6M
  mp_digit *out = malloc(num * sizeof(mp_digit));
1971
19.6M
1972
19.6M
  assert(out != NULL); /* for debugging */
1973
19.6M
#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
19.6M
1983
19.6M
  return out;
1984
19.6M
}
1985
1986
STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
1987
642k
{
1988
642k
#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
642k
  mp_digit *new = realloc(old, nsize * sizeof(mp_digit));
1998
642k
1999
642k
  assert(new != NULL); /* for debugging */
2000
642k
#endif
2001
642k
  return new;
2002
642k
}
2003
2004
STATIC void s_free(void *ptr)
2005
19.6M
{
2006
19.6M
  free(ptr);
2007
19.6M
}
2008
2009
STATIC int      s_pad(mp_int z, mp_size min)
2010
181M
{
2011
181M
  if (
MP_ALLOC181M
(z) < min181M
)
{13.2M
2012
13.2M
    mp_size nsize = ROUND_PREC(min);
2013
13.2M
    mp_digit *tmp;
2014
13.2M
2015
13.2M
    if (
(void *)z->digits == (void *)z13.2M
)
{12.5M
2016
12.5M
      if ((tmp = s_alloc(nsize)) == NULL)
2017
0
        return 0;
2018
12.5M
2019
12.5M
      
COPY12.5M
(MP_DIGITS(z), tmp, MP_USED(z));12.5M
2020
12.5M
    }
2021
642k
    else 
if (642k
(tmp = s_realloc(642k
MP_DIGITS642k
(z),
MP_ALLOC642k
(z), nsize)) == NULL)
2022
0
      return 0;
2023
13.2M
2024
13.2M
    
MP_DIGITS13.2M
(z) = tmp;13.2M
2025
13.2M
    MP_ALLOC(z) = nsize;
2026
13.2M
  }
2027
181M
2028
181M
  return 1;
2029
181M
}
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
915k
{
2034
910k
  mp_usmall uv = (mp_usmall) (value < 0) ? 
-value4.90k
:
value910k
;
2035
915k
  s_ufake(z, uv, vbuf);
2036
915k
  if (value < 0)
2037
4.90k
    z->sign = MP_NEG;
2038
915k
}
2039
2040
STATIC void      s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[])
2041
36.0M
{
2042
36.0M
  mp_size ndig = (mp_size) s_uvpack(value, vbuf);
2043
36.0M
2044
36.0M
  z->used = ndig;
2045
36.0M
  z->alloc = MP_VALUE_DIGITS(value);
2046
36.0M
  z->sign = MP_ZPOS;
2047
36.0M
  z->digits = vbuf;
2048
36.0M
}
2049
2050
STATIC int      s_cdig(mp_digit *da, mp_digit *db, mp_size len)
2051
50.4M
{
2052
50.4M
  mp_digit *dat = da + len - 1, *dbt = db + len - 1;
2053
50.4M
2054
69.2M
  for (/* */; 
len != 069.2M
;
--len, --dat, --dbt18.7M
)
{61.5M
2055
61.5M
    if (*dat > *dbt)
2056
11.2M
      return 1;
2057
50.3M
    else 
if (50.3M
*dat < *dbt50.3M
)
2058
31.5M
      return -1;
2059
61.5M
  }
2060
50.4M
2061
7.66M
  return 0;
2062
50.4M
}
2063
2064
STATIC int       s_uvpack(mp_usmall uv, mp_digit t[])
2065
36.0M
{
2066
36.0M
  int ndig = 0;
2067
36.0M
2068
36.0M
  if (uv == 0)
2069
787k
    t[ndig++] = 0;
2070
35.2M
  else {
2071
102M
    while (
uv != 0102M
)
{67.0M
2072
67.0M
      t[ndig++] = (mp_digit) uv;
2073
67.0M
      uv >>= MP_DIGIT_BIT/2;
2074
67.0M
      uv >>= MP_DIGIT_BIT/2;
2075
67.0M
    }
2076
35.2M
  }
2077
36.0M
2078
36.0M
  return ndig;
2079
36.0M
}
2080
2081
STATIC int      s_ucmp(mp_int a, mp_int b)
2082
122M
{
2083
122M
  mp_size  ua = 
MP_USED122M
(a), ub =
MP_USED122M
(b);
2084
122M
2085
122M
  if (ua > ub)
2086
16.6M
    return 1;
2087
105M
  else 
if (105M
ub > ua105M
)
2088
54.9M
    return -1;
2089
105M
  else
2090
50.4M
    
return s_cdig(50.4M
MP_DIGITS50.4M
(a),
MP_DIGITS50.4M
(b), ua);
2091
122M
}
2092
2093
STATIC int      s_vcmp(mp_int a, mp_small v)
2094
34.9M
{
2095
28.5M
  mp_usmall uv = (v < 0) ? 
-(mp_usmall) v6.38M
:
(mp_usmall) v28.5M
;
2096
34.9M
  return s_uvcmp(a, uv);
2097
34.9M
}
2098
2099
STATIC int      s_uvcmp(mp_int a, mp_usmall uv)
2100
34.9M
{
2101
34.9M
  mpz_t vtmp;
2102
34.9M
  mp_digit vdig[MP_VALUE_DIGITS(uv)];
2103
34.9M
2104
34.9M
  s_ufake(&vtmp, uv, vdig);
2105
34.9M
  return s_ucmp(a, &vtmp);
2106
34.9M
}
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
5.72M
{
2111
5.72M
  mp_size pos;
2112
5.72M
  mp_word w = 0;
2113
5.72M
2114
5.72M
  /* Insure that da is the longer of the two to simplify later code */
2115
5.72M
  if (
size_b > size_a5.72M
)
{743k
2116
743k
    SWAP(mp_digit *, da, db);
2117
743k
    SWAP(mp_size, size_a, size_b);
2118
743k
  }
2119
5.72M
2120
5.72M
  /* Add corresponding digits until the shorter number runs out */
2121
12.7M
  for (pos = 0; 
pos < size_b12.7M
;
++pos, ++da, ++db, ++dc7.05M
)
{7.05M
2122
7.05M
    w = w + (mp_word) *da + (mp_word) *db;
2123
7.05M
    *dc = LOWER_HALF(w);
2124
7.05M
    w = UPPER_HALF(w);
2125
7.05M
  }
2126
5.72M
2127
5.72M
  /* Propagate carries as far as necessary */
2128
11.0M
  for (/* */; 
pos < size_a11.0M
;
++pos, ++da, ++dc5.33M
)
{5.33M
2129
5.33M
    w = w + *da;
2130
5.33M
2131
5.33M
    *dc = LOWER_HALF(w);
2132
5.33M
    w = UPPER_HALF(w);
2133
5.33M
  }
2134
5.72M
2135
5.72M
  /* Return carry out */
2136
5.72M
  return (mp_digit)w;
2137
5.72M
}
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
79.3M
{
2142
79.3M
  mp_size pos;
2143
79.3M
  mp_word w = 0;
2144
79.3M
2145
79.3M
  /* We assume that |a| >= |b| so this should definitely hold */
2146
79.3M
  assert(size_a >= size_b);
2147
79.3M
2148
79.3M
  /* Subtract corresponding digits and propagate borrow */
2149
213M
  for (pos = 0; 
pos < size_b213M
;
++pos, ++da, ++db, ++dc134M
)
{134M
2150
134M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2151
134M
   (mp_word)*da) - w - (mp_word)*db;
2152
134M
2153
134M
    *dc = LOWER_HALF(w);
2154
134M
    w = (UPPER_HALF(w) == 0);
2155
134M
  }
2156
79.3M
2157
79.3M
  /* Finish the subtraction for remaining upper digits of da */
2158
164M
  for (/* */; 
pos < size_a164M
;
++pos, ++da, ++dc84.8M
)
{84.8M
2159
84.8M
    w = ((mp_word)MP_DIGIT_MAX + 1 +  /* MP_RADIX */
2160
84.8M
   (mp_word)*da) - w;
2161
84.8M
2162
84.8M
    *dc = LOWER_HALF(w);
2163
84.8M
    w = (UPPER_HALF(w) == 0);
2164
84.8M
  }
2165
79.3M
2166
79.3M
  /* If there is a borrow out at the end, it violates the precondition */
2167
79.3M
  assert(w == 0);
2168
79.3M
}
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
6.78M
{
2173
6.78M
  mp_size  bot_size;
2174
6.78M
2175
6.78M
  /* Make sure b is the smaller of the two input values */
2176
6.78M
  if (
size_b > size_a6.78M
)
{2.94M
2177
2.94M
    SWAP(mp_digit *, da, db);
2178
2.94M
    SWAP(mp_size, size_a, size_b);
2179
2.94M
  }
2180
6.78M
2181
6.78M
  /* Insure that the bottom is the larger half in an odd-length split; the code
2182
6.78M
     below relies on this being true.
2183
6.78M
   */
2184
6.78M
  bot_size = (size_a + 1) / 2;
2185
6.78M
2186
6.78M
  /* If the values are big enough to bother with recursion, use the Karatsuba
2187
6.78M
     algorithm to compute the product; otherwise use the normal multiplication
2188
6.78M
     algorithm
2189
6.78M
   */
2190
6.78M
  if (multiply_threshold &&
2191
6.78M
      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
6.78M
  else {
2249
6.78M
    s_umul(da, db, dc, size_a, size_b);
2250
6.78M
  }
2251
6.78M
2252
6.78M
  return 1;
2253
6.78M
}
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
6.78M
{
2258
6.78M
  mp_size a, b;
2259
6.78M
  mp_word w;
2260
6.78M
2261
24.1M
  for (a = 0; 
a < size_a24.1M
;
++a, ++dc, ++da17.3M
)
{17.3M
2262
17.3M
    mp_digit *dct = dc;
2263
17.3M
    mp_digit *dbt = db;
2264
17.3M
2265
17.3M
    if (*da == 0)
2266
1.26M
      continue;
2267
17.3M
2268
16.1M
    w = 0;
2269
38.6M
    for (b = 0; 
b < size_b38.6M
;
++b, ++dbt, ++dct22.4M
)
{22.4M
2270
22.4M
      w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
2271
22.4M
2272
22.4M
      *dct = LOWER_HALF(w);
2273
22.4M
      w = UPPER_HALF(w);
2274
22.4M
    }
2275
16.1M
2276
16.1M
    *dct = (mp_digit)w;
2277
16.1M
  }
2278
6.78M
}
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
960
{
2387
960
  mp_word w = 0;
2388
960
  mp_digit *da = MP_DIGITS(a);
2389
960
  mp_size ua = MP_USED(a);
2390
960
2391
960
  w = (mp_word)*da + b;
2392
960
  *da++ = LOWER_HALF(w);
2393
960
  w = UPPER_HALF(w);
2394
960
2395
1.16k
  for (ua -= 1; 
ua > 01.16k
;
--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
960
2402
960
  if (
w960
)
{2
2403
2
    *da = (mp_digit)w;
2404
2
    MP_USED(a) += 1;
2405
2
  }
2406
960
}
2407
2408
STATIC void      s_dmul(mp_int a, mp_digit b)
2409
960
{
2410
960
  mp_word w = 0;
2411
960
  mp_digit *da = MP_DIGITS(a);
2412
960
  mp_size ua = MP_USED(a);
2413
960
2414
2.10k
  while (
ua > 02.10k
)
{1.14k
2415
1.14k
    w = (mp_word)*da * b + w;
2416
1.14k
    *da++ = LOWER_HALF(w);
2417
1.14k
    w = UPPER_HALF(w);
2418
1.14k
    --ua;
2419
1.14k
  }
2420
960
2421
960
  if (
w960
)
{20
2422
20
    *da = (mp_digit)w;
2423
20
    MP_USED(a) += 1;
2424
20
  }
2425
960
}
2426
2427
STATIC void      s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
2428
2.08M
{
2429
2.08M
  mp_word  w = 0;
2430
2.08M
2431
10.6M
  while (
size_a > 010.6M
)
{8.53M
2432
8.53M
    w = (mp_word)*da++ * (mp_word)b + w;
2433
8.53M
2434
8.53M
    *dc++ = LOWER_HALF(w);
2435
8.53M
    w = UPPER_HALF(w);
2436
8.53M
    --size_a;
2437
8.53M
  }
2438
2.08M
2439
2.08M
  if (w)
2440
0
    
*dc = 0
LOWER_HALF0
(w);
2441
2.08M
}
2442
2443
STATIC mp_digit  s_ddiv(mp_int a, mp_digit b)
2444
826k
{
2445
826k
  mp_word w = 0, qdigit;
2446
826k
  mp_size ua = MP_USED(a);
2447
826k
  mp_digit *da = MP_DIGITS(a) + ua - 1;
2448
826k
2449
2.84M
  for (/* */; 
ua > 02.84M
;
--ua, --da2.02M
)
{2.02M
2450
2.02M
    w = (w << MP_DIGIT_BIT) | *da;
2451
2.02M
2452
2.02M
    if (
w >= b2.02M
)
{1.54M
2453
1.54M
      qdigit = w / b;
2454
1.54M
      w = w % b;
2455
1.54M
    }
2456
474k
    else {
2457
474k
      qdigit = 0;
2458
474k
    }
2459
2.02M
2460
2.02M
    *da = (mp_digit)qdigit;
2461
2.02M
  }
2462
826k
2463
826k
  CLAMP(a);
2464
826k
  return (mp_digit)w;
2465
826k
}
2466
2467
STATIC void     s_qdiv(mp_int z, mp_size p2)
2468
80.5M
{
2469
80.5M
  mp_size ndig = p2 / 
MP_DIGIT_BIT80.5M
, nbits = p2 %
MP_DIGIT_BIT80.5M
;
2470
80.5M
  mp_size uz = MP_USED(z);
2471
80.5M
2472
80.5M
  if (
ndig80.5M
)
{238k
2473
238k
    mp_size  mark;
2474
238k
    mp_digit *to, *from;
2475
238k
2476
238k
    if (
ndig >= uz238k
)
{0
2477
0
      mp_int_zero(z);
2478
0
      return;
2479
0
    }
2480
238k
2481
238k
    
to = 238k
MP_DIGITS238k
(z); from = to + ndig;
2482
238k
2483
557k
    for (mark = ndig; 
mark < uz557k
;
++mark318k
)
2484
318k
      *to++ = *from++;
2485
238k
2486
238k
    MP_USED(z) = uz - ndig;
2487
238k
  }
2488
80.5M
2489
80.5M
  
if (80.5M
nbits80.5M
)
{76.1M
2490
76.1M
    mp_digit d = 0, *dz, save;
2491
76.1M
    mp_size  up = MP_DIGIT_BIT - nbits;
2492
76.1M
2493
76.1M
    uz = MP_USED(z);
2494
76.1M
    dz = MP_DIGITS(z) + uz - 1;
2495
76.1M
2496
285M
    for (/* */; 
uz > 0285M
;
--uz, --dz208M
)
{208M
2497
208M
      save = *dz;
2498
208M
2499
208M
      *dz = (*dz >> nbits) | (d << up);
2500
208M
      d = save;
2501
208M
    }
2502
76.1M
2503
76.1M
    CLAMP(z);
2504
76.1M
  }
2505
80.5M
2506
80.5M
  if (
MP_USED80.5M
(z) == 1 && 80.5M
z->digits[0] == 028.0M
)
2507
776k
    
MP_SIGN776k
(z) = MP_ZPOS776k
;
2508
80.5M
}
2509
2510
STATIC void     s_qmod(mp_int z, mp_size p2)
2511
122k
{
2512
122k
  mp_size start = p2 / 
MP_DIGIT_BIT122k
+ 1, rest = p2 %
MP_DIGIT_BIT122k
;
2513
122k
  mp_size uz = MP_USED(z);
2514
122k
  mp_digit mask = (1 << rest) - 1;
2515
122k
2516
122k
  if (
start <= uz122k
)
{122k
2517
122k
    MP_USED(z) = start;
2518
122k
    z->digits[start - 1] &= mask;
2519
122k
    CLAMP(z);
2520
122k
  }
2521
122k
}
2522
2523
STATIC int      s_qmul(mp_int z, mp_size p2)
2524
4.14M
{
2525
4.14M
  mp_size   uz, need, rest, extra, i;
2526
4.14M
  mp_digit *from, *to, d;
2527
4.14M
2528
4.14M
  if (p2 == 0)
2529
1.37M
    return 1;
2530
4.14M
2531
2.76M
  
uz = 2.76M
MP_USED2.76M
(z);
2532
2.76M
  need = p2 / 
MP_DIGIT_BIT2.76M
; rest = p2 %
MP_DIGIT_BIT2.76M
;
2533
2.76M
2534
2.76M
  /* Figure out if we need an extra digit at the top end; this occurs if the
2535
2.76M
     topmost `rest' bits of the high-order digit of z are not zero, meaning
2536
2.76M
     they will be shifted off the end if not preserved */
2537
2.76M
  extra = 0;
2538
2.76M
  if (
rest != 02.76M
)
{2.71M
2539
2.71M
    mp_digit *dz = MP_DIGITS(z) + uz - 1;
2540
2.71M
2541
2.71M
    if (
(*dz >> (2.71M
MP_DIGIT_BIT2.71M
- rest)) != 0)
2542
578k
      extra = 1;
2543
2.71M
  }
2544
2.76M
2545
2.76M
  if (!s_pad(z, uz + need + extra))
2546
0
    return 0;
2547
2.76M
2548
2.76M
  /* If we need to shift by whole digits, do that in one pass, then
2549
2.76M
     to back and shift by partial digits.
2550
2.76M
   */
2551
2.76M
  
if (2.76M
need > 02.76M
)
{79.8k
2552
79.8k
    from = MP_DIGITS(z) + uz - 1;
2553
79.8k
    to = from + need;
2554
79.8k
2555
159k
    for (i = 0; 
i < uz159k
;
++i80.0k
)
2556
80.0k
      *to-- = *from--;
2557
79.8k
2558
79.8k
    ZERO(MP_DIGITS(z), need);
2559
79.8k
    uz += need;
2560
79.8k
  }
2561
2.76M
2562
2.76M
  if (
rest2.76M
)
{2.71M
2563
2.71M
    d = 0;
2564
10.0M
    for (i = need, from = 
MP_DIGITS2.71M
(z) + need;
i < uz10.0M
;
++i, ++from7.37M
)
{7.37M
2565
7.37M
      mp_digit save = *from;
2566
7.37M
      
2567
7.37M
      *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
2568
7.37M
      d = save;
2569
7.37M
    }
2570
2.71M
2571
2.71M
    d >>= (MP_DIGIT_BIT - rest);
2572
2.71M
    if (
d != 02.71M
)
{578k
2573
578k
      *from = d;
2574
578k
      uz += extra;
2575
578k
    }
2576
2.71M
  }
2577
2.76M
2578
2.76M
  MP_USED(z) = uz;
2579
2.76M
  CLAMP(z);
2580
2.76M
2581
2.76M
  return 1;
2582
2.76M
}
2583
2584
/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
2585
   The sign of the result is always zero/positive.
2586
 */
2587
STATIC int       s_qsub(mp_int z, mp_size p2)
2588
0
{
2589
0
  mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
2590
0
  mp_size  tdig = (p2 / MP_DIGIT_BIT), pos;
2591
0
  mp_word  w = 0;
2592
0
2593
0
  if (!s_pad(z, tdig + 1))
2594
0
    return 0;
2595
0
2596
0
  
for (pos = 0, zp = 0
MP_DIGITS0
(z);
pos < tdig0
;
++pos, ++zp0
)
{0
2597
0
    w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp;
2598
0
2599
0
    *zp = LOWER_HALF(w);
2600
0
    w = 
UPPER_HALF0
(w) ?
00
:
10
;
2601
0
  }
2602
0
2603
0
  w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp;
2604
0
  *zp = LOWER_HALF(w);
2605
0
2606
0
  assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
2607
0
2608
0
  MP_SIGN(z) = MP_ZPOS;
2609
0
  CLAMP(z);
2610
0
2611
0
  return 1;
2612
0
}
2613
2614
STATIC int      s_dp2k(mp_int z)
2615
78.9M
{
2616
78.9M
  int       k = 0;
2617
78.9M
  mp_digit *dp = MP_DIGITS(z), d;
2618
78.9M
2619
78.9M
  if (
MP_USED78.9M
(z) == 1 && 78.9M
*dp == 025.2M
)
2620
0
    return 1;
2621
78.9M
2622
79.3M
  
while (78.9M
*dp == 079.3M
)
{401k
2623
401k
    k += MP_DIGIT_BIT;
2624
401k
    ++dp;
2625
401k
  }
2626
78.9M
2627
78.9M
  d = *dp;
2628
231M
  while (
(d & 1) == 0231M
)
{152M
2629
152M
    d >>= 1;
2630
152M
    ++k;
2631
152M
  }
2632
78.9M
2633
78.9M
  return k;
2634
78.9M
}
2635
2636
STATIC int       s_isp2(mp_int z)
2637
2.44M
{
2638
2.44M
  mp_size uz = MP_USED(z), k = 0;
2639
2.44M
  mp_digit *dz = MP_DIGITS(z), d;
2640
2.44M
2641
2.49M
  while (
uz > 12.49M
)
{890k
2642
890k
    if (*dz++ != 0)
2643
838k
      return -1;
2644
52.3k
    
k += 52.3k
MP_DIGIT_BIT52.3k
;
2645
52.3k
    --uz;
2646
52.3k
  }
2647
2.44M
2648
1.60M
  d = *dz;
2649
3.80M
  while (
d > 13.80M
)
{3.00M
2650
3.00M
    if (d & 1)
2651
803k
      return -1;
2652
2.20M
    ++k; d >>= 1;
2653
2.20M
  }
2654
1.60M
2655
799k
  return (int) k;
2656
1.60M
}
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
839k
{
2679
839k
  mp_digit d = b->digits[MP_USED(b) - 1];
2680
839k
  int k = 0;
2681
839k
2682
14.9M
  while (
d < (mp_digit) (1 << (14.9M
MP_DIGIT_BIT14.9M
- 1)))
{ /* d < (MP_RADIX / 2) */14.1M
2683
14.1M
    d <<= 1;
2684
14.1M
    ++k;
2685
14.1M
  }
2686
839k
2687
839k
  /* These multiplications can't fail */
2688
839k
  if (
k != 0839k
)
{824k
2689
824k
    (void) s_qmul(a, (mp_size) k);
2690
824k
    (void) s_qmul(b, (mp_size) k);
2691
824k
  }
2692
839k
2693
839k
  return k;
2694
839k
}
2695
2696
STATIC mp_result s_brmu(mp_int z, mp_int m)
2697
0
{
2698
0
  mp_size um = MP_USED(m) * 2;
2699
0
2700
0
  if (!s_pad(z, um))
2701
0
    return MP_MEMORY;
2702
0
2703
0
  
s_2expt(z, 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
1.64M
STATIC mp_result s_udiv_knuth(mp_int u, mp_int v) {
2940
1.64M
  mpz_t q, r, t;
2941
1.64M
  mp_result
2942
1.64M
  res = MP_OK;
2943
1.64M
  int k,j;
2944
1.64M
  mp_size m,n;
2945
1.64M
2946
1.64M
  /* Force signs to positive */
2947
1.64M
  MP_SIGN(u) = MP_ZPOS;
2948
1.64M
  MP_SIGN(v) = MP_ZPOS;
2949
1.64M
2950
1.64M
  /* Use simple division algorithm when v is only one digit long */
2951
1.64M
  if (
MP_USED1.64M
(v) == 11.64M
)
{801k
2952
801k
    mp_digit d, rem;
2953
801k
    d   = v->digits[0];
2954
801k
    rem = s_ddiv(u, d);
2955
801k
    mp_int_set_value(v, rem);
2956
801k
    return MP_OK;
2957
801k
  }
2958
1.64M
2959
1.64M
  /************************************************************/
2960
1.64M
  /* Algorithm D */
2961
1.64M
  /************************************************************/
2962
1.64M
  /* The n and m variables are defined as used by Knuth.
2963
1.64M
     u is an n digit number with digits u_{n-1}..u_0.
2964
1.64M
     v is an n+m digit number with digits from v_{m+n-1}..v_0.
2965
1.64M
     We require that n > 1 and m >= 0 */
2966
839k
  
n = 839k
MP_USED839k
(v);
2967
839k
  m = MP_USED(u) - n;
2968
839k
  assert(n >  1);
2969
839k
  assert(m >= 0);
2970
839k
2971
839k
  /************************************************************/
2972
839k
  /* D1: Normalize.
2973
839k
     The normalization step provides the necessary condition for Theorem B,
2974
839k
     which states that the quotient estimate for q_j, call it qhat
2975
839k
2976
839k
       qhat = u_{j+n}u_{j+n-1} / v_{n-1}
2977
839k
2978
839k
     is bounded by
2979
839k
2980
839k
      qhat - 2 <= q_j <= qhat.
2981
839k
2982
839k
     That is, qhat is always greater than the actual quotient digit q,
2983
839k
     and it is never more than two larger than the actual quotient digit.  */
2984
839k
  k = s_norm(u, v);
2985
839k
2986
839k
  /* Extend size of u by one if needed.
2987
839k
2988
839k
     The algorithm begins with a value of u that has one more digit of input.
2989
839k
     The normalization step sets u_{m+n}..u_0 = 2^k * u_{m+n-1}..u_0. If the
2990
839k
     multiplication did not increase the number of digits of u, we need to add
2991
839k
     a leading zero here.
2992
839k
   */
2993
839k
  if (
k == 0 || 839k
MP_USED824k
(u) != m + n + 1824k
)
{318k
2994
318k
    if (!s_pad(u, m+n+1))
2995
0
      return MP_MEMORY;
2996
318k
    u->digits[m+n] = 0;
2997
318k
    u->used = m+n+1;
2998
318k
  }
2999
839k
3000
839k
  /* Add a leading 0 to v.
3001
839k
3002
839k
     The multiplication in step D4 multiplies qhat * 0v_{n-1}..v_0.  We need to
3003
839k
     add the leading zero to v here to ensure that the multiplication will
3004
839k
     produce the full n+1 digit result.  */
3005
839k
  
if (839k
!s_pad(v, n+1)839k
)
return MP_MEMORY0
;
v->digits[n] = 0;839k
3006
839k
3007
839k
  /* Initialize temporary variables q and t.
3008
839k
     q allocates space for m+1 digits to store the quotient digits
3009
839k
     t allocates space for n+1 digits to hold the result of q_j*v */
3010
839k
  if (
(res = mp_int_init_size(&q, m + 1)) != MP_OK839k
)
return res0
;
3011
839k
  
if (839k
(res = mp_int_init_size(&t, n + 1)) != MP_OK839k
)
goto CLEANUP0
;
3012
839k
3013
839k
  /************************************************************/
3014
839k
  /* D2: Initialize j */
3015
839k
  j = m;
3016
839k
  r.digits = MP_DIGITS(u) + j;  /* The contents of r are shared with u */
3017
839k
  r.used   = n + 1;
3018
839k
  r.sign   = MP_ZPOS;
3019
839k
  r.alloc  = MP_ALLOC(u);
3020
839k
  ZERO(t.digits, t.alloc);
3021
839k
3022
839k
  /* Calculate the m+1 digits of the quotient result */
3023
2.75M
  for (; 
j >= 02.75M
;
j--1.91M
)
{1.91M
3024
1.91M
    /************************************************************/
3025
1.91M
    /* D3: Calculate q' */
3026
1.91M
    /* r->digits is aligned to position j of the number u */
3027
1.91M
    mp_word pfx, qhat;
3028
1.91M
    pfx   = r.digits[n];
3029
1.91M
    pfx <<= MP_DIGIT_BIT / 2;
3030
1.91M
    pfx <<= MP_DIGIT_BIT / 2;
3031
1.91M
    pfx |= r.digits[n-1]; /* pfx = u_{j+n}{j+n-1} */
3032
1.91M
3033
1.91M
    qhat = pfx / v->digits[n-1];
3034
1.91M
    /* Check to see if qhat > b, and decrease qhat if so.
3035
1.91M
       Theorem B guarantess that qhat is at most 2 larger than the
3036
1.91M
       actual value, so it is possible that qhat is greater than
3037
1.91M
       the maximum value that will fit in a digit */
3038
1.91M
    if (
qhat > 1.91M
MP_DIGIT_MAX1.91M
)
3039
86
      
qhat = 86
MP_DIGIT_MAX86
;
3040
1.91M
3041
1.91M
    /************************************************************/
3042
1.91M
    /* D4,D5,D6: Multiply qhat * v and test for a correct value of q
3043
1.91M
3044
1.91M
       We proceed a bit different than the way described by Knuth. This way is
3045
1.91M
       simpler but less efficent. Instead of doing the multiply and subtract
3046
1.91M
       then checking for underflow, we first do the multiply of qhat * v and
3047
1.91M
       see if it is larger than the current remainder r. If it is larger, we
3048
1.91M
       decrease qhat by one and try again. We may need to decrease qhat one
3049
1.91M
       more time before we get a value that is smaller than r.
3050
1.91M
3051
1.91M
       This way is less efficent than Knuth becuase we do more multiplies, but
3052
1.91M
       we do not need to worry about underflow this way.  */
3053
1.91M
    /* t = qhat * v */
3054
1.91M
    s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1); t.used = n + 1;
3055
1.91M
    CLAMP(&t);
3056
1.91M
3057
1.91M
    /* Clamp r for the comparison. Comparisons do not like leading zeros. */
3058
1.91M
    CLAMP(&r);
3059
1.91M
    if (
s_ucmp(&t, &r) > 01.91M
)
{ /* would the remainder be negative? */157k
3060
157k
      qhat -= 1;   /* try a smaller q */
3061
157k
      s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1);
3062
157k
      t.used = n + 1; CLAMP(&t);
3063
157k
      if (
s_ucmp(&t, &r) > 0157k
)
{ /* would the remainder be negative? */4.39k
3064
4.39k
        assert(qhat > 0);
3065
4.39k
        qhat -= 1; /* try a smaller q */
3066
4.39k
        s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1);
3067
4.39k
        t.used = n + 1; CLAMP(&t);
3068
4.39k
      }
3069
157k
      assert(s_ucmp(&t, &r) <=  0 && "The mathematics failed us.");
3070
157k
    }
3071
1.91M
    /* Unclamp r. The D algorithm expects r = u_{j+n}..u_j to always be n+1
3072
1.91M
       digits long. */
3073
1.91M
    r.used = n + 1;
3074
1.91M
3075
1.91M
    /************************************************************/
3076
1.91M
    /* D4: Multiply and subtract */
3077
1.91M
    /* note: The multiply was completed above so we only need to subtract here.
3078
1.91M
     **/
3079
1.91M
    s_usub(r.digits, t.digits, r.digits, r.used, t.used);
3080
1.91M
3081
1.91M
    /************************************************************/
3082
1.91M
    /* D5: Test remainder */
3083
1.91M
    /* note: Not needed because we always check that qhat is the correct value
3084
1.91M
     *       before performing the subtract.
3085
1.91M
     *       Value cast to mp_digit to prevent warning, qhat has been clamped to MP_DIGIT_MAX */
3086
1.91M
    q.digits[j] = (mp_digit)qhat;
3087
1.91M
3088
1.91M
    /************************************************************/
3089
1.91M
    /* D6: Add back */
3090
1.91M
    /* note: Not needed because we always check that qhat is the correct value
3091
1.91M
     *       before performing the subtract. */
3092
1.91M
3093
1.91M
    /************************************************************/
3094
1.91M
    /* D7: Loop on j */
3095
1.91M
    r.digits--;
3096
1.91M
    ZERO(t.digits, t.alloc);
3097
1.91M
  }
3098
839k
3099
839k
  /* Get rid of leading zeros in q */
3100
839k
  q.used = m + 1;
3101
839k
  CLAMP(&q);
3102
839k
3103
839k
  /* Denormalize the remainder */
3104
839k
  CLAMP(u); /* use u here because the r.digits pointer is off-by-one */
3105
839k
  if (k != 0)
3106
824k
    s_qdiv(u, k);
3107
839k
3108
839k
  mp_int_copy(u, v);  /* ok:  0 <= r < v */
3109
839k
  mp_int_copy(&q, u); /* ok:  q <= u     */
3110
839k
3111
839k
  mp_int_clear(&t);
3112
839k
 CLEANUP:
3113
839k
  mp_int_clear(&q);
3114
839k
  return res;
3115
839k
}
3116
3117
STATIC int       s_outlen(mp_int z, mp_size r)
3118
1.43k
{
3119
1.43k
  mp_result bits;
3120
1.43k
  double raw;
3121
1.43k
3122
1.43k
  assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
3123
1.43k
3124
1.43k
  bits = mp_int_count_bits(z);
3125
1.43k
  raw = (double)bits * s_log2[r];
3126
1.43k
3127
1.43k
  return (int)(raw + 0.999999);
3128
1.43k
}
3129
3130
STATIC mp_size   s_inlen(int len, mp_size r)
3131
78
{
3132
78
  double  raw = (double)len / s_log2[r];
3133
78
  mp_size bits = (mp_size)(raw + 0.5);
3134
78
3135
78
  return (mp_size)((bits + (
MP_DIGIT_BIT78
- 1)) /
MP_DIGIT_BIT78
) + 1;
3136
78
}
3137
3138
STATIC int       s_ch2val(char c, int r)
3139
1.03k
{
3140
1.03k
  int out;
3141
1.03k
3142
1.03k
  if (isdigit((unsigned char) c))
3143
1.03k
    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.03k
3149
1.03k
  
return (out >= r) ? 1.03k
-10
:
out1.03k
;
3150
1.03k
}
3151
3152
STATIC char      s_val2ch(int v, int caps)
3153
24.9k
{
3154
24.9k
  assert(v >= 0);
3155
24.9k
3156
24.9k
  if (v < 10)
3157
24.9k
    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
24.9k
}
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 */