crypto/external/bsd/heimdal/dist/lib/hcrypto/libtommath/mtest/mpi.c

   1 /*      $NetBSD: mpi.c,v 1.1.1.3 2014/04/24 12:45:39 pettai Exp $       */
   2
   3 /*
   4     mpi.c
   5
   6     by Michael J. Fromberger <sting@linguist.dartmouth.edu>
   7     Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved
   8
   9     Arbitrary precision integer arithmetic library
  10
  11     Id: mpi.c,v 1.2 2005/05/05 14:38:47 tom Exp
  12  */
  13
  14 #include "mpi.h"
  15 #include <stdlib.h>
  16 #include <string.h>
  17 #include <ctype.h>
  18
  19 #if MP_DEBUG
  20 #include <stdio.h>
  21
  22 #define DIAG(T,V) {fprintf(stderr,T);mp_print(V,stderr);fputc('\n',stderr);}
  23 #else
  24 #define DIAG(T,V)
  25 #endif
  26
  27 /*
  28    If MP_LOGTAB is not defined, use the math library to compute the
  29    logarithms on the fly.  Otherwise, use the static table below.
  30    Pick which works best for your system.
  31  */
  32 #if MP_LOGTAB
  33
  34 /* {{{ s_logv_2[] - log table for 2 in various bases */
  35
  36 /*
  37   A table of the logs of 2 for various bases (the 0 and 1 entries of
  38   this table are meaningless and should not be referenced).
  39
  40   This table is used to compute output lengths for the mp_toradix()
  41   function.  Since a number n in radix r takes up about log_r(n)
  42   digits, we estimate the output size by taking the least integer
  43   greater than log_r(n), where:
  44
  45   log_r(n) = log_2(n) * log_r(2)
  46
  47   This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
  48   which are the output bases supported.
  49  */
  50
  51 #include "logtab.h"
  52
  53 /* }}} */
  54 #define LOG_V_2(R)  s_logv_2[(R)]
  55
  56 #else
  57
  58 #include <math.h>
  59 #define LOG_V_2(R)  (log(2.0)/log(R))
  60
  61 #endif
  62
  63 /* Default precision for newly created mp_int's      */
  64 static unsigned int s_mp_defprec = MP_DEFPREC;
  65
  66 /* {{{ Digit arithmetic macros */
  67
  68 /*
  69   When adding and multiplying digits, the results can be larger than
  70   can be contained in an mp_digit.  Thus, an mp_word is used.  These
  71   macros mask off the upper and lower digits of the mp_word (the
  72   mp_word may be more than 2 mp_digits wide, but we only concern
  73   ourselves with the low-order 2 mp_digits)
  74
  75   If your mp_word DOES have more than 2 mp_digits, you need to
  76   uncomment the first line, and comment out the second.
  77  */
  78
  79 /* #define  CARRYOUT(W)  (((W)>>DIGIT_BIT)&MP_DIGIT_MAX) */
  80 #define  CARRYOUT(W)  ((W)>>DIGIT_BIT)
  81 #define  ACCUM(W)     ((W)&MP_DIGIT_MAX)
  82
  83 /* }}} */
  84
  85 /* {{{ Comparison constants */
  86
  87 #define  MP_LT       -1
  88 #define  MP_EQ        0
  89 #define  MP_GT        1
  90
  91 /* }}} */
  92
  93 /* {{{ Constant strings */
  94
  95 /* Constant strings returned by mp_strerror() */
  96 static const char *mp_err_string[] = {
  97   "unknown result code",     /* say what?            */
  98   "boolean true",            /* MP_OKAY, MP_YES      */
  99   "boolean false",           /* MP_NO                */
 100   "out of memory",           /* MP_MEM               */
 101   "argument out of range",   /* MP_RANGE             */
 102   "invalid input parameter", /* MP_BADARG            */
 103   "result is undefined"      /* MP_UNDEF             */
 104 };
 105
 106 /* Value to digit maps for radix conversion   */
 107
 108 /* s_dmap_1 - standard digits and letters */
 109 static const char *s_dmap_1 =
 110   "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
 111
 112 #if 0
 113 /* s_dmap_2 - base64 ordering for digits  */
 114 static const char *s_dmap_2 =
 115   "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
 116 #endif
 117
 118 /* }}} */
 119
 120 /* {{{ Static function declarations */
 121
 122 /*
 123    If MP_MACRO is false, these will be defined as actual functions;
 124    otherwise, suitable macro definitions will be used.  This works
 125    around the fact that ANSI C89 doesn't support an 'inline' keyword
 126    (although I hear C9x will ... about bloody time).  At present, the
 127    macro definitions are identical to the function bodies, but they'll
 128    expand in place, instead of generating a function call.
 129
 130    I chose these particular functions to be made into macros because
 131    some profiling showed they are called a lot on a typical workload,
 132    and yet they are primarily housekeeping.
 133  */
 134 #if MP_MACRO == 0
 135  void     s_mp_setz(mp_digit *dp, mp_size count); /* zero digits           */
 136  void     s_mp_copy(mp_digit *sp, mp_digit *dp, mp_size count); /* copy    */
 137  void    *s_mp_alloc(size_t nb, size_t ni);       /* general allocator     */
 138  void     s_mp_free(void *ptr);                   /* general free function */
 139 #else
 140
 141  /* Even if these are defined as macros, we need to respect the settings
 142     of the MP_MEMSET and MP_MEMCPY configuration options...
 143   */
 144  #if MP_MEMSET == 0
 145   #define  s_mp_setz(dp, count) \
 146        {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=0;}
 147  #else
 148   #define  s_mp_setz(dp, count) memset(dp, 0, (count) * sizeof(mp_digit))
 149  #endif /* MP_MEMSET */
 150
 151  #if MP_MEMCPY == 0
 152   #define  s_mp_copy(sp, dp, count) \
 153        {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=(sp)[ix];}
 154  #else
 155   #define  s_mp_copy(sp, dp, count) memcpy(dp, sp, (count) * sizeof(mp_digit))
 156  #endif /* MP_MEMCPY */
 157
 158  #define  s_mp_alloc(nb, ni)  calloc(nb, ni)
 159  #define  s_mp_free(ptr) {if(ptr) free(ptr);}
 160 #endif /* MP_MACRO */
 161
 162 mp_err   s_mp_grow(mp_int *mp, mp_size min);   /* increase allocated size */
 163 mp_err   s_mp_pad(mp_int *mp, mp_size min);    /* left pad with zeroes    */
 164
 165 void     s_mp_clamp(mp_int *mp);               /* clip leading zeroes     */
 166
 167 void     s_mp_exch(mp_int *a, mp_int *b);      /* swap a and b in place   */
 168
 169 mp_err   s_mp_lshd(mp_int *mp, mp_size p);     /* left-shift by p digits  */
 170 void     s_mp_rshd(mp_int *mp, mp_size p);     /* right-shift by p digits */
 171 void     s_mp_div_2d(mp_int *mp, mp_digit d);  /* divide by 2^d in place  */
 172 void     s_mp_mod_2d(mp_int *mp, mp_digit d);  /* modulo 2^d in place     */
 173 mp_err   s_mp_mul_2d(mp_int *mp, mp_digit d);  /* multiply by 2^d in place*/
 174 void     s_mp_div_2(mp_int *mp);               /* divide by 2 in place    */
 175 mp_err   s_mp_mul_2(mp_int *mp);               /* multiply by 2 in place  */
 176 mp_digit s_mp_norm(mp_int *a, mp_int *b);      /* normalize for division  */
 177 mp_err   s_mp_add_d(mp_int *mp, mp_digit d);   /* unsigned digit addition */
 178 mp_err   s_mp_sub_d(mp_int *mp, mp_digit d);   /* unsigned digit subtract */
 179 mp_err   s_mp_mul_d(mp_int *mp, mp_digit d);   /* unsigned digit multiply */
 180 mp_err   s_mp_div_d(mp_int *mp, mp_digit d, mp_digit *r);
 181                                                /* unsigned digit divide   */
 182 mp_err   s_mp_reduce(mp_int *x, mp_int *m, mp_int *mu);
 183                                                /* Barrett reduction       */
 184 mp_err   s_mp_add(mp_int *a, mp_int *b);       /* magnitude addition      */
 185 mp_err   s_mp_sub(mp_int *a, mp_int *b);       /* magnitude subtract      */
 186 mp_err   s_mp_mul(mp_int *a, mp_int *b);       /* magnitude multiply      */
 187 #if 0
 188 void     s_mp_kmul(mp_digit *a, mp_digit *b, mp_digit *out, mp_size len);
 189                                                /* multiply buffers in place */
 190 #endif
 191 #if MP_SQUARE
 192 mp_err   s_mp_sqr(mp_int *a);                  /* magnitude square        */
 193 #else
 194 #define  s_mp_sqr(a) s_mp_mul(a, a)
 195 #endif
 196 mp_err   s_mp_div(mp_int *a, mp_int *b);       /* magnitude divide        */
 197 mp_err   s_mp_2expt(mp_int *a, mp_digit k);    /* a = 2^k                 */
 198 int      s_mp_cmp(mp_int *a, mp_int *b);       /* magnitude comparison    */
 199 int      s_mp_cmp_d(mp_int *a, mp_digit d);    /* magnitude digit compare */
 200 int      s_mp_ispow2(mp_int *v);               /* is v a power of 2?      */
 201 int      s_mp_ispow2d(mp_digit d);             /* is d a power of 2?      */
 202
 203 int      s_mp_tovalue(char ch, int r);          /* convert ch to value    */
 204 char     s_mp_todigit(int val, int r, int low); /* convert val to digit   */
 205 int      s_mp_outlen(int bits, int r);          /* output length in bytes */
 206
 207 /* }}} */
 208
 209 /* {{{ Default precision manipulation */
 210
 211 unsigned int mp_get_prec(void)
 212 {
 213   return s_mp_defprec;
 214
 215 } /* end mp_get_prec() */
 216
 217 void         mp_set_prec(unsigned int prec)
 218 {
 219   if(prec == 0)
 220     s_mp_defprec = MP_DEFPREC;
 221   else
 222     s_mp_defprec = prec;
 223
 224 } /* end mp_set_prec() */
 225
 226 /* }}} */
 227
 228 /*------------------------------------------------------------------------*/
 229 /* {{{ mp_init(mp) */
 230
 231 /*
 232   mp_init(mp)
 233
 234   Initialize a new zero-valued mp_int.  Returns MP_OKAY if successful,
 235   MP_MEM if memory could not be allocated for the structure.
 236  */
 237
 238 mp_err mp_init(mp_int *mp)
 239 {
 240   return mp_init_size(mp, s_mp_defprec);
 241
 242 } /* end mp_init() */
 243
 244 /* }}} */
 245
 246 /* {{{ mp_init_array(mp[], count) */
 247
 248 mp_err mp_init_array(mp_int mp[], int count)
 249 {
 250   mp_err  res;
 251   int     pos;
 252
 253   ARGCHK(mp !=NULL && count > 0, MP_BADARG);
 254
 255   for(pos = 0; pos < count; ++pos) {
 256     if((res = mp_init(&mp[pos])) != MP_OKAY)
 257       goto CLEANUP;
 258   }
 259
 260   return MP_OKAY;
 261
 262  CLEANUP:
 263   while(--pos >= 0)
 264     mp_clear(&mp[pos]);
 265
 266   return res;
 267
 268 } /* end mp_init_array() */
 269
 270 /* }}} */
 271
 272 /* {{{ mp_init_size(mp, prec) */
 273
 274 /*
 275   mp_init_size(mp, prec)
 276
 277   Initialize a new zero-valued mp_int with at least the given
 278   precision; returns MP_OKAY if successful, or MP_MEM if memory could
 279   not be allocated for the structure.
 280  */
 281
 282 mp_err mp_init_size(mp_int *mp, mp_size prec)
 283 {
 284   ARGCHK(mp != NULL && prec > 0, MP_BADARG);
 285
 286   if((DIGITS(mp) = s_mp_alloc(prec, sizeof(mp_digit))) == NULL)
 287     return MP_MEM;
 288
 289   SIGN(mp) = MP_ZPOS;
 290   USED(mp) = 1;
 291   ALLOC(mp) = prec;
 292
 293   return MP_OKAY;
 294
 295 } /* end mp_init_size() */
 296
 297 /* }}} */
 298
 299 /* {{{ mp_init_copy(mp, from) */
 300
 301 /*
 302   mp_init_copy(mp, from)
 303
 304   Initialize mp as an exact copy of from.  Returns MP_OKAY if
 305   successful, MP_MEM if memory could not be allocated for the new
 306   structure.
 307  */
 308
 309 mp_err mp_init_copy(mp_int *mp, mp_int *from)
 310 {
 311   ARGCHK(mp != NULL && from != NULL, MP_BADARG);
 312
 313   if(mp == from)
 314     return MP_OKAY;
 315
 316   if((DIGITS(mp) = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL)
 317     return MP_MEM;
 318
 319   s_mp_copy(DIGITS(from), DIGITS(mp), USED(from));
 320   USED(mp) = USED(from);
 321   ALLOC(mp) = USED(from);
 322   SIGN(mp) = SIGN(from);
 323
 324   return MP_OKAY;
 325
 326 } /* end mp_init_copy() */
 327
 328 /* }}} */
 329
 330 /* {{{ mp_copy(from, to) */
 331
 332 /*
 333   mp_copy(from, to)
 334
 335   Copies the mp_int 'from' to the mp_int 'to'.  It is presumed that
 336   'to' has already been initialized (if not, use mp_init_copy()
 337   instead). If 'from' and 'to' are identical, nothing happens.
 338  */
 339
 340 mp_err mp_copy(mp_int *from, mp_int *to)
 341 {
 342   ARGCHK(from != NULL && to != NULL, MP_BADARG);
 343
 344   if(from == to)
 345     return MP_OKAY;
 346
 347   { /* copy */
 348     mp_digit   *tmp;
 349
 350     /*
 351       If the allocated buffer in 'to' already has enough space to hold
 352       all the used digits of 'from', we'll re-use it to avoid hitting
 353       the memory allocater more than necessary; otherwise, we'd have
 354       to grow anyway, so we just allocate a hunk and make the copy as
 355       usual
 356      */
 357     if(ALLOC(to) >= USED(from)) {
 358       s_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));
 359       s_mp_copy(DIGITS(from), DIGITS(to), USED(from));
 360
 361     } else {
 362       if((tmp = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL)
 363         return MP_MEM;
 364
 365       s_mp_copy(DIGITS(from), tmp, USED(from));
 366
 367       if(DIGITS(to) != NULL) {
 368 #if MP_CRYPTO
 369         s_mp_setz(DIGITS(to), ALLOC(to));
 370 #endif
 371         s_mp_free(DIGITS(to));
 372       }
 373
 374       DIGITS(to) = tmp;
 375       ALLOC(to) = USED(from);
 376     }
 377
 378     /* Copy the precision and sign from the original */
 379     USED(to) = USED(from);
 380     SIGN(to) = SIGN(from);
 381   } /* end copy */
 382
 383   return MP_OKAY;
 384
 385 } /* end mp_copy() */
 386
 387 /* }}} */
 388
 389 /* {{{ mp_exch(mp1, mp2) */
 390
 391 /*
 392   mp_exch(mp1, mp2)
 393
 394   Exchange mp1 and mp2 without allocating any intermediate memory
 395   (well, unless you count the stack space needed for this call and the
 396   locals it creates...).  This cannot fail.
 397  */
 398
 399 void mp_exch(mp_int *mp1, mp_int *mp2)
 400 {
 401 #if MP_ARGCHK == 2
 402   assert(mp1 != NULL && mp2 != NULL);
 403 #else
 404   if(mp1 == NULL || mp2 == NULL)
 405     return;
 406 #endif
 407
 408   s_mp_exch(mp1, mp2);
 409
 410 } /* end mp_exch() */
 411
 412 /* }}} */
 413
 414 /* {{{ mp_clear(mp) */
 415
 416 /*
 417   mp_clear(mp)
 418
 419   Release the storage used by an mp_int, and void its fields so that
 420   if someone calls mp_clear() again for the same int later, we won't
 421   get tollchocked.
 422  */
 423
 424 void   mp_clear(mp_int *mp)
 425 {
 426   if(mp == NULL)
 427     return;
 428
 429   if(DIGITS(mp) != NULL) {
 430 #if MP_CRYPTO
 431     s_mp_setz(DIGITS(mp), ALLOC(mp));
 432 #endif
 433     s_mp_free(DIGITS(mp));
 434     DIGITS(mp) = NULL;
 435   }
 436
 437   USED(mp) = 0;
 438   ALLOC(mp) = 0;
 439
 440 } /* end mp_clear() */
 441
 442 /* }}} */
 443
 444 /* {{{ mp_clear_array(mp[], count) */
 445
 446 void   mp_clear_array(mp_int mp[], int count)
 447 {
 448   ARGCHK(mp != NULL && count > 0, MP_BADARG);
 449
 450   while(--count >= 0)
 451     mp_clear(&mp[count]);
 452
 453 } /* end mp_clear_array() */
 454
 455 /* }}} */
 456
 457 /* {{{ mp_zero(mp) */
 458
 459 /*
 460   mp_zero(mp)
 461
 462   Set mp to zero.  Does not change the allocated size of the structure,
 463   and therefore cannot fail (except on a bad argument, which we ignore)
 464  */
 465 void   mp_zero(mp_int *mp)
 466 {
 467   if(mp == NULL)
 468     return;
 469
 470   s_mp_setz(DIGITS(mp), ALLOC(mp));
 471   USED(mp) = 1;
 472   SIGN(mp) = MP_ZPOS;
 473
 474 } /* end mp_zero() */
 475
 476 /* }}} */
 477
 478 /* {{{ mp_set(mp, d) */
 479
 480 void   mp_set(mp_int *mp, mp_digit d)
 481 {
 482   if(mp == NULL)
 483     return;
 484
 485   mp_zero(mp);
 486   DIGIT(mp, 0) = d;
 487
 488 } /* end mp_set() */
 489
 490 /* }}} */
 491
 492 /* {{{ mp_set_int(mp, z) */
 493
 494 mp_err mp_set_int(mp_int *mp, long z)
 495 {
 496   int            ix;
 497   unsigned long  v = abs(z);
 498   mp_err         res;
 499
 500   ARGCHK(mp != NULL, MP_BADARG);
 501
 502   mp_zero(mp);
 503   if(z == 0)
 504     return MP_OKAY;  /* shortcut for zero */
 505
 506   for(ix = sizeof(long) - 1; ix >= 0; ix--) {
 507
 508     if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
 509       return res;
 510
 511     res = s_mp_add_d(mp,
 512                      (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));
 513     if(res != MP_OKAY)
 514       return res;
 515
 516   }
 517
 518   if(z < 0)
 519     SIGN(mp) = MP_NEG;
 520
 521   return MP_OKAY;
 522
 523 } /* end mp_set_int() */
 524
 525 /* }}} */
 526
 527 /*------------------------------------------------------------------------*/
 528 /* {{{ Digit arithmetic */
 529
 530 /* {{{ mp_add_d(a, d, b) */
 531
 532 /*
 533   mp_add_d(a, d, b)
 534
 535   Compute the sum b = a + d, for a single digit d.  Respects the sign of
 536   its primary addend (single digits are unsigned anyway).
 537  */
 538
 539 mp_err mp_add_d(mp_int *a, mp_digit d, mp_int *b)
 540 {
 541   mp_err   res = MP_OKAY;
 542
 543   ARGCHK(a != NULL && b != NULL, MP_BADARG);
 544
 545   if((res = mp_copy(a, b)) != MP_OKAY)
 546     return res;
 547
 548   if(SIGN(b) == MP_ZPOS) {
 549     res = s_mp_add_d(b, d);
 550   } else if(s_mp_cmp_d(b, d) >= 0) {
 551     res = s_mp_sub_d(b, d);
 552   } else {
 553     SIGN(b) = MP_ZPOS;
 554
 555     DIGIT(b, 0) = d - DIGIT(b, 0);
 556   }
 557
 558   return res;
 559
 560 } /* end mp_add_d() */
 561
 562 /* }}} */
 563
 564 /* {{{ mp_sub_d(a, d, b) */
 565
 566 /*
 567   mp_sub_d(a, d, b)
 568
 569   Compute the difference b = a - d, for a single digit d.  Respects the
 570   sign of its subtrahend (single digits are unsigned anyway).
 571  */
 572
 573 mp_err mp_sub_d(mp_int *a, mp_digit d, mp_int *b)
 574 {
 575   mp_err   res;
 576
 577   ARGCHK(a != NULL && b != NULL, MP_BADARG);
 578
 579   if((res = mp_copy(a, b)) != MP_OKAY)
 580     return res;
 581
 582   if(SIGN(b) == MP_NEG) {
 583     if((res = s_mp_add_d(b, d)) != MP_OKAY)
 584       return res;
 585
 586   } else if(s_mp_cmp_d(b, d) >= 0) {
 587     if((res = s_mp_sub_d(b, d)) != MP_OKAY)
 588       return res;
 589
 590   } else {
 591     mp_neg(b, b);
 592
 593     DIGIT(b, 0) = d - DIGIT(b, 0);
 594     SIGN(b) = MP_NEG;
 595   }
 596
 597   if(s_mp_cmp_d(b, 0) == 0)
 598     SIGN(b) = MP_ZPOS;
 599
 600   return MP_OKAY;
 601
 602 } /* end mp_sub_d() */
 603
 604 /* }}} */
 605
 606 /* {{{ mp_mul_d(a, d, b) */
 607
 608 /*
 609   mp_mul_d(a, d, b)
 610
 611   Compute the product b = a * d, for a single digit d.  Respects the sign
 612   of its multiplicand (single digits are unsigned anyway)
 613  */
 614
 615 mp_err mp_mul_d(mp_int *a, mp_digit d, mp_int *b)
 616 {
 617   mp_err  res;
 618
 619   ARGCHK(a != NULL && b != NULL, MP_BADARG);
 620
 621   if(d == 0) {
 622     mp_zero(b);
 623     return MP_OKAY;
 624   }
 625
 626   if((res = mp_copy(a, b)) != MP_OKAY)
 627     return res;
 628
 629   res = s_mp_mul_d(b, d);
 630
 631   return res;
 632
 633 } /* end mp_mul_d() */
 634
 635 /* }}} */
 636
 637 /* {{{ mp_mul_2(a, c) */
 638
 639 mp_err mp_mul_2(mp_int *a, mp_int *c)
 640 {
 641   mp_err  res;
 642
 643   ARGCHK(a != NULL && c != NULL, MP_BADARG);
 644
 645   if((res = mp_copy(a, c)) != MP_OKAY)
 646     return res;
 647
 648   return s_mp_mul_2(c);
 649
 650 } /* end mp_mul_2() */
 651
 652 /* }}} */
 653
 654 /* {{{ mp_div_d(a, d, q, r) */
 655
 656 /*
 657   mp_div_d(a, d, q, r)
 658
 659   Compute the quotient q = a / d and remainder r = a mod d, for a
 660   single digit d.  Respects the sign of its divisor (single digits are
 661   unsigned anyway).
 662  */
 663
 664 mp_err mp_div_d(mp_int *a, mp_digit d, mp_int *q, mp_digit *r)
 665 {
 666   mp_err   res;
 667   mp_digit rem;
 668   int      pow;
 669
 670   ARGCHK(a != NULL, MP_BADARG);
 671
 672   if(d == 0)
 673     return MP_RANGE;
 674
 675   /* Shortcut for powers of two ... */
 676   if((pow = s_mp_ispow2d(d)) >= 0) {
 677     mp_digit  mask;
 678
 679     mask = (1 << pow) - 1;
 680     rem = DIGIT(a, 0) & mask;
 681
 682     if(q) {
 683       mp_copy(a, q);
 684       s_mp_div_2d(q, pow);
 685     }
 686
 687     if(r)
 688       *r = rem;
 689
 690     return MP_OKAY;
 691   }
 692
 693   /*
 694     If the quotient is actually going to be returned, we'll try to
 695     avoid hitting the memory allocator by copying the dividend into it
 696     and doing the division there.  This can't be any _worse_ than
 697     always copying, and will sometimes be better (since it won't make
 698     another copy)
 699
 700     If it's not going to be returned, we need to allocate a temporary
 701     to hold the quotient, which will just be discarded.
 702    */
 703   if(q) {
 704     if((res = mp_copy(a, q)) != MP_OKAY)
 705       return res;
 706
 707     res = s_mp_div_d(q, d, &rem);
 708     if(s_mp_cmp_d(q, 0) == MP_EQ)
 709       SIGN(q) = MP_ZPOS;
 710
 711   } else {
 712     mp_int  qp;
 713
 714     if((res = mp_init_copy(&qp, a)) != MP_OKAY)
 715       return res;
 716
 717     res = s_mp_div_d(&qp, d, &rem);
 718     if(s_mp_cmp_d(&qp, 0) == 0)
 719       SIGN(&qp) = MP_ZPOS;
 720
 721     mp_clear(&qp);
 722   }
 723
 724   if(r)
 725     *r = rem;
 726
 727   return res;
 728
 729 } /* end mp_div_d() */
 730
 731 /* }}} */
 732
 733 /* {{{ mp_div_2(a, c) */
 734
 735 /*
 736   mp_div_2(a, c)
 737
 738   Compute c = a / 2, disregarding the remainder.
 739  */
 740
 741 mp_err mp_div_2(mp_int *a, mp_int *c)
 742 {
 743   mp_err  res;
 744
 745   ARGCHK(a != NULL && c != NULL, MP_BADARG);
 746
 747   if((res = mp_copy(a, c)) != MP_OKAY)
 748     return res;
 749
 750   s_mp_div_2(c);
 751
 752   return MP_OKAY;
 753
 754 } /* end mp_div_2() */
 755
 756 /* }}} */
 757
 758 /* {{{ mp_expt_d(a, d, b) */
 759
 760 mp_err mp_expt_d(mp_int *a, mp_digit d, mp_int *c)
 761 {
 762   mp_int   s, x;
 763   mp_err   res;
 764
 765   ARGCHK(a != NULL && c != NULL, MP_BADARG);
 766
 767   if((res = mp_init(&s)) != MP_OKAY)
 768     return res;
 769   if((res = mp_init_copy(&x, a)) != MP_OKAY)
 770     goto X;
 771
 772   DIGIT(&s, 0) = 1;
 773
 774   while(d != 0) {
 775     if(d & 1) {
 776       if((res = s_mp_mul(&s, &x)) != MP_OKAY)
 777         goto CLEANUP;
 778     }
 779
 780     d >>= 1;
 781
 782     if((res = s_mp_sqr(&x)) != MP_OKAY)
 783       goto CLEANUP;
 784   }
 785
 786   s_mp_exch(&s, c);
 787
 788 CLEANUP:
 789   mp_clear(&x);
 790 X:
 791   mp_clear(&s);
 792
 793   return res;
 794
 795 } /* end mp_expt_d() */
 796
 797 /* }}} */
 798
 799 /* }}} */
 800
 801 /*------------------------------------------------------------------------*/
 802 /* {{{ Full arithmetic */
 803
 804 /* {{{ mp_abs(a, b) */
 805
 806 /*
 807   mp_abs(a, b)
 808
 809   Compute b = |a|.  'a' and 'b' may be identical.
 810  */
 811
 812 mp_err mp_abs(mp_int *a, mp_int *b)
 813 {
 814   mp_err   res;
 815
 816   ARGCHK(a != NULL && b != NULL, MP_BADARG);
 817
 818   if((res = mp_copy(a, b)) != MP_OKAY)
 819     return res;
 820
 821   SIGN(b) = MP_ZPOS;
 822
 823   return MP_OKAY;
 824
 825 } /* end mp_abs() */
 826
 827 /* }}} */
 828
 829 /* {{{ mp_neg(a, b) */
 830
 831 /*
 832   mp_neg(a, b)
 833
 834   Compute b = -a.  'a' and 'b' may be identical.
 835  */
 836
 837 mp_err mp_neg(mp_int *a, mp_int *b)
 838 {
 839   mp_err   res;
 840
 841   ARGCHK(a != NULL && b != NULL, MP_BADARG);
 842
 843   if((res = mp_copy(a, b)) != MP_OKAY)
 844     return res;
 845
 846   if(s_mp_cmp_d(b, 0) == MP_EQ)
 847     SIGN(b) = MP_ZPOS;
 848   else
 849     SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG;
 850
 851   return MP_OKAY;
 852
 853 } /* end mp_neg() */
 854
 855 /* }}} */
 856
 857 /* {{{ mp_add(a, b, c) */
 858
 859 /*
 860   mp_add(a, b, c)
 861
 862   Compute c = a + b.  All parameters may be identical.
 863  */
 864
 865 mp_err mp_add(mp_int *a, mp_int *b, mp_int *c)
 866 {
 867   mp_err  res;
 868   int     cmp;
 869
 870   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
 871
 872   if(SIGN(a) == SIGN(b)) { /* same sign:  add values, keep sign */
 873
 874     /* Commutativity of addition lets us do this in either order,
 875        so we avoid having to use a temporary even if the result
 876        is supposed to replace the output
 877      */
 878     if(c == b) {
 879       if((res = s_mp_add(c, a)) != MP_OKAY)
 880         return res;
 881     } else {
 882       if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
 883         return res;
 884
 885       if((res = s_mp_add(c, b)) != MP_OKAY)
 886         return res;
 887     }
 888
 889   } else if((cmp = s_mp_cmp(a, b)) > 0) {  /* different sign: a > b   */
 890
 891     /* If the output is going to be clobbered, we will use a temporary
 892        variable; otherwise, we'll do it without touching the memory
 893        allocator at all, if possible
 894      */
 895     if(c == b) {
 896       mp_int  tmp;
 897
 898       if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
 899         return res;
 900       if((res = s_mp_sub(&tmp, b)) != MP_OKAY) {
 901         mp_clear(&tmp);
 902         return res;
 903       }
 904
 905       s_mp_exch(&tmp, c);
 906       mp_clear(&tmp);
 907
 908     } else {
 909
 910       if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
 911         return res;
 912       if((res = s_mp_sub(c, b)) != MP_OKAY)
 913         return res;
 914
 915     }
 916
 917   } else if(cmp == 0) {             /* different sign, a == b   */
 918
 919     mp_zero(c);
 920     return MP_OKAY;
 921
 922   } else {                          /* different sign: a < b    */
 923
 924     /* See above... */
 925     if(c == a) {
 926       mp_int  tmp;
 927
 928       if((res = mp_init_copy(&tmp, b)) != MP_OKAY)
 929         return res;
 930       if((res = s_mp_sub(&tmp, a)) != MP_OKAY) {
 931         mp_clear(&tmp);
 932         return res;
 933       }
 934
 935       s_mp_exch(&tmp, c);
 936       mp_clear(&tmp);
 937
 938     } else {
 939
 940       if(c != b && (res = mp_copy(b, c)) != MP_OKAY)
 941         return res;
 942       if((res = s_mp_sub(c, a)) != MP_OKAY)
 943         return res;
 944
 945     }
 946   }
 947
 948   if(USED(c) == 1 && DIGIT(c, 0) == 0)
 949     SIGN(c) = MP_ZPOS;
 950
 951   return MP_OKAY;
 952
 953 } /* end mp_add() */
 954
 955 /* }}} */
 956
 957 /* {{{ mp_sub(a, b, c) */
 958
 959 /*
 960   mp_sub(a, b, c)
 961
 962   Compute c = a - b.  All parameters may be identical.
 963  */
 964
 965 mp_err mp_sub(mp_int *a, mp_int *b, mp_int *c)
 966 {
 967   mp_err  res;
 968   int     cmp;
 969
 970   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
 971
 972   if(SIGN(a) != SIGN(b)) {
 973     if(c == a) {
 974       if((res = s_mp_add(c, b)) != MP_OKAY)
 975         return res;
 976     } else {
 977       if(c != b && ((res = mp_copy(b, c)) != MP_OKAY))
 978         return res;
 979       if((res = s_mp_add(c, a)) != MP_OKAY)
 980         return res;
 981       SIGN(c) = SIGN(a);
 982     }
 983
 984   } else if((cmp = s_mp_cmp(a, b)) > 0) { /* Same sign, a > b */
 985     if(c == b) {
 986       mp_int  tmp;
 987
 988       if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
 989         return res;
 990       if((res = s_mp_sub(&tmp, b)) != MP_OKAY) {
 991         mp_clear(&tmp);
 992         return res;
 993       }
 994       s_mp_exch(&tmp, c);
 995       mp_clear(&tmp);
 996
 997     } else {
 998       if(c != a && ((res = mp_copy(a, c)) != MP_OKAY))
 999         return res;
1000
1001       if((res = s_mp_sub(c, b)) != MP_OKAY)
1002         return res;
1003     }
1004
1005   } else if(cmp == 0) {  /* Same sign, equal magnitude */
1006     mp_zero(c);
1007     return MP_OKAY;
1008
1009   } else {               /* Same sign, b > a */
1010     if(c == a) {
1011       mp_int  tmp;
1012
1013       if((res = mp_init_copy(&tmp, b)) != MP_OKAY)
1014         return res;
1015
1016       if((res = s_mp_sub(&tmp, a)) != MP_OKAY) {
1017         mp_clear(&tmp);
1018         return res;
1019       }
1020       s_mp_exch(&tmp, c);
1021       mp_clear(&tmp);
1022
1023     } else {
1024       if(c != b && ((res = mp_copy(b, c)) != MP_OKAY))
1025         return res;
1026
1027       if((res = s_mp_sub(c, a)) != MP_OKAY)
1028         return res;
1029     }
1030
1031     SIGN(c) = !SIGN(b);
1032   }
1033
1034   if(USED(c) == 1 && DIGIT(c, 0) == 0)
1035     SIGN(c) = MP_ZPOS;
1036
1037   return MP_OKAY;
1038
1039 } /* end mp_sub() */
1040
1041 /* }}} */
1042
1043 /* {{{ mp_mul(a, b, c) */
1044
1045 /*
1046   mp_mul(a, b, c)
1047
1048   Compute c = a * b.  All parameters may be identical.
1049  */
1050
1051 mp_err mp_mul(mp_int *a, mp_int *b, mp_int *c)
1052 {
1053   mp_err   res;
1054   mp_sign  sgn;
1055
1056   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1057
1058   sgn = (SIGN(a) == SIGN(b)) ? MP_ZPOS : MP_NEG;
1059
1060   if(c == b) {
1061     if((res = s_mp_mul(c, a)) != MP_OKAY)
1062       return res;
1063
1064   } else {
1065     if((res = mp_copy(a, c)) != MP_OKAY)
1066       return res;
1067
1068     if((res = s_mp_mul(c, b)) != MP_OKAY)
1069       return res;
1070   }
1071
1072   if(sgn == MP_ZPOS || s_mp_cmp_d(c, 0) == MP_EQ)
1073     SIGN(c) = MP_ZPOS;
1074   else
1075     SIGN(c) = sgn;
1076
1077   return MP_OKAY;
1078
1079 } /* end mp_mul() */
1080
1081 /* }}} */
1082
1083 /* {{{ mp_mul_2d(a, d, c) */
1084
1085 /*
1086   mp_mul_2d(a, d, c)
1087
1088   Compute c = a * 2^d.  a may be the same as c.
1089  */
1090
1091 mp_err mp_mul_2d(mp_int *a, mp_digit d, mp_int *c)
1092 {
1093   mp_err   res;
1094
1095   ARGCHK(a != NULL && c != NULL, MP_BADARG);
1096
1097   if((res = mp_copy(a, c)) != MP_OKAY)
1098     return res;
1099
1100   if(d == 0)
1101     return MP_OKAY;
1102
1103   return s_mp_mul_2d(c, d);
1104
1105 } /* end mp_mul() */
1106
1107 /* }}} */
1108
1109 /* {{{ mp_sqr(a, b) */
1110
1111 #if MP_SQUARE
1112 mp_err mp_sqr(mp_int *a, mp_int *b)
1113 {
1114   mp_err   res;
1115
1116   ARGCHK(a != NULL && b != NULL, MP_BADARG);
1117
1118   if((res = mp_copy(a, b)) != MP_OKAY)
1119     return res;
1120
1121   if((res = s_mp_sqr(b)) != MP_OKAY)
1122     return res;
1123
1124   SIGN(b) = MP_ZPOS;
1125
1126   return MP_OKAY;
1127
1128 } /* end mp_sqr() */
1129 #endif
1130
1131 /* }}} */
1132
1133 /* {{{ mp_div(a, b, q, r) */
1134
1135 /*
1136   mp_div(a, b, q, r)
1137
1138   Compute q = a / b and r = a mod b.  Input parameters may be re-used
1139   as output parameters.  If q or r is NULL, that portion of the
1140   computation will be discarded (although it will still be computed)
1141
1142   Pay no attention to the hacker behind the curtain.
1143  */
1144
1145 mp_err mp_div(mp_int *a, mp_int *b, mp_int *q, mp_int *r)
1146 {
1147   mp_err   res;
1148   mp_int   qtmp, rtmp;
1149   int      cmp;
1150
1151   ARGCHK(a != NULL && b != NULL, MP_BADARG);
1152
1153   if(mp_cmp_z(b) == MP_EQ)
1154     return MP_RANGE;
1155
1156   /* If a <= b, we can compute the solution without division, and
1157      avoid any memory allocation
1158    */
1159   if((cmp = s_mp_cmp(a, b)) < 0) {
1160     if(r) {
1161       if((res = mp_copy(a, r)) != MP_OKAY)
1162         return res;
1163     }
1164
1165     if(q)
1166       mp_zero(q);
1167
1168     return MP_OKAY;
1169
1170   } else if(cmp == 0) {
1171
1172     /* Set quotient to 1, with appropriate sign */
1173     if(q) {
1174       int qneg = (SIGN(a) != SIGN(b));
1175
1176       mp_set(q, 1);
1177       if(qneg)
1178         SIGN(q) = MP_NEG;
1179     }
1180
1181     if(r)
1182       mp_zero(r);
1183
1184     return MP_OKAY;
1185   }
1186
1187   /* If we get here, it means we actually have to do some division */
1188
1189   /* Set up some temporaries... */
1190   if((res = mp_init_copy(&qtmp, a)) != MP_OKAY)
1191     return res;
1192   if((res = mp_init_copy(&rtmp, b)) != MP_OKAY)
1193     goto CLEANUP;
1194
1195   if((res = s_mp_div(&qtmp, &rtmp)) != MP_OKAY)
1196     goto CLEANUP;
1197
1198   /* Compute the signs for the output  */
1199   SIGN(&rtmp) = SIGN(a); /* Sr = Sa              */
1200   if(SIGN(a) == SIGN(b))
1201     SIGN(&qtmp) = MP_ZPOS;  /* Sq = MP_ZPOS if Sa = Sb */
1202   else
1203     SIGN(&qtmp) = MP_NEG;   /* Sq = MP_NEG if Sa != Sb */
1204
1205   if(s_mp_cmp_d(&qtmp, 0) == MP_EQ)
1206     SIGN(&qtmp) = MP_ZPOS;
1207   if(s_mp_cmp_d(&rtmp, 0) == MP_EQ)
1208     SIGN(&rtmp) = MP_ZPOS;
1209
1210   /* Copy output, if it is needed      */
1211   if(q)
1212     s_mp_exch(&qtmp, q);
1213
1214   if(r)
1215     s_mp_exch(&rtmp, r);
1216
1217 CLEANUP:
1218   mp_clear(&rtmp);
1219   mp_clear(&qtmp);
1220
1221   return res;
1222
1223 } /* end mp_div() */
1224
1225 /* }}} */
1226
1227 /* {{{ mp_div_2d(a, d, q, r) */
1228
1229 mp_err mp_div_2d(mp_int *a, mp_digit d, mp_int *q, mp_int *r)
1230 {
1231   mp_err  res;
1232
1233   ARGCHK(a != NULL, MP_BADARG);
1234
1235   if(q) {
1236     if((res = mp_copy(a, q)) != MP_OKAY)
1237       return res;
1238
1239     s_mp_div_2d(q, d);
1240   }
1241
1242   if(r) {
1243     if((res = mp_copy(a, r)) != MP_OKAY)
1244       return res;
1245
1246     s_mp_mod_2d(r, d);
1247   }
1248
1249   return MP_OKAY;
1250
1251 } /* end mp_div_2d() */
1252
1253 /* }}} */
1254
1255 /* {{{ mp_expt(a, b, c) */
1256
1257 /*
1258   mp_expt(a, b, c)
1259
1260   Compute c = a ** b, that is, raise a to the b power.  Uses a
1261   standard iterative square-and-multiply technique.
1262  */
1263
1264 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c)
1265 {
1266   mp_int   s, x;
1267   mp_err   res;
1268   mp_digit d;
1269   int      dig, bit;
1270
1271   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1272
1273   if(mp_cmp_z(b) < 0)
1274     return MP_RANGE;
1275
1276   if((res = mp_init(&s)) != MP_OKAY)
1277     return res;
1278
1279   mp_set(&s, 1);
1280
1281   if((res = mp_init_copy(&x, a)) != MP_OKAY)
1282     goto X;
1283
1284   /* Loop over low-order digits in ascending order */
1285   for(dig = 0; dig < (USED(b) - 1); dig++) {
1286     d = DIGIT(b, dig);
1287
1288     /* Loop over bits of each non-maximal digit */
1289     for(bit = 0; bit < DIGIT_BIT; bit++) {
1290       if(d & 1) {
1291         if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1292           goto CLEANUP;
1293       }
1294
1295       d >>= 1;
1296
1297       if((res = s_mp_sqr(&x)) != MP_OKAY)
1298         goto CLEANUP;
1299     }
1300   }
1301
1302   /* Consider now the last digit... */
1303   d = DIGIT(b, dig);
1304
1305   while(d) {
1306     if(d & 1) {
1307       if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1308         goto CLEANUP;
1309     }
1310
1311     d >>= 1;
1312
1313     if((res = s_mp_sqr(&x)) != MP_OKAY)
1314       goto CLEANUP;
1315   }
1316
1317   if(mp_iseven(b))
1318     SIGN(&s) = SIGN(a);
1319
1320   res = mp_copy(&s, c);
1321
1322 CLEANUP:
1323   mp_clear(&x);
1324 X:
1325   mp_clear(&s);
1326
1327   return res;
1328
1329 } /* end mp_expt() */
1330
1331 /* }}} */
1332
1333 /* {{{ mp_2expt(a, k) */
1334
1335 /* Compute a = 2^k */
1336
1337 mp_err mp_2expt(mp_int *a, mp_digit k)
1338 {
1339   ARGCHK(a != NULL, MP_BADARG);
1340
1341   return s_mp_2expt(a, k);
1342
1343 } /* end mp_2expt() */
1344
1345 /* }}} */
1346
1347 /* {{{ mp_mod(a, m, c) */
1348
1349 /*
1350   mp_mod(a, m, c)
1351
1352   Compute c = a (mod m).  Result will always be 0 <= c < m.
1353  */
1354
1355 mp_err mp_mod(mp_int *a, mp_int *m, mp_int *c)
1356 {
1357   mp_err  res;
1358   int     mag;
1359
1360   ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
1361
1362   if(SIGN(m) == MP_NEG)
1363     return MP_RANGE;
1364
1365   /*
1366      If |a| > m, we need to divide to get the remainder and take the
1367      absolute value.
1368
1369      If |a| < m, we don't need to do any division, just copy and adjust
1370      the sign (if a is negative).
1371
1372      If |a| == m, we can simply set the result to zero.
1373
1374      This order is intended to minimize the average path length of the
1375      comparison chain on common workloads -- the most frequent cases are
1376      that |a| != m, so we do those first.
1377    */
1378   if((mag = s_mp_cmp(a, m)) > 0) {
1379     if((res = mp_div(a, m, NULL, c)) != MP_OKAY)
1380       return res;
1381
1382     if(SIGN(c) == MP_NEG) {
1383       if((res = mp_add(c, m, c)) != MP_OKAY)
1384         return res;
1385     }
1386
1387   } else if(mag < 0) {
1388     if((res = mp_copy(a, c)) != MP_OKAY)
1389       return res;
1390
1391     if(mp_cmp_z(a) < 0) {
1392       if((res = mp_add(c, m, c)) != MP_OKAY)
1393         return res;
1394
1395     }
1396
1397   } else {
1398     mp_zero(c);
1399
1400   }
1401
1402   return MP_OKAY;
1403
1404 } /* end mp_mod() */
1405
1406 /* }}} */
1407
1408 /* {{{ mp_mod_d(a, d, c) */
1409
1410 /*
1411   mp_mod_d(a, d, c)
1412
1413   Compute c = a (mod d).  Result will always be 0 <= c < d
1414  */
1415 mp_err mp_mod_d(mp_int *a, mp_digit d, mp_digit *c)
1416 {
1417   mp_err   res;
1418   mp_digit rem;
1419
1420   ARGCHK(a != NULL && c != NULL, MP_BADARG);
1421
1422   if(s_mp_cmp_d(a, d) > 0) {
1423     if((res = mp_div_d(a, d, NULL, &rem)) != MP_OKAY)
1424       return res;
1425
1426   } else {
1427     if(SIGN(a) == MP_NEG)
1428       rem = d - DIGIT(a, 0);
1429     else
1430       rem = DIGIT(a, 0);
1431   }
1432
1433   if(c)
1434     *c = rem;
1435
1436   return MP_OKAY;
1437
1438 } /* end mp_mod_d() */
1439
1440 /* }}} */
1441
1442 /* {{{ mp_sqrt(a, b) */
1443
1444 /*
1445   mp_sqrt(a, b)
1446
1447   Compute the integer square root of a, and store the result in b.
1448   Uses an integer-arithmetic version of Newton's iterative linear
1449   approximation technique to determine this value; the result has the
1450   following two properties:
1451
1452      b^2 <= a
1453      (b+1)^2 >= a
1454
1455   It is a range error to pass a negative value.
1456  */
1457 mp_err mp_sqrt(mp_int *a, mp_int *b)
1458 {
1459   mp_int   x, t;
1460   mp_err   res;
1461
1462   ARGCHK(a != NULL && b != NULL, MP_BADARG);
1463
1464   /* Cannot take square root of a negative value */
1465   if(SIGN(a) == MP_NEG)
1466     return MP_RANGE;
1467
1468   /* Special cases for zero and one, trivial     */
1469   if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ)
1470     return mp_copy(a, b);
1471
1472   /* Initialize the temporaries we'll use below  */
1473   if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
1474     return res;
1475
1476   /* Compute an initial guess for the iteration as a itself */
1477   if((res = mp_init_copy(&x, a)) != MP_OKAY)
1478     goto X;
1479
1480 s_mp_rshd(&x, (USED(&x)/2)+1);
1481 mp_add_d(&x, 1, &x);
1482
1483   for(;;) {
1484     /* t = (x * x) - a */
1485     mp_copy(&x, &t);      /* can't fail, t is big enough for original x */
1486     if((res = mp_sqr(&t, &t)) != MP_OKAY ||
1487        (res = mp_sub(&t, a, &t)) != MP_OKAY)
1488       goto CLEANUP;
1489
1490     /* t = t / 2x       */
1491     s_mp_mul_2(&x);
1492     if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)
1493       goto CLEANUP;
1494     s_mp_div_2(&x);
1495
1496     /* Terminate the loop, if the quotient is zero */
1497     if(mp_cmp_z(&t) == MP_EQ)
1498       break;
1499
1500     /* x = x - t       */
1501     if((res = mp_sub(&x, &t, &x)) != MP_OKAY)
1502       goto CLEANUP;
1503
1504   }
1505
1506   /* Copy result to output parameter */
1507   mp_sub_d(&x, 1, &x);
1508   s_mp_exch(&x, b);
1509
1510  CLEANUP:
1511   mp_clear(&x);
1512  X:
1513   mp_clear(&t);
1514
1515   return res;
1516
1517 } /* end mp_sqrt() */
1518
1519 /* }}} */
1520
1521 /* }}} */
1522
1523 /*------------------------------------------------------------------------*/
1524 /* {{{ Modular arithmetic */
1525
1526 #if MP_MODARITH
1527 /* {{{ mp_addmod(a, b, m, c) */
1528
1529 /*
1530   mp_addmod(a, b, m, c)
1531
1532   Compute c = (a + b) mod m
1533  */
1534
1535 mp_err mp_addmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1536 {
1537   mp_err  res;
1538
1539   ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1540
1541   if((res = mp_add(a, b, c)) != MP_OKAY)
1542     return res;
1543   if((res = mp_mod(c, m, c)) != MP_OKAY)
1544     return res;
1545
1546   return MP_OKAY;
1547
1548 }
1549
1550 /* }}} */
1551
1552 /* {{{ mp_submod(a, b, m, c) */
1553
1554 /*
1555   mp_submod(a, b, m, c)
1556
1557   Compute c = (a - b) mod m
1558  */
1559
1560 mp_err mp_submod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1561 {
1562   mp_err  res;
1563
1564   ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1565
1566   if((res = mp_sub(a, b, c)) != MP_OKAY)
1567     return res;
1568   if((res = mp_mod(c, m, c)) != MP_OKAY)
1569     return res;
1570
1571   return MP_OKAY;
1572
1573 }
1574
1575 /* }}} */
1576
1577 /* {{{ mp_mulmod(a, b, m, c) */
1578
1579 /*
1580   mp_mulmod(a, b, m, c)
1581
1582   Compute c = (a * b) mod m
1583  */
1584
1585 mp_err mp_mulmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1586 {
1587   mp_err  res;
1588
1589   ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1590
1591   if((res = mp_mul(a, b, c)) != MP_OKAY)
1592     return res;
1593   if((res = mp_mod(c, m, c)) != MP_OKAY)
1594     return res;
1595
1596   return MP_OKAY;
1597
1598 }
1599
1600 /* }}} */
1601
1602 /* {{{ mp_sqrmod(a, m, c) */
1603
1604 #if MP_SQUARE
1605 mp_err mp_sqrmod(mp_int *a, mp_int *m, mp_int *c)
1606 {
1607   mp_err  res;
1608
1609   ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
1610
1611   if((res = mp_sqr(a, c)) != MP_OKAY)
1612     return res;
1613   if((res = mp_mod(c, m, c)) != MP_OKAY)
1614     return res;
1615
1616   return MP_OKAY;
1617
1618 } /* end mp_sqrmod() */
1619 #endif
1620
1621 /* }}} */
1622
1623 /* {{{ mp_exptmod(a, b, m, c) */
1624
1625 /*
1626   mp_exptmod(a, b, m, c)
1627
1628   Compute c = (a ** b) mod m.  Uses a standard square-and-multiply
1629   method with modular reductions at each step. (This is basically the
1630   same code as mp_expt(), except for the addition of the reductions)
1631
1632   The modular reductions are done using Barrett's algorithm (see
1633   s_mp_reduce() below for details)
1634  */
1635
1636 mp_err mp_exptmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1637 {
1638   mp_int   s, x, mu;
1639   mp_err   res;
1640   mp_digit d, *db = DIGITS(b);
1641   mp_size  ub = USED(b);
1642   int      dig, bit;
1643
1644   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1645
1646   if(mp_cmp_z(b) < 0 || mp_cmp_z(m) <= 0)
1647     return MP_RANGE;
1648
1649   if((res = mp_init(&s)) != MP_OKAY)
1650     return res;
1651   if((res = mp_init_copy(&x, a)) != MP_OKAY)
1652     goto X;
1653   if((res = mp_mod(&x, m, &x)) != MP_OKAY ||
1654      (res = mp_init(&mu)) != MP_OKAY)
1655     goto MU;
1656
1657   mp_set(&s, 1);
1658
1659   /* mu = b^2k / m */
1660   s_mp_add_d(&mu, 1);
1661   s_mp_lshd(&mu, 2 * USED(m));
1662   if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY)
1663     goto CLEANUP;
1664
1665   /* Loop over digits of b in ascending order, except highest order */
1666   for(dig = 0; dig < (ub - 1); dig++) {
1667     d = *db++;
1668
1669     /* Loop over the bits of the lower-order digits */
1670     for(bit = 0; bit < DIGIT_BIT; bit++) {
1671       if(d & 1) {
1672         if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1673           goto CLEANUP;
1674         if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)
1675           goto CLEANUP;
1676       }
1677
1678       d >>= 1;
1679
1680       if((res = s_mp_sqr(&x)) != MP_OKAY)
1681         goto CLEANUP;
1682       if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)
1683         goto CLEANUP;
1684     }
1685   }
1686
1687   /* Now do the last digit... */
1688   d = *db;
1689
1690   while(d) {
1691     if(d & 1) {
1692       if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1693         goto CLEANUP;
1694       if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)
1695         goto CLEANUP;
1696     }
1697
1698     d >>= 1;
1699
1700     if((res = s_mp_sqr(&x)) != MP_OKAY)
1701       goto CLEANUP;
1702     if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)
1703       goto CLEANUP;
1704   }
1705
1706   s_mp_exch(&s, c);
1707
1708  CLEANUP:
1709   mp_clear(&mu);
1710  MU:
1711   mp_clear(&x);
1712  X:
1713   mp_clear(&s);
1714
1715   return res;
1716
1717 } /* end mp_exptmod() */
1718
1719 /* }}} */
1720
1721 /* {{{ mp_exptmod_d(a, d, m, c) */
1722
1723 mp_err mp_exptmod_d(mp_int *a, mp_digit d, mp_int *m, mp_int *c)
1724 {
1725   mp_int   s, x;
1726   mp_err   res;
1727
1728   ARGCHK(a != NULL && c != NULL, MP_BADARG);
1729
1730   if((res = mp_init(&s)) != MP_OKAY)
1731     return res;
1732   if((res = mp_init_copy(&x, a)) != MP_OKAY)
1733     goto X;
1734
1735   mp_set(&s, 1);
1736
1737   while(d != 0) {
1738     if(d & 1) {
1739       if((res = s_mp_mul(&s, &x)) != MP_OKAY ||
1740          (res = mp_mod(&s, m, &s)) != MP_OKAY)
1741         goto CLEANUP;
1742     }
1743
1744     d /= 2;
1745
1746     if((res = s_mp_sqr(&x)) != MP_OKAY ||
1747        (res = mp_mod(&x, m, &x)) != MP_OKAY)
1748       goto CLEANUP;
1749   }
1750
1751   s_mp_exch(&s, c);
1752
1753 CLEANUP:
1754   mp_clear(&x);
1755 X:
1756   mp_clear(&s);
1757
1758   return res;
1759
1760 } /* end mp_exptmod_d() */
1761
1762 /* }}} */
1763 #endif /* if MP_MODARITH */
1764
1765 /* }}} */
1766
1767 /*------------------------------------------------------------------------*/
1768 /* {{{ Comparison functions */
1769
1770 /* {{{ mp_cmp_z(a) */
1771
1772 /*
1773   mp_cmp_z(a)
1774
1775   Compare a <=> 0.  Returns <0 if a<0, 0 if a=0, >0 if a>0.
1776  */
1777
1778 int    mp_cmp_z(mp_int *a)
1779 {
1780   if(SIGN(a) == MP_NEG)
1781     return MP_LT;
1782   else if(USED(a) == 1 && DIGIT(a, 0) == 0)
1783     return MP_EQ;
1784   else
1785     return MP_GT;
1786
1787 } /* end mp_cmp_z() */
1788
1789 /* }}} */
1790
1791 /* {{{ mp_cmp_d(a, d) */
1792
1793 /*
1794   mp_cmp_d(a, d)
1795
1796   Compare a <=> d.  Returns <0 if a<d, 0 if a=d, >0 if a>d
1797  */
1798
1799 int    mp_cmp_d(mp_int *a, mp_digit d)
1800 {
1801   ARGCHK(a != NULL, MP_EQ);
1802
1803   if(SIGN(a) == MP_NEG)
1804     return MP_LT;
1805
1806   return s_mp_cmp_d(a, d);
1807
1808 } /* end mp_cmp_d() */
1809
1810 /* }}} */
1811
1812 /* {{{ mp_cmp(a, b) */
1813
1814 int    mp_cmp(mp_int *a, mp_int *b)
1815 {
1816   ARGCHK(a != NULL && b != NULL, MP_EQ);
1817
1818   if(SIGN(a) == SIGN(b)) {
1819     int  mag;
1820
1821     if((mag = s_mp_cmp(a, b)) == MP_EQ)
1822       return MP_EQ;
1823
1824     if(SIGN(a) == MP_ZPOS)
1825       return mag;
1826     else
1827       return -mag;
1828
1829   } else if(SIGN(a) == MP_ZPOS) {
1830     return MP_GT;
1831   } else {
1832     return MP_LT;
1833   }
1834
1835 } /* end mp_cmp() */
1836
1837 /* }}} */
1838
1839 /* {{{ mp_cmp_mag(a, b) */
1840
1841 /*
1842   mp_cmp_mag(a, b)
1843
1844   Compares |a| <=> |b|, and returns an appropriate comparison result
1845  */
1846
1847 int    mp_cmp_mag(mp_int *a, mp_int *b)
1848 {
1849   ARGCHK(a != NULL && b != NULL, MP_EQ);
1850
1851   return s_mp_cmp(a, b);
1852
1853 } /* end mp_cmp_mag() */
1854
1855 /* }}} */
1856
1857 /* {{{ mp_cmp_int(a, z) */
1858
1859 /*
1860   This just converts z to an mp_int, and uses the existing comparison
1861   routines.  This is sort of inefficient, but it's not clear to me how
1862   frequently this wil get used anyway.  For small positive constants,
1863   you can always use mp_cmp_d(), and for zero, there is mp_cmp_z().
1864  */
1865 int    mp_cmp_int(mp_int *a, long z)
1866 {
1867   mp_int  tmp;
1868   int     out;
1869
1870   ARGCHK(a != NULL, MP_EQ);
1871
1872   mp_init(&tmp); mp_set_int(&tmp, z);
1873   out = mp_cmp(a, &tmp);
1874   mp_clear(&tmp);
1875
1876   return out;
1877
1878 } /* end mp_cmp_int() */
1879
1880 /* }}} */
1881
1882 /* {{{ mp_isodd(a) */
1883
1884 /*
1885   mp_isodd(a)
1886
1887   Returns a true (non-zero) value if a is odd, false (zero) otherwise.
1888  */
1889 int    mp_isodd(mp_int *a)
1890 {
1891   ARGCHK(a != NULL, 0);
1892
1893   return (DIGIT(a, 0) & 1);
1894
1895 } /* end mp_isodd() */
1896
1897 /* }}} */
1898
1899 /* {{{ mp_iseven(a) */
1900
1901 int    mp_iseven(mp_int *a)
1902 {
1903   return !mp_isodd(a);
1904
1905 } /* end mp_iseven() */
1906
1907 /* }}} */
1908
1909 /* }}} */
1910
1911 /*------------------------------------------------------------------------*/
1912 /* {{{ Number theoretic functions */
1913
1914 #if MP_NUMTH
1915 /* {{{ mp_gcd(a, b, c) */
1916
1917 /*
1918   Like the old mp_gcd() function, except computes the GCD using the
1919   binary algorithm due to Josef Stein in 1961 (via Knuth).
1920  */
1921 mp_err mp_gcd(mp_int *a, mp_int *b, mp_int *c)
1922 {
1923   mp_err   res;
1924   mp_int   u, v, t;
1925   mp_size  k = 0;
1926
1927   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1928
1929   if(mp_cmp_z(a) == MP_EQ && mp_cmp_z(b) == MP_EQ)
1930       return MP_RANGE;
1931   if(mp_cmp_z(a) == MP_EQ) {
1932     return mp_copy(b, c);
1933   } else if(mp_cmp_z(b) == MP_EQ) {
1934     return mp_copy(a, c);
1935   }
1936
1937   if((res = mp_init(&t)) != MP_OKAY)
1938     return res;
1939   if((res = mp_init_copy(&u, a)) != MP_OKAY)
1940     goto U;
1941   if((res = mp_init_copy(&v, b)) != MP_OKAY)
1942     goto V;
1943
1944   SIGN(&u) = MP_ZPOS;
1945   SIGN(&v) = MP_ZPOS;
1946
1947   /* Divide out common factors of 2 until at least 1 of a, b is even */
1948   while(mp_iseven(&u) && mp_iseven(&v)) {
1949     s_mp_div_2(&u);
1950     s_mp_div_2(&v);
1951     ++k;
1952   }
1953
1954   /* Initialize t */
1955   if(mp_isodd(&u)) {
1956     if((res = mp_copy(&v, &t)) != MP_OKAY)
1957       goto CLEANUP;
1958
1959     /* t = -v */
1960     if(SIGN(&v) == MP_ZPOS)
1961       SIGN(&t) = MP_NEG;
1962     else
1963       SIGN(&t) = MP_ZPOS;
1964
1965   } else {
1966     if((res = mp_copy(&u, &t)) != MP_OKAY)
1967       goto CLEANUP;
1968
1969   }
1970
1971   for(;;) {
1972     while(mp_iseven(&t)) {
1973       s_mp_div_2(&t);
1974     }
1975
1976     if(mp_cmp_z(&t) == MP_GT) {
1977       if((res = mp_copy(&t, &u)) != MP_OKAY)
1978         goto CLEANUP;
1979
1980     } else {
1981       if((res = mp_copy(&t, &v)) != MP_OKAY)
1982         goto CLEANUP;
1983
1984       /* v = -t */
1985       if(SIGN(&t) == MP_ZPOS)
1986         SIGN(&v) = MP_NEG;
1987       else
1988         SIGN(&v) = MP_ZPOS;
1989     }
1990
1991     if((res = mp_sub(&u, &v, &t)) != MP_OKAY)
1992       goto CLEANUP;
1993
1994     if(s_mp_cmp_d(&t, 0) == MP_EQ)
1995       break;
1996   }
1997
1998   s_mp_2expt(&v, k);       /* v = 2^k   */
1999   res = mp_mul(&u, &v, c); /* c = u * v */
2000
2001  CLEANUP:
2002   mp_clear(&v);
2003  V:
2004   mp_clear(&u);
2005  U:
2006   mp_clear(&t);
2007
2008   return res;
2009
2010 } /* end mp_bgcd() */
2011
2012 /* }}} */
2013
2014 /* {{{ mp_lcm(a, b, c) */
2015
2016 /* We compute the least common multiple using the rule:
2017
2018    ab = [a, b](a, b)
2019
2020    ... by computing the product, and dividing out the gcd.
2021  */
2022
2023 mp_err mp_lcm(mp_int *a, mp_int *b, mp_int *c)
2024 {
2025   mp_int  gcd, prod;
2026   mp_err  res;
2027
2028   ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
2029
2030   /* Set up temporaries */
2031   if((res = mp_init(&gcd)) != MP_OKAY)
2032     return res;
2033   if((res = mp_init(&prod)) != MP_OKAY)
2034     goto GCD;
2035
2036   if((res = mp_mul(a, b, &prod)) != MP_OKAY)
2037     goto CLEANUP;
2038   if((res = mp_gcd(a, b, &gcd)) != MP_OKAY)
2039     goto CLEANUP;
2040
2041   res = mp_div(&prod, &gcd, c, NULL);
2042
2043  CLEANUP:
2044   mp_clear(&prod);
2045  GCD:
2046   mp_clear(&gcd);
2047
2048   return res;
2049
2050 } /* end mp_lcm() */
2051
2052 /* }}} */
2053
2054 /* {{{ mp_xgcd(a, b, g, x, y) */
2055
2056 /*
2057   mp_xgcd(a, b, g, x, y)
2058
2059   Compute g = (a, b) and values x and y satisfying Bezout's identity
2060   (that is, ax + by = g).  This uses the extended binary GCD algorithm
2061   based on the Stein algorithm used for mp_gcd()
2062  */
2063
2064 mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
2065 {
2066   mp_int   gx, xc, yc, u, v, A, B, C, D;
2067   mp_int  *clean[9];
2068   mp_err   res;
2069   int      last = -1;
2070
2071   if(mp_cmp_z(b) == 0)
2072     return MP_RANGE;
2073
2074   /* Initialize all these variables we need */
2075   if((res = mp_init(&u)) != MP_OKAY) goto CLEANUP;
2076   clean[++last] = &u;
2077   if((res = mp_init(&v)) != MP_OKAY) goto CLEANUP;
2078   clean[++last] = &v;
2079   if((res = mp_init(&gx)) != MP_OKAY) goto CLEANUP;
2080   clean[++last] = &gx;
2081   if((res = mp_init(&A)) != MP_OKAY) goto CLEANUP;
2082   clean[++last] = &A;
2083   if((res = mp_init(&B)) != MP_OKAY) goto CLEANUP;
2084   clean[++last] = &B;
2085   if((res = mp_init(&C)) != MP_OKAY) goto CLEANUP;
2086   clean[++last] = &C;
2087   if((res = mp_init(&D)) != MP_OKAY) goto CLEANUP;
2088   clean[++last] = &D;
2089   if((res = mp_init_copy(&xc, a)) != MP_OKAY) goto CLEANUP;
2090   clean[++last] = &xc;
2091   mp_abs(&xc, &xc);
2092   if((res = mp_init_copy(&yc, b)) != MP_OKAY) goto CLEANUP;
2093   clean[++last] = &yc;
2094   mp_abs(&yc, &yc);
2095
2096   mp_set(&gx, 1);
2097
2098   /* Divide by two until at least one of them is even */
2099   while(mp_iseven(&xc) && mp_iseven(&yc)) {
2100     s_mp_div_2(&xc);
2101     s_mp_div_2(&yc);
2102     if((res = s_mp_mul_2(&gx)) != MP_OKAY)
2103       goto CLEANUP;
2104   }
2105
2106   mp_copy(&xc, &u);
2107   mp_copy(&yc, &v);
2108   mp_set(&A, 1); mp_set(&D, 1);
2109
2110   /* Loop through binary GCD algorithm */
2111   for(;;) {
2112     while(mp_iseven(&u)) {
2113       s_mp_div_2(&u);
2114
2115       if(mp_iseven(&A) && mp_iseven(&B)) {
2116         s_mp_div_2(&A); s_mp_div_2(&B);
2117       } else {
2118         if((res = mp_add(&A, &yc, &A)) != MP_OKAY) goto CLEANUP;
2119         s_mp_div_2(&A);
2120         if((res = mp_sub(&B, &xc, &B)) != MP_OKAY) goto CLEANUP;
2121         s_mp_div_2(&B);
2122       }
2123     }
2124
2125     while(mp_iseven(&v)) {
2126       s_mp_div_2(&v);
2127
2128       if(mp_iseven(&C) && mp_iseven(&D)) {
2129         s_mp_div_2(&C); s_mp_div_2(&D);
2130       } else {
2131         if((res = mp_add(&C, &yc, &C)) != MP_OKAY) goto CLEANUP;
2132         s_mp_div_2(&C);
2133         if((res = mp_sub(&D, &xc, &D)) != MP_OKAY) goto CLEANUP;
2134         s_mp_div_2(&D);
2135       }
2136     }
2137
2138     if(mp_cmp(&u, &v) >= 0) {
2139       if((res = mp_sub(&u, &v, &u)) != MP_OKAY) goto CLEANUP;
2140       if((res = mp_sub(&A, &C, &A)) != MP_OKAY) goto CLEANUP;
2141       if((res = mp_sub(&B, &D, &B)) != MP_OKAY) goto CLEANUP;
2142
2143     } else {
2144       if((res = mp_sub(&v, &u, &v)) != MP_OKAY) goto CLEANUP;
2145       if((res = mp_sub(&C, &A, &C)) != MP_OKAY) goto CLEANUP;
2146       if((res = mp_sub(&D, &B, &D)) != MP_OKAY) goto CLEANUP;
2147
2148     }
2149
2150     /* If we're done, copy results to output */
2151     if(mp_cmp_z(&u) == 0) {
2152       if(x)
2153         if((res = mp_copy(&C, x)) != MP_OKAY) goto CLEANUP;
2154
2155       if(y)
2156         if((res = mp_copy(&D, y)) != MP_OKAY) goto CLEANUP;
2157
2158       if(g)
2159         if((res = mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP;
2160
2161       break;
2162     }
2163   }
2164
2165  CLEANUP:
2166   while(last >= 0)
2167     mp_clear(clean[last--]);
2168
2169   return res;
2170
2171 } /* end mp_xgcd() */
2172
2173 /* }}} */
2174
2175 /* {{{ mp_invmod(a, m, c) */
2176
2177 /*
2178   mp_invmod(a, m, c)
2179
2180   Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
2181   This is equivalent to the question of whether (a, m) = 1.  If not,
2182   MP_UNDEF is returned, and there is no inverse.
2183  */
2184
2185 mp_err mp_invmod(mp_int *a, mp_int *m, mp_int *c)
2186 {
2187   mp_int  g, x;
2188   mp_err  res;
2189
2190   ARGCHK(a && m && c, MP_BADARG);
2191
2192   if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
2193     return MP_RANGE;
2194
2195   if((res = mp_init(&g)) != MP_OKAY)
2196     return res;
2197   if((res = mp_init(&x)) != MP_OKAY)
2198     goto X;
2199
2200   if((res = mp_xgcd(a, m, &g, &x, NULL)) != MP_OKAY)
2201     goto CLEANUP;
2202
2203   if(mp_cmp_d(&g, 1) != MP_EQ) {
2204     res = MP_UNDEF;
2205     goto CLEANUP;
2206   }
2207
2208   res = mp_mod(&x, m, c);
2209   SIGN(c) = SIGN(a);
2210
2211 CLEANUP:
2212   mp_clear(&x);
2213 X:
2214   mp_clear(&g);
2215
2216   return res;
2217
2218 } /* end mp_invmod() */
2219
2220 /* }}} */
2221 #endif /* if MP_NUMTH */
2222
2223 /* }}} */
2224
2225 /*------------------------------------------------------------------------*/
2226 /* {{{ mp_print(mp, ofp) */
2227
2228 #if MP_IOFUNC
2229 /*
2230   mp_print(mp, ofp)
2231
2232   Print a textual representation of the given mp_int on the output
2233   stream 'ofp'.  Output is generated using the internal radix.
2234  */
2235
2236 void   mp_print(mp_int *mp, FILE *ofp)
2237 {
2238   int   ix;
2239
2240   if(mp == NULL || ofp == NULL)
2241     return;
2242
2243   fputc((SIGN(mp) == MP_NEG) ? '-' : '+', ofp);
2244
2245   for(ix = USED(mp) - 1; ix >= 0; ix--) {
2246     fprintf(ofp, DIGIT_FMT, DIGIT(mp, ix));
2247   }
2248
2249 } /* end mp_print() */
2250
2251 #endif /* if MP_IOFUNC */
2252
2253 /* }}} */
2254
2255 /*------------------------------------------------------------------------*/
2256 /* {{{ More I/O Functions */
2257
2258 /* {{{ mp_read_signed_bin(mp, str, len) */
2259
2260 /*
2261    mp_read_signed_bin(mp, str, len)
2262
2263    Read in a raw value (base 256) into the given mp_int
2264  */
2265
2266 mp_err  mp_read_signed_bin(mp_int *mp, unsigned char *str, int len)
2267 {
2268   mp_err         res;
2269
2270   ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
2271
2272   if((res = mp_read_unsigned_bin(mp, str + 1, len - 1)) == MP_OKAY) {
2273     /* Get sign from first byte */
2274     if(str[0])
2275       SIGN(mp) = MP_NEG;
2276     else
2277       SIGN(mp) = MP_ZPOS;
2278   }
2279
2280   return res;
2281
2282 } /* end mp_read_signed_bin() */
2283
2284 /* }}} */
2285
2286 /* {{{ mp_signed_bin_size(mp) */
2287
2288 int    mp_signed_bin_size(mp_int *mp)
2289 {
2290   ARGCHK(mp != NULL, 0);
2291
2292   return mp_unsigned_bin_size(mp) + 1;
2293
2294 } /* end mp_signed_bin_size() */
2295
2296 /* }}} */
2297
2298 /* {{{ mp_to_signed_bin(mp, str) */
2299
2300 mp_err mp_to_signed_bin(mp_int *mp, unsigned char *str)
2301 {
2302   ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2303
2304   /* Caller responsible for allocating enough memory (use mp_raw_size(mp)) */
2305   str[0] = (char)SIGN(mp);
2306
2307   return mp_to_unsigned_bin(mp, str + 1);
2308
2309 } /* end mp_to_signed_bin() */
2310
2311 /* }}} */
2312
2313 /* {{{ mp_read_unsigned_bin(mp, str, len) */
2314
2315 /*
2316   mp_read_unsigned_bin(mp, str, len)
2317
2318   Read in an unsigned value (base 256) into the given mp_int
2319  */
2320
2321 mp_err  mp_read_unsigned_bin(mp_int *mp, unsigned char *str, int len)
2322 {
2323   int     ix;
2324   mp_err  res;
2325
2326   ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
2327
2328   mp_zero(mp);
2329
2330   for(ix = 0; ix < len; ix++) {
2331     if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
2332       return res;
2333
2334     if((res = mp_add_d(mp, str[ix], mp)) != MP_OKAY)
2335       return res;
2336   }
2337
2338   return MP_OKAY;
2339
2340 } /* end mp_read_unsigned_bin() */
2341
2342 /* }}} */
2343
2344 /* {{{ mp_unsigned_bin_size(mp) */
2345
2346 int     mp_unsigned_bin_size(mp_int *mp)
2347 {
2348   mp_digit   topdig;
2349   int        count;
2350
2351   ARGCHK(mp != NULL, 0);
2352
2353   /* Special case for the value zero */
2354   if(USED(mp) == 1 && DIGIT(mp, 0) == 0)
2355     return 1;
2356
2357   count = (USED(mp) - 1) * sizeof(mp_digit);
2358   topdig = DIGIT(mp, USED(mp) - 1);
2359
2360   while(topdig != 0) {
2361     ++count;
2362     topdig >>= CHAR_BIT;
2363   }
2364
2365   return count;
2366
2367 } /* end mp_unsigned_bin_size() */
2368
2369 /* }}} */
2370
2371 /* {{{ mp_to_unsigned_bin(mp, str) */
2372
2373 mp_err mp_to_unsigned_bin(mp_int *mp, unsigned char *str)
2374 {
2375   mp_digit      *dp, *end, d;
2376   unsigned char *spos;
2377
2378   ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2379
2380   dp = DIGITS(mp);
2381   end = dp + USED(mp) - 1;
2382   spos = str;
2383
2384   /* Special case for zero, quick test */
2385   if(dp == end && *dp == 0) {
2386     *str = '\0';
2387     return MP_OKAY;
2388   }
2389
2390   /* Generate digits in reverse order */
2391   while(dp < end) {
2392     int      ix;
2393
2394     d = *dp;
2395     for(ix = 0; ix < sizeof(mp_digit); ++ix) {
2396       *spos = d & UCHAR_MAX;
2397       d >>= CHAR_BIT;
2398       ++spos;
2399     }
2400
2401     ++dp;
2402   }
2403
2404   /* Now handle last digit specially, high order zeroes are not written */
2405   d = *end;
2406   while(d != 0) {
2407     *spos = d & UCHAR_MAX;
2408     d >>= CHAR_BIT;
2409     ++spos;
2410   }
2411
2412   /* Reverse everything to get digits in the correct order */
2413   while(--spos > str) {
2414     unsigned char t = *str;
2415     *str = *spos;
2416     *spos = t;
2417
2418     ++str;
2419   }
2420
2421   return MP_OKAY;
2422
2423 } /* end mp_to_unsigned_bin() */
2424
2425 /* }}} */
2426
2427 /* {{{ mp_count_bits(mp) */
2428
2429 int    mp_count_bits(mp_int *mp)
2430 {
2431   int      len;
2432   mp_digit d;
2433
2434   ARGCHK(mp != NULL, MP_BADARG);
2435
2436   len = DIGIT_BIT * (USED(mp) - 1);
2437   d = DIGIT(mp, USED(mp) - 1);
2438
2439   while(d != 0) {
2440     ++len;
2441     d >>= 1;
2442   }
2443
2444   return len;
2445
2446 } /* end mp_count_bits() */
2447
2448 /* }}} */
2449
2450 /* {{{ mp_read_radix(mp, str, radix) */
2451
2452 /*
2453   mp_read_radix(mp, str, radix)
2454
2455   Read an integer from the given string, and set mp to the resulting
2456   value.  The input is presumed to be in base 10.  Leading non-digit
2457   characters are ignored, and the function reads until a non-digit
2458   character or the end of the string.
2459  */
2460
2461 mp_err  mp_read_radix(mp_int *mp, unsigned char *str, int radix)
2462 {
2463   int     ix = 0, val = 0;
2464   mp_err  res;
2465   mp_sign sig = MP_ZPOS;
2466
2467   ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX,
2468          MP_BADARG);
2469
2470   mp_zero(mp);
2471
2472   /* Skip leading non-digit characters until a digit or '-' or '+' */
2473   while(str[ix] &&
2474         (s_mp_tovalue(str[ix], radix) < 0) &&
2475         str[ix] != '-' &&
2476         str[ix] != '+') {
2477     ++ix;
2478   }
2479
2480   if(str[ix] == '-') {
2481     sig = MP_NEG;
2482     ++ix;
2483   } else if(str[ix] == '+') {
2484     sig = MP_ZPOS; /* this is the default anyway... */
2485     ++ix;
2486   }
2487
2488   while((val = s_mp_tovalue(str[ix], radix)) >= 0) {
2489     if((res = s_mp_mul_d(mp, radix)) != MP_OKAY)
2490       return res;
2491     if((res = s_mp_add_d(mp, val)) != MP_OKAY)
2492       return res;
2493     ++ix;
2494   }
2495
2496   if(s_mp_cmp_d(mp, 0) == MP_EQ)
2497     SIGN(mp) = MP_ZPOS;
2498   else
2499     SIGN(mp) = sig;
2500
2501   return MP_OKAY;
2502
2503 } /* end mp_read_radix() */
2504
2505 /* }}} */
2506
2507 /* {{{ mp_radix_size(mp, radix) */
2508
2509 int    mp_radix_size(mp_int *mp, int radix)
2510 {
2511   int  len;
2512   ARGCHK(mp != NULL, 0);
2513
2514   len = s_mp_outlen(mp_count_bits(mp), radix) + 1; /* for NUL terminator */
2515
2516   if(mp_cmp_z(mp) < 0)
2517     ++len; /* for sign */
2518
2519   return len;
2520
2521 } /* end mp_radix_size() */
2522
2523 /* }}} */
2524
2525 /* {{{ mp_value_radix_size(num, qty, radix) */
2526
2527 /* num = number of digits
2528    qty = number of bits per digit
2529    radix = target base
2530
2531    Return the number of digits in the specified radix that would be
2532    needed to express 'num' digits of 'qty' bits each.
2533  */
2534 int    mp_value_radix_size(int num, int qty, int radix)
2535 {
2536   ARGCHK(num >= 0 && qty > 0 && radix >= 2 && radix <= MAX_RADIX, 0);
2537
2538   return s_mp_outlen(num * qty, radix);
2539
2540 } /* end mp_value_radix_size() */
2541
2542 /* }}} */
2543
2544 /* {{{ mp_toradix(mp, str, radix) */
2545
2546 mp_err mp_toradix(mp_int *mp, unsigned char *str, int radix)
2547 {
2548   int  ix, pos = 0;
2549
2550   ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2551   ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);
2552
2553   if(mp_cmp_z(mp) == MP_EQ) {
2554     str[0] = '0';
2555     str[1] = '\0';
2556   } else {
2557     mp_err   res;
2558     mp_int   tmp;
2559     mp_sign  sgn;
2560     mp_digit rem, rdx = (mp_digit)radix;
2561     char     ch;
2562
2563     if((res = mp_init_copy(&tmp, mp)) != MP_OKAY)
2564       return res;
2565
2566     /* Save sign for later, and take absolute value */
2567     sgn = SIGN(&tmp); SIGN(&tmp) = MP_ZPOS;
2568
2569     /* Generate output digits in reverse order      */
2570     while(mp_cmp_z(&tmp) != 0) {
2571       if((res = s_mp_div_d(&tmp, rdx, &rem)) != MP_OKAY) {
2572         mp_clear(&tmp);
2573         return res;
2574       }
2575
2576       /* Generate digits, use capital letters */
2577       ch = s_mp_todigit(rem, radix, 0);
2578
2579       str[pos++] = ch;
2580     }
2581
2582     /* Add - sign if original value was negative */
2583     if(sgn == MP_NEG)
2584       str[pos++] = '-';
2585
2586     /* Add trailing NUL to end the string        */
2587     str[pos--] = '\0';
2588
2589     /* Reverse the digits and sign indicator     */
2590     ix = 0;
2591     while(ix < pos) {
2592       char tmp = str[ix];
2593
2594       str[ix] = str[pos];
2595       str[pos] = tmp;
2596       ++ix;
2597       --pos;
2598     }
2599
2600     mp_clear(&tmp);
2601   }
2602
2603   return MP_OKAY;
2604
2605 } /* end mp_toradix() */
2606
2607 /* }}} */
2608
2609 /* {{{ mp_char2value(ch, r) */
2610
2611 int    mp_char2value(char ch, int r)
2612 {
2613   return s_mp_tovalue(ch, r);
2614
2615 } /* end mp_tovalue() */
2616
2617 /* }}} */
2618
2619 /* }}} */
2620
2621 /* {{{ mp_strerror(ec) */
2622
2623 /*
2624   mp_strerror(ec)
2625
2626   Return a string describing the meaning of error code 'ec'.  The
2627   string returned is allocated in static memory, so the caller should
2628   not attempt to modify or free the memory associated with this
2629   string.
2630  */
2631 const char  *mp_strerror(mp_err ec)
2632 {
2633   int   aec = (ec < 0) ? -ec : ec;
2634
2635   /* Code values are negative, so the senses of these comparisons
2636      are accurate */
2637   if(ec < MP_LAST_CODE || ec > MP_OKAY) {
2638     return mp_err_string[0];  /* unknown error code */
2639   } else {
2640     return mp_err_string[aec + 1];
2641   }
2642
2643 } /* end mp_strerror() */
2644
2645 /* }}} */
2646
2647 /*========================================================================*/
2648 /*------------------------------------------------------------------------*/
2649 /* Static function definitions (internal use only)                        */
2650
2651 /* {{{ Memory management */
2652
2653 /* {{{ s_mp_grow(mp, min) */
2654
2655 /* Make sure there are at least 'min' digits allocated to mp              */
2656 mp_err   s_mp_grow(mp_int *mp, mp_size min)
2657 {
2658   if(min > ALLOC(mp)) {
2659     mp_digit   *tmp;
2660
2661     /* Set min to next nearest default precision block size */
2662     min = ((min + (s_mp_defprec - 1)) / s_mp_defprec) * s_mp_defprec;
2663
2664     if((tmp = s_mp_alloc(min, sizeof(mp_digit))) == NULL)
2665       return MP_MEM;
2666
2667     s_mp_copy(DIGITS(mp), tmp, USED(mp));
2668
2669 #if MP_CRYPTO
2670     s_mp_setz(DIGITS(mp), ALLOC(mp));
2671 #endif
2672     s_mp_free(DIGITS(mp));
2673     DIGITS(mp) = tmp;
2674     ALLOC(mp) = min;
2675   }
2676
2677   return MP_OKAY;
2678
2679 } /* end s_mp_grow() */
2680
2681 /* }}} */
2682
2683 /* {{{ s_mp_pad(mp, min) */
2684
2685 /* Make sure the used size of mp is at least 'min', growing if needed     */
2686 mp_err   s_mp_pad(mp_int *mp, mp_size min)
2687 {
2688   if(min > USED(mp)) {
2689     mp_err  res;
2690
2691     /* Make sure there is room to increase precision  */
2692     if(min > ALLOC(mp) && (res = s_mp_grow(mp, min)) != MP_OKAY)
2693       return res;
2694
2695     /* Increase precision; should already be 0-filled */
2696     USED(mp) = min;
2697   }
2698
2699   return MP_OKAY;
2700
2701 } /* end s_mp_pad() */
2702
2703 /* }}} */
2704
2705 /* {{{ s_mp_setz(dp, count) */
2706
2707 #if MP_MACRO == 0
2708 /* Set 'count' digits pointed to by dp to be zeroes                       */
2709 void s_mp_setz(mp_digit *dp, mp_size count)
2710 {
2711 #if MP_MEMSET == 0
2712   int  ix;
2713
2714   for(ix = 0; ix < count; ix++)
2715     dp[ix] = 0;
2716 #else
2717   memset(dp, 0, count * sizeof(mp_digit));
2718 #endif
2719
2720 } /* end s_mp_setz() */
2721 #endif
2722
2723 /* }}} */
2724
2725 /* {{{ s_mp_copy(sp, dp, count) */
2726
2727 #if MP_MACRO == 0
2728 /* Copy 'count' digits from sp to dp                                      */
2729 void s_mp_copy(mp_digit *sp, mp_digit *dp, mp_size count)
2730 {
2731 #if MP_MEMCPY == 0
2732   int  ix;
2733
2734   for(ix = 0; ix < count; ix++)
2735     dp[ix] = sp[ix];
2736 #else
2737   memcpy(dp, sp, count * sizeof(mp_digit));
2738 #endif
2739
2740 } /* end s_mp_copy() */
2741 #endif
2742
2743 /* }}} */
2744
2745 /* {{{ s_mp_alloc(nb, ni) */
2746
2747 #if MP_MACRO == 0
2748 /* Allocate ni records of nb bytes each, and return a pointer to that     */
2749 void    *s_mp_alloc(size_t nb, size_t ni)
2750 {
2751   return calloc(nb, ni);
2752
2753 } /* end s_mp_alloc() */
2754 #endif
2755
2756 /* }}} */
2757
2758 /* {{{ s_mp_free(ptr) */
2759
2760 #if MP_MACRO == 0
2761 /* Free the memory pointed to by ptr                                      */
2762 void     s_mp_free(void *ptr)
2763 {
2764   if(ptr)
2765     free(ptr);
2766
2767 } /* end s_mp_free() */
2768 #endif
2769
2770 /* }}} */
2771
2772 /* {{{ s_mp_clamp(mp) */
2773
2774 /* Remove leading zeroes from the given value                             */
2775 void     s_mp_clamp(mp_int *mp)
2776 {
2777   mp_size   du = USED(mp);
2778   mp_digit *zp = DIGITS(mp) + du - 1;
2779
2780   while(du > 1 && !*zp--)
2781     --du;
2782
2783   USED(mp) = du;
2784
2785 } /* end s_mp_clamp() */
2786
2787
2788 /* }}} */
2789
2790 /* {{{ s_mp_exch(a, b) */
2791
2792 /* Exchange the data for a and b; (b, a) = (a, b)                         */
2793 void     s_mp_exch(mp_int *a, mp_int *b)
2794 {
2795   mp_int   tmp;
2796
2797   tmp = *a;
2798   *a = *b;
2799   *b = tmp;
2800
2801 } /* end s_mp_exch() */
2802
2803 /* }}} */
2804
2805 /* }}} */
2806
2807 /* {{{ Arithmetic helpers */
2808
2809 /* {{{ s_mp_lshd(mp, p) */
2810
2811 /*
2812    Shift mp leftward by p digits, growing if needed, and zero-filling
2813    the in-shifted digits at the right end.  This is a convenient
2814    alternative to multiplication by powers of the radix
2815  */
2816
2817 mp_err   s_mp_lshd(mp_int *mp, mp_size p)
2818 {
2819   mp_err   res;
2820   mp_size  pos;
2821   mp_digit *dp;
2822   int     ix;
2823
2824   if(p == 0)
2825     return MP_OKAY;
2826
2827   if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY)
2828     return res;
2829
2830   pos = USED(mp) - 1;
2831   dp = DIGITS(mp);
2832
2833   /* Shift all the significant figures over as needed */
2834   for(ix = pos - p; ix >= 0; ix--)
2835     dp[ix + p] = dp[ix];
2836
2837   /* Fill the bottom digits with zeroes */
2838   for(ix = 0; ix < p; ix++)
2839     dp[ix] = 0;
2840
2841   return MP_OKAY;
2842
2843 } /* end s_mp_lshd() */
2844
2845 /* }}} */
2846
2847 /* {{{ s_mp_rshd(mp, p) */
2848
2849 /*
2850    Shift mp rightward by p digits.  Maintains the invariant that
2851    digits above the precision are all zero.  Digits shifted off the
2852    end are lost.  Cannot fail.
2853  */
2854
2855 void     s_mp_rshd(mp_int *mp, mp_size p)
2856 {
2857   mp_size  ix;
2858   mp_digit *dp;
2859
2860   if(p == 0)
2861     return;
2862
2863   /* Shortcut when all digits are to be shifted off */
2864   if(p >= USED(mp)) {
2865     s_mp_setz(DIGITS(mp), ALLOC(mp));
2866     USED(mp) = 1;
2867     SIGN(mp) = MP_ZPOS;
2868     return;
2869   }
2870
2871   /* Shift all the significant figures over as needed */
2872   dp = DIGITS(mp);
2873   for(ix = p; ix < USED(mp); ix++)
2874     dp[ix - p] = dp[ix];
2875
2876   /* Fill the top digits with zeroes */
2877   ix -= p;
2878   while(ix < USED(mp))
2879     dp[ix++] = 0;
2880
2881   /* Strip off any leading zeroes    */
2882   s_mp_clamp(mp);
2883
2884 } /* end s_mp_rshd() */
2885
2886 /* }}} */
2887
2888 /* {{{ s_mp_div_2(mp) */
2889
2890 /* Divide by two -- take advantage of radix properties to do it fast      */
2891 void     s_mp_div_2(mp_int *mp)
2892 {
2893   s_mp_div_2d(mp, 1);
2894
2895 } /* end s_mp_div_2() */
2896
2897 /* }}} */
2898
2899 /* {{{ s_mp_mul_2(mp) */
2900
2901 mp_err s_mp_mul_2(mp_int *mp)
2902 {
2903   int      ix;
2904   mp_digit kin = 0, kout, *dp = DIGITS(mp);
2905   mp_err   res;
2906
2907   /* Shift digits leftward by 1 bit */
2908   for(ix = 0; ix < USED(mp); ix++) {
2909     kout = (dp[ix] >> (DIGIT_BIT - 1)) & 1;
2910     dp[ix] = (dp[ix] << 1) | kin;
2911
2912     kin = kout;
2913   }
2914
2915   /* Deal with rollover from last digit */
2916   if(kin) {
2917     if(ix >= ALLOC(mp)) {
2918       if((res = s_mp_grow(mp, ALLOC(mp) + 1)) != MP_OKAY)
2919         return res;
2920       dp = DIGITS(mp);
2921     }
2922
2923     dp[ix] = kin;
2924     USED(mp) += 1;
2925   }
2926
2927   return MP_OKAY;
2928
2929 } /* end s_mp_mul_2() */
2930
2931 /* }}} */
2932
2933 /* {{{ s_mp_mod_2d(mp, d) */
2934
2935 /*
2936   Remainder the integer by 2^d, where d is a number of bits.  This
2937   amounts to a bitwise AND of the value, and does not require the full
2938   division code
2939  */
2940 void     s_mp_mod_2d(mp_int *mp, mp_digit d)
2941 {
2942   unsigned int  ndig = (d / DIGIT_BIT), nbit = (d % DIGIT_BIT);
2943   unsigned int  ix;
2944   mp_digit      dmask, *dp = DIGITS(mp);
2945
2946   if(ndig >= USED(mp))
2947     return;
2948
2949   /* Flush all the bits above 2^d in its digit */
2950   dmask = (1 << nbit) - 1;
2951   dp[ndig] &= dmask;
2952
2953   /* Flush all digits above the one with 2^d in it */
2954   for(ix = ndig + 1; ix < USED(mp); ix++)
2955     dp[ix] = 0;
2956
2957   s_mp_clamp(mp);
2958
2959 } /* end s_mp_mod_2d() */
2960
2961 /* }}} */
2962
2963 /* {{{ s_mp_mul_2d(mp, d) */
2964
2965 /*
2966   Multiply by the integer 2^d, where d is a number of bits.  This
2967   amounts to a bitwise shift of the value, and does not require the
2968   full multiplication code.
2969  */
2970 mp_err    s_mp_mul_2d(mp_int *mp, mp_digit d)
2971 {
2972   mp_err   res;
2973   mp_digit save, next, mask, *dp;
2974   mp_size  used;
2975   int      ix;
2976
2977   if((res = s_mp_lshd(mp, d / DIGIT_BIT)) != MP_OKAY)
2978     return res;
2979
2980   dp = DIGITS(mp); used = USED(mp);
2981   d %= DIGIT_BIT;
2982
2983   mask = (1 << d) - 1;
2984
2985   /* If the shift requires another digit, make sure we've got one to
2986      work with */
2987   if((dp[used - 1] >> (DIGIT_BIT - d)) & mask) {
2988     if((res = s_mp_grow(mp, used + 1)) != MP_OKAY)
2989       return res;
2990     dp = DIGITS(mp);
2991   }
2992
2993   /* Do the shifting... */
2994   save = 0;
2995   for(ix = 0; ix < used; ix++) {
2996     next = (dp[ix] >> (DIGIT_BIT - d)) & mask;
2997     dp[ix] = (dp[ix] << d) | save;
2998     save = next;
2999   }
3000
3001   /* If, at this point, we have a nonzero carryout into the next
3002      digit, we'll increase the size by one digit, and store it...
3003    */
3004   if(save) {
3005     dp[used] = save;
3006     USED(mp) += 1;
3007   }
3008
3009   s_mp_clamp(mp);
3010   return MP_OKAY;
3011
3012 } /* end s_mp_mul_2d() */
3013
3014 /* }}} */
3015
3016 /* {{{ s_mp_div_2d(mp, d) */
3017
3018 /*
3019   Divide the integer by 2^d, where d is a number of bits.  This
3020   amounts to a bitwise shift of the value, and does not require the
3021   full division code (used in Barrett reduction, see below)
3022  */
3023 void     s_mp_div_2d(mp_int *mp, mp_digit d)
3024 {
3025   int       ix;
3026   mp_digit  save, next, mask, *dp = DIGITS(mp);
3027
3028   s_mp_rshd(mp, d / DIGIT_BIT);
3029   d %= DIGIT_BIT;
3030
3031   mask = (1 << d) - 1;
3032
3033   save = 0;
3034   for(ix = USED(mp) - 1; ix >= 0; ix--) {
3035     next = dp[ix] & mask;
3036     dp[ix] = (dp[ix] >> d) | (save << (DIGIT_BIT - d));
3037     save = next;
3038   }
3039
3040   s_mp_clamp(mp);
3041
3042 } /* end s_mp_div_2d() */
3043
3044 /* }}} */
3045
3046 /* {{{ s_mp_norm(a, b) */
3047
3048 /*
3049   s_mp_norm(a, b)
3050
3051   Normalize a and b for division, where b is the divisor.  In order
3052   that we might make good guesses for quotient digits, we want the
3053   leading digit of b to be at least half the radix, which we
3054   accomplish by multiplying a and b by a constant.  This constant is
3055   returned (so that it can be divided back out of the remainder at the
3056   end of the division process).
3057
3058   We multiply by the smallest power of 2 that gives us a leading digit
3059   at least half the radix.  By choosing a power of 2, we simplify the
3060   multiplication and division steps to simple shifts.
3061  */
3062 mp_digit s_mp_norm(mp_int *a, mp_int *b)
3063 {
3064   mp_digit  t, d = 0;
3065
3066   t = DIGIT(b, USED(b) - 1);
3067   while(t < (RADIX / 2)) {
3068     t <<= 1;
3069     ++d;
3070   }
3071
3072   if(d != 0) {
3073     s_mp_mul_2d(a, d);
3074     s_mp_mul_2d(b, d);
3075   }
3076
3077   return d;
3078
3079 } /* end s_mp_norm() */
3080
3081 /* }}} */
3082
3083 /* }}} */
3084
3085 /* {{{ Primitive digit arithmetic */
3086
3087 /* {{{ s_mp_add_d(mp, d) */
3088
3089 /* Add d to |mp| in place                                                 */
3090 mp_err   s_mp_add_d(mp_int *mp, mp_digit d)    /* unsigned digit addition */
3091 {
3092   mp_word   w, k = 0;
3093   mp_size   ix = 1, used = USED(mp);
3094   mp_digit *dp = DIGITS(mp);
3095
3096   w = dp[0] + d;
3097   dp[0] = ACCUM(w);
3098   k = CARRYOUT(w);
3099
3100   while(ix < used && k) {
3101     w = dp[ix] + k;
3102     dp[ix] = ACCUM(w);
3103     k = CARRYOUT(w);
3104     ++ix;
3105   }
3106
3107   if(k != 0) {
3108     mp_err  res;
3109
3110     if((res = s_mp_pad(mp, USED(mp) + 1)) != MP_OKAY)
3111       return res;
3112
3113     DIGIT(mp, ix) = k;
3114   }
3115
3116   return MP_OKAY;
3117
3118 } /* end s_mp_add_d() */
3119
3120 /* }}} */
3121
3122 /* {{{ s_mp_sub_d(mp, d) */
3123
3124 /* Subtract d from |mp| in place, assumes |mp| > d                        */
3125 mp_err   s_mp_sub_d(mp_int *mp, mp_digit d)    /* unsigned digit subtract */
3126 {
3127   mp_word   w, b = 0;
3128   mp_size   ix = 1, used = USED(mp);
3129   mp_digit *dp = DIGITS(mp);
3130
3131   /* Compute initial subtraction    */
3132   w = (RADIX + dp[0]) - d;
3133   b = CARRYOUT(w) ? 0 : 1;
3134   dp[0] = ACCUM(w);
3135
3136   /* Propagate borrows leftward     */
3137   while(b && ix < used) {
3138     w = (RADIX + dp[ix]) - b;
3139     b = CARRYOUT(w) ? 0 : 1;
3140     dp[ix] = ACCUM(w);
3141     ++ix;
3142   }
3143
3144   /* Remove leading zeroes          */
3145   s_mp_clamp(mp);
3146
3147   /* If we have a borrow out, it's a violation of the input invariant */
3148   if(b)
3149     return MP_RANGE;
3150   else
3151     return MP_OKAY;
3152
3153 } /* end s_mp_sub_d() */
3154
3155 /* }}} */
3156
3157 /* {{{ s_mp_mul_d(a, d) */
3158
3159 /* Compute a = a * d, single digit multiplication                         */
3160 mp_err   s_mp_mul_d(mp_int *a, mp_digit d)
3161 {
3162   mp_word w, k = 0;
3163   mp_size ix, max;
3164   mp_err  res;
3165   mp_digit *dp = DIGITS(a);
3166
3167   /*
3168     Single-digit multiplication will increase the precision of the
3169     output by at most one digit.  However, we can detect when this
3170     will happen -- if the high-order digit of a, times d, gives a
3171     two-digit result, then the precision of the result will increase;
3172     otherwise it won't.  We use this fact to avoid calling s_mp_pad()
3173     unless absolutely necessary.
3174    */
3175   max = USED(a);
3176   w = dp[max - 1] * d;
3177   if(CARRYOUT(w) != 0) {
3178     if((res = s_mp_pad(a, max + 1)) != MP_OKAY)
3179       return res;
3180     dp = DIGITS(a);
3181   }
3182
3183   for(ix = 0; ix < max; ix++) {
3184     w = (dp[ix] * d) + k;
3185     dp[ix] = ACCUM(w);
3186     k = CARRYOUT(w);
3187   }
3188
3189   /* If there is a precision increase, take care of it here; the above
3190      test guarantees we have enough storage to do this safely.
3191    */
3192   if(k) {
3193     dp[max] = k;
3194     USED(a) = max + 1;
3195   }
3196
3197   s_mp_clamp(a);
3198
3199   return MP_OKAY;
3200
3201 } /* end s_mp_mul_d() */
3202
3203 /* }}} */
3204
3205 /* {{{ s_mp_div_d(mp, d, r) */
3206
3207 /*
3208   s_mp_div_d(mp, d, r)
3209
3210   Compute the quotient mp = mp / d and remainder r = mp mod d, for a
3211   single digit d.  If r is null, the remainder will be discarded.
3212  */
3213
3214 mp_err   s_mp_div_d(mp_int *mp, mp_digit d, mp_digit *r)
3215 {
3216   mp_word   w = 0, t;
3217   mp_int    quot;
3218   mp_err    res;
3219   mp_digit *dp = DIGITS(mp), *qp;
3220   int       ix;
3221
3222   if(d == 0)
3223     return MP_RANGE;
3224
3225   /* Make room for the quotient */
3226   if((res = mp_init_size(&quot, USED(mp))) != MP_OKAY)
3227     return res;
3228
3229   USED(&quot) = USED(mp); /* so clamping will work below */
3230   qp = DIGITS(&quot);
3231
3232   /* Divide without subtraction */
3233   for(ix = USED(mp) - 1; ix >= 0; ix--) {
3234     w = (w << DIGIT_BIT) | dp[ix];
3235
3236     if(w >= d) {
3237       t = w / d;
3238       w = w % d;
3239     } else {
3240       t = 0;
3241     }
3242
3243     qp[ix] = t;
3244   }
3245
3246   /* Deliver the remainder, if desired */
3247   if(r)
3248     *r = w;
3249
3250   s_mp_clamp(&quot);
3251   mp_exch(&quot, mp);
3252   mp_clear(&quot);
3253
3254   return MP_OKAY;
3255
3256 } /* end s_mp_div_d() */
3257
3258 /* }}} */
3259
3260 /* }}} */
3261
3262 /* {{{ Primitive full arithmetic */
3263
3264 /* {{{ s_mp_add(a, b) */
3265
3266 /* Compute a = |a| + |b|                                                  */
3267 mp_err   s_mp_add(mp_int *a, mp_int *b)        /* magnitude addition      */
3268 {
3269   mp_word   w = 0;
3270   mp_digit *pa, *pb;
3271   mp_size   ix, used = USED(b);
3272   mp_err    res;
3273
3274   /* Make sure a has enough precision for the output value */
3275   if((used > USED(a)) && (res = s_mp_pad(a, used)) != MP_OKAY)
3276     return res;
3277
3278   /*
3279     Add up all digits up to the precision of b.  If b had initially
3280     the same precision as a, or greater, we took care of it by the
3281     padding step above, so there is no problem.  If b had initially
3282     less precision, we'll have to make sure the carry out is duly
3283     propagated upward among the higher-order digits of the sum.
3284    */
3285   pa = DIGITS(a);
3286   pb = DIGITS(b);
3287   for(ix = 0; ix < used; ++ix) {
3288     w += *pa + *pb++;
3289     *pa++ = ACCUM(w);
3290     w = CARRYOUT(w);
3291   }
3292
3293   /* If we run out of 'b' digits before we're actually done, make
3294      sure the carries get propagated upward...
3295    */
3296   used = USED(a);
3297   while(w && ix < used) {
3298     w += *pa;
3299     *pa++ = ACCUM(w);
3300     w = CARRYOUT(w);
3301     ++ix;
3302   }
3303
3304   /* If there's an overall carry out, increase precision and include
3305      it.  We could have done this initially, but why touch the memory
3306      allocator unless we're sure we have to?
3307    */
3308   if(w) {
3309     if((res = s_mp_pad(a, used + 1)) != MP_OKAY)
3310       return res;
3311
3312     DIGIT(a, ix) = w;  /* pa may not be valid after s_mp_pad() call */
3313   }
3314
3315   return MP_OKAY;
3316
3317 } /* end s_mp_add() */
3318
3319 /* }}} */
3320
3321 /* {{{ s_mp_sub(a, b) */
3322
3323 /* Compute a = |a| - |b|, assumes |a| >= |b|                              */
3324 mp_err   s_mp_sub(mp_int *a, mp_int *b)        /* magnitude subtract      */
3325 {
3326   mp_word   w = 0;
3327   mp_digit *pa, *pb;
3328   mp_size   ix, used = USED(b);
3329
3330   /*
3331     Subtract and propagate borrow.  Up to the precision of b, this
3332     accounts for the digits of b; after that, we just make sure the
3333     carries get to the right place.  This saves having to pad b out to
3334     the precision of a just to make the loops work right...
3335    */
3336   pa = DIGITS(a);
3337   pb = DIGITS(b);
3338
3339   for(ix = 0; ix < used; ++ix) {
3340     w = (RADIX + *pa) - w - *pb++;
3341     *pa++ = ACCUM(w);
3342     w = CARRYOUT(w) ? 0 : 1;
3343   }
3344
3345   used = USED(a);
3346   while(ix < used) {
3347     w = RADIX + *pa - w;
3348     *pa++ = ACCUM(w);
3349     w = CARRYOUT(w) ? 0 : 1;
3350     ++ix;
3351   }
3352
3353   /* Clobber any leading zeroes we created    */
3354   s_mp_clamp(a);
3355
3356   /*
3357      If there was a borrow out, then |b| > |a| in violation
3358      of our input invariant.  We've already done the work,
3359      but we'll at least complain about it...
3360    */
3361   if(w)
3362     return MP_RANGE;
3363   else
3364     return MP_OKAY;
3365
3366 } /* end s_mp_sub() */
3367
3368 /* }}} */
3369
3370 mp_err   s_mp_reduce(mp_int *x, mp_int *m, mp_int *mu)
3371 {
3372   mp_int   q;
3373   mp_err   res;
3374   mp_size  um = USED(m);
3375
3376   if((res = mp_init_copy(&q, x)) != MP_OKAY)
3377     return res;
3378
3379   s_mp_rshd(&q, um - 1);       /* q1 = x / b^(k-1)  */
3380   s_mp_mul(&q, mu);            /* q2 = q1 * mu      */
3381   s_mp_rshd(&q, um + 1);       /* q3 = q2 / b^(k+1) */
3382
3383   /* x = x mod b^(k+1), quick (no division) */
3384   s_mp_mod_2d(x, (mp_digit)(DIGIT_BIT * (um + 1)));
3385
3386   /* q = q * m mod b^(k+1), quick (no division), uses the short multiplier */
3387 #ifndef SHRT_MUL
3388   s_mp_mul(&q, m);
3389   s_mp_mod_2d(&q, (mp_digit)(DIGIT_BIT * (um + 1)));
3390 #else
3391   s_mp_mul_dig(&q, m, um + 1);
3392 #endif
3393
3394   /* x = x - q */
3395   if((res = mp_sub(x, &q, x)) != MP_OKAY)
3396     goto CLEANUP;
3397
3398   /* If x < 0, add b^(k+1) to it */
3399   if(mp_cmp_z(x) < 0) {
3400     mp_set(&q, 1);
3401     if((res = s_mp_lshd(&q, um + 1)) != MP_OKAY)
3402       goto CLEANUP;
3403     if((res = mp_add(x, &q, x)) != MP_OKAY)
3404       goto CLEANUP;
3405   }
3406
3407   /* Back off if it's too big */
3408   while(mp_cmp(x, m) >= 0) {
3409     if((res = s_mp_sub(x, m)) != MP_OKAY)
3410       break;
3411   }
3412
3413  CLEANUP:
3414   mp_clear(&q);
3415
3416   return res;
3417
3418 } /* end s_mp_reduce() */
3419
3420
3421
3422 /* {{{ s_mp_mul(a, b) */
3423
3424 /* Compute a = |a| * |b|                                                  */
3425 mp_err   s_mp_mul(mp_int *a, mp_int *b)
3426 {
3427   mp_word   w, k = 0;
3428   mp_int    tmp;
3429   mp_err    res;
3430   mp_size   ix, jx, ua = USED(a), ub = USED(b);
3431   mp_digit *pa, *pb, *pt, *pbt;
3432
3433   if((res = mp_init_size(&tmp, ua + ub)) != MP_OKAY)
3434     return res;
3435
3436   /* This has the effect of left-padding with zeroes... */
3437   USED(&tmp) = ua + ub;
3438
3439   /* We're going to need the base value each iteration */
3440   pbt = DIGITS(&tmp);
3441
3442   /* Outer loop:  Digits of b */
3443
3444   pb = DIGITS(b);
3445   for(ix = 0; ix < ub; ++ix, ++pb) {
3446     if(*pb == 0)
3447       continue;
3448
3449     /* Inner product:  Digits of a */
3450     pa = DIGITS(a);
3451     for(jx = 0; jx < ua; ++jx, ++pa) {
3452       pt = pbt + ix + jx;
3453       w = *pb * *pa + k + *pt;
3454       *pt = ACCUM(w);
3455       k = CARRYOUT(w);
3456     }
3457
3458     pbt[ix + jx] = k;
3459     k = 0;
3460   }
3461
3462   s_mp_clamp(&tmp);
3463   s_mp_exch(&tmp, a);
3464
3465   mp_clear(&tmp);
3466
3467   return MP_OKAY;
3468
3469 } /* end s_mp_mul() */
3470
3471 /* }}} */
3472
3473 /* {{{ s_mp_kmul(a, b, out, len) */
3474
3475 #if 0
3476 void   s_mp_kmul(mp_digit *a, mp_digit *b, mp_digit *out, mp_size len)
3477 {
3478   mp_word   w, k = 0;
3479   mp_size   ix, jx;
3480   mp_digit *pa, *pt;
3481
3482   for(ix = 0; ix < len; ++ix, ++b) {
3483     if(*b == 0)
3484       continue;
3485
3486     pa = a;
3487     for(jx = 0; jx < len; ++jx, ++pa) {
3488       pt = out + ix + jx;
3489       w = *b * *pa + k + *pt;
3490       *pt = ACCUM(w);
3491       k = CARRYOUT(w);
3492     }
3493
3494     out[ix + jx] = k;
3495     k = 0;
3496   }
3497
3498 } /* end s_mp_kmul() */
3499 #endif
3500
3501 /* }}} */
3502
3503 /* {{{ s_mp_sqr(a) */
3504
3505 /*
3506   Computes the square of a, in place.  This can be done more
3507   efficiently than a general multiplication, because many of the
3508   computation steps are redundant when squaring.  The inner product
3509   step is a bit more complicated, but we save a fair number of
3510   iterations of the multiplication loop.
3511  */
3512 #if MP_SQUARE
3513 mp_err   s_mp_sqr(mp_int *a)
3514 {
3515   mp_word  w, k = 0;
3516   mp_int   tmp;
3517   mp_err   res;
3518   mp_size  ix, jx, kx, used = USED(a);
3519   mp_digit *pa1, *pa2, *pt, *pbt;
3520
3521   if((res = mp_init_size(&tmp, 2 * used)) != MP_OKAY)
3522     return res;
3523
3524   /* Left-pad with zeroes */
3525   USED(&tmp) = 2 * used;
3526
3527   /* We need the base value each time through the loop */
3528   pbt = DIGITS(&tmp);
3529
3530   pa1 = DIGITS(a);
3531   for(ix = 0; ix < used; ++ix, ++pa1) {
3532     if(*pa1 == 0)
3533       continue;
3534
3535     w = DIGIT(&tmp, ix + ix) + (*pa1 * *pa1);
3536
3537     pbt[ix + ix] = ACCUM(w);
3538     k = CARRYOUT(w);
3539
3540     /*
3541       The inner product is computed as:
3542
3543          (C, S) = t[i,j] + 2 a[i] a[j] + C
3544
3545       This can overflow what can be represented in an mp_word, and
3546       since C arithmetic does not provide any way to check for
3547       overflow, we have to check explicitly for overflow conditions
3548       before they happen.
3549      */
3550     for(jx = ix + 1, pa2 = DIGITS(a) + jx; jx < used; ++jx, ++pa2) {
3551       mp_word  u = 0, v;
3552
3553       /* Store this in a temporary to avoid indirections later */
3554       pt = pbt + ix + jx;
3555
3556       /* Compute the multiplicative step */
3557       w = *pa1 * *pa2;
3558
3559       /* If w is more than half MP_WORD_MAX, the doubling will
3560          overflow, and we need to record a carry out into the next
3561          word */
3562       u = (w >> (MP_WORD_BIT - 1)) & 1;
3563
3564       /* Double what we've got, overflow will be ignored as defined
3565          for C arithmetic (we've already noted if it is to occur)
3566        */
3567       w *= 2;
3568
3569       /* Compute the additive step */
3570       v = *pt + k;
3571
3572       /* If we do not already have an overflow carry, check to see
3573          if the addition will cause one, and set the carry out if so
3574        */
3575       u |= ((MP_WORD_MAX - v) < w);
3576
3577       /* Add in the rest, again ignoring overflow */
3578       w += v;
3579
3580       /* Set the i,j digit of the output */
3581       *pt = ACCUM(w);
3582
3583       /* Save carry information for the next iteration of the loop.
3584          This is why k must be an mp_word, instead of an mp_digit */
3585       k = CARRYOUT(w) | (u << DIGIT_BIT);
3586
3587     } /* for(jx ...) */
3588
3589     /* Set the last digit in the cycle and reset the carry */
3590     k = DIGIT(&tmp, ix + jx) + k;
3591     pbt[ix + jx] = ACCUM(k);
3592     k = CARRYOUT(k);
3593
3594     /* If we are carrying out, propagate the carry to the next digit
3595        in the output.  This may cascade, so we have to be somewhat
3596        circumspect -- but we will have enough precision in the output
3597        that we won't overflow
3598      */
3599     kx = 1;
3600     while(k) {
3601       k = pbt[ix + jx + kx] + 1;
3602       pbt[ix + jx + kx] = ACCUM(k);
3603       k = CARRYOUT(k);
3604       ++kx;
3605     }
3606   } /* for(ix ...) */
3607
3608   s_mp_clamp(&tmp);
3609   s_mp_exch(&tmp, a);
3610
3611   mp_clear(&tmp);
3612
3613   return MP_OKAY;
3614
3615 } /* end s_mp_sqr() */
3616 #endif
3617
3618 /* }}} */
3619
3620 /* {{{ s_mp_div(a, b) */
3621
3622 /*
3623   s_mp_div(a, b)
3624
3625   Compute a = a / b and b = a mod b.  Assumes b > a.
3626  */
3627
3628 mp_err   s_mp_div(mp_int *a, mp_int *b)
3629 {
3630   mp_int   quot, rem, t;
3631   mp_word  q;
3632   mp_err   res;
3633   mp_digit d;
3634   int      ix;
3635
3636   if(mp_cmp_z(b) == 0)
3637     return MP_RANGE;
3638
3639   /* Shortcut if b is power of two */
3640   if((ix = s_mp_ispow2(b)) >= 0) {
3641     mp_copy(a, b);  /* need this for remainder */
3642     s_mp_div_2d(a, (mp_digit)ix);
3643     s_mp_mod_2d(b, (mp_digit)ix);
3644
3645     return MP_OKAY;
3646   }
3647
3648   /* Allocate space to store the quotient */
3649   if((res = mp_init_size(&quot, USED(a))) != MP_OKAY)
3650     return res;
3651
3652   /* A working temporary for division     */
3653   if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
3654     goto T;
3655
3656   /* Allocate space for the remainder     */
3657   if((res = mp_init_size(&rem, USED(a))) != MP_OKAY)
3658     goto REM;
3659
3660   /* Normalize to optimize guessing       */
3661   d = s_mp_norm(a, b);
3662
3663   /* Perform the division itself...woo!   */
3664   ix = USED(a) - 1;
3665
3666   while(ix >= 0) {
3667     /* Find a partial substring of a which is at least b */
3668     while(s_mp_cmp(&rem, b) < 0 && ix >= 0) {
3669       if((res = s_mp_lshd(&rem, 1)) != MP_OKAY)
3670         goto CLEANUP;
3671
3672       if((res = s_mp_lshd(&quot, 1)) != MP_OKAY)
3673         goto CLEANUP;
3674
3675       DIGIT(&rem, 0) = DIGIT(a, ix);
3676       s_mp_clamp(&rem);
3677       --ix;
3678     }
3679
3680     /* If we didn't find one, we're finished dividing    */
3681     if(s_mp_cmp(&rem, b) < 0)
3682       break;
3683
3684     /* Compute a guess for the next quotient digit       */
3685     q = DIGIT(&rem, USED(&rem) - 1);
3686     if(q <= DIGIT(b, USED(b) - 1) && USED(&rem) > 1)
3687       q = (q << DIGIT_BIT) | DIGIT(&rem, USED(&rem) - 2);
3688
3689     q /= DIGIT(b, USED(b) - 1);
3690
3691     /* The guess can be as much as RADIX + 1 */
3692     if(q >= RADIX)
3693       q = RADIX - 1;
3694
3695     /* See what that multiplies out to                   */
3696     mp_copy(b, &t);
3697     if((res = s_mp_mul_d(&t, q)) != MP_OKAY)
3698       goto CLEANUP;
3699
3700     /*
3701        If it's too big, back it off.  We should not have to do this
3702        more than once, or, in rare cases, twice.  Knuth describes a
3703        method by which this could be reduced to a maximum of once, but
3704        I didn't implement that here.
3705      */
3706     while(s_mp_cmp(&t, &rem) > 0) {
3707       --q;
3708       s_mp_sub(&t, b);
3709     }
3710
3711     /* At this point, q should be the right next digit   */
3712     if((res = s_mp_sub(&rem, &t)) != MP_OKAY)
3713       goto CLEANUP;
3714
3715     /*
3716       Include the digit in the quotient.  We allocated enough memory
3717       for any quotient we could ever possibly get, so we should not
3718       have to check for failures here
3719      */
3720     DIGIT(&quot, 0) = q;
3721   }
3722
3723   /* Denormalize remainder                */
3724   if(d != 0)
3725     s_mp_div_2d(&rem, d);
3726
3727   s_mp_clamp(&quot);
3728   s_mp_clamp(&rem);
3729
3730   /* Copy quotient back to output         */
3731   s_mp_exch(&quot, a);
3732
3733   /* Copy remainder back to output        */
3734   s_mp_exch(&rem, b);
3735
3736 CLEANUP:
3737   mp_clear(&rem);
3738 REM:
3739   mp_clear(&t);
3740 T:
3741   mp_clear(&quot);
3742
3743   return res;
3744
3745 } /* end s_mp_div() */
3746
3747 /* }}} */
3748
3749 /* {{{ s_mp_2expt(a, k) */
3750
3751 mp_err   s_mp_2expt(mp_int *a, mp_digit k)
3752 {
3753   mp_err    res;
3754   mp_size   dig, bit;
3755
3756   dig = k / DIGIT_BIT;
3757   bit = k % DIGIT_BIT;
3758
3759   mp_zero(a);
3760   if((res = s_mp_pad(a, dig + 1)) != MP_OKAY)
3761     return res;
3762
3763   DIGIT(a, dig) |= (1 << bit);
3764
3765   return MP_OKAY;
3766
3767 } /* end s_mp_2expt() */
3768
3769 /* }}} */
3770
3771
3772 /* }}} */
3773
3774 /* }}} */
3775
3776 /* {{{ Primitive comparisons */
3777
3778 /* {{{ s_mp_cmp(a, b) */
3779
3780 /* Compare |a| <=> |b|, return 0 if equal, <0 if a<b, >0 if a>b           */
3781 int      s_mp_cmp(mp_int *a, mp_int *b)
3782 {
3783   mp_size   ua = USED(a), ub = USED(b);
3784
3785   if(ua > ub)
3786     return MP_GT;
3787   else if(ua < ub)
3788     return MP_LT;
3789   else {
3790     int      ix = ua - 1;
3791     mp_digit *ap = DIGITS(a) + ix, *bp = DIGITS(b) + ix;
3792
3793     while(ix >= 0) {
3794       if(*ap > *bp)
3795         return MP_GT;
3796       else if(*ap < *bp)
3797         return MP_LT;
3798
3799       --ap; --bp; --ix;
3800     }
3801
3802     return MP_EQ;
3803   }
3804
3805 } /* end s_mp_cmp() */
3806
3807 /* }}} */
3808
3809 /* {{{ s_mp_cmp_d(a, d) */
3810
3811 /* Compare |a| <=> d, return 0 if equal, <0 if a<d, >0 if a>d             */
3812 int      s_mp_cmp_d(mp_int *a, mp_digit d)
3813 {
3814   mp_size  ua = USED(a);
3815   mp_digit *ap = DIGITS(a);
3816
3817   if(ua > 1)
3818     return MP_GT;
3819
3820   if(*ap < d)
3821     return MP_LT;
3822   else if(*ap > d)
3823     return MP_GT;
3824   else
3825     return MP_EQ;
3826
3827 } /* end s_mp_cmp_d() */
3828
3829 /* }}} */
3830
3831 /* {{{ s_mp_ispow2(v) */
3832
3833 /*
3834   Returns -1 if the value is not a power of two; otherwise, it returns
3835   k such that v = 2^k, i.e. lg(v).
3836  */
3837 int      s_mp_ispow2(mp_int *v)
3838 {
3839   mp_digit d, *dp;
3840   mp_size  uv = USED(v);
3841   int      extra = 0, ix;
3842
3843   d = DIGIT(v, uv - 1); /* most significant digit of v */
3844
3845   while(d && ((d & 1) == 0)) {
3846     d >>= 1;
3847     ++extra;
3848   }
3849
3850   if(d == 1) {
3851     ix = uv - 2;
3852     dp = DIGITS(v) + ix;
3853
3854     while(ix >= 0) {
3855       if(*dp)
3856         return -1; /* not a power of two */
3857
3858       --dp; --ix;
3859     }
3860
3861     return ((uv - 1) * DIGIT_BIT) + extra;
3862   }
3863
3864   return -1;
3865
3866 } /* end s_mp_ispow2() */
3867
3868 /* }}} */
3869
3870 /* {{{ s_mp_ispow2d(d) */
3871
3872 int      s_mp_ispow2d(mp_digit d)
3873 {
3874   int   pow = 0;
3875
3876   while((d & 1) == 0) {
3877     ++pow; d >>= 1;
3878   }
3879
3880   if(d == 1)
3881     return pow;
3882
3883   return -1;
3884
3885 } /* end s_mp_ispow2d() */
3886
3887 /* }}} */
3888
3889 /* }}} */
3890
3891 /* {{{ Primitive I/O helpers */
3892
3893 /* {{{ s_mp_tovalue(ch, r) */
3894
3895 /*
3896   Convert the given character to its digit value, in the given radix.
3897   If the given character is not understood in the given radix, -1 is
3898   returned.  Otherwise the digit's numeric value is returned.
3899
3900   The results will be odd if you use a radix < 2 or > 62, you are
3901   expected to know what you're up to.
3902  */
3903 int      s_mp_tovalue(char ch, int r)
3904 {
3905   int    val, xch;
3906
3907   if(r > 36)
3908     xch = ch;
3909   else
3910     xch = toupper(ch);
3911
3912   if(isdigit(xch))
3913     val = xch - '0';
3914   else if(isupper(xch))
3915     val = xch - 'A' + 10;
3916   else if(islower(xch))
3917     val = xch - 'a' + 36;
3918   else if(xch == '+')
3919     val = 62;
3920   else if(xch == '/')
3921     val = 63;
3922   else
3923     return -1;
3924
3925   if(val < 0 || val >= r)
3926     return -1;
3927
3928   return val;
3929
3930 } /* end s_mp_tovalue() */
3931
3932 /* }}} */
3933
3934 /* {{{ s_mp_todigit(val, r, low) */
3935
3936 /*
3937   Convert val to a radix-r digit, if possible.  If val is out of range
3938   for r, returns zero.  Otherwise, returns an ASCII character denoting
3939   the value in the given radix.
3940
3941   The results may be odd if you use a radix < 2 or > 64, you are
3942   expected to know what you're doing.
3943  */
3944
3945 char     s_mp_todigit(int val, int r, int low)
3946 {
3947   char   ch;
3948
3949   if(val < 0 || val >= r)
3950     return 0;
3951
3952   ch = s_dmap_1[val];
3953
3954   if(r <= 36 && low)
3955     ch = tolower(ch);
3956
3957   return ch;
3958
3959 } /* end s_mp_todigit() */
3960
3961 /* }}} */
3962
3963 /* {{{ s_mp_outlen(bits, radix) */
3964
3965 /*
3966    Return an estimate for how long a string is needed to hold a radix
3967    r representation of a number with 'bits' significant bits.
3968
3969    Does not include space for a sign or a NUL terminator.
3970  */
3971 int      s_mp_outlen(int bits, int r)
3972 {
3973   return (int)((double)bits * LOG_V_2(r));
3974
3975 } /* end s_mp_outlen() */
3976
3977 /* }}} */
3978
3979 /* }}} */
3980
3981 /*------------------------------------------------------------------------*/
3982 /* HERE THERE BE DRAGONS                                                  */
3983 /* crc==4242132123, version==2, Sat Feb 02 06:43:52 2002 */
3984
3985 /* Source: /cvs/libtom/libtommath/mtest/mpi.c,v  */
3986 /* Revision: 1.2  */
3987 /* Date: 2005/05/05 14:38:47  */