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