source4/heimdal/lib/des/imath/imath.c

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