PostgreSQL Source Code  git master
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
imath.c
Go to the documentation of this file.
1 /* imath version 1.3 */
2 /*
3  Name: imath.c
4  Purpose: Arbitrary precision integer arithmetic routines.
5  Author: M. J. Fromberger <http://spinning-yarns.org/michael/sw/>
6  Info: Id: imath.c 21 2006-04-02 18:58:36Z sting
7 
8  Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
9 
10  Permission is hereby granted, free of charge, to any person
11  obtaining a copy of this software and associated documentation files
12  (the "Software"), to deal in the Software without restriction,
13  including without limitation the rights to use, copy, modify, merge,
14  publish, distribute, sublicense, and/or sell copies of the Software,
15  and to permit persons to whom the Software is furnished to do so,
16  subject to the following conditions:
17 
18  The above copyright notice and this permission notice shall be
19  included in all copies or substantial portions of the Software.
20 
21  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
22  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
23  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
24  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
25  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
26  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
27  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
28  SOFTWARE.
29  */
30 /* contrib/pgcrypto/imath.c */
31 
32 #include "postgres.h"
33 #include "px.h"
34 #include "imath.h"
35 
36 #undef assert
37 #define assert(TEST) Assert(TEST)
38 #define TRACEABLE_CLAMP 0
39 #define TRACEABLE_FREE 0
40 
41 /* {{{ Constants */
42 
43 const mp_result MP_OK = 0; /* no error, all is well */
44 const mp_result MP_FALSE = 0; /* boolean false */
45 const mp_result MP_TRUE = -1; /* boolean true */
46 const mp_result MP_MEMORY = -2; /* out of memory */
47 const mp_result MP_RANGE = -3; /* argument out of range */
48 const mp_result MP_UNDEF = -4; /* result undefined */
49 const mp_result MP_TRUNC = -5; /* output truncated */
50 const mp_result MP_BADARG = -6; /* invalid null argument */
51 
52 const mp_sign MP_NEG = 1; /* value is strictly negative */
53 const mp_sign MP_ZPOS = 0; /* value is non-negative */
54 
55 static const char *s_unknown_err = "unknown result code";
56 static const char *s_error_msg[] = {
57  "error code 0",
58  "boolean true",
59  "out of memory",
60  "argument out of range",
61  "result undefined",
62  "output truncated",
63  "invalid null argument",
64  NULL
65 };
66 
67 /* }}} */
68 
69 /* Optional library flags */
70 #define MP_CAP_DIGITS 1 /* flag bit to capitalize letter digits */
71 
72 /* Argument checking macros
73  Use CHECK() where a return value is required; NRCHECK() elsewhere */
74 #define CHECK(TEST) assert(TEST)
75 #define NRCHECK(TEST) assert(TEST)
76 
77 /* {{{ Logarithm table for computing output sizes */
78 
79 /* The ith entry of this table gives the value of log_i(2).
80 
81  An integer value n requires ceil(log_i(n)) digits to be represented
82  in base i. Since it is easy to compute lg(n), by counting bits, we
83  can compute log_i(n) = lg(n) * log_i(2).
84  */
85 static const double s_log2[] = {
86  0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
87  0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
88  0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
89  0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */
90  0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */
91  0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */
92  0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */
93  0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */
94  0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */
95  0.193426404, 0.191958720, 0.190551412, 0.189200360, /* 36 37 38 39 */
96  0.187901825, 0.186652411, 0.185449023, 0.184288833, /* 40 41 42 43 */
97  0.183169251, 0.182087900, 0.181042597, 0.180031327, /* 44 45 46 47 */
98  0.179052232, 0.178103594, 0.177183820, 0.176291434, /* 48 49 50 51 */
99  0.175425064, 0.174583430, 0.173765343, 0.172969690, /* 52 53 54 55 */
100  0.172195434, 0.171441601, 0.170707280, 0.169991616, /* 56 57 58 59 */
101  0.169293808, 0.168613099, 0.167948779, 0.167300179, /* 60 61 62 63 */
102  0.166666667
103 };
104 
105 /* }}} */
106 /* {{{ Various macros */
107 
108 /* Return the number of digits needed to represent a static value */
109 #define MP_VALUE_DIGITS(V) \
110 ((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
111 
112 /* Round precision P to nearest word boundary */
113 #define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
114 
115 /* Set array P of S digits to zero */
116 #define ZERO(P, S) \
117 do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P);memset(p__,0,i__);}while(0)
118 
119 /* Copy S digits from array P to array Q */
120 #define COPY(P, Q, S) \
121 do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P),*q__=(Q);\
122 memcpy(q__,p__,i__);}while(0)
123 
124 /* Reverse N elements of type T in array A */
125 #define REV(T, A, N) \
126 do{T *u_=(A),*v_=u_+(N)-1;while(u_<v_){T xch=*u_;*u_++=*v_;*v_--=xch;}}while(0)
127 
128 #if TRACEABLE_CLAMP
129 #define CLAMP(Z) s_clamp(Z)
130 #else
131 #define CLAMP(Z) \
132 do{mp_int z_=(Z);mp_size uz_=MP_USED(z_);mp_digit *dz_=MP_DIGITS(z_)+uz_-1;\
133 while(uz_ > 1 && (*dz_-- == 0)) --uz_;MP_USED(z_)=uz_;}while(0)
134 #endif
135 
136 #undef MIN
137 #undef MAX
138 #define MIN(A, B) ((B)<(A)?(B):(A))
139 #define MAX(A, B) ((B)>(A)?(B):(A))
140 #define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
141 
142 #define TEMP(K) (temp + (K))
143 #define SETUP(E, C) \
144 do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
145 
146 #define CMPZ(Z) \
147 (((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
148 
149 #define UMUL(X, Y, Z) \
150 do{mp_size ua_=MP_USED(X),ub_=MP_USED(Y);mp_size o_=ua_+ub_;\
151 ZERO(MP_DIGITS(Z),o_);\
152 (void) s_kmul(MP_DIGITS(X),MP_DIGITS(Y),MP_DIGITS(Z),ua_,ub_);\
153 MP_USED(Z)=o_;CLAMP(Z);}while(0)
154 
155 #define USQR(X, Z) \
156 do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
157 (void) s_ksqr(MP_DIGITS(X),MP_DIGITS(Z),ua_);MP_USED(Z)=o_;CLAMP(Z);}while(0)
158 
159 #define UPPER_HALF(W) ((mp_word)((W) >> MP_DIGIT_BIT))
160 #define LOWER_HALF(W) ((mp_digit)(W))
161 #define HIGH_BIT_SET(W) ((W) >> (MP_WORD_BIT - 1))
162 #define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
163 
164 /* }}} */
165 
166 /* Default number of digits allocated to a new mp_int */
168 
169 /* Minimum number of digits to invoke recursive multiply */
171 
172 /* Default library configuration flags */
174 
175 /* Allocate a buffer of (at least) num digits, or return
176  NULL if that couldn't be done. */
177 static mp_digit *s_alloc(mp_size num);
178 
179 #if TRACEABLE_FREE
180 static void s_free(void *ptr);
181 #else
182 #define s_free(P) px_free(P)
183 #endif
184 
185 /* Insure that z has at least min digits allocated, resizing if
186  necessary. Returns true if successful, false if out of memory. */
187 static int s_pad(mp_int z, mp_size min);
188 
189 /* Normalize by removing leading zeroes (except when z = 0) */
190 #if TRACEABLE_CLAMP
191 static void s_clamp(mp_int z);
192 #endif
193 
194 /* Fill in a "fake" mp_int on the stack with a given value */
195 static void s_fake(mp_int z, int value, mp_digit vbuf[]);
196 
197 /* Compare two runs of digits of given length, returns <0, 0, >0 */
198 static int s_cdig(mp_digit *da, mp_digit *db, mp_size len);
199 
200 /* Pack the unsigned digits of v into array t */
201 static int s_vpack(int v, mp_digit t[]);
202 
203 /* Compare magnitudes of a and b, returns <0, 0, >0 */
204 static int s_ucmp(mp_int a, mp_int b);
205 
206 /* Compare magnitudes of a and v, returns <0, 0, >0 */
207 static int s_vcmp(mp_int a, int v);
208 
209 /* Unsigned magnitude addition; assumes dc is big enough.
210  Carry out is returned (no memory allocated). */
211 static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
212  mp_size size_a, mp_size size_b);
213 
214 /* Unsigned magnitude subtraction. Assumes dc is big enough. */
215 static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
216  mp_size size_a, mp_size size_b);
217 
218 /* Unsigned recursive multiplication. Assumes dc is big enough. */
219 static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
220  mp_size size_a, mp_size size_b);
221 
222 /* Unsigned magnitude multiplication. Assumes dc is big enough. */
223 static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
224  mp_size size_a, mp_size size_b);
225 
226 /* Unsigned recursive squaring. Assumes dc is big enough. */
227 static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
228 
229 /* Unsigned magnitude squaring. Assumes dc is big enough. */
230 static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
231 
232 /* Single digit addition. Assumes a is big enough. */
233 static void s_dadd(mp_int a, mp_digit b);
234 
235 /* Single digit multiplication. Assumes a is big enough. */
236 static void s_dmul(mp_int a, mp_digit b);
237 
238 /* Single digit multiplication on buffers; assumes dc is big enough. */
239 static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
240  mp_size size_a);
241 
242 /* Single digit division. Replaces a with the quotient,
243  returns the remainder. */
244 static mp_digit s_ddiv(mp_int a, mp_digit b);
245 
246 /* Quick division by a power of 2, replaces z (no allocation) */
247 static void s_qdiv(mp_int z, mp_size p2);
248 
249 /* Quick remainder by a power of 2, replaces z (no allocation) */
250 static void s_qmod(mp_int z, mp_size p2);
251 
252 /* Quick multiplication by a power of 2, replaces z.
253  Allocates if necessary; returns false in case this fails. */
254 static int s_qmul(mp_int z, mp_size p2);
255 
256 /* Quick subtraction from a power of 2, replaces z.
257  Allocates if necessary; returns false in case this fails. */
258 static int s_qsub(mp_int z, mp_size p2);
259 
260 /* Return maximum k such that 2^k divides z. */
261 static int s_dp2k(mp_int z);
262 
263 /* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
264 static int s_isp2(mp_int z);
265 
266 /* Set z to 2^k. May allocate; returns false in case this fails. */
267 static int s_2expt(mp_int z, int k);
268 
269 /* Normalize a and b for division, returns normalization constant */
270 static int s_norm(mp_int a, mp_int b);
271 
272 /* Compute constant mu for Barrett reduction, given modulus m, result
273  replaces z, m is untouched. */
274 static mp_result s_brmu(mp_int z, mp_int m);
275 
276 /* Reduce a modulo m, using Barrett's algorithm. */
277 static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
278 
279 /* Modular exponentiation, using Barrett reduction */
280 static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
281 
282 /* Unsigned magnitude division. Assumes |a| > |b|. Allocates
283  temporaries; overwrites a with quotient, b with remainder. */
284 static mp_result s_udiv(mp_int a, mp_int b);
285 
286 /* Compute the number of digits in radix r required to represent the
287  given value. Does not account for sign flags, terminators, etc. */
288 static int s_outlen(mp_int z, mp_size r);
289 
290 /* Guess how many digits of precision will be needed to represent a
291  radix r value of the specified number of digits. Returns a value
292  guaranteed to be no smaller than the actual number required. */
293 static mp_size s_inlen(int len, mp_size r);
294 
295 /* Convert a character to a digit value in radix r, or
296  -1 if out of range */
297 static int s_ch2val(char c, int r);
298 
299 /* Convert a digit value to a character */
300 static char s_val2ch(int v, int caps);
301 
302 /* Take 2's complement of a buffer in place */
303 static void s_2comp(unsigned char *buf, int len);
304 
305 /* Convert a value to binary, ignoring sign. On input, *limpos is the
306  bound on how many bytes should be written to buf; on output, *limpos
307  is set to the number of bytes actually written. */
308 static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad);
309 
310 #if 0
311 /* Dump a representation of the mp_int to standard output */
312 void s_print(char *tag, mp_int z);
313 void s_print_buf(char *tag, mp_digit *buf, mp_size num);
314 #endif
315 
316 /* {{{ get_default_precision() */
317 
318 mp_size
320 {
321  return default_precision;
322 }
323 
324 /* }}} */
325 
326 /* {{{ mp_set_default_precision(s) */
327 
328 void
330 {
331  NRCHECK(s > 0);
332 
334 }
335 
336 /* }}} */
337 
338 /* {{{ mp_get_multiply_threshold() */
339 
340 mp_size
342 {
343  return multiply_threshold;
344 }
345 
346 /* }}} */
347 
348 /* {{{ mp_set_multiply_threshold(s) */
349 
350 void
352 {
353  multiply_threshold = s;
354 }
355 
356 /* }}} */
357 
358 /* {{{ mp_int_init(z) */
359 
360 mp_result
362 {
364 }
365 
366 /* }}} */
367 
368 /* {{{ mp_int_alloc() */
369 
370 mp_int
372 {
373  mp_int out = px_alloc(sizeof(mpz_t));
374 
375  assert(out != NULL);
376  out->digits = NULL;
377  out->used = 0;
378  out->alloc = 0;
379  out->sign = 0;
380 
381  return out;
382 }
383 
384 /* }}} */
385 
386 /* {{{ mp_int_init_size(z, prec) */
387 
388 mp_result
390 {
391  CHECK(z != NULL);
392 
393  prec = (mp_size) ROUND_PREC(prec);
394  prec = MAX(prec, default_precision);
395 
396  if ((MP_DIGITS(z) = s_alloc(prec)) == NULL)
397  return MP_MEMORY;
398 
399  z->digits[0] = 0;
400  MP_USED(z) = 1;
401  MP_ALLOC(z) = prec;
402  MP_SIGN(z) = MP_ZPOS;
403 
404  return MP_OK;
405 }
406 
407 /* }}} */
408 
409 /* {{{ mp_int_init_copy(z, old) */
410 
411 mp_result
413 {
414  mp_result res;
415  mp_size uold,
416  target;
417 
418  CHECK(z != NULL && old != NULL);
419 
420  uold = MP_USED(old);
421  target = MAX(uold, default_precision);
422 
423  if ((res = mp_int_init_size(z, target)) != MP_OK)
424  return res;
425 
426  MP_USED(z) = uold;
427  MP_SIGN(z) = MP_SIGN(old);
428  COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
429 
430  return MP_OK;
431 }
432 
433 /* }}} */
434 
435 /* {{{ mp_int_init_value(z, value) */
436 
437 mp_result
439 {
440  mp_result res;
441 
442  CHECK(z != NULL);
443 
444  if ((res = mp_int_init(z)) != MP_OK)
445  return res;
446 
447  return mp_int_set_value(z, value);
448 }
449 
450 /* }}} */
451 
452 /* {{{ mp_int_set_value(z, value) */
453 
454 mp_result
456 {
457  mp_size ndig;
458 
459  CHECK(z != NULL);
460 
461  /* How many digits to copy */
462  ndig = (mp_size) MP_VALUE_DIGITS(value);
463 
464  if (!s_pad(z, ndig))
465  return MP_MEMORY;
466 
467  MP_USED(z) = (mp_size) s_vpack(value, MP_DIGITS(z));
468  MP_SIGN(z) = (value < 0) ? MP_NEG : MP_ZPOS;
469 
470  return MP_OK;
471 }
472 
473 /* }}} */
474 
475 /* {{{ mp_int_clear(z) */
476 
477 void
479 {
480  if (z == NULL)
481  return;
482 
483  if (MP_DIGITS(z) != NULL)
484  {
485  s_free(MP_DIGITS(z));
486  MP_DIGITS(z) = NULL;
487  }
488 }
489 
490 /* }}} */
491 
492 /* {{{ mp_int_free(z) */
493 
494 void
496 {
497  NRCHECK(z != NULL);
498 
499  if (z->digits != NULL)
500  mp_int_clear(z);
501 
502  px_free(z);
503 }
504 
505 /* }}} */
506 
507 /* {{{ mp_int_copy(a, c) */
508 
509 mp_result
511 {
512  CHECK(a != NULL && c != NULL);
513 
514  if (a != c)
515  {
516  mp_size ua = MP_USED(a);
517  mp_digit *da,
518  *dc;
519 
520  if (!s_pad(c, ua))
521  return MP_MEMORY;
522 
523  da = MP_DIGITS(a);
524  dc = MP_DIGITS(c);
525  COPY(da, dc, ua);
526 
527  MP_USED(c) = ua;
528  MP_SIGN(c) = MP_SIGN(a);
529  }
530 
531  return MP_OK;
532 }
533 
534 /* }}} */
535 
536 /* {{{ mp_int_swap(a, c) */
537 
538 void
540 {
541  if (a != c)
542  {
543  mpz_t tmp = *a;
544 
545  *a = *c;
546  *c = tmp;
547  }
548 }
549 
550 /* }}} */
551 
552 /* {{{ mp_int_zero(z) */
553 
554 void
556 {
557  NRCHECK(z != NULL);
558 
559  z->digits[0] = 0;
560  MP_USED(z) = 1;
561  MP_SIGN(z) = MP_ZPOS;
562 }
563 
564 /* }}} */
565 
566 /* {{{ mp_int_abs(a, c) */
567 
568 mp_result
570 {
571  mp_result res;
572 
573  CHECK(a != NULL && c != NULL);
574 
575  if ((res = mp_int_copy(a, c)) != MP_OK)
576  return res;
577 
578  MP_SIGN(c) = MP_ZPOS;
579  return MP_OK;
580 }
581 
582 /* }}} */
583 
584 /* {{{ mp_int_neg(a, c) */
585 
586 mp_result
588 {
589  mp_result res;
590 
591  CHECK(a != NULL && c != NULL);
592 
593  if ((res = mp_int_copy(a, c)) != MP_OK)
594  return res;
595 
596  if (CMPZ(c) != 0)
597  MP_SIGN(c) = 1 - MP_SIGN(a);
598 
599  return MP_OK;
600 }
601 
602 /* }}} */
603 
604 /* {{{ mp_int_add(a, b, c) */
605 
606 mp_result
608 {
609  mp_size ua,
610  ub,
611  uc,
612  max;
613 
614  CHECK(a != NULL && b != NULL && c != NULL);
615 
616  ua = MP_USED(a);
617  ub = MP_USED(b);
618  uc = MP_USED(c);
619  max = MAX(ua, ub);
620 
621  if (MP_SIGN(a) == MP_SIGN(b))
622  {
623  /* Same sign -- add magnitudes, preserve sign of addends */
624  mp_digit carry;
625 
626  if (!s_pad(c, max))
627  return MP_MEMORY;
628 
629  carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
630  uc = max;
631 
632  if (carry)
633  {
634  if (!s_pad(c, max + 1))
635  return MP_MEMORY;
636 
637  c->digits[max] = carry;
638  ++uc;
639  }
640 
641  MP_USED(c) = uc;
642  MP_SIGN(c) = MP_SIGN(a);
643 
644  }
645  else
646  {
647  /* Different signs -- subtract magnitudes, preserve sign of greater */
648  mp_int x,
649  y;
650  int cmp = s_ucmp(a, b); /* magnitude comparison, sign ignored */
651 
652  /* Set x to max(a, b), y to min(a, b) to simplify later code */
653  if (cmp >= 0)
654  {
655  x = a;
656  y = b;
657  }
658  else
659  {
660  x = b;
661  y = a;
662  }
663 
664  if (!s_pad(c, MP_USED(x)))
665  return MP_MEMORY;
666 
667  /* Subtract smaller from larger */
668  s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
669  MP_USED(c) = MP_USED(x);
670  CLAMP(c);
671 
672  /* Give result the sign of the larger */
673  MP_SIGN(c) = MP_SIGN(x);
674  }
675 
676  return MP_OK;
677 }
678 
679 /* }}} */
680 
681 /* {{{ mp_int_add_value(a, value, c) */
682 
683 mp_result
685 {
686  mpz_t vtmp;
687  mp_digit vbuf[MP_VALUE_DIGITS(value)];
688 
689  s_fake(&vtmp, value, vbuf);
690 
691  return mp_int_add(a, &vtmp, c);
692 }
693 
694 /* }}} */
695 
696 /* {{{ mp_int_sub(a, b, c) */
697 
698 mp_result
700 {
701  mp_size ua,
702  ub,
703  uc,
704  max;
705 
706  CHECK(a != NULL && b != NULL && c != NULL);
707 
708  ua = MP_USED(a);
709  ub = MP_USED(b);
710  uc = MP_USED(c);
711  max = MAX(ua, ub);
712 
713  if (MP_SIGN(a) != MP_SIGN(b))
714  {
715  /* Different signs -- add magnitudes and keep sign of a */
716  mp_digit carry;
717 
718  if (!s_pad(c, max))
719  return MP_MEMORY;
720 
721  carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
722  uc = max;
723 
724  if (carry)
725  {
726  if (!s_pad(c, max + 1))
727  return MP_MEMORY;
728 
729  c->digits[max] = carry;
730  ++uc;
731  }
732 
733  MP_USED(c) = uc;
734  MP_SIGN(c) = MP_SIGN(a);
735 
736  }
737  else
738  {
739  /* Same signs -- subtract magnitudes */
740  mp_int x,
741  y;
742  mp_sign osign;
743  int cmp = s_ucmp(a, b);
744 
745  if (!s_pad(c, max))
746  return MP_MEMORY;
747 
748  if (cmp >= 0)
749  {
750  x = a;
751  y = b;
752  osign = MP_ZPOS;
753  }
754  else
755  {
756  x = b;
757  y = a;
758  osign = MP_NEG;
759  }
760 
761  if (MP_SIGN(a) == MP_NEG && cmp != 0)
762  osign = 1 - osign;
763 
764  s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
765  MP_USED(c) = MP_USED(x);
766  CLAMP(c);
767 
768  MP_SIGN(c) = osign;
769  }
770 
771  return MP_OK;
772 }
773 
774 /* }}} */
775 
776 /* {{{ mp_int_sub_value(a, value, c) */
777 
778 mp_result
780 {
781  mpz_t vtmp;
782  mp_digit vbuf[MP_VALUE_DIGITS(value)];
783 
784  s_fake(&vtmp, value, vbuf);
785 
786  return mp_int_sub(a, &vtmp, c);
787 }
788 
789 /* }}} */
790 
791 /* {{{ mp_int_mul(a, b, c) */
792 
793 mp_result
795 {
796  mp_digit *out;
797  mp_size osize,
798  ua,
799  ub,
800  p = 0;
801  mp_sign osign;
802 
803  CHECK(a != NULL && b != NULL && c != NULL);
804 
805  /* If either input is zero, we can shortcut multiplication */
806  if (mp_int_compare_zero(a) == 0 || mp_int_compare_zero(b) == 0)
807  {
808  mp_int_zero(c);
809  return MP_OK;
810  }
811 
812  /* Output is positive if inputs have same sign, otherwise negative */
813  osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG;
814 
815  /*
816  * If the output is not equal to any of the inputs, we'll write the
817  * results there directly; otherwise, allocate a temporary space.
818  */
819  ua = MP_USED(a);
820  ub = MP_USED(b);
821  osize = MAX(ua, ub);
822  osize = 4 * ((osize + 1) / 2);
823 
824  if (c == a || c == b)
825  {
826  p = ROUND_PREC(osize);
827  p = MAX(p, default_precision);
828 
829  if ((out = s_alloc(p)) == NULL)
830  return MP_MEMORY;
831  }
832  else
833  {
834  if (!s_pad(c, osize))
835  return MP_MEMORY;
836 
837  out = MP_DIGITS(c);
838  }
839  ZERO(out, osize);
840 
841  if (!s_kmul(MP_DIGITS(a), MP_DIGITS(b), out, ua, ub))
842  return MP_MEMORY;
843 
844  /*
845  * If we allocated a new buffer, get rid of whatever memory c was already
846  * using, and fix up its fields to reflect that.
847  */
848  if (out != MP_DIGITS(c))
849  {
850  s_free(MP_DIGITS(c));
851  MP_DIGITS(c) = out;
852  MP_ALLOC(c) = p;
853  }
854 
855  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
856  CLAMP(c); /* ... right here */
857  MP_SIGN(c) = osign;
858 
859  return MP_OK;
860 }
861 
862 /* }}} */
863 
864 /* {{{ mp_int_mul_value(a, value, c) */
865 
866 mp_result
868 {
869  mpz_t vtmp;
870  mp_digit vbuf[MP_VALUE_DIGITS(value)];
871 
872  s_fake(&vtmp, value, vbuf);
873 
874  return mp_int_mul(a, &vtmp, c);
875 }
876 
877 /* }}} */
878 
879 /* {{{ mp_int_mul_pow2(a, p2, c) */
880 
881 mp_result
883 {
884  mp_result res;
885 
886  CHECK(a != NULL && c != NULL && p2 >= 0);
887 
888  if ((res = mp_int_copy(a, c)) != MP_OK)
889  return res;
890 
891  if (s_qmul(c, (mp_size) p2))
892  return MP_OK;
893  else
894  return MP_MEMORY;
895 }
896 
897 /* }}} */
898 
899 /* {{{ mp_int_sqr(a, c) */
900 
901 mp_result
903 {
904  mp_digit *out;
905  mp_size osize,
906  p = 0;
907 
908  CHECK(a != NULL && c != NULL);
909 
910  /* Get a temporary buffer big enough to hold the result */
911  osize = (mp_size) 4 *((MP_USED(a) + 1) / 2);
912 
913  if (a == c)
914  {
915  p = ROUND_PREC(osize);
916  p = MAX(p, default_precision);
917 
918  if ((out = s_alloc(p)) == NULL)
919  return MP_MEMORY;
920  }
921  else
922  {
923  if (!s_pad(c, osize))
924  return MP_MEMORY;
925 
926  out = MP_DIGITS(c);
927  }
928  ZERO(out, osize);
929 
930  s_ksqr(MP_DIGITS(a), out, MP_USED(a));
931 
932  /*
933  * Get rid of whatever memory c was already using, and fix up its fields
934  * to reflect the new digit array it's using
935  */
936  if (out != MP_DIGITS(c))
937  {
938  s_free(MP_DIGITS(c));
939  MP_DIGITS(c) = out;
940  MP_ALLOC(c) = p;
941  }
942 
943  MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
944  CLAMP(c); /* ... right here */
945  MP_SIGN(c) = MP_ZPOS;
946 
947  return MP_OK;
948 }
949 
950 /* }}} */
951 
952 /* {{{ mp_int_div(a, b, q, r) */
953 
954 mp_result
956 {
957  int cmp,
958  last = 0,
959  lg;
960  mp_result res = MP_OK;
961  mpz_t temp[2];
962  mp_int qout,
963  rout;
964  mp_sign sa = MP_SIGN(a),
965  sb = MP_SIGN(b);
966 
967  CHECK(a != NULL && b != NULL && q != r);
968 
969  if (CMPZ(b) == 0)
970  return MP_UNDEF;
971  else if ((cmp = s_ucmp(a, b)) < 0)
972  {
973  /*
974  * If |a| < |b|, no division is required: q = 0, r = a
975  */
976  if (r && (res = mp_int_copy(a, r)) != MP_OK)
977  return res;
978 
979  if (q)
980  mp_int_zero(q);
981 
982  return MP_OK;
983  }
984  else if (cmp == 0)
985  {
986  /*
987  * If |a| = |b|, no division is required: q = 1 or -1, r = 0
988  */
989  if (r)
990  mp_int_zero(r);
991 
992  if (q)
993  {
994  mp_int_zero(q);
995  q->digits[0] = 1;
996 
997  if (sa != sb)
998  MP_SIGN(q) = MP_NEG;
999  }
1000 
1001  return MP_OK;
1002  }
1003 
1004  /*
1005  * When |a| > |b|, real division is required. We need someplace to store
1006  * quotient and remainder, but q and r are allowed to be NULL or to
1007  * overlap with the inputs.
1008  */
1009  if ((lg = s_isp2(b)) < 0)
1010  {
1011  if (q && b != q && (res = mp_int_copy(a, q)) == MP_OK)
1012  {
1013  qout = q;
1014  }
1015  else
1016  {
1017  qout = TEMP(last);
1018  SETUP(mp_int_init_copy(TEMP(last), a), last);
1019  }
1020 
1021  if (r && a != r && (res = mp_int_copy(b, r)) == MP_OK)
1022  {
1023  rout = r;
1024  }
1025  else
1026  {
1027  rout = TEMP(last);
1028  SETUP(mp_int_init_copy(TEMP(last), b), last);
1029  }
1030 
1031  if ((res = s_udiv(qout, rout)) != MP_OK)
1032  goto CLEANUP;
1033  }
1034  else
1035  {
1036  if (q && (res = mp_int_copy(a, q)) != MP_OK)
1037  goto CLEANUP;
1038  if (r && (res = mp_int_copy(a, r)) != MP_OK)
1039  goto CLEANUP;
1040 
1041  if (q)
1042  s_qdiv(q, (mp_size) lg);
1043  qout = q;
1044  if (r)
1045  s_qmod(r, (mp_size) lg);
1046  rout = r;
1047  }
1048 
1049  /* Recompute signs for output */
1050  if (rout)
1051  {
1052  MP_SIGN(rout) = sa;
1053  if (CMPZ(rout) == 0)
1054  MP_SIGN(rout) = MP_ZPOS;
1055  }
1056  if (qout)
1057  {
1058  MP_SIGN(qout) = (sa == sb) ? MP_ZPOS : MP_NEG;
1059  if (CMPZ(qout) == 0)
1060  MP_SIGN(qout) = MP_ZPOS;
1061  }
1062 
1063  if (q && (res = mp_int_copy(qout, q)) != MP_OK)
1064  goto CLEANUP;
1065  if (r && (res = mp_int_copy(rout, r)) != MP_OK)
1066  goto CLEANUP;
1067 
1068 CLEANUP:
1069  while (--last >= 0)
1070  mp_int_clear(TEMP(last));
1071 
1072  return res;
1073 }
1074 
1075 /* }}} */
1076 
1077 /* {{{ mp_int_mod(a, m, c) */
1078 
1079 mp_result
1081 {
1082  mp_result res;
1083  mpz_t tmp;
1084  mp_int out;
1085 
1086  if (m == c)
1087  {
1088  if ((res = mp_int_init(&tmp)) != MP_OK)
1089  return res;
1090 
1091  out = &tmp;
1092  }
1093  else
1094  {
1095  out = c;
1096  }
1097 
1098  if ((res = mp_int_div(a, m, NULL, out)) != MP_OK)
1099  goto CLEANUP;
1100 
1101  if (CMPZ(out) < 0)
1102  res = mp_int_add(out, m, c);
1103  else
1104  res = mp_int_copy(out, c);
1105 
1106 CLEANUP:
1107  if (out != c)
1108  mp_int_clear(&tmp);
1109 
1110  return res;
1111 }
1112 
1113 /* }}} */
1114 
1115 
1116 /* {{{ mp_int_div_value(a, value, q, r) */
1117 
1118 mp_result
1120 {
1121  mpz_t vtmp,
1122  rtmp;
1123  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1124  mp_result res;
1125 
1126  if ((res = mp_int_init(&rtmp)) != MP_OK)
1127  return res;
1128  s_fake(&vtmp, value, vbuf);
1129 
1130  if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
1131  goto CLEANUP;
1132 
1133  if (r)
1134  (void) mp_int_to_int(&rtmp, r); /* can't fail */
1135 
1136 CLEANUP:
1137  mp_int_clear(&rtmp);
1138  return res;
1139 }
1140 
1141 /* }}} */
1142 
1143 /* {{{ mp_int_div_pow2(a, p2, q, r) */
1144 
1145 mp_result
1147 {
1148  mp_result res = MP_OK;
1149 
1150  CHECK(a != NULL && p2 >= 0 && q != r);
1151 
1152  if (q != NULL && (res = mp_int_copy(a, q)) == MP_OK)
1153  s_qdiv(q, (mp_size) p2);
1154 
1155  if (res == MP_OK && r != NULL && (res = mp_int_copy(a, r)) == MP_OK)
1156  s_qmod(r, (mp_size) p2);
1157 
1158  return res;
1159 }
1160 
1161 /* }}} */
1162 
1163 /* {{{ mp_int_expt(a, b, c) */
1164 
1165 mp_result
1167 {
1168  mpz_t t;
1169  mp_result res;
1170  unsigned int v = abs(b);
1171 
1172  CHECK(b >= 0 && c != NULL);
1173 
1174  if ((res = mp_int_init_copy(&t, a)) != MP_OK)
1175  return res;
1176 
1177  (void) mp_int_set_value(c, 1);
1178  while (v != 0)
1179  {
1180  if (v & 1)
1181  {
1182  if ((res = mp_int_mul(c, &t, c)) != MP_OK)
1183  goto CLEANUP;
1184  }
1185 
1186  v >>= 1;
1187  if (v == 0)
1188  break;
1189 
1190  if ((res = mp_int_sqr(&t, &t)) != MP_OK)
1191  goto CLEANUP;
1192  }
1193 
1194 CLEANUP:
1195  mp_int_clear(&t);
1196  return res;
1197 }
1198 
1199 /* }}} */
1200 
1201 /* {{{ mp_int_expt_value(a, b, c) */
1202 
1203 mp_result
1205 {
1206  mpz_t t;
1207  mp_result res;
1208  unsigned int v = abs(b);
1209 
1210  CHECK(b >= 0 && c != NULL);
1211 
1212  if ((res = mp_int_init_value(&t, a)) != MP_OK)
1213  return res;
1214 
1215  (void) mp_int_set_value(c, 1);
1216  while (v != 0)
1217  {
1218  if (v & 1)
1219  {
1220  if ((res = mp_int_mul(c, &t, c)) != MP_OK)
1221  goto CLEANUP;
1222  }
1223 
1224  v >>= 1;
1225  if (v == 0)
1226  break;
1227 
1228  if ((res = mp_int_sqr(&t, &t)) != MP_OK)
1229  goto CLEANUP;
1230  }
1231 
1232 CLEANUP:
1233  mp_int_clear(&t);
1234  return res;
1235 }
1236 
1237 /* }}} */
1238 
1239 /* {{{ mp_int_compare(a, b) */
1240 
1241 int
1243 {
1244  mp_sign sa;
1245 
1246  CHECK(a != NULL && b != NULL);
1247 
1248  sa = MP_SIGN(a);
1249  if (sa == MP_SIGN(b))
1250  {
1251  int cmp = s_ucmp(a, b);
1252 
1253  /*
1254  * If they're both zero or positive, the normal comparison applies; if
1255  * both negative, the sense is reversed.
1256  */
1257  if (sa == MP_ZPOS)
1258  return cmp;
1259  else
1260  return -cmp;
1261 
1262  }
1263  else
1264  {
1265  if (sa == MP_ZPOS)
1266  return 1;
1267  else
1268  return -1;
1269  }
1270 }
1271 
1272 /* }}} */
1273 
1274 /* {{{ mp_int_compare_unsigned(a, b) */
1275 
1276 int
1278 {
1279  NRCHECK(a != NULL && b != NULL);
1280 
1281  return s_ucmp(a, b);
1282 }
1283 
1284 /* }}} */
1285 
1286 /* {{{ mp_int_compare_zero(z) */
1287 
1288 int
1290 {
1291  NRCHECK(z != NULL);
1292 
1293  if (MP_USED(z) == 1 && z->digits[0] == 0)
1294  return 0;
1295  else if (MP_SIGN(z) == MP_ZPOS)
1296  return 1;
1297  else
1298  return -1;
1299 }
1300 
1301 /* }}} */
1302 
1303 /* {{{ mp_int_compare_value(z, value) */
1304 
1305 int
1307 {
1308  mp_sign vsign = (value < 0) ? MP_NEG : MP_ZPOS;
1309  int cmp;
1310 
1311  CHECK(z != NULL);
1312 
1313  if (vsign == MP_SIGN(z))
1314  {
1315  cmp = s_vcmp(z, value);
1316 
1317  if (vsign == MP_ZPOS)
1318  return cmp;
1319  else
1320  return -cmp;
1321  }
1322  else
1323  {
1324  if (value < 0)
1325  return 1;
1326  else
1327  return -1;
1328  }
1329 }
1330 
1331 /* }}} */
1332 
1333 /* {{{ mp_int_exptmod(a, b, m, c) */
1334 
1335 mp_result
1337 {
1338  mp_result res;
1339  mp_size um;
1340  mpz_t temp[3];
1341  mp_int s;
1342  int last = 0;
1343 
1344  CHECK(a != NULL && b != NULL && c != NULL && m != NULL);
1345 
1346  /* Zero moduli and negative exponents are not considered. */
1347  if (CMPZ(m) == 0)
1348  return MP_UNDEF;
1349  if (CMPZ(b) < 0)
1350  return MP_RANGE;
1351 
1352  um = MP_USED(m);
1353  SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
1354  SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
1355 
1356  if (c == b || c == m)
1357  {
1358  SETUP(mp_int_init_size(TEMP(2), 2 * um), last);
1359  s = TEMP(2);
1360  }
1361  else
1362  {
1363  s = c;
1364  }
1365 
1366  if ((res = mp_int_mod(a, m, TEMP(0))) != MP_OK)
1367  goto CLEANUP;
1368 
1369  if ((res = s_brmu(TEMP(1), m)) != MP_OK)
1370  goto CLEANUP;
1371 
1372  if ((res = s_embar(TEMP(0), b, m, TEMP(1), s)) != MP_OK)
1373  goto CLEANUP;
1374 
1375  res = mp_int_copy(s, c);
1376 
1377 CLEANUP:
1378  while (--last >= 0)
1379  mp_int_clear(TEMP(last));
1380 
1381  return res;
1382 }
1383 
1384 /* }}} */
1385 
1386 /* {{{ mp_int_exptmod_evalue(a, value, m, c) */
1387 
1388 mp_result
1390 {
1391  mpz_t vtmp;
1392  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1393 
1394  s_fake(&vtmp, value, vbuf);
1395 
1396  return mp_int_exptmod(a, &vtmp, m, c);
1397 }
1398 
1399 /* }}} */
1400 
1401 /* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
1402 
1403 mp_result
1405  mp_int m, mp_int c)
1406 {
1407  mpz_t vtmp;
1408  mp_digit vbuf[MP_VALUE_DIGITS(value)];
1409 
1410  s_fake(&vtmp, value, vbuf);
1411 
1412  return mp_int_exptmod(&vtmp, b, m, c);
1413 }
1414 
1415 /* }}} */
1416 
1417 /* {{{ mp_int_exptmod_known(a, b, m, mu, c) */
1418 
1419 mp_result
1421 {
1422  mp_result res;
1423  mp_size um;
1424  mpz_t temp[2];
1425  mp_int s;
1426  int last = 0;
1427 
1428  CHECK(a && b && m && c);
1429 
1430  /* Zero moduli and negative exponents are not considered. */
1431  if (CMPZ(m) == 0)
1432  return MP_UNDEF;
1433  if (CMPZ(b) < 0)
1434  return MP_RANGE;
1435 
1436  um = MP_USED(m);
1437  SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
1438 
1439  if (c == b || c == m)
1440  {
1441  SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
1442  s = TEMP(1);
1443  }
1444  else
1445  {
1446  s = c;
1447  }
1448 
1449  if ((res = mp_int_mod(a, m, TEMP(0))) != MP_OK)
1450  goto CLEANUP;
1451 
1452  if ((res = s_embar(TEMP(0), b, m, mu, s)) != MP_OK)
1453  goto CLEANUP;
1454 
1455  res = mp_int_copy(s, c);
1456 
1457 CLEANUP:
1458  while (--last >= 0)
1459  mp_int_clear(TEMP(last));
1460 
1461  return res;
1462 }
1463 
1464 /* }}} */
1465 
1466 /* {{{ mp_int_redux_const(m, c) */
1467 
1468 mp_result
1470 {
1471  CHECK(m != NULL && c != NULL && m != c);
1472 
1473  return s_brmu(c, m);
1474 }
1475 
1476 /* }}} */
1477 
1478 /* {{{ mp_int_invmod(a, m, c) */
1479 
1480 mp_result
1482 {
1483  mp_result res;
1484  mp_sign sa;
1485  int last = 0;
1486  mpz_t temp[2];
1487 
1488  CHECK(a != NULL && m != NULL && c != NULL);
1489 
1490  if (CMPZ(a) == 0 || CMPZ(m) <= 0)
1491  return MP_RANGE;
1492 
1493  sa = MP_SIGN(a); /* need this for the result later */
1494 
1495  for (last = 0; last < 2; ++last)
1496  if ((res = mp_int_init(TEMP(last))) != MP_OK)
1497  goto CLEANUP;
1498 
1499  if ((res = mp_int_egcd(a, m, TEMP(0), TEMP(1), NULL)) != MP_OK)
1500  goto CLEANUP;
1501 
1502  if (mp_int_compare_value(TEMP(0), 1) != 0)
1503  {
1504  res = MP_UNDEF;
1505  goto CLEANUP;
1506  }
1507 
1508  /* It is first necessary to constrain the value to the proper range */
1509  if ((res = mp_int_mod(TEMP(1), m, TEMP(1))) != MP_OK)
1510  goto CLEANUP;
1511 
1512  /*
1513  * Now, if 'a' was originally negative, the value we have is actually the
1514  * magnitude of the negative representative; to get the positive value we
1515  * have to subtract from the modulus. Otherwise, the value is okay as it
1516  * stands.
1517  */
1518  if (sa == MP_NEG)
1519  res = mp_int_sub(m, TEMP(1), c);
1520  else
1521  res = mp_int_copy(TEMP(1), c);
1522 
1523 CLEANUP:
1524  while (--last >= 0)
1525  mp_int_clear(TEMP(last));
1526 
1527  return res;
1528 }
1529 
1530 /* }}} */
1531 
1532 /* {{{ mp_int_gcd(a, b, c) */
1533 
1534 /* Binary GCD algorithm due to Josef Stein, 1961 */
1535 mp_result
1537 {
1538  int ca,
1539  cb,
1540  k = 0;
1541  mpz_t u,
1542  v,
1543  t;
1544  mp_result res;
1545 
1546  CHECK(a != NULL && b != NULL && c != NULL);
1547 
1548  ca = CMPZ(a);
1549  cb = CMPZ(b);
1550  if (ca == 0 && cb == 0)
1551  return MP_UNDEF;
1552  else if (ca == 0)
1553  return mp_int_abs(b, c);
1554  else if (cb == 0)
1555  return mp_int_abs(a, c);
1556 
1557  if ((res = mp_int_init(&t)) != MP_OK)
1558  return res;
1559  if ((res = mp_int_init_copy(&u, a)) != MP_OK)
1560  goto U;
1561  if ((res = mp_int_init_copy(&v, b)) != MP_OK)
1562  goto V;
1563 
1564  MP_SIGN(&u) = MP_ZPOS;
1565  MP_SIGN(&v) = MP_ZPOS;
1566 
1567  { /* Divide out common factors of 2 from u and v */
1568  int div2_u = s_dp2k(&u),
1569  div2_v = s_dp2k(&v);
1570 
1571  k = MIN(div2_u, div2_v);
1572  s_qdiv(&u, (mp_size) k);
1573  s_qdiv(&v, (mp_size) k);
1574  }
1575 
1576  if (mp_int_is_odd(&u))
1577  {
1578  if ((res = mp_int_neg(&v, &t)) != MP_OK)
1579  goto CLEANUP;
1580  }
1581  else
1582  {
1583  if ((res = mp_int_copy(&u, &t)) != MP_OK)
1584  goto CLEANUP;
1585  }
1586 
1587  for (;;)
1588  {
1589  s_qdiv(&t, s_dp2k(&t));
1590 
1591  if (CMPZ(&t) > 0)
1592  {
1593  if ((res = mp_int_copy(&t, &u)) != MP_OK)
1594  goto CLEANUP;
1595  }
1596  else
1597  {
1598  if ((res = mp_int_neg(&t, &v)) != MP_OK)
1599  goto CLEANUP;
1600  }
1601 
1602  if ((res = mp_int_sub(&u, &v, &t)) != MP_OK)
1603  goto CLEANUP;
1604 
1605  if (CMPZ(&t) == 0)
1606  break;
1607  }
1608 
1609  if ((res = mp_int_abs(&u, c)) != MP_OK)
1610  goto CLEANUP;
1611  if (!s_qmul(c, (mp_size) k))
1612  res = MP_MEMORY;
1613 
1614 CLEANUP:
1615  mp_int_clear(&v);
1616 V: mp_int_clear(&u);
1617 U: mp_int_clear(&t);
1618 
1619  return res;
1620 }
1621 
1622 /* }}} */
1623 
1624 /* {{{ mp_int_egcd(a, b, c, x, y) */
1625 
1626 /* This is the binary GCD algorithm again, but this time we keep track
1627  of the elementary matrix operations as we go, so we can get values
1628  x and y satisfying c = ax + by.
1629  */
1630 mp_result
1632  mp_int x, mp_int y)
1633 {
1634  int k,
1635  last = 0,
1636  ca,
1637  cb;
1638  mpz_t temp[8];
1639  mp_result res;
1640 
1641  CHECK(a != NULL && b != NULL && c != NULL &&
1642  (x != NULL || y != NULL));
1643 
1644  ca = CMPZ(a);
1645  cb = CMPZ(b);
1646  if (ca == 0 && cb == 0)
1647  return MP_UNDEF;
1648  else if (ca == 0)
1649  {
1650  if ((res = mp_int_abs(b, c)) != MP_OK)
1651  return res;
1652  mp_int_zero(x);
1653  (void) mp_int_set_value(y, 1);
1654  return MP_OK;
1655  }
1656  else if (cb == 0)
1657  {
1658  if ((res = mp_int_abs(a, c)) != MP_OK)
1659  return res;
1660  (void) mp_int_set_value(x, 1);
1661  mp_int_zero(y);
1662  return MP_OK;
1663  }
1664 
1665  /*
1666  * Initialize temporaries: A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7
1667  */
1668  for (last = 0; last < 4; ++last)
1669  {
1670  if ((res = mp_int_init(TEMP(last))) != MP_OK)
1671  goto CLEANUP;
1672  }
1673  TEMP(0)->digits[0] = 1;
1674  TEMP(3)->digits[0] = 1;
1675 
1676  SETUP(mp_int_init_copy(TEMP(4), a), last);
1677  SETUP(mp_int_init_copy(TEMP(5), b), last);
1678 
1679  /* We will work with absolute values here */
1680  MP_SIGN(TEMP(4)) = MP_ZPOS;
1681  MP_SIGN(TEMP(5)) = MP_ZPOS;
1682 
1683  { /* Divide out common factors of 2 from u and v */
1684  int div2_u = s_dp2k(TEMP(4)),
1685  div2_v = s_dp2k(TEMP(5));
1686 
1687  k = MIN(div2_u, div2_v);
1688  s_qdiv(TEMP(4), k);
1689  s_qdiv(TEMP(5), k);
1690  }
1691 
1692  SETUP(mp_int_init_copy(TEMP(6), TEMP(4)), last);
1693  SETUP(mp_int_init_copy(TEMP(7), TEMP(5)), last);
1694 
1695  for (;;)
1696  {
1697  while (mp_int_is_even(TEMP(4)))
1698  {
1699  s_qdiv(TEMP(4), 1);
1700 
1701  if (mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1)))
1702  {
1703  if ((res = mp_int_add(TEMP(0), TEMP(7), TEMP(0))) != MP_OK)
1704  goto CLEANUP;
1705  if ((res = mp_int_sub(TEMP(1), TEMP(6), TEMP(1))) != MP_OK)
1706  goto CLEANUP;
1707  }
1708 
1709  s_qdiv(TEMP(0), 1);
1710  s_qdiv(TEMP(1), 1);
1711  }
1712 
1713  while (mp_int_is_even(TEMP(5)))
1714  {
1715  s_qdiv(TEMP(5), 1);
1716 
1717  if (mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3)))
1718  {
1719  if ((res = mp_int_add(TEMP(2), TEMP(7), TEMP(2))) != MP_OK)
1720  goto CLEANUP;
1721  if ((res = mp_int_sub(TEMP(3), TEMP(6), TEMP(3))) != MP_OK)
1722  goto CLEANUP;
1723  }
1724 
1725  s_qdiv(TEMP(2), 1);
1726  s_qdiv(TEMP(3), 1);
1727  }
1728 
1729  if (mp_int_compare(TEMP(4), TEMP(5)) >= 0)
1730  {
1731  if ((res = mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != MP_OK)
1732  goto CLEANUP;
1733  if ((res = mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != MP_OK)
1734  goto CLEANUP;
1735  if ((res = mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != MP_OK)
1736  goto CLEANUP;
1737  }
1738  else
1739  {
1740  if ((res = mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != MP_OK)
1741  goto CLEANUP;
1742  if ((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK)
1743  goto CLEANUP;
1744  if ((res = mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != MP_OK)
1745  goto CLEANUP;
1746  }
1747 
1748  if (CMPZ(TEMP(4)) == 0)
1749  {
1750  if (x && (res = mp_int_copy(TEMP(2), x)) != MP_OK)
1751  goto CLEANUP;
1752  if (y && (res = mp_int_copy(TEMP(3), y)) != MP_OK)
1753  goto CLEANUP;
1754  if (c)
1755  {
1756  if (!s_qmul(TEMP(5), k))
1757  {
1758  res = MP_MEMORY;
1759  goto CLEANUP;
1760  }
1761 
1762  res = mp_int_copy(TEMP(5), c);
1763  }
1764 
1765  break;
1766  }
1767  }
1768 
1769 CLEANUP:
1770  while (--last >= 0)
1771  mp_int_clear(TEMP(last));
1772 
1773  return res;
1774 }
1775 
1776 /* }}} */
1777 
1778 /* {{{ mp_int_divisible_value(a, v) */
1779 
1780 int
1782 {
1783  int rem = 0;
1784 
1785  if (mp_int_div_value(a, v, NULL, &rem) != MP_OK)
1786  return 0;
1787 
1788  return rem == 0;
1789 }
1790 
1791 /* }}} */
1792 
1793 /* {{{ mp_int_is_pow2(z) */
1794 
1795 int
1797 {
1798  CHECK(z != NULL);
1799 
1800  return s_isp2(z);
1801 }
1802 
1803 /* }}} */
1804 
1805 /* {{{ mp_int_sqrt(a, c) */
1806 
1807 mp_result
1809 {
1810  mp_result res = MP_OK;
1811  mpz_t temp[2];
1812  int last = 0;
1813 
1814  CHECK(a != NULL && c != NULL);
1815 
1816  /* The square root of a negative value does not exist in the integers. */
1817  if (MP_SIGN(a) == MP_NEG)
1818  return MP_UNDEF;
1819 
1820  SETUP(mp_int_init_copy(TEMP(last), a), last);
1821  SETUP(mp_int_init(TEMP(last)), last);
1822 
1823  for (;;)
1824  {
1825  if ((res = mp_int_sqr(TEMP(0), TEMP(1))) != MP_OK)
1826  goto CLEANUP;
1827 
1828  if (mp_int_compare_unsigned(a, TEMP(1)) == 0)
1829  break;
1830 
1831  if ((res = mp_int_copy(a, TEMP(1))) != MP_OK)
1832  goto CLEANUP;
1833  if ((res = mp_int_div(TEMP(1), TEMP(0), TEMP(1), NULL)) != MP_OK)
1834  goto CLEANUP;
1835  if ((res = mp_int_add(TEMP(0), TEMP(1), TEMP(1))) != MP_OK)
1836  goto CLEANUP;
1837  if ((res = mp_int_div_pow2(TEMP(1), 1, TEMP(1), NULL)) != MP_OK)
1838  goto CLEANUP;
1839 
1840  if (mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0)
1841  break;
1842  if ((res = mp_int_sub_value(TEMP(0), 1, TEMP(0))) != MP_OK)
1843  goto CLEANUP;
1844  if (mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0)
1845  break;
1846 
1847  if ((res = mp_int_copy(TEMP(1), TEMP(0))) != MP_OK)
1848  goto CLEANUP;
1849  }
1850 
1851  res = mp_int_copy(TEMP(0), c);
1852 
1853 CLEANUP:
1854  while (--last >= 0)
1855  mp_int_clear(TEMP(last));
1856 
1857  return res;
1858 }
1859 
1860 /* }}} */
1861 
1862 /* {{{ mp_int_to_int(z, out) */
1863 
1864 mp_result
1865 mp_int_to_int(mp_int z, int *out)
1866 {
1867  unsigned int uv = 0;
1868  mp_size uz;
1869  mp_digit *dz;
1870  mp_sign sz;
1871 
1872  CHECK(z != NULL);
1873 
1874  /* Make sure the value is representable as an int */
1875  sz = MP_SIGN(z);
1876  if ((sz == MP_ZPOS && mp_int_compare_value(z, INT_MAX) > 0) ||
1877  mp_int_compare_value(z, INT_MIN) < 0)
1878  return MP_RANGE;
1879 
1880  uz = MP_USED(z);
1881  dz = MP_DIGITS(z) + uz - 1;
1882 
1883  while (uz > 0)
1884  {
1885  uv <<= MP_DIGIT_BIT / 2;
1886  uv = (uv << (MP_DIGIT_BIT / 2)) | *dz--;
1887  --uz;
1888  }
1889 
1890  if (out)
1891  *out = (sz == MP_NEG) ? -(int) uv : (int) uv;
1892 
1893  return MP_OK;
1894 }
1895 
1896 /* }}} */
1897 
1898 /* {{{ mp_int_to_string(z, radix, str, limit) */
1899 
1900 mp_result
1902  char *str, int limit)
1903 {
1904  mp_result res;
1905  int cmp = 0;
1906 
1907  CHECK(z != NULL && str != NULL && limit >= 2);
1908 
1909  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1910  return MP_RANGE;
1911 
1912  if (CMPZ(z) == 0)
1913  {
1914  *str++ = s_val2ch(0, mp_flags & MP_CAP_DIGITS);
1915  }
1916  else
1917  {
1918  mpz_t tmp;
1919  char *h,
1920  *t;
1921 
1922  if ((res = mp_int_init_copy(&tmp, z)) != MP_OK)
1923  return res;
1924 
1925  if (MP_SIGN(z) == MP_NEG)
1926  {
1927  *str++ = '-';
1928  --limit;
1929  }
1930  h = str;
1931 
1932  /* Generate digits in reverse order until finished or limit reached */
1933  for ( /* */ ; limit > 0; --limit)
1934  {
1935  mp_digit d;
1936 
1937  if ((cmp = CMPZ(&tmp)) == 0)
1938  break;
1939 
1940  d = s_ddiv(&tmp, (mp_digit) radix);
1941  *str++ = s_val2ch(d, mp_flags & MP_CAP_DIGITS);
1942  }
1943  t = str - 1;
1944 
1945  /* Put digits back in correct output order */
1946  while (h < t)
1947  {
1948  char tc = *h;
1949 
1950  *h++ = *t;
1951  *t-- = tc;
1952  }
1953 
1954  mp_int_clear(&tmp);
1955  }
1956 
1957  *str = '\0';
1958  if (cmp == 0)
1959  return MP_OK;
1960  else
1961  return MP_TRUNC;
1962 }
1963 
1964 /* }}} */
1965 
1966 /* {{{ mp_int_string_len(z, radix) */
1967 
1968 mp_result
1970 {
1971  int len;
1972 
1973  CHECK(z != NULL);
1974 
1975  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1976  return MP_RANGE;
1977 
1978  len = s_outlen(z, radix) + 1; /* for terminator */
1979 
1980  /* Allow for sign marker on negatives */
1981  if (MP_SIGN(z) == MP_NEG)
1982  len += 1;
1983 
1984  return len;
1985 }
1986 
1987 /* }}} */
1988 
1989 /* {{{ mp_int_read_string(z, radix, *str) */
1990 
1991 /* Read zero-terminated string into z */
1992 mp_result
1993 mp_int_read_string(mp_int z, mp_size radix, const char *str)
1994 {
1995  return mp_int_read_cstring(z, radix, str, NULL);
1996 
1997 }
1998 
1999 /* }}} */
2000 
2001 /* {{{ mp_int_read_cstring(z, radix, *str, **end) */
2002 
2003 mp_result
2004 mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
2005 {
2006  int ch;
2007 
2008  CHECK(z != NULL && str != NULL);
2009 
2010  if (radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
2011  return MP_RANGE;
2012 
2013  /* Skip leading whitespace */
2014  while (isspace((unsigned char) *str))
2015  ++str;
2016 
2017  /* Handle leading sign tag (+/-, positive default) */
2018  switch (*str)
2019  {
2020  case '-':
2021  MP_SIGN(z) = MP_NEG;
2022  ++str;
2023  break;
2024  case '+':
2025  ++str; /* fallthrough */
2026  default:
2027  MP_SIGN(z) = MP_ZPOS;
2028  break;
2029  }
2030 
2031  /* Skip leading zeroes */
2032  while ((ch = s_ch2val(*str, radix)) == 0)
2033  ++str;
2034 
2035  /* Make sure there is enough space for the value */
2036  if (!s_pad(z, s_inlen(strlen(str), radix)))
2037  return MP_MEMORY;
2038 
2039  MP_USED(z) = 1;
2040  z->digits[0] = 0;
2041 
2042  while (*str != '\0' && ((ch = s_ch2val(*str, radix)) >= 0))
2043  {
2044  s_dmul(z, (mp_digit) radix);
2045  s_dadd(z, (mp_digit) ch);
2046  ++str;
2047  }
2048 
2049  CLAMP(z);
2050 
2051  /* Override sign for zero, even if negative specified. */
2052  if (CMPZ(z) == 0)
2053  MP_SIGN(z) = MP_ZPOS;
2054 
2055  if (end != NULL)
2056  *end = (char *) str;
2057 
2058  /*
2059  * Return a truncation error if the string has unprocessed characters
2060  * remaining, so the caller can tell if the whole string was done
2061  */
2062  if (*str != '\0')
2063  return MP_TRUNC;
2064  else
2065  return MP_OK;
2066 }
2067 
2068 /* }}} */
2069 
2070 /* {{{ mp_int_count_bits(z) */
2071 
2072 mp_result
2074 {
2075  mp_size nbits = 0,
2076  uz;
2077  mp_digit d;
2078 
2079  CHECK(z != NULL);
2080 
2081  uz = MP_USED(z);
2082  if (uz == 1 && z->digits[0] == 0)
2083  return 1;
2084 
2085  --uz;
2086  nbits = uz * MP_DIGIT_BIT;
2087  d = z->digits[uz];
2088 
2089  while (d != 0)
2090  {
2091  d >>= 1;
2092  ++nbits;
2093  }
2094 
2095  return nbits;
2096 }
2097 
2098 /* }}} */
2099 
2100 /* {{{ mp_int_to_binary(z, buf, limit) */
2101 
2102 mp_result
2103 mp_int_to_binary(mp_int z, unsigned char *buf, int limit)
2104 {
2105  static const int PAD_FOR_2C = 1;
2106 
2107  mp_result res;
2108  int limpos = limit;
2109 
2110  CHECK(z != NULL && buf != NULL);
2111 
2112  res = s_tobin(z, buf, &limpos, PAD_FOR_2C);
2113 
2114  if (MP_SIGN(z) == MP_NEG)
2115  s_2comp(buf, limpos);
2116 
2117  return res;
2118 }
2119 
2120 /* }}} */
2121 
2122 /* {{{ mp_int_read_binary(z, buf, len) */
2123 
2124 mp_result
2125 mp_int_read_binary(mp_int z, unsigned char *buf, int len)
2126 {
2127  mp_size need,
2128  i;
2129  unsigned char *tmp;
2130  mp_digit *dz;
2131 
2132  CHECK(z != NULL && buf != NULL && len > 0);
2133 
2134  /* Figure out how many digits are needed to represent this value */
2135  need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
2136  if (!s_pad(z, need))
2137  return MP_MEMORY;
2138 
2139  mp_int_zero(z);
2140 
2141  /*
2142  * If the high-order bit is set, take the 2's complement before reading
2143  * the value (it will be restored afterward)
2144  */
2145  if (buf[0] >> (CHAR_BIT - 1))
2146  {
2147  MP_SIGN(z) = MP_NEG;
2148  s_2comp(buf, len);
2149  }
2150 
2151  dz = MP_DIGITS(z);
2152  for (tmp = buf, i = len; i > 0; --i, ++tmp)
2153  {
2154  s_qmul(z, (mp_size) CHAR_BIT);
2155  *dz |= *tmp;
2156  }
2157 
2158  /* Restore 2's complement if we took it before */
2159  if (MP_SIGN(z) == MP_NEG)
2160  s_2comp(buf, len);
2161 
2162  return MP_OK;
2163 }
2164 
2165 /* }}} */
2166 
2167 /* {{{ mp_int_binary_len(z) */
2168 
2169 mp_result
2171 {
2172  mp_result res = mp_int_count_bits(z);
2173  int bytes = mp_int_unsigned_len(z);
2174 
2175  if (res <= 0)
2176  return res;
2177 
2178  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
2179 
2180  /*
2181  * If the highest-order bit falls exactly on a byte boundary, we need to
2182  * pad with an extra byte so that the sign will be read correctly when
2183  * reading it back in.
2184  */
2185  if (bytes * CHAR_BIT == res)
2186  ++bytes;
2187 
2188  return bytes;
2189 }
2190 
2191 /* }}} */
2192 
2193 /* {{{ mp_int_to_unsigned(z, buf, limit) */
2194 
2195 mp_result
2196 mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit)
2197 {
2198  static const int NO_PADDING = 0;
2199 
2200  CHECK(z != NULL && buf != NULL);
2201 
2202  return s_tobin(z, buf, &limit, NO_PADDING);
2203 }
2204 
2205 /* }}} */
2206 
2207 /* {{{ mp_int_read_unsigned(z, buf, len) */
2208 
2209 mp_result
2210 mp_int_read_unsigned(mp_int z, unsigned char *buf, int len)
2211 {
2212  mp_size need,
2213  i;
2214  unsigned char *tmp;
2215  mp_digit *dz;
2216 
2217  CHECK(z != NULL && buf != NULL && len > 0);
2218 
2219  /* Figure out how many digits are needed to represent this value */
2220  need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
2221  if (!s_pad(z, need))
2222  return MP_MEMORY;
2223 
2224  mp_int_zero(z);
2225 
2226  dz = MP_DIGITS(z);
2227  for (tmp = buf, i = len; i > 0; --i, ++tmp)
2228  {
2229  (void) s_qmul(z, CHAR_BIT);
2230  *dz |= *tmp;
2231  }
2232 
2233  return MP_OK;
2234 }
2235 
2236 /* }}} */
2237 
2238 /* {{{ mp_int_unsigned_len(z) */
2239 
2240 mp_result
2242 {
2243  mp_result res = mp_int_count_bits(z);
2244  int bytes;
2245 
2246  if (res <= 0)
2247  return res;
2248 
2249  bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
2250 
2251  return bytes;
2252 }
2253 
2254 /* }}} */
2255 
2256 /* {{{ mp_error_string(res) */
2257 
2258 const char *
2260 {
2261  int ix;
2262 
2263  if (res > 0)
2264  return s_unknown_err;
2265 
2266  res = -res;
2267  for (ix = 0; ix < res && s_error_msg[ix] != NULL; ++ix)
2268  ;
2269 
2270  if (s_error_msg[ix] != NULL)
2271  return s_error_msg[ix];
2272  else
2273  return s_unknown_err;
2274 }
2275 
2276 /* }}} */
2277 
2278 /*------------------------------------------------------------------------*/
2279 /* Private functions for internal use. These make assumptions. */
2280 
2281 /* {{{ s_alloc(num) */
2282 
2283 static mp_digit *
2285 {
2286  mp_digit *out = px_alloc(num * sizeof(mp_digit));
2287 
2288  assert(out != NULL); /* for debugging */
2289 
2290  return out;
2291 }
2292 
2293 /* }}} */
2294 
2295 /* {{{ s_realloc(old, num) */
2296 
2297 static mp_digit *
2299 {
2300  mp_digit *new = px_realloc(old, num * sizeof(mp_digit));
2301 
2302  assert(new != NULL); /* for debugging */
2303 
2304  return new;
2305 }
2306 
2307 /* }}} */
2308 
2309 /* {{{ s_free(ptr) */
2310 
2311 #if TRACEABLE_FREE
2312 static void
2313 s_free(void *ptr)
2314 {
2315  px_free(ptr);
2316 }
2317 #endif
2318 
2319 /* }}} */
2320 
2321 /* {{{ s_pad(z, min) */
2322 
2323 static int
2325 {
2326  if (MP_ALLOC(z) < min)
2327  {
2328  mp_size nsize = ROUND_PREC(min);
2329  mp_digit *tmp = s_realloc(MP_DIGITS(z), nsize);
2330 
2331  if (tmp == NULL)
2332  return 0;
2333 
2334  MP_DIGITS(z) = tmp;
2335  MP_ALLOC(z) = nsize;
2336  }
2337 
2338  return 1;
2339 }
2340 
2341 /* }}} */
2342 
2343 /* {{{ s_clamp(z) */
2344 
2345 #if TRACEABLE_CLAMP
2346 static void
2347 s_clamp(mp_int z)
2348 {
2349  mp_size uz = MP_USED(z);
2350  mp_digit *zd = MP_DIGITS(z) + uz - 1;
2351 
2352  while (uz > 1 && (*zd-- == 0))
2353  --uz;
2354 
2355  MP_USED(z) = uz;
2356 }
2357 #endif
2358 
2359 /* }}} */
2360 
2361 /* {{{ s_fake(z, value, vbuf) */
2362 
2363 static void
2364 s_fake(mp_int z, int value, mp_digit vbuf[])
2365 {
2366  mp_size uv = (mp_size) s_vpack(value, vbuf);
2367 
2368  z->used = uv;
2369  z->alloc = MP_VALUE_DIGITS(value);
2370  z->sign = (value < 0) ? MP_NEG : MP_ZPOS;
2371  z->digits = vbuf;
2372 }
2373 
2374 /* }}} */
2375 
2376 /* {{{ s_cdig(da, db, len) */
2377 
2378 static int
2380 {
2381  mp_digit *dat = da + len - 1,
2382  *dbt = db + len - 1;
2383 
2384  for ( /* */ ; len != 0; --len, --dat, --dbt)
2385  {
2386  if (*dat > *dbt)
2387  return 1;
2388  else if (*dat < *dbt)
2389  return -1;
2390  }
2391 
2392  return 0;
2393 }
2394 
2395 /* }}} */
2396 
2397 /* {{{ s_vpack(v, t[]) */
2398 
2399 static int
2400 s_vpack(int v, mp_digit t[])
2401 {
2402  unsigned int uv = (unsigned int) ((v < 0) ? -v : v);
2403  int ndig = 0;
2404 
2405  if (uv == 0)
2406  t[ndig++] = 0;
2407  else
2408  {
2409  while (uv != 0)
2410  {
2411  t[ndig++] = (mp_digit) uv;
2412  uv >>= MP_DIGIT_BIT / 2;
2413  uv >>= MP_DIGIT_BIT / 2;
2414  }
2415  }
2416 
2417  return ndig;
2418 }
2419 
2420 /* }}} */
2421 
2422 /* {{{ s_ucmp(a, b) */
2423 
2424 static int
2426 {
2427  mp_size ua = MP_USED(a),
2428  ub = MP_USED(b);
2429 
2430  if (ua > ub)
2431  return 1;
2432  else if (ub > ua)
2433  return -1;
2434  else
2435  return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua);
2436 }
2437 
2438 /* }}} */
2439 
2440 /* {{{ s_vcmp(a, v) */
2441 
2442 static int
2443 s_vcmp(mp_int a, int v)
2444 {
2445  mp_digit vdig[MP_VALUE_DIGITS(v)];
2446  int ndig = 0;
2447  mp_size ua = MP_USED(a);
2448 
2449  ndig = s_vpack(v, vdig);
2450 
2451  if (ua > ndig)
2452  return 1;
2453  else if (ua < ndig)
2454  return -1;
2455  else
2456  return s_cdig(MP_DIGITS(a), vdig, ndig);
2457 }
2458 
2459 /* }}} */
2460 
2461 /* {{{ s_uadd(da, db, dc, size_a, size_b) */
2462 
2463 static mp_digit
2465  mp_size size_a, mp_size size_b)
2466 {
2467  mp_size pos;
2468  mp_word w = 0;
2469 
2470  /* Insure that da is the longer of the two to simplify later code */
2471  if (size_b > size_a)
2472  {
2473  SWAP(mp_digit *, da, db);
2474  SWAP(mp_size, size_a, size_b);
2475  }
2476 
2477  /* Add corresponding digits until the shorter number runs out */
2478  for (pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc)
2479  {
2480  w = w + (mp_word) *da + (mp_word) *db;
2481  *dc = LOWER_HALF(w);
2482  w = UPPER_HALF(w);
2483  }
2484 
2485  /* Propagate carries as far as necessary */
2486  for ( /* */ ; pos < size_a; ++pos, ++da, ++dc)
2487  {
2488  w = w + *da;
2489 
2490  *dc = LOWER_HALF(w);
2491  w = UPPER_HALF(w);
2492  }
2493 
2494  /* Return carry out */
2495  return (mp_digit) w;
2496 }
2497 
2498 /* }}} */
2499 
2500 /* {{{ s_usub(da, db, dc, size_a, size_b) */
2501 
2502 static void
2504  mp_size size_a, mp_size size_b)
2505 {
2506  mp_size pos;
2507  mp_word w = 0;
2508 
2509  /* We assume that |a| >= |b| so this should definitely hold */
2510  assert(size_a >= size_b);
2511 
2512  /* Subtract corresponding digits and propagate borrow */
2513  for (pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc)
2514  {
2515  w = ((mp_word) MP_DIGIT_MAX + 1 + /* MP_RADIX */
2516  (mp_word) *da) - w - (mp_word) *db;
2517 
2518  *dc = LOWER_HALF(w);
2519  w = (UPPER_HALF(w) == 0);
2520  }
2521 
2522  /* Finish the subtraction for remaining upper digits of da */
2523  for ( /* */ ; pos < size_a; ++pos, ++da, ++dc)
2524  {
2525  w = ((mp_word) MP_DIGIT_MAX + 1 + /* MP_RADIX */
2526  (mp_word) *da) - w;
2527 
2528  *dc = LOWER_HALF(w);
2529  w = (UPPER_HALF(w) == 0);
2530  }
2531 
2532  /* If there is a borrow out at the end, it violates the precondition */
2533  assert(w == 0);
2534 }
2535 
2536 /* }}} */
2537 
2538 /* {{{ s_kmul(da, db, dc, size_a, size_b) */
2539 
2540 static int
2542  mp_size size_a, mp_size size_b)
2543 {
2544  mp_size bot_size;
2545 
2546  /* Make sure b is the smaller of the two input values */
2547  if (size_b > size_a)
2548  {
2549  SWAP(mp_digit *, da, db);
2550  SWAP(mp_size, size_a, size_b);
2551  }
2552 
2553  /*
2554  * Insure that the bottom is the larger half in an odd-length split; the
2555  * code below relies on this being true.
2556  */
2557  bot_size = (size_a + 1) / 2;
2558 
2559  /*
2560  * If the values are big enough to bother with recursion, use the
2561  * Karatsuba algorithm to compute the product; otherwise use the normal
2562  * multiplication algorithm
2563  */
2564  if (multiply_threshold &&
2565  size_a >= multiply_threshold &&
2566  size_b > bot_size)
2567  {
2568 
2569  mp_digit *t1,
2570  *t2,
2571  *t3,
2572  carry;
2573 
2574  mp_digit *a_top = da + bot_size;
2575  mp_digit *b_top = db + bot_size;
2576 
2577  mp_size at_size = size_a - bot_size;
2578  mp_size bt_size = size_b - bot_size;
2579  mp_size buf_size = 2 * bot_size;
2580 
2581  /*
2582  * Do a single allocation for all three temporary buffers needed; each
2583  * buffer must be big enough to hold the product of two bottom halves,
2584  * and one buffer needs space for the completed product; twice the
2585  * space is plenty.
2586  */
2587  if ((t1 = s_alloc(4 * buf_size)) == NULL)
2588  return 0;
2589  t2 = t1 + buf_size;
2590  t3 = t2 + buf_size;
2591  ZERO(t1, 4 * buf_size);
2592 
2593  /*
2594  * t1 and t2 are initially used as temporaries to compute the inner
2595  * product (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
2596  */
2597  carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */
2598  t1[bot_size] = carry;
2599 
2600  carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */
2601  t2[bot_size] = carry;
2602 
2603  (void) s_kmul(t1, t2, t3, bot_size + 1, bot_size + 1); /* t3 = t1 * t2 */
2604 
2605  /*
2606  * Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so
2607  * that we're left with only the pieces we want: t3 = a1b0 + a0b1
2608  */
2609  ZERO(t1, buf_size);
2610  ZERO(t2, buf_size);
2611  (void) s_kmul(da, db, t1, bot_size, bot_size); /* t1 = a0 * b0 */
2612  (void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */
2613 
2614  /* Subtract out t1 and t2 to get the inner product */
2615  s_usub(t3, t1, t3, buf_size + 2, buf_size);
2616  s_usub(t3, t2, t3, buf_size + 2, buf_size);
2617 
2618  /* Assemble the output value */
2619  COPY(t1, dc, buf_size);
2620  carry = s_uadd(t3, dc + bot_size, dc + bot_size,
2621  buf_size + 1, buf_size);
2622  assert(carry == 0);
2623 
2624  carry = s_uadd(t2, dc + 2 * bot_size, dc + 2 * bot_size,
2625  buf_size, buf_size);
2626  assert(carry == 0);
2627 
2628  s_free(t1); /* note t2 and t3 are just internal pointers
2629  * to t1 */
2630  }
2631  else
2632  {
2633  s_umul(da, db, dc, size_a, size_b);
2634  }
2635 
2636  return 1;
2637 }
2638 
2639 /* }}} */
2640 
2641 /* {{{ s_umul(da, db, dc, size_a, size_b) */
2642 
2643 static void
2645  mp_size size_a, mp_size size_b)
2646 {
2647  mp_size a,
2648  b;
2649  mp_word w;
2650 
2651  for (a = 0; a < size_a; ++a, ++dc, ++da)
2652  {
2653  mp_digit *dct = dc;
2654  mp_digit *dbt = db;
2655 
2656  if (*da == 0)
2657  continue;
2658 
2659  w = 0;
2660  for (b = 0; b < size_b; ++b, ++dbt, ++dct)
2661  {
2662  w = (mp_word) *da * (mp_word) *dbt + w + (mp_word) *dct;
2663 
2664  *dct = LOWER_HALF(w);
2665  w = UPPER_HALF(w);
2666  }
2667 
2668  *dct = (mp_digit) w;
2669  }
2670 }
2671 
2672 /* }}} */
2673 
2674 /* {{{ s_ksqr(da, dc, size_a) */
2675 
2676 static int
2677 s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2678 {
2679  if (multiply_threshold && size_a > multiply_threshold)
2680  {
2681  mp_size bot_size = (size_a + 1) / 2;
2682  mp_digit *a_top = da + bot_size;
2683  mp_digit *t1,
2684  *t2,
2685  *t3;
2686  mp_size at_size = size_a - bot_size;
2687  mp_size buf_size = 2 * bot_size;
2688 
2689  if ((t1 = s_alloc(4 * buf_size)) == NULL)
2690  return 0;
2691  t2 = t1 + buf_size;
2692  t3 = t2 + buf_size;
2693  ZERO(t1, 4 * buf_size);
2694 
2695  (void) s_ksqr(da, t1, bot_size); /* t1 = a0 ^ 2 */
2696  (void) s_ksqr(a_top, t2, at_size); /* t2 = a1 ^ 2 */
2697 
2698  (void) s_kmul(da, a_top, t3, bot_size, at_size); /* t3 = a0 * a1 */
2699 
2700  /* Quick multiply t3 by 2, shifting left (can't overflow) */
2701  {
2702  int i,
2703  top = bot_size + at_size;
2704  mp_word w,
2705  save = 0;
2706 
2707  for (i = 0; i < top; ++i)
2708  {
2709  w = t3[i];
2710  w = (w << 1) | save;
2711  t3[i] = LOWER_HALF(w);
2712  save = UPPER_HALF(w);
2713  }
2714  t3[i] = LOWER_HALF(save);
2715  }
2716 
2717  /* Assemble the output value */
2718  COPY(t1, dc, 2 * bot_size);
2719  (void) s_uadd(t3, dc + bot_size, dc + bot_size,
2720  buf_size + 1, buf_size + 1);
2721 
2722  (void) s_uadd(t2, dc + 2 * bot_size, dc + 2 * bot_size,
2723  buf_size, buf_size);
2724 
2725  px_free(t1); /* note that t2 and t2 are internal pointers
2726  * only */
2727 
2728  }
2729  else
2730  {
2731  s_usqr(da, dc, size_a);
2732  }
2733 
2734  return 1;
2735 }
2736 
2737 /* }}} */
2738 
2739 /* {{{ s_usqr(da, dc, size_a) */
2740 
2741 static void
2742 s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2743 {
2744  mp_size i,
2745  j;
2746  mp_word w;
2747 
2748  for (i = 0; i < size_a; ++i, dc += 2, ++da)
2749  {
2750  mp_digit *dct = dc,
2751  *dat = da;
2752 
2753  if (*da == 0)
2754  continue;
2755 
2756  /* Take care of the first digit, no rollover */
2757  w = (mp_word) *dat * (mp_word) *dat + (mp_word) *dct;
2758  *dct = LOWER_HALF(w);
2759  w = UPPER_HALF(w);
2760  ++dat;
2761  ++dct;
2762 
2763  for (j = i + 1; j < size_a; ++j, ++dat, ++dct)
2764  {
2765  mp_word t = (mp_word) *da * (mp_word) *dat;
2766  mp_word u = w + (mp_word) *dct,
2767  ov = 0;
2768 
2769  /* Check if doubling t will overflow a word */
2770  if (HIGH_BIT_SET(t))
2771  ov = 1;
2772 
2773  w = t + t;
2774 
2775  /* Check if adding u to w will overflow a word */
2776  if (ADD_WILL_OVERFLOW(w, u))
2777  ov = 1;
2778 
2779  w += u;
2780 
2781  *dct = LOWER_HALF(w);
2782  w = UPPER_HALF(w);
2783  if (ov)
2784  {
2785  w += MP_DIGIT_MAX; /* MP_RADIX */
2786  ++w;
2787  }
2788  }
2789 
2790  w = w + *dct;
2791  *dct = (mp_digit) w;
2792  while ((w = UPPER_HALF(w)) != 0)
2793  {
2794  ++dct;
2795  w = w + *dct;
2796  *dct = LOWER_HALF(w);
2797  }
2798 
2799  assert(w == 0);
2800  }
2801 }
2802 
2803 /* }}} */
2804 
2805 /* {{{ s_dadd(a, b) */
2806 
2807 static void
2809 {
2810  mp_word w = 0;
2811  mp_digit *da = MP_DIGITS(a);
2812  mp_size ua = MP_USED(a);
2813 
2814  w = (mp_word) *da + b;
2815  *da++ = LOWER_HALF(w);
2816  w = UPPER_HALF(w);
2817 
2818  for (ua -= 1; ua > 0; --ua, ++da)
2819  {
2820  w = (mp_word) *da + w;
2821 
2822  *da = LOWER_HALF(w);
2823  w = UPPER_HALF(w);
2824  }
2825 
2826  if (w)
2827  {
2828  *da = (mp_digit) w;
2829  MP_USED(a) += 1;
2830  }
2831 }
2832 
2833 /* }}} */
2834 
2835 /* {{{ s_dmul(a, b) */
2836 
2837 static void
2839 {
2840  mp_word w = 0;
2841  mp_digit *da = MP_DIGITS(a);
2842  mp_size ua = MP_USED(a);
2843 
2844  while (ua > 0)
2845  {
2846  w = (mp_word) *da * b + w;
2847  *da++ = LOWER_HALF(w);
2848  w = UPPER_HALF(w);
2849  --ua;
2850  }
2851 
2852  if (w)
2853  {
2854  *da = (mp_digit) w;
2855  MP_USED(a) += 1;
2856  }
2857 }
2858 
2859 /* }}} */
2860 
2861 /* {{{ s_dbmul(da, b, dc, size_a) */
2862 
2863 static void
2865 {
2866  mp_word w = 0;
2867 
2868  while (size_a > 0)
2869  {
2870  w = (mp_word) *da++ * (mp_word) b + w;
2871 
2872  *dc++ = LOWER_HALF(w);
2873  w = UPPER_HALF(w);
2874  --size_a;
2875  }
2876 
2877  if (w)
2878  *dc = LOWER_HALF(w);
2879 }
2880 
2881 /* }}} */
2882 
2883 /* {{{ s_ddiv(da, d, dc, size_a) */
2884 
2885 static mp_digit
2887 {
2888  mp_word w = 0,
2889  qdigit;
2890  mp_size ua = MP_USED(a);
2891  mp_digit *da = MP_DIGITS(a) + ua - 1;
2892 
2893  for ( /* */ ; ua > 0; --ua, --da)
2894  {
2895  w = (w << MP_DIGIT_BIT) | *da;
2896 
2897  if (w >= b)
2898  {
2899  qdigit = w / b;
2900  w = w % b;
2901  }
2902  else
2903  {
2904  qdigit = 0;
2905  }
2906 
2907  *da = (mp_digit) qdigit;
2908  }
2909 
2910  CLAMP(a);
2911  return (mp_digit) w;
2912 }
2913 
2914 /* }}} */
2915 
2916 /* {{{ s_qdiv(z, p2) */
2917 
2918 static void
2920 {
2921  mp_size ndig = p2 / MP_DIGIT_BIT,
2922  nbits = p2 % MP_DIGIT_BIT;
2923  mp_size uz = MP_USED(z);
2924 
2925  if (ndig)
2926  {
2927  mp_size mark;
2928  mp_digit *to,
2929  *from;
2930 
2931  if (ndig >= uz)
2932  {
2933  mp_int_zero(z);
2934  return;
2935  }
2936 
2937  to = MP_DIGITS(z);
2938  from = to + ndig;
2939 
2940  for (mark = ndig; mark < uz; ++mark)
2941  *to++ = *from++;
2942 
2943  MP_USED(z) = uz - ndig;
2944  }
2945 
2946  if (nbits)
2947  {
2948  mp_digit d = 0,
2949  *dz,
2950  save;
2951  mp_size up = MP_DIGIT_BIT - nbits;
2952 
2953  uz = MP_USED(z);
2954  dz = MP_DIGITS(z) + uz - 1;
2955 
2956  for ( /* */ ; uz > 0; --uz, --dz)
2957  {
2958  save = *dz;
2959 
2960  *dz = (*dz >> nbits) | (d << up);
2961  d = save;
2962  }
2963 
2964  CLAMP(z);
2965  }
2966 
2967  if (MP_USED(z) == 1 && z->digits[0] == 0)
2968  MP_SIGN(z) = MP_ZPOS;
2969 }
2970 
2971 /* }}} */
2972 
2973 /* {{{ s_qmod(z, p2) */
2974 
2975 static void
2977 {
2978  mp_size start = p2 / MP_DIGIT_BIT + 1,
2979  rest = p2 % MP_DIGIT_BIT;
2980  mp_size uz = MP_USED(z);
2981  mp_digit mask = (1 << rest) - 1;
2982 
2983  if (start <= uz)
2984  {
2985  MP_USED(z) = start;
2986  z->digits[start - 1] &= mask;
2987  CLAMP(z);
2988  }
2989 }
2990 
2991 /* }}} */
2992 
2993 /* {{{ s_qmul(z, p2) */
2994 
2995 static int
2997 {
2998  mp_size uz,
2999  need,
3000  rest,
3001  extra,
3002  i;
3003  mp_digit *from,
3004  *to,
3005  d;
3006 
3007  if (p2 == 0)
3008  return 1;
3009 
3010  uz = MP_USED(z);
3011  need = p2 / MP_DIGIT_BIT;
3012  rest = p2 % MP_DIGIT_BIT;
3013 
3014  /*
3015  * Figure out if we need an extra digit at the top end; this occurs if the
3016  * topmost `rest' bits of the high-order digit of z are not zero, meaning
3017  * they will be shifted off the end if not preserved
3018  */
3019  extra = 0;
3020  if (rest != 0)
3021  {
3022  mp_digit *dz = MP_DIGITS(z) + uz - 1;
3023 
3024  if ((*dz >> (MP_DIGIT_BIT - rest)) != 0)
3025  extra = 1;
3026  }
3027 
3028  if (!s_pad(z, uz + need + extra))
3029  return 0;
3030 
3031  /*
3032  * If we need to shift by whole digits, do that in one pass, then to back
3033  * and shift by partial digits.
3034  */
3035  if (need > 0)
3036  {
3037  from = MP_DIGITS(z) + uz - 1;
3038  to = from + need;
3039 
3040  for (i = 0; i < uz; ++i)
3041  *to-- = *from--;
3042 
3043  ZERO(MP_DIGITS(z), need);
3044  uz += need;
3045  }
3046 
3047  if (rest)
3048  {
3049  d = 0;
3050  for (i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from)
3051  {
3052  mp_digit save = *from;
3053 
3054  *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
3055  d = save;
3056  }
3057 
3058  d >>= (MP_DIGIT_BIT - rest);
3059  if (d != 0)
3060  {
3061  *from = d;
3062  uz += extra;
3063  }
3064  }
3065 
3066  MP_USED(z) = uz;
3067  CLAMP(z);
3068 
3069  return 1;
3070 }
3071 
3072 /* }}} */
3073 
3074 /* {{{ s_qsub(z, p2) */
3075 
3076 /* Subtract |z| from 2^p2, assuming 2^p2 > |z|, and set z to be positive */
3077 static int
3079 {
3080  mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)),
3081  *zp;
3082  mp_size tdig = (p2 / MP_DIGIT_BIT),
3083  pos;
3084  mp_word w = 0;
3085 
3086  if (!s_pad(z, tdig + 1))
3087  return 0;
3088 
3089  for (pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp)
3090  {
3091  w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word) *zp;
3092 
3093  *zp = LOWER_HALF(w);
3094  w = UPPER_HALF(w) ? 0 : 1;
3095  }
3096 
3097  w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word) *zp;
3098  *zp = LOWER_HALF(w);
3099 
3100  assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
3101 
3102  MP_SIGN(z) = MP_ZPOS;
3103  CLAMP(z);
3104 
3105  return 1;
3106 }
3107 
3108 /* }}} */
3109 
3110 /* {{{ s_dp2k(z) */
3111 
3112 static int
3114 {
3115  int k = 0;
3116  mp_digit *dp = MP_DIGITS(z),
3117  d;
3118 
3119  if (MP_USED(z) == 1 && *dp == 0)
3120  return 1;
3121 
3122  while (*dp == 0)
3123  {
3124  k += MP_DIGIT_BIT;
3125  ++dp;
3126  }
3127 
3128  d = *dp;
3129  while ((d & 1) == 0)
3130  {
3131  d >>= 1;
3132  ++k;
3133  }
3134 
3135  return k;
3136 }
3137 
3138 /* }}} */
3139 
3140 /* {{{ s_isp2(z) */
3141 
3142 static int
3144 {
3145  mp_size uz = MP_USED(z),
3146  k = 0;
3147  mp_digit *dz = MP_DIGITS(z),
3148  d;
3149 
3150  while (uz > 1)
3151  {
3152  if (*dz++ != 0)
3153  return -1;
3154  k += MP_DIGIT_BIT;
3155  --uz;
3156  }
3157 
3158  d = *dz;
3159  while (d > 1)
3160  {
3161  if (d & 1)
3162  return -1;
3163  ++k;
3164  d >>= 1;
3165  }
3166 
3167  return (int) k;
3168 }
3169 
3170 /* }}} */
3171 
3172 /* {{{ s_2expt(z, k) */
3173 
3174 static int
3175 s_2expt(mp_int z, int k)
3176 {
3177  mp_size ndig,
3178  rest;
3179  mp_digit *dz;
3180 
3181  ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT;
3182  rest = k % MP_DIGIT_BIT;
3183 
3184  if (!s_pad(z, ndig))
3185  return 0;
3186 
3187  dz = MP_DIGITS(z);
3188  ZERO(dz, ndig);
3189  *(dz + ndig - 1) = (1 << rest);
3190  MP_USED(z) = ndig;
3191 
3192  return 1;
3193 }
3194 
3195 /* }}} */
3196 
3197 /* {{{ s_norm(a, b) */
3198 
3199 static int
3201 {
3202  mp_digit d = b->digits[MP_USED(b) - 1];
3203  int k = 0;
3204 
3205  while (d < (mp_digit) ((mp_digit) 1 << (MP_DIGIT_BIT - 1)))
3206  { /* d < (MP_RADIX / 2) */
3207  d <<= 1;
3208  ++k;
3209  }
3210 
3211  /* These multiplications can't fail */
3212  if (k != 0)
3213  {
3214  (void) s_qmul(a, (mp_size) k);
3215  (void) s_qmul(b, (mp_size) k);
3216  }
3217 
3218  return k;
3219 }
3220 
3221 /* }}} */
3222 
3223 /* {{{ s_brmu(z, m) */
3224 
3225 static mp_result
3227 {
3228  mp_size um = MP_USED(m) * 2;
3229 
3230  if (!s_pad(z, um))
3231  return MP_MEMORY;
3232 
3233  s_2expt(z, MP_DIGIT_BIT * um);
3234  return mp_int_div(z, m, z, NULL);
3235 }
3236 
3237 /* }}} */
3238 
3239 /* {{{ s_reduce(x, m, mu, q1, q2) */
3240 
3241 static int
3243 {
3244  mp_size um = MP_USED(m),
3245  umb_p1,
3246  umb_m1;
3247 
3248  umb_p1 = (um + 1) * MP_DIGIT_BIT;
3249  umb_m1 = (um - 1) * MP_DIGIT_BIT;
3250 
3251  if (mp_int_copy(x, q1) != MP_OK)
3252  return 0;
3253 
3254  /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
3255  s_qdiv(q1, umb_m1);
3256  UMUL(q1, mu, q2);
3257  s_qdiv(q2, umb_p1);
3258 
3259  /* Set x = x mod b^(k+1) */
3260  s_qmod(x, umb_p1);
3261 
3262  /*
3263  * Now, q is a guess for the quotient a / m. Compute x - q * m mod
3264  * b^(k+1), replacing x. This may be off by a factor of 2m, but no more
3265  * than that.
3266  */
3267  UMUL(q2, m, q1);
3268  s_qmod(q1, umb_p1);
3269  (void) mp_int_sub(x, q1, x); /* can't fail */
3270 
3271  /*
3272  * The result may be < 0; if it is, add b^(k+1) to pin it in the proper
3273  * range.
3274  */
3275  if ((CMPZ(x) < 0) && !s_qsub(x, umb_p1))
3276  return 0;
3277 
3278  /*
3279  * If x > m, we need to back it off until it is in range. This will be
3280  * required at most twice.
3281  */
3282  if (mp_int_compare(x, m) >= 0)
3283  (void) mp_int_sub(x, m, x);
3284  if (mp_int_compare(x, m) >= 0)
3285  (void) mp_int_sub(x, m, x);
3286 
3287  /* At this point, x has been properly reduced. */
3288  return 1;
3289 }
3290 
3291 /* }}} */
3292 
3293 /* {{{ s_embar(a, b, m, mu, c) */
3294 
3295 /* Perform modular exponentiation using Barrett's method, where mu is
3296  the reduction constant for m. Assumes a < m, b > 0. */
3297 static mp_result
3299 {
3300  mp_digit *db,
3301  *dbt,
3302  umu,
3303  d;
3304  mpz_t temp[3];
3305  mp_result res;
3306  int last = 0;
3307 
3308  umu = MP_USED(mu);
3309  db = MP_DIGITS(b);
3310  dbt = db + MP_USED(b) - 1;
3311 
3312  while (last < 3)
3313  {
3314  SETUP(mp_int_init_size(TEMP(last), 4 * umu), last);
3315  ZERO(MP_DIGITS(TEMP(last - 1)), MP_ALLOC(TEMP(last - 1)));
3316  }
3317 
3318  (void) mp_int_set_value(c, 1);
3319 
3320  /* Take care of low-order digits */
3321  while (db < dbt)
3322  {
3323  int i;
3324 
3325  for (d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1)
3326  {
3327  if (d & 1)
3328  {
3329  /* The use of a second temporary avoids allocation */
3330  UMUL(c, a, TEMP(0));
3331  if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2)))
3332  {
3333  res = MP_MEMORY;
3334  goto CLEANUP;
3335  }
3336  mp_int_copy(TEMP(0), c);
3337  }
3338 
3339 
3340  USQR(a, TEMP(0));
3341  assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
3342  if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2)))
3343  {
3344  res = MP_MEMORY;
3345  goto CLEANUP;
3346  }
3347  assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
3348  mp_int_copy(TEMP(0), a);
3349 
3350 
3351  }
3352 
3353  ++db;
3354  }
3355 
3356  /* Take care of highest-order digit */
3357  d = *dbt;
3358  for (;;)
3359  {
3360  if (d & 1)
3361  {
3362  UMUL(c, a, TEMP(0));
3363  if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2)))
3364  {
3365  res = MP_MEMORY;
3366  goto CLEANUP;
3367  }
3368  mp_int_copy(TEMP(0), c);
3369  }
3370 
3371  d >>= 1;
3372  if (!d)
3373  break;
3374 
3375  USQR(a, TEMP(0));
3376  if (!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2)))
3377  {
3378  res = MP_MEMORY;
3379  goto CLEANUP;
3380  }
3381  (void) mp_int_copy(TEMP(0), a);
3382  }
3383 
3384 CLEANUP:
3385  while (--last >= 0)
3386  mp_int_clear(TEMP(last));
3387 
3388  return res;
3389 }
3390 
3391 /* }}} */
3392 
3393 /* {{{ s_udiv(a, b) */
3394 
3395 /* Precondition: a >= b and b > 0
3396  Postcondition: a' = a / b, b' = a % b
3397  */
3398 static mp_result
3400 {
3401  mpz_t q,
3402  r,
3403  t;
3404  mp_size ua,
3405  ub,
3406  qpos = 0;
3407  mp_digit *da,
3408  btop;
3409  mp_result res = MP_OK;
3410  int k,
3411  skip = 0;
3412 
3413  /* Force signs to positive */
3414  MP_SIGN(a) = MP_ZPOS;
3415  MP_SIGN(b) = MP_ZPOS;
3416 
3417  /* Normalize, per Knuth */
3418  k = s_norm(a, b);
3419 
3420  ua = MP_USED(a);
3421  ub = MP_USED(b);
3422  btop = b->digits[ub - 1];
3423  if ((res = mp_int_init_size(&q, ua)) != MP_OK)
3424  return res;
3425  if ((res = mp_int_init_size(&t, ua + 1)) != MP_OK)
3426  goto CLEANUP;
3427 
3428  da = MP_DIGITS(a);
3429  r.digits = da + ua - 1; /* The contents of r are shared with a */
3430  r.used = 1;
3431  r.sign = MP_ZPOS;
3432  r.alloc = MP_ALLOC(a);
3433  ZERO(t.digits, t.alloc);
3434 
3435  /* Solve for quotient digits, store in q.digits in reverse order */
3436  while (r.digits >= da)
3437  {
3438  assert(qpos <= q.alloc);
3439 
3440  if (s_ucmp(b, &r) > 0)
3441  {
3442  r.digits -= 1;
3443  r.used += 1;
3444 
3445  if (++skip > 1)
3446  q.digits[qpos++] = 0;
3447 
3448  CLAMP(&r);
3449  }
3450  else
3451  {
3452  mp_word pfx = r.digits[r.used - 1];
3453  mp_word qdigit;
3454 
3455  if (r.used > 1 && (pfx < btop || r.digits[r.used - 2] == 0))
3456  {
3457  pfx <<= MP_DIGIT_BIT / 2;
3458  pfx <<= MP_DIGIT_BIT / 2;
3459  pfx |= r.digits[r.used - 2];
3460  }
3461 
3462  qdigit = pfx / btop;
3463  if (qdigit > MP_DIGIT_MAX)
3464  qdigit = 1;
3465 
3466  s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub);
3467  t.used = ub + 1;
3468  CLAMP(&t);
3469  while (s_ucmp(&t, &r) > 0)
3470  {
3471  --qdigit;
3472  (void) mp_int_sub(&t, b, &t); /* cannot fail */
3473  }
3474 
3475  s_usub(r.digits, t.digits, r.digits, r.used, t.used);
3476  CLAMP(&r);
3477 
3478  q.digits[qpos++] = (mp_digit) qdigit;
3479  ZERO(t.digits, t.used);
3480  skip = 0;
3481  }
3482  }
3483 
3484  /* Put quotient digits in the correct order, and discard extra zeroes */
3485  q.used = qpos;
3486  REV(mp_digit, q.digits, qpos);
3487  CLAMP(&q);
3488 
3489  /* Denormalize the remainder */
3490  CLAMP(a);
3491  if (k != 0)
3492  s_qdiv(a, k);
3493 
3494  mp_int_copy(a, b); /* ok: 0 <= r < b */
3495  mp_int_copy(&q, a); /* ok: q <= a */
3496 
3497  mp_int_clear(&t);
3498 CLEANUP:
3499  mp_int_clear(&q);
3500  return res;
3501 }
3502 
3503 /* }}} */
3504 
3505 /* {{{ s_outlen(z, r) */
3506 
3507 /* Precondition: 2 <= r < 64 */
3508 static int
3510 {
3511  mp_result bits;
3512  double raw;
3513 
3514  bits = mp_int_count_bits(z);
3515  raw = (double) bits *s_log2[r];
3516 
3517  return (int) (raw + 0.999999);
3518 }
3519 
3520 /* }}} */
3521 
3522 /* {{{ s_inlen(len, r) */
3523 
3524 static mp_size
3525 s_inlen(int len, mp_size r)
3526 {
3527  double raw = (double) len / s_log2[r];
3528  mp_size bits = (mp_size) (raw + 0.5);
3529 
3530  return (mp_size) ((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT);
3531 }
3532 
3533 /* }}} */
3534 
3535 /* {{{ s_ch2val(c, r) */
3536 
3537 static int
3538 s_ch2val(char c, int r)
3539 {
3540  int out;
3541 
3542  if (isdigit((unsigned char) c))
3543  out = c - '0';
3544  else if (r > 10 && isalpha((unsigned char) c))
3545  out = toupper((unsigned char) c) - 'A' + 10;
3546  else
3547  return -1;
3548 
3549  return (out >= r) ? -1 : out;
3550 }
3551 
3552 /* }}} */
3553 
3554 /* {{{ s_val2ch(v, caps) */
3555 
3556 static char
3557 s_val2ch(int v, int caps)
3558 {
3559  assert(v >= 0);
3560 
3561  if (v < 10)
3562  return v + '0';
3563  else
3564  {
3565  char out = (v - 10) + 'a';
3566 
3567  if (caps)
3568  return toupper((unsigned char) out);
3569  else
3570  return out;
3571  }
3572 }
3573 
3574 /* }}} */
3575 
3576 /* {{{ s_2comp(buf, len) */
3577 
3578 static void
3579 s_2comp(unsigned char *buf, int len)
3580 {
3581  int i;
3582  unsigned short s = 1;
3583 
3584  for (i = len - 1; i >= 0; --i)
3585  {
3586  unsigned char c = ~buf[i];
3587 
3588  s = c + s;
3589  c = s & UCHAR_MAX;
3590  s >>= CHAR_BIT;
3591 
3592  buf[i] = c;
3593  }
3594 
3595  /* last carry out is ignored */
3596 }
3597 
3598 /* }}} */
3599 
3600 /* {{{ s_tobin(z, buf, *limpos) */
3601 
3602 static mp_result
3603 s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
3604 {
3605  mp_size uz;
3606  mp_digit *dz;
3607  int pos = 0,
3608  limit = *limpos;
3609 
3610  uz = MP_USED(z);
3611  dz = MP_DIGITS(z);
3612  while (uz > 0 && pos < limit)
3613  {
3614  mp_digit d = *dz++;
3615  int i;
3616 
3617  for (i = sizeof(mp_digit); i > 0 && pos < limit; --i)
3618  {
3619  buf[pos++] = (unsigned char) d;
3620  d >>= CHAR_BIT;
3621 
3622  /* Don't write leading zeroes */
3623  if (d == 0 && uz == 1)
3624  i = 0; /* exit loop without signaling truncation */
3625  }
3626 
3627  /* Detect truncation (loop exited with pos >= limit) */
3628  if (i > 0)
3629  break;
3630 
3631  --uz;
3632  }
3633 
3634  if (pad != 0 && (buf[pos - 1] >> (CHAR_BIT - 1)))
3635  {
3636  if (pos < limit)
3637  buf[pos++] = 0;
3638  else
3639  uz = 1;
3640  }
3641 
3642  /* Digits are in reverse order, fix that */
3643  REV(unsigned char, buf, pos);
3644 
3645  /* Return the number of bytes actually written */
3646  *limpos = pos;
3647 
3648  return (uz == 0) ? MP_OK : MP_TRUNC;
3649 }
3650 
3651 /* }}} */
3652 
3653 /* {{{ s_print(tag, z) */
3654 
3655 #if 0
3656 void
3657 s_print(char *tag, mp_int z)
3658 {
3659  int i;
3660 
3661  fprintf(stderr, "%s: %c ", tag,
3662  (MP_SIGN(z) == MP_NEG) ? '-' : '+');
3663 
3664  for (i = MP_USED(z) - 1; i >= 0; --i)
3665  fprintf(stderr, "%0*X", (int) (MP_DIGIT_BIT / 4), z->digits[i]);
3666 
3667  fputc('\n', stderr);
3668 
3669 }
3670 
3671 void
3672 s_print_buf(char *tag, mp_digit *buf, mp_size num)
3673 {
3674  int i;
3675 
3676  fprintf(stderr, "%s: ", tag);
3677 
3678  for (i = num - 1; i >= 0; --i)
3679  fprintf(stderr, "%0*X", (int) (MP_DIGIT_BIT / 4), buf[i]);
3680 
3681  fputc('\n', stderr);
3682 }
3683 #endif
3684 
3685 /* }}} */
3686 
3687 /* HERE THERE BE DRAGONS */
mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
Definition: imath.c:794
static const char * s_error_msg[]
Definition: imath.c:56
static int s_pad(mp_int z, mp_size min)
Definition: imath.c:2324
int mp_int_compare_unsigned(mp_int a, mp_int b)
Definition: imath.c:1277
static mp_size multiply_threshold
Definition: imath.c:170
static char s_val2ch(int v, int caps)
Definition: imath.c:3557
static int s_dp2k(mp_int z)
Definition: imath.c:3113
static void skip(struct vars *v)
Definition: regc_lex.c:1109
static mp_size default_precision
Definition: imath.c:167
static struct @76 value
#define MP_DIGIT_MAX
Definition: imath.h:47
mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c)
Definition: imath.c:1481
static int s_outlen(mp_int z, mp_size r)
Definition: imath.c:3509
#define UPPER_HALF(W)
Definition: imath.c:159
const mp_result MP_TRUE
Definition: imath.c:45
#define MP_ALLOC(Z)
Definition: imath.h:67
mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit)
Definition: imath.c:2103
mp_result mp_int_to_string(mp_int z, mp_size radix, char *str, int limit)
Definition: imath.c:1901
mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len)
Definition: imath.c:2210
#define MAX(A, B)
Definition: imath.c:139
#define USQR(X, Z)
Definition: imath.c:155
const mp_sign MP_NEG
Definition: imath.c:52
static int s_qmul(mp_int z, mp_size p2)
Definition: imath.c:2996
#define s_free(P)
Definition: imath.c:182
mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str)
Definition: imath.c:1993
mp_digit * digits
Definition: imath.h:59
mp_size used
Definition: imath.h:61
mp_result mp_int_sub_value(mp_int a, int value, mp_int c)
Definition: imath.c:779
mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c)
Definition: imath.c:1536
mp_result mp_int_div_pow2(mp_int a, int p2, mp_int q, mp_int r)
Definition: imath.c:1146
#define MP_DIGIT_BIT
Definition: imath.h:80
static mp_digit * s_realloc(mp_digit *old, mp_size num)
Definition: imath.c:2298
#define MP_CAP_DIGITS
Definition: imath.c:70
#define px_free(p)
Definition: px.h:46
#define NRCHECK(TEST)
Definition: imath.c:75
mp_size mp_get_default_precision(void)
Definition: imath.c:319
static void s_qmod(mp_int z, mp_size p2)
Definition: imath.c:2976
mp_result mp_int_copy(mp_int a, mp_int c)
Definition: imath.c:510
static int s_2expt(mp_int z, int k)
Definition: imath.c:3175
mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit)
Definition: imath.c:2196
static const double s_log2[]
Definition: imath.c:85
static mp_word mp_flags
Definition: imath.c:173
mp_result mp_int_mod(mp_int a, mp_int m, mp_int c)
Definition: imath.c:1080
#define ZERO(P, S)
Definition: imath.c:116
#define MP_VALUE_DIGITS(V)
Definition: imath.c:109
static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
Definition: imath.c:2742
#define mp_int_is_even(Z)
Definition: imath.h:90
#define HIGH_BIT_SET(W)
Definition: imath.c:161
mp_result mp_int_count_bits(mp_int z)
Definition: imath.c:2073
const mp_result MP_RANGE
Definition: imath.c:47
mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
Definition: imath.c:1420
mp_result mp_int_expt_value(int a, int b, mp_int c)
Definition: imath.c:1204
mp_result mp_int_div_value(mp_int a, int value, mp_int q, int *r)
Definition: imath.c:1119
void mp_int_swap(mp_int a, mp_int c)
Definition: imath.c:539
#define CMPZ(Z)
Definition: imath.c:146
static void s_2comp(unsigned char *buf, int len)
Definition: imath.c:3579
mp_result mp_int_binary_len(mp_int z)
Definition: imath.c:2170
static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
Definition: imath.c:3298
const mp_sign MP_ZPOS
Definition: imath.c:53
static mp_size s_inlen(int len, mp_size r)
Definition: imath.c:3525
mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
Definition: imath.c:699
mp_result mp_int_init_copy(mp_int z, mp_int old)
Definition: imath.c:412
const mp_result MP_FALSE
Definition: imath.c:44
const mp_result MP_TRUNC
Definition: imath.c:49
#define SETUP(E, C)
Definition: imath.c:143
#define ROUND_PREC(P)
Definition: imath.c:113
mp_result mp_int_to_int(mp_int z, int *out)
Definition: imath.c:1865
mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c)
Definition: imath.c:1336
mp_result mp_int_init(mp_int z)
Definition: imath.c:361
int mp_int_compare_zero(mp_int z)
Definition: imath.c:1289
static mp_digit * s_alloc(mp_size num)
Definition: imath.c:2284
void mp_int_free(mp_int z)
Definition: imath.c:495
char * c
static char * buf
Definition: pg_test_fsync.c:65
mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
Definition: imath.c:607
static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
Definition: imath.c:2864
mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r)
Definition: imath.c:955
#define mp_int_is_odd(Z)
Definition: imath.h:89
#define assert(TEST)
Definition: imath.c:37
mp_result mp_int_init_value(mp_int z, int value)
Definition: imath.c:438
void mp_int_zero(mp_int z)
Definition: imath.c:555
mp_result mp_int_sqrt(mp_int a, mp_int c)
Definition: imath.c:1808
mp_result mp_int_mul_pow2(mp_int a, int p2, mp_int c)
Definition: imath.c:882
mp_result mp_int_init_size(mp_int z, mp_size prec)
Definition: imath.c:389
mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len)
Definition: imath.c:2125
mp_sign sign
Definition: imath.h:62
static void s_dadd(mp_int a, mp_digit b)
Definition: imath.c:2808
#define TEMP(K)
Definition: imath.c:142
static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
Definition: imath.c:3242
#define px_realloc(p, s)
Definition: px.h:45
#define UMUL(X, Y, Z)
Definition: imath.c:149
mp_result mp_int_add_value(mp_int a, int value, mp_int c)
Definition: imath.c:684
#define MP_DIGITS(Z)
Definition: imath.h:66
static int s_isp2(mp_int z)
Definition: imath.c:3143
mp_result mp_int_unsigned_len(mp_int z)
Definition: imath.c:2241
static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b)
Definition: imath.c:2464
#define LOWER_HALF(W)
Definition: imath.c:160
mp_result mp_int_exptmod_bvalue(int value, mp_int b, mp_int m, mp_int c)
Definition: imath.c:1404
mp_result mp_int_set_value(mp_int z, int value)
Definition: imath.c:455
const mp_result MP_OK
Definition: imath.c:43
const mp_result MP_BADARG
Definition: imath.c:50
int mp_int_divisible_value(mp_int a, int v)
Definition: imath.c:1781
mp_result mp_int_sqr(mp_int a, mp_int c)
Definition: imath.c:902
void mp_set_default_precision(mp_size s)
Definition: imath.c:329
const char * mp_error_string(mp_result res)
Definition: imath.c:2259
const mp_result MP_UNDEF
Definition: imath.c:48
static void s_qdiv(mp_int z, mp_size p2)
Definition: imath.c:2919
mp_result mp_int_exptmod_evalue(mp_int a, int value, mp_int m, mp_int c)
Definition: imath.c:1389
static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b)
Definition: imath.c:2503
#define MP_USED(Z)
Definition: imath.h:68
#define MP_SIGN(Z)
Definition: imath.h:69
static int s_vpack(int v, mp_digit t[])
Definition: imath.c:2400
#define NULL
Definition: c.h:229
static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
Definition: imath.c:3603
static int s_vcmp(mp_int a, int v)
Definition: imath.c:2443
mp_result mp_int_redux_const(mp_int m, mp_int c)
Definition: imath.c:1469
mp_result mp_int_neg(mp_int a, mp_int c)
Definition: imath.c:587
static int s_norm(mp_int a, mp_int b)
Definition: imath.c:3200
static mp_result s_brmu(mp_int z, mp_int m)
Definition: imath.c:3226
#define ADD_WILL_OVERFLOW(W, V)
Definition: imath.c:162
unsigned char mp_sign
Definition: imath.h:39
static const char * s_unknown_err
Definition: imath.c:55
mp_result mp_int_abs(mp_int a, mp_int c)
Definition: imath.c:569
void mp_int_clear(mp_int z)
Definition: imath.c:478
mp_result mp_int_expt(mp_int a, int b, mp_int c)
Definition: imath.c:1166
static void s_fake(mp_int z, int value, mp_digit vbuf[])
Definition: imath.c:2364
mp_size mp_get_multiply_threshold(void)
Definition: imath.c:341
static mp_result s_udiv(mp_int a, mp_int b)
Definition: imath.c:3399
static int s_cdig(mp_digit *da, mp_digit *db, mp_size len)
Definition: imath.c:2379
static int s_ch2val(char c, int r)
Definition: imath.c:3538
#define px_alloc(s)
Definition: px.h:44
static int s_ucmp(mp_int a, mp_int b)
Definition: imath.c:2425
int mp_result
Definition: imath.h:41
unsigned int mp_size
Definition: imath.h:40
mp_size alloc
Definition: imath.h:60
mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
Definition: imath.c:2004
#define REV(T, A, N)
Definition: imath.c:125
#define CHECK(TEST)
Definition: imath.c:74
mp_int mp_int_alloc(void)
Definition: imath.c:371
mp_result mp_int_mul_value(mp_int a, int value, mp_int c)
Definition: imath.c:867
static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b)
Definition: imath.c:2644
static void s_dmul(mp_int a, mp_digit b)
Definition: imath.c:2838
static mp_digit s_ddiv(mp_int a, mp_digit b)
Definition: imath.c:2886
int i
static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
Definition: imath.c:2677
#define CLAMP(Z)
Definition: imath.c:131
#define SWAP(T, A, B)
Definition: imath.c:140
const mp_result MP_MEMORY
Definition: imath.c:46
uint32 mp_digit
Definition: imath.h:44
#define MIN(A, B)
Definition: imath.c:138
int mp_int_compare_value(mp_int z, int value)
Definition: imath.c:1306
mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c, mp_int x, mp_int y)
Definition: imath.c:1631
Definition: imath.h:57
#define MP_MAX_RADIX
Definition: imath.h:84
mp_result mp_int_string_len(mp_int z, mp_size radix)
Definition: imath.c:1969
int mp_int_is_pow2(mp_int z)
Definition: imath.c:1796
static int s_qsub(mp_int z, mp_size p2)
Definition: imath.c:3078
uint64 mp_word
Definition: imath.h:45
void mp_set_multiply_threshold(mp_size s)
Definition: imath.c:351
#define COPY(P, Q, S)
Definition: imath.c:120
static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc, mp_size size_a, mp_size size_b)
Definition: imath.c:2541
static int cmp(const chr *x, const chr *y, size_t len)
Definition: regc_locale.c:742
int mp_int_compare(mp_int a, mp_int b)
Definition: imath.c:1242