package/libs/nettle/patches/0001-Updated-mini-gmp.patch

   1 From 57700c26d73cf7fa6f5cfaec1145eccf388acab9 Mon Sep 17 00:00:00 2001
   2 From: Nikos Mavrogiannopoulos <nmav@gnutls.org>
   3 Date: Sun, 9 Mar 2014 11:27:42 +0100
   4 Subject: [PATCH 1/5] Updated mini-gmp
   5
   6 ---
   7  mini-gmp.c | 890 +++++++++++++++++++++++++++++++++++++++----------------------
   8  mini-gmp.h |  55 +++-
   9  2 files changed, 618 insertions(+), 327 deletions(-)
  10
  11 diff --git a/mini-gmp.c b/mini-gmp.c
  12 index 8b6f070..766df30 100644
  13 --- a/mini-gmp.c
  14 +++ b/mini-gmp.c
  15 @@ -2,24 +2,33 @@
  16
  17     Contributed to the GNU project by Niels Möller
  18
  19 -Copyright 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001,
  20 -2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013
  21 -Free Software Foundation, Inc.
  22 +Copyright 1991-1997, 1999-2014 Free Software Foundation, Inc.
  23
  24  This file is part of the GNU MP Library.
  25
  26  The GNU MP Library is free software; you can redistribute it and/or modify
  27 -it under the terms of the GNU Lesser General Public License as published by
  28 -the Free Software Foundation; either version 3 of the License, or (at your
  29 -option) any later version.
  30 +it under the terms of either:
  31 +
  32 +  * the GNU Lesser General Public License as published by the Free
  33 +    Software Foundation; either version 3 of the License, or (at your
  34 +    option) any later version.
  35 +
  36 +or
  37 +
  38 +  * the GNU General Public License as published by the Free Software
  39 +    Foundation; either version 2 of the License, or (at your option) any
  40 +    later version.
  41 +
  42 +or both in parallel, as here.
  43
  44  The GNU MP Library is distributed in the hope that it will be useful, but
  45  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  46 -or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  47 -License for more details.
  48 +or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  49 +for more details.
  50
  51 -You should have received a copy of the GNU Lesser General Public License
  52 -along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  53 +You should have received copies of the GNU General Public License and the
  54 +GNU Lesser General Public License along with the GNU MP Library.  If not,
  55 +see https://www.gnu.org/licenses/.  */
  56
  57  /* NOTE: All functions in this file which are not declared in
  58     mini-gmp.h are internal, and are not intended to be compatible
  59 @@ -222,11 +231,13 @@ along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  60    } while (0)
  61  #define MPZ_SRCPTR_SWAP(x, y)                                          \
  62    do {                                                                 \
  63 -    mpz_srcptr __mpz_srcptr_swap__tmp = (x);                   \
  64 +    mpz_srcptr __mpz_srcptr_swap__tmp = (x);                           \
  65      (x) = (y);                                                         \
  66      (y) = __mpz_srcptr_swap__tmp;                                      \
  67    } while (0)
  68
  69 +const int mp_bits_per_limb = GMP_LIMB_BITS;
  70 +
  71  \f
  72  /* Memory allocation and other helper functions. */
  73  static void
  74 @@ -342,12 +353,10 @@ mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
  75  int
  76  mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
  77  {
  78 -  for (; n > 0; n--)
  79 +  while (--n >= 0)
  80      {
  81 -      if (ap[n-1] < bp[n-1])
  82 -       return -1;
  83 -      else if (ap[n-1] > bp[n-1])
  84 -       return 1;
  85 +      if (ap[n] != bp[n])
  86 +       return ap[n] > bp[n] ? 1 : -1;
  87      }
  88    return 0;
  89  }
  90 @@ -355,10 +364,8 @@ mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
  91  static int
  92  mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
  93  {
  94 -  if (an > bn)
  95 -    return 1;
  96 -  else if (an < bn)
  97 -    return -1;
  98 +  if (an != bn)
  99 +    return an < bn ? -1 : 1;
 100    else
 101      return mpn_cmp (ap, bp, an);
 102  }
 103 @@ -373,20 +380,31 @@ mpn_normalized_size (mp_srcptr xp, mp_size_t n)
 104
 105  #define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0)
 106
 107 +void
 108 +mpn_zero (mp_ptr rp, mp_size_t n)
 109 +{
 110 +  mp_size_t i;
 111 +
 112 +  for (i = 0; i < n; i++)
 113 +    rp[i] = 0;
 114 +}
 115 +
 116  mp_limb_t
 117  mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 118  {
 119    mp_size_t i;
 120
 121    assert (n > 0);
 122 -
 123 -  for (i = 0; i < n; i++)
 124 +  i = 0;
 125 +  do
 126      {
 127        mp_limb_t r = ap[i] + b;
 128        /* Carry out */
 129        b = (r < b);
 130        rp[i] = r;
 131      }
 132 +  while (++i < n);
 133 +
 134    return b;
 135  }
 136
 137 @@ -429,7 +447,8 @@ mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 138
 139    assert (n > 0);
 140
 141 -  for (i = 0; i < n; i++)
 142 +  i = 0;
 143 +  do
 144      {
 145        mp_limb_t a = ap[i];
 146        /* Carry out */
 147 @@ -437,6 +456,8 @@ mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 148        rp[i] = a - b;
 149        b = cy;
 150      }
 151 +  while (++i < n);
 152 +
 153    return b;
 154  }
 155
 156 @@ -602,7 +623,7 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
 157    retval = low_limb >> tnc;
 158    high_limb = (low_limb << cnt);
 159
 160 -  for (i = n - 1; i != 0; i--)
 161 +  for (i = n; --i != 0;)
 162      {
 163        low_limb = *--up;
 164        *--rp = high_limb | (low_limb >> tnc);
 165 @@ -630,7 +651,7 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
 166    retval = (high_limb << tnc);
 167    low_limb = high_limb >> cnt;
 168
 169 -  for (i = n - 1; i != 0; i--)
 170 +  for (i = n; --i != 0;)
 171      {
 172        high_limb = *up++;
 173        *rp++ = low_limb | (high_limb << tnc);
 174 @@ -641,6 +662,46 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
 175    return retval;
 176  }
 177
 178 +static mp_bitcnt_t
 179 +mpn_common_scan (mp_limb_t limb, mp_size_t i, mp_srcptr up, mp_size_t un,
 180 +                mp_limb_t ux)
 181 +{
 182 +  unsigned cnt;
 183 +
 184 +  assert (ux == 0 || ux == GMP_LIMB_MAX);
 185 +  assert (0 <= i && i <= un );
 186 +
 187 +  while (limb == 0)
 188 +    {
 189 +      i++;
 190 +      if (i == un)
 191 +       return (ux == 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
 192 +      limb = ux ^ up[i];
 193 +    }
 194 +  gmp_ctz (cnt, limb);
 195 +  return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
 196 +}
 197 +
 198 +mp_bitcnt_t
 199 +mpn_scan1 (mp_srcptr ptr, mp_bitcnt_t bit)
 200 +{
 201 +  mp_size_t i;
 202 +  i = bit / GMP_LIMB_BITS;
 203 +
 204 +  return mpn_common_scan ( ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
 205 +                         i, ptr, i, 0);
 206 +}
 207 +
 208 +mp_bitcnt_t
 209 +mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
 210 +{
 211 +  mp_size_t i;
 212 +  i = bit / GMP_LIMB_BITS;
 213 +
 214 +  return mpn_common_scan (~ptr[i] & (GMP_LIMB_MAX << (bit % GMP_LIMB_BITS)),
 215 +                         i, ptr, i, GMP_LIMB_MAX);
 216 +}
 217 +
 218  \f
 219  /* MPN division interface. */
 220  mp_limb_t
 221 @@ -715,8 +776,7 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
 222        if (r < th)
 223         {
 224           m--;
 225 -         if (r > u1 || (r == u1 && tl > u0))
 226 -           m--;
 227 +         m -= ((r > u1) | ((r == u1) & (tl > u0)));
 228         }
 229      }
 230
 231 @@ -836,14 +896,20 @@ mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
 232    assert (d > 0);
 233
 234    /* Special case for powers of two. */
 235 -  if (d > 1 && (d & (d-1)) == 0)
 236 +  if ((d & (d-1)) == 0)
 237      {
 238 -      unsigned shift;
 239        mp_limb_t r = np[0] & (d-1);
 240 -      gmp_ctz (shift, d);
 241        if (qp)
 242 -       mpn_rshift (qp, np, nn, shift);
 243 -
 244 +       {
 245 +         if (d <= 1)
 246 +           mpn_copyi (qp, np, nn);
 247 +         else
 248 +           {
 249 +             unsigned shift;
 250 +             gmp_ctz (shift, d);
 251 +             mpn_rshift (qp, np, nn, shift);
 252 +           }
 253 +       }
 254        return r;
 255      }
 256    else
 257 @@ -880,7 +946,8 @@ mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
 258
 259    r0 = np[nn - 1];
 260
 261 -  for (i = nn - 2; i >= 0; i--)
 262 +  i = nn - 2;
 263 +  do
 264      {
 265        mp_limb_t n0, q;
 266        n0 = np[i];
 267 @@ -889,6 +956,7 @@ mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
 268        if (qp)
 269         qp[i] = q;
 270      }
 271 +  while (--i >= 0);
 272
 273    if (shift > 0)
 274      {
 275 @@ -930,18 +998,19 @@ mpn_div_qr_pi1 (mp_ptr qp,
 276
 277    assert (dn > 2);
 278    assert (nn >= dn);
 279 -  assert ((dp[dn-1] & GMP_LIMB_HIGHBIT) != 0);
 280
 281    d1 = dp[dn - 1];
 282    d0 = dp[dn - 2];
 283
 284 +  assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
 285    /* Iteration variable is the index of the q limb.
 286     *
 287     * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
 288     * by            <d1,          d0,        dp[dn-3],  ..., dp[0] >
 289     */
 290
 291 -  for (i = nn - dn; i >= 0; i--)
 292 +  i = nn - dn;
 293 +  do
 294      {
 295        mp_limb_t n0 = np[dn-1+i];
 296
 297 @@ -973,6 +1042,7 @@ mpn_div_qr_pi1 (mp_ptr qp,
 298        if (qp)
 299         qp[i] = q;
 300      }
 301 +  while (--i >= 0);
 302
 303    np[dn - 1] = n1;
 304  }
 305 @@ -994,7 +1064,9 @@ mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
 306        mp_limb_t nh;
 307        unsigned shift;
 308
 309 -      assert (dp[dn-1] & GMP_LIMB_HIGHBIT);
 310 +      assert (inv->d1 == dp[dn-1]);
 311 +      assert (inv->d0 == dp[dn-2]);
 312 +      assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
 313
 314        shift = inv->shift;
 315        if (shift > 0)
 316 @@ -1002,9 +1074,6 @@ mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
 317        else
 318         nh = 0;
 319
 320 -      assert (inv->d1 == dp[dn-1]);
 321 -      assert (inv->d0 == dp[dn-2]);
 322 -
 323        mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
 324
 325        if (shift > 0)
 326 @@ -1238,15 +1307,14 @@ mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
 327  {
 328    mp_size_t rn;
 329    mp_limb_t w;
 330 -  unsigned first;
 331    unsigned k;
 332    size_t j;
 333
 334 -  first = 1 + (sn - 1) % info->exp;
 335 +  k = 1 + (sn - 1) % info->exp;
 336
 337    j = 0;
 338    w = sp[j++];
 339 -  for (k = 1; k < first; k++)
 340 +  for (; --k > 0; )
 341      w = w * b + sp[j++];
 342
 343    rp[0] = w;
 344 @@ -1300,7 +1368,7 @@ mpz_init (mpz_t r)
 345  }
 346
 347  /* The utility of this function is a bit limited, since many functions
 348 -   assings the result variable using mpz_swap. */
 349 +   assigns the result variable using mpz_swap. */
 350  void
 351  mpz_init2 (mpz_t r, mp_bitcnt_t bits)
 352  {
 353 @@ -1422,7 +1490,7 @@ mpz_fits_ulong_p (const mpz_t u)
 354  {
 355    mp_size_t us = u->_mp_size;
 356
 357 -  return us == 0 || us == 1;
 358 +  return (us == (us > 0));
 359  }
 360
 361  long int
 362 @@ -1459,6 +1527,48 @@ mpz_getlimbn (const mpz_t u, mp_size_t n)
 363      return 0;
 364  }
 365
 366 +void
 367 +mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
 368 +{
 369 +  mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
 370 +}
 371 +
 372 +mp_srcptr
 373 +mpz_limbs_read (mpz_srcptr x)
 374 +{
 375 +  return x->_mp_d;;
 376 +}
 377 +
 378 +mp_ptr
 379 +mpz_limbs_modify (mpz_t x, mp_size_t n)
 380 +{
 381 +  assert (n > 0);
 382 +  return MPZ_REALLOC (x, n);
 383 +}
 384 +
 385 +mp_ptr
 386 +mpz_limbs_write (mpz_t x, mp_size_t n)
 387 +{
 388 +  return mpz_limbs_modify (x, n);
 389 +}
 390 +
 391 +void
 392 +mpz_limbs_finish (mpz_t x, mp_size_t xs)
 393 +{
 394 +  mp_size_t xn;
 395 +  xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
 396 +  x->_mp_size = xs < 0 ? -xn : xn;
 397 +}
 398 +
 399 +mpz_srcptr
 400 +mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
 401 +{
 402 +  x->_mp_alloc = 0;
 403 +  x->_mp_d = (mp_ptr) xp;
 404 +  mpz_limbs_finish (x, xs);
 405 +  return x;
 406 +}
 407 +
 408  \f
 409  /* Conversions and comparison to double. */
 410  void
 411 @@ -1473,19 +1583,15 @@ mpz_set_d (mpz_t r, double x)
 412
 413    /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
 414       zero or infinity. */
 415 -  if (x == 0.0 || x != x || x == x * 0.5)
 416 +  if (x != x || x == x * 0.5)
 417      {
 418        r->_mp_size = 0;
 419        return;
 420      }
 421
 422 -  if (x < 0.0)
 423 -    {
 424 -      x = - x;
 425 -      sign = 1;
 426 -    }
 427 -  else
 428 -    sign = 0;
 429 +  sign = x < 0.0 ;
 430 +  if (sign)
 431 +    x = - x;
 432
 433    if (x < 1.0)
 434      {
 435 @@ -1502,8 +1608,9 @@ mpz_set_d (mpz_t r, double x)
 436    f = (mp_limb_t) x;
 437    x -= f;
 438    assert (x < 1.0);
 439 -  rp[rn-1] = f;
 440 -  for (i = rn-1; i-- > 0; )
 441 +  i = rn-1;
 442 +  rp[i] = f;
 443 +  while (--i >= 0)
 444      {
 445        x = B * x;
 446        f = (mp_limb_t) x;
 447 @@ -1611,12 +1718,7 @@ mpz_sgn (const mpz_t u)
 448  {
 449    mp_size_t usize = u->_mp_size;
 450
 451 -  if (usize > 0)
 452 -    return 1;
 453 -  else if (usize < 0)
 454 -    return -1;
 455 -  else
 456 -    return 0;
 457 +  return (usize > 0) - (usize < 0);
 458  }
 459
 460  int
 461 @@ -1635,10 +1737,9 @@ mpz_cmp_si (const mpz_t u, long v)
 462        mp_limb_t ul = u->_mp_d[0];
 463        if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
 464         return -1;
 465 -      else if ( (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul)
 466 -       return 1;
 467 +      else
 468 +       return (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul;
 469      }
 470 -  return 0;
 471  }
 472
 473  int
 474 @@ -1653,12 +1754,8 @@ mpz_cmp_ui (const mpz_t u, unsigned long v)
 475    else
 476      {
 477        mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
 478 -      if (ul > v)
 479 -       return 1;
 480 -      else if (ul < v)
 481 -       return -1;
 482 +      return (ul > v) - (ul < v);
 483      }
 484 -  return 0;
 485  }
 486
 487  int
 488 @@ -1667,16 +1764,12 @@ mpz_cmp (const mpz_t a, const mpz_t b)
 489    mp_size_t asize = a->_mp_size;
 490    mp_size_t bsize = b->_mp_size;
 491
 492 -  if (asize > bsize)
 493 -    return 1;
 494 -  else if (asize < bsize)
 495 -    return -1;
 496 -  else if (asize > 0)
 497 +  if (asize != bsize)
 498 +    return (asize < bsize) ? -1 : 1;
 499 +  else if (asize >= 0)
 500      return mpn_cmp (a->_mp_d, b->_mp_d, asize);
 501 -  else if (asize < 0)
 502 -    return -mpn_cmp (a->_mp_d, b->_mp_d, -asize);
 503    else
 504 -    return 0;
 505 +    return mpn_cmp (b->_mp_d, a->_mp_d, -asize);
 506  }
 507
 508  int
 509 @@ -1690,12 +1783,7 @@ mpz_cmpabs_ui (const mpz_t u, unsigned long v)
 510
 511    ul = (un == 1) ? u->_mp_d[0] : 0;
 512
 513 -  if (ul > v)
 514 -    return 1;
 515 -  else if (ul < v)
 516 -    return -1;
 517 -  else
 518 -    return 0;
 519 +  return (ul > v) - (ul < v);
 520  }
 521
 522  int
 523 @@ -1753,7 +1841,7 @@ mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
 524
 525    cy = mpn_add_1 (rp, a->_mp_d, an, b);
 526    rp[an] = cy;
 527 -  an += (cy > 0);
 528 +  an += cy;
 529
 530    return an;
 531  }
 532 @@ -1815,20 +1903,21 @@ mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
 533  {
 534    mp_size_t an = GMP_ABS (a->_mp_size);
 535    mp_size_t bn = GMP_ABS (b->_mp_size);
 536 -  mp_size_t rn;
 537    mp_ptr rp;
 538    mp_limb_t cy;
 539
 540 -  rn = GMP_MAX (an, bn);
 541 -  rp = MPZ_REALLOC (r, rn + 1);
 542 -  if (an >= bn)
 543 -    cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
 544 -  else
 545 -    cy = mpn_add (rp, b->_mp_d, bn, a->_mp_d, an);
 546 +  if (an < bn)
 547 +    {
 548 +      MPZ_SRCPTR_SWAP (a, b);
 549 +      MP_SIZE_T_SWAP (an, bn);
 550 +    }
 551
 552 -  rp[rn] = cy;
 553 +  rp = MPZ_REALLOC (r, an + 1);
 554 +  cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
 555 +
 556 +  rp[an] = cy;
 557
 558 -  return rn + (cy > 0);
 559 +  return an + cy;
 560  }
 561
 562  static mp_size_t
 563 @@ -1899,31 +1988,26 @@ mpz_mul_si (mpz_t r, const mpz_t u, long int v)
 564  void
 565  mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
 566  {
 567 -  mp_size_t un;
 568 -  mpz_t t;
 569 +  mp_size_t un, us;
 570    mp_ptr tp;
 571    mp_limb_t cy;
 572
 573 -  un = GMP_ABS (u->_mp_size);
 574 +  us = u->_mp_size;
 575
 576 -  if (un == 0 || v == 0)
 577 +  if (us == 0 || v == 0)
 578      {
 579        r->_mp_size = 0;
 580        return;
 581      }
 582
 583 -  mpz_init2 (t, (un + 1) * GMP_LIMB_BITS);
 584 +  un = GMP_ABS (us);
 585
 586 -  tp = t->_mp_d;
 587 +  tp = MPZ_REALLOC (r, un + 1);
 588    cy = mpn_mul_1 (tp, u->_mp_d, un, v);
 589    tp[un] = cy;
 590
 591 -  t->_mp_size = un + (cy > 0);
 592 -  if (u->_mp_size < 0)
 593 -    t->_mp_size = - t->_mp_size;
 594 -
 595 -  mpz_swap (r, t);
 596 -  mpz_clear (t);
 597 +  un += (cy > 0);
 598 +  r->_mp_size = (us < 0) ? - un : un;
 599  }
 600
 601  void
 602 @@ -1934,8 +2018,8 @@ mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
 603    mpz_t t;
 604    mp_ptr tp;
 605
 606 -  un = GMP_ABS (u->_mp_size);
 607 -  vn = GMP_ABS (v->_mp_size);
 608 +  un = u->_mp_size;
 609 +  vn = v->_mp_size;
 610
 611    if (un == 0 || vn == 0)
 612      {
 613 @@ -1943,7 +2027,10 @@ mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
 614        return;
 615      }
 616
 617 -  sign = (u->_mp_size ^ v->_mp_size) < 0;
 618 +  sign = (un ^ vn) < 0;
 619 +
 620 +  un = GMP_ABS (un);
 621 +  vn = GMP_ABS (vn);
 622
 623    mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
 624
 625 @@ -1996,6 +2083,46 @@ mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
 626    r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
 627  }
 628
 629 +void
 630 +mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
 631 +{
 632 +  mpz_t t;
 633 +  mpz_init (t);
 634 +  mpz_mul_ui (t, u, v);
 635 +  mpz_add (r, r, t);
 636 +  mpz_clear (t);
 637 +}
 638 +
 639 +void
 640 +mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
 641 +{
 642 +  mpz_t t;
 643 +  mpz_init (t);
 644 +  mpz_mul_ui (t, u, v);
 645 +  mpz_sub (r, r, t);
 646 +  mpz_clear (t);
 647 +}
 648 +
 649 +void
 650 +mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
 651 +{
 652 +  mpz_t t;
 653 +  mpz_init (t);
 654 +  mpz_mul (t, u, v);
 655 +  mpz_add (r, r, t);
 656 +  mpz_clear (t);
 657 +}
 658 +
 659 +void
 660 +mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
 661 +{
 662 +  mpz_t t;
 663 +  mpz_init (t);
 664 +  mpz_mul (t, u, v);
 665 +  mpz_sub (r, r, t);
 666 +  mpz_clear (t);
 667 +}
 668 +
 669  \f
 670  /* MPZ division */
 671  enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
 672 @@ -2060,8 +2187,7 @@ mpz_div_qr (mpz_t q, mpz_t r,
 673        mp_size_t qn, rn;
 674        mpz_t tq, tr;
 675
 676 -      mpz_init (tr);
 677 -      mpz_set (tr, n);
 678 +      mpz_init_set (tr, n);
 679        np = tr->_mp_d;
 680
 681        qn = nn - dn + 1;
 682 @@ -2171,10 +2297,7 @@ mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
 683  void
 684  mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
 685  {
 686 -  if (d->_mp_size >= 0)
 687 -    mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
 688 -  else
 689 -    mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
 690 +  mpz_div_qr (NULL, r, n, d, d->_mp_size >= 0 ? GMP_DIV_FLOOR : GMP_DIV_CEIL);
 691  }
 692
 693  static void
 694 @@ -2184,7 +2307,7 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
 695    mp_size_t un, qn;
 696    mp_size_t limb_cnt;
 697    mp_ptr qp;
 698 -  mp_limb_t adjust;
 699 +  int adjust;
 700
 701    un = u->_mp_size;
 702    if (un == 0)
 703 @@ -2226,7 +2349,8 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
 704
 705    q->_mp_size = qn;
 706
 707 -  mpz_add_ui (q, q, adjust);
 708 +  if (adjust)
 709 +    mpz_add_ui (q, q, 1);
 710    if (un < 0)
 711      mpz_neg (q, q);
 712  }
 713 @@ -2303,7 +2427,7 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
 714             {
 715               /* r > 0, need to flip sign. */
 716               rp[i] = ~rp[i] + 1;
 717 -             for (i++; i < rn; i++)
 718 +             while (++i < rn)
 719                 rp[i] = ~rp[i];
 720
 721               rp[rn-1] &= mask;
 722 @@ -2366,6 +2490,24 @@ mpz_divisible_p (const mpz_t n, const mpz_t d)
 723    return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
 724  }
 725
 726 +int
 727 +mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
 728 +{
 729 +  mpz_t t;
 730 +  int res;
 731 +
 732 +  /* a == b (mod 0) iff a == b */
 733 +  if (mpz_sgn (m) == 0)
 734 +    return (mpz_cmp (a, b) == 0);
 735 +
 736 +  mpz_init (t);
 737 +  mpz_sub (t, a, b);
 738 +  res = mpz_divisible_p (t, m);
 739 +  mpz_clear (t);
 740 +
 741 +  return res;
 742 +}
 743 +
 744  static unsigned long
 745  mpz_div_qr_ui (mpz_t q, mpz_t r,
 746                const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
 747 @@ -2579,32 +2721,16 @@ mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
 748  }
 749
 750  static mp_bitcnt_t
 751 -mpz_make_odd (mpz_t r, const mpz_t u)
 752 +mpz_make_odd (mpz_t r)
 753  {
 754 -  mp_size_t un, rn, i;
 755 -  mp_ptr rp;
 756 -  unsigned shift;
 757 -
 758 -  un = GMP_ABS (u->_mp_size);
 759 -  assert (un > 0);
 760 +  mp_bitcnt_t shift;
 761
 762 -  for (i = 0; u->_mp_d[i] == 0; i++)
 763 -    ;
 764 -
 765 -  gmp_ctz (shift, u->_mp_d[i]);
 766 -
 767 -  rn = un - i;
 768 -  rp = MPZ_REALLOC (r, rn);
 769 -  if (shift > 0)
 770 -    {
 771 -      mpn_rshift (rp, u->_mp_d + i, rn, shift);
 772 -      rn -= (rp[rn-1] == 0);
 773 -    }
 774 -  else
 775 -    mpn_copyi (rp, u->_mp_d + i, rn);
 776 +  assert (r->_mp_size > 0);
 777 +  /* Count trailing zeros, equivalent to mpn_scan1, because we know that there is a 1 */
 778 +  shift = mpn_common_scan (r->_mp_d[0], 0, r->_mp_d, 0, 0);
 779 +  mpz_tdiv_q_2exp (r, r, shift);
 780
 781 -  r->_mp_size = rn;
 782 -  return i * GMP_LIMB_BITS + shift;
 783 +  return shift;
 784  }
 785
 786  void
 787 @@ -2627,8 +2753,10 @@ mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
 788    mpz_init (tu);
 789    mpz_init (tv);
 790
 791 -  uz = mpz_make_odd (tu, u);
 792 -  vz = mpz_make_odd (tv, v);
 793 +  mpz_abs (tu, u);
 794 +  uz = mpz_make_odd (tu);
 795 +  mpz_abs (tv, v);
 796 +  vz = mpz_make_odd (tv);
 797    gz = GMP_MIN (uz, vz);
 798
 799    if (tu->_mp_size < tv->_mp_size)
 800 @@ -2644,7 +2772,7 @@ mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
 801        {
 802         int c;
 803
 804 -       mpz_make_odd (tu, tu);
 805 +       mpz_make_odd (tu);
 806         c = mpz_cmp (tu, tv);
 807         if (c == 0)
 808           {
 809 @@ -2706,8 +2834,10 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
 810    mpz_init (t0);
 811    mpz_init (t1);
 812
 813 -  uz = mpz_make_odd (tu, u);
 814 -  vz = mpz_make_odd (tv, v);
 815 +  mpz_abs (tu, u);
 816 +  uz = mpz_make_odd (tu);
 817 +  mpz_abs (tv, v);
 818 +  vz = mpz_make_odd (tv);
 819    gz = GMP_MIN (uz, vz);
 820
 821    uz -= gz;
 822 @@ -2755,7 +2885,7 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
 823    if (tu->_mp_size > 0)
 824      {
 825        mp_bitcnt_t shift;
 826 -      shift = mpz_make_odd (tu, tu);
 827 +      shift = mpz_make_odd (tu);
 828        mpz_mul_2exp (t0, t0, shift);
 829        mpz_mul_2exp (s0, s0, shift);
 830        power += shift;
 831 @@ -2778,7 +2908,7 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
 832               mpz_add (t0, t0, t1);
 833               mpz_add (s0, s0, s1);
 834
 835 -             shift = mpz_make_odd (tv, tv);
 836 +             shift = mpz_make_odd (tv);
 837               mpz_mul_2exp (t1, t1, shift);
 838               mpz_mul_2exp (s1, s1, shift);
 839             }
 840 @@ -2788,7 +2918,7 @@ mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
 841               mpz_add (t1, t0, t1);
 842               mpz_add (s1, s0, s1);
 843
 844 -             shift = mpz_make_odd (tu, tu);
 845 +             shift = mpz_make_odd (tu);
 846               mpz_mul_2exp (t0, t0, shift);
 847               mpz_mul_2exp (s0, s0, shift);
 848             }
 849 @@ -2926,12 +3056,16 @@ mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
 850    mpz_t tr;
 851    mpz_init_set_ui (tr, 1);
 852
 853 -  for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1)
 854 +  bit = GMP_ULONG_HIGHBIT;
 855 +  do
 856      {
 857        mpz_mul (tr, tr, tr);
 858        if (e & bit)
 859         mpz_mul (tr, tr, b);
 860 +      bit >>= 1;
 861      }
 862 +  while (bit > 0);
 863 +
 864    mpz_swap (r, tr);
 865    mpz_clear (tr);
 866  }
 867 @@ -2987,7 +3121,7 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
 868    if (e->_mp_size < 0)
 869      {
 870        if (!mpz_invert (base, b, m))
 871 -       gmp_die ("mpz_powm: Negative exponent and non-invertibe base.");
 872 +       gmp_die ("mpz_powm: Negative exponent and non-invertible base.");
 873      }
 874    else
 875      {
 876 @@ -3019,7 +3153,8 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
 877        mp_limb_t w = e->_mp_d[en];
 878        mp_limb_t bit;
 879
 880 -      for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1)
 881 +      bit = GMP_LIMB_HIGHBIT;
 882 +      do
 883         {
 884           mpz_mul (tr, tr, tr);
 885           if (w & bit)
 886 @@ -3029,7 +3164,9 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
 887               mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
 888               tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
 889             }
 890 +         bit >>= 1;
 891         }
 892 +      while (bit > 0);
 893      }
 894
 895    /* Final reduction */
 896 @@ -3064,21 +3201,26 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
 897    mpz_t t, u;
 898
 899    sgn = y->_mp_size < 0;
 900 -  if (sgn && (z & 1) == 0)
 901 +  if ((~z & sgn) != 0)
 902      gmp_die ("mpz_rootrem: Negative argument, with even root.");
 903    if (z == 0)
 904      gmp_die ("mpz_rootrem: Zeroth root.");
 905
 906    if (mpz_cmpabs_ui (y, 1) <= 0) {
 907 -    mpz_set (x, y);
 908 +    if (x)
 909 +      mpz_set (x, y);
 910      if (r)
 911        r->_mp_size = 0;
 912      return;
 913    }
 914
 915 -  mpz_init (t);
 916    mpz_init (u);
 917 -  mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
 918 +  {
 919 +    mp_bitcnt_t tb;
 920 +    tb = mpz_sizeinbase (y, 2) / z + 1;
 921 +    mpz_init2 (t, tb);
 922 +    mpz_setbit (t, tb);
 923 +  }
 924
 925    if (z == 2) /* simplify sqrt loop: z-1 == 1 */
 926      do {
 927 @@ -3110,7 +3252,8 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
 928      mpz_pow_ui (t, u, z);
 929      mpz_sub (r, y, t);
 930    }
 931 -  mpz_swap (x, u);
 932 +  if (x)
 933 +    mpz_swap (x, u);
 934    mpz_clear (u);
 935    mpz_clear (t);
 936  }
 937 @@ -3142,19 +3285,56 @@ mpz_sqrt (mpz_t s, const mpz_t u)
 938    mpz_rootrem (s, NULL, u, 2);
 939  }
 940
 941 +int
 942 +mpz_perfect_square_p (const mpz_t u)
 943 +{
 944 +  if (u->_mp_size <= 0)
 945 +    return (u->_mp_size == 0);
 946 +  else
 947 +    return mpz_root (NULL, u, 2);
 948 +}
 949 +
 950 +int
 951 +mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
 952 +{
 953 +  mpz_t t;
 954 +
 955 +  assert (n > 0);
 956 +  assert (p [n-1] != 0);
 957 +  return mpz_root (NULL, mpz_roinit_n (t, p, n), 2);
 958 +}
 959 +
 960 +mp_size_t
 961 +mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
 962 +{
 963 +  mpz_t s, r, u;
 964 +  mp_size_t res;
 965 +
 966 +  assert (n > 0);
 967 +  assert (p [n-1] != 0);
 968 +
 969 +  mpz_init (r);
 970 +  mpz_init (s);
 971 +  mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2);
 972 +
 973 +  assert (s->_mp_size == (n+1)/2);
 974 +  mpn_copyd (sp, s->_mp_d, s->_mp_size);
 975 +  mpz_clear (s);
 976 +  res = r->_mp_size;
 977 +  if (rp)
 978 +    mpn_copyd (rp, r->_mp_d, res);
 979 +  mpz_clear (r);
 980 +  return res;
 981 +}
 982  \f
 983  /* Combinatorics */
 984
 985  void
 986  mpz_fac_ui (mpz_t x, unsigned long n)
 987  {
 988 -  if (n < 2) {
 989 -    mpz_set_ui (x, 1);
 990 -    return;
 991 -  }
 992 -  mpz_set_ui (x, n);
 993 -  for (;--n > 1;)
 994 -    mpz_mul_ui (x, x, n);
 995 +  mpz_set_ui (x, n + (n == 0));
 996 +  for (;n > 2;)
 997 +    mpz_mul_ui (x, x, --n);
 998  }
 999
1000  void
1001 @@ -3162,25 +3342,120 @@ mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
1002  {
1003    mpz_t t;
1004
1005 -  if (k > n) {
1006 -    r->_mp_size = 0;
1007 -    return;
1008 -  }
1009 -  mpz_fac_ui (r, n);
1010 +  mpz_set_ui (r, k <= n);
1011 +
1012 +  if (k > (n >> 1))
1013 +    k = (k <= n) ? n - k : 0;
1014 +
1015    mpz_init (t);
1016    mpz_fac_ui (t, k);
1017 -  mpz_divexact (r, r, t);
1018 -  mpz_fac_ui (t, n - k);
1019 +
1020 +  for (; k > 0; k--)
1021 +      mpz_mul_ui (r, r, n--);
1022 +
1023    mpz_divexact (r, r, t);
1024    mpz_clear (t);
1025  }
1026
1027  \f
1028 +/* Primality testing */
1029 +static int
1030 +gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
1031 +                const mpz_t q, mp_bitcnt_t k)
1032 +{
1033 +  mp_bitcnt_t i;
1034 +
1035 +  /* Caller must initialize y to the base. */
1036 +  mpz_powm (y, y, q, n);
1037 +
1038 +  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
1039 +    return 1;
1040 +
1041 +  for (i = 1; i < k; i++)
1042 +    {
1043 +      mpz_powm_ui (y, y, 2, n);
1044 +      if (mpz_cmp (y, nm1) == 0)
1045 +       return 1;
1046 +      if (mpz_cmp_ui (y, 1) == 0)
1047 +       return 0;
1048 +    }
1049 +  return 0;
1050 +}
1051 +
1052 +/* This product is 0xc0cfd797, and fits in 32 bits. */
1053 +#define GMP_PRIME_PRODUCT \
1054 +  (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
1055 +
1056 +/* Bit (p+1)/2 is set, for each odd prime <= 61 */
1057 +#define GMP_PRIME_MASK 0xc96996dcUL
1058 +
1059 +int
1060 +mpz_probab_prime_p (const mpz_t n, int reps)
1061 +{
1062 +  mpz_t nm1;
1063 +  mpz_t q;
1064 +  mpz_t y;
1065 +  mp_bitcnt_t k;
1066 +  int is_prime;
1067 +  int j;
1068 +
1069 +  /* Note that we use the absolute value of n only, for compatibility
1070 +     with the real GMP. */
1071 +  if (mpz_even_p (n))
1072 +    return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
1073 +
1074 +  /* Above test excludes n == 0 */
1075 +  assert (n->_mp_size != 0);
1076 +
1077 +  if (mpz_cmpabs_ui (n, 64) < 0)
1078 +    return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
1079 +
1080 +  if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
1081 +    return 0;
1082 +
1083 +  /* All prime factors are >= 31. */
1084 +  if (mpz_cmpabs_ui (n, 31*31) < 0)
1085 +    return 2;
1086 +
1087 +  /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
1088 +     j^2 + j + 41 using Euler's polynomial. We potentially stop early,
1089 +     if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
1090 +     30 (a[30] == 971 > 31*31 == 961). */
1091 +
1092 +  mpz_init (nm1);
1093 +  mpz_init (q);
1094 +  mpz_init (y);
1095 +
1096 +  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
1097 +  mpz_abs (nm1, n);
1098 +  mpz_sub_ui (nm1, nm1, 1);
1099 +  k = mpz_scan1 (nm1, 0);
1100 +  mpz_tdiv_q_2exp (q, nm1, k);
1101 +
1102 +  for (j = 0, is_prime = 1; is_prime && j < reps; j++)
1103 +    {
1104 +      mpz_set_ui (y, (unsigned long) j*j+j+41);
1105 +      if (mpz_cmp (y, nm1) >= 0)
1106 +       {
1107 +         /* Don't try any further bases. */
1108 +         assert (j >= 30);
1109 +         break;
1110 +       }
1111 +      is_prime &= gmp_millerrabin (n, nm1, y, q, k);
1112 +    }
1113 +  mpz_clear (nm1);
1114 +  mpz_clear (q);
1115 +  mpz_clear (y);
1116 +
1117 +  return is_prime;
1118 +}
1119 +
1120 +\f
1121  /* Logical operations and bit manipulation. */
1122
1123  /* Numbers are treated as if represented in two's complement (and
1124     infinitely sign extended). For a negative values we get the two's
1125 -   complement from -x = ~x + 1, where ~ is bitwise complementt.
1126 +   complement from -x = ~x + 1, where ~ is bitwise complement.
1127     Negation transforms
1128
1129       xxxx10...0
1130 @@ -3374,7 +3649,8 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
1131    up = u->_mp_d;
1132    vp = v->_mp_d;
1133
1134 -  for (i = 0; i < vn; i++)
1135 +  i = 0;
1136 +  do
1137      {
1138        ul = (up[i] ^ ux) + uc;
1139        uc = ul < uc;
1140 @@ -3386,6 +3662,7 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
1141        rc = rl < rc;
1142        rp[i] = rl;
1143      }
1144 +  while (++i < vn);
1145    assert (vc == 0);
1146
1147    for (; i < rn; i++)
1148 @@ -3445,7 +3722,8 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
1149    up = u->_mp_d;
1150    vp = v->_mp_d;
1151
1152 -  for (i = 0; i < vn; i++)
1153 +  i = 0;
1154 +  do
1155      {
1156        ul = (up[i] ^ ux) + uc;
1157        uc = ul < uc;
1158 @@ -3457,6 +3735,7 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
1159        rc = rl < rc;
1160        rp[i] = rl;
1161      }
1162 +  while (++i < vn);
1163    assert (vc == 0);
1164
1165    for (; i < rn; i++)
1166 @@ -3512,7 +3791,8 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
1167    up = u->_mp_d;
1168    vp = v->_mp_d;
1169
1170 -  for (i = 0; i < vn; i++)
1171 +  i = 0;
1172 +  do
1173      {
1174        ul = (up[i] ^ ux) + uc;
1175        uc = ul < uc;
1176 @@ -3524,6 +3804,7 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
1177        rc = rl < rc;
1178        rp[i] = rl;
1179      }
1180 +  while (++i < vn);
1181    assert (vc == 0);
1182
1183    for (; i < un; i++)
1184 @@ -3561,20 +3842,28 @@ gmp_popcount_limb (mp_limb_t x)
1185  }
1186
1187  mp_bitcnt_t
1188 -mpz_popcount (const mpz_t u)
1189 +mpn_popcount (mp_srcptr p, mp_size_t n)
1190  {
1191 -  mp_size_t un, i;
1192 +  mp_size_t i;
1193    mp_bitcnt_t c;
1194
1195 +  for (c = 0, i = 0; i < n; i++)
1196 +    c += gmp_popcount_limb (p[i]);
1197 +
1198 +  return c;
1199 +}
1200 +
1201 +mp_bitcnt_t
1202 +mpz_popcount (const mpz_t u)
1203 +{
1204 +  mp_size_t un;
1205 +
1206    un = u->_mp_size;
1207
1208    if (un < 0)
1209      return ~(mp_bitcnt_t) 0;
1210
1211 -  for (c = 0, i = 0; i < un; i++)
1212 -    c += gmp_popcount_limb (u->_mp_d[i]);
1213 -
1214 -  return c;
1215 +  return mpn_popcount (u->_mp_d, un);
1216  }
1217
1218  mp_bitcnt_t
1219 @@ -3591,16 +3880,13 @@ mpz_hamdist (const mpz_t u, const mpz_t v)
1220    if ( (un ^ vn) < 0)
1221      return ~(mp_bitcnt_t) 0;
1222
1223 -  if (un < 0)
1224 +  comp = - (uc = vc = (un < 0));
1225 +  if (uc)
1226      {
1227        assert (vn < 0);
1228        un = -un;
1229        vn = -vn;
1230 -      uc = vc = 1;
1231 -      comp = - (mp_limb_t) 1;
1232      }
1233 -  else
1234 -    uc = vc = comp = 0;
1235
1236    up = u->_mp_d;
1237    vp = v->_mp_d;
1238 @@ -3636,10 +3922,8 @@ mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
1239  {
1240    mp_ptr up;
1241    mp_size_t us, un, i;
1242 -  mp_limb_t limb, ux, uc;
1243 -  unsigned cnt;
1244 +  mp_limb_t limb, ux;
1245
1246 -  up = u->_mp_d;
1247    us = u->_mp_size;
1248    un = GMP_ABS (us);
1249    i = starting_bit / GMP_LIMB_BITS;
1250 @@ -3649,36 +3933,24 @@ mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
1251    if (i >= un)
1252      return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
1253
1254 -  if (us < 0)
1255 -    {
1256 -      ux = GMP_LIMB_MAX;
1257 -      uc = mpn_zero_p (up, i);
1258 -    }
1259 -  else
1260 -    ux = uc = 0;
1261 -
1262 -  limb = (ux ^ up[i]) + uc;
1263 -  uc = limb < uc;
1264 -
1265 -  /* Mask to 0 all bits before starting_bit, thus ignoring them. */
1266 -  limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
1267 +  up = u->_mp_d;
1268 +  ux = 0;
1269 +  limb = up[i];
1270
1271 -  while (limb == 0)
1272 +  if (starting_bit != 0)
1273      {
1274 -      i++;
1275 -      if (i == un)
1276 +      if (us < 0)
1277         {
1278 -         assert (uc == 0);
1279 -         /* For the u > 0 case, this can happen only for the first
1280 -            masked limb. For the u < 0 case, it happens when the
1281 -            highest limbs of the absolute value are all ones. */
1282 -         return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
1283 +         ux = mpn_zero_p (up, i);
1284 +         limb = ~ limb + ux;
1285 +         ux = - (mp_limb_t) (limb >= ux);
1286         }
1287 -      limb = (ux ^ up[i]) + uc;
1288 -      uc = limb < uc;
1289 +
1290 +      /* Mask to 0 all bits before starting_bit, thus ignoring them. */
1291 +      limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
1292      }
1293 -  gmp_ctz (cnt, limb);
1294 -  return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
1295 +
1296 +  return mpn_common_scan (limb, i, up, un, ux);
1297  }
1298
1299  mp_bitcnt_t
1300 @@ -3686,46 +3958,28 @@ mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
1301  {
1302    mp_ptr up;
1303    mp_size_t us, un, i;
1304 -  mp_limb_t limb, ux, uc;
1305 -  unsigned cnt;
1306 +  mp_limb_t limb, ux;
1307
1308 -  up = u->_mp_d;
1309    us = u->_mp_size;
1310 +  ux = - (mp_limb_t) (us >= 0);
1311    un = GMP_ABS (us);
1312    i = starting_bit / GMP_LIMB_BITS;
1313
1314    /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
1315       u<0.  Notice this test picks up all cases of u==0 too. */
1316    if (i >= un)
1317 -    return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0);
1318 +    return (ux ? starting_bit : ~(mp_bitcnt_t) 0);
1319
1320 -  if (us < 0)
1321 -    {
1322 -      ux = GMP_LIMB_MAX;
1323 -      uc = mpn_zero_p (up, i);
1324 -    }
1325 -  else
1326 -    ux = uc = 0;
1327 +  up = u->_mp_d;
1328 +  limb = up[i] ^ ux;
1329
1330 -  limb = (ux ^ up[i]) + uc;
1331 -  uc = limb < uc;
1332 +  if (ux == 0)
1333 +    limb -= mpn_zero_p (up, i); /* limb = ~(~limb + zero_p) */
1334
1335 -  /* Mask to 1 all bits before starting_bit, thus ignoring them. */
1336 -  limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1;
1337 +  /* Mask all bits before starting_bit, thus ignoring them. */
1338 +  limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
1339
1340 -  while (limb == GMP_LIMB_MAX)
1341 -    {
1342 -      i++;
1343 -      if (i == un)
1344 -       {
1345 -         assert (uc == 0);
1346 -         return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0);
1347 -       }
1348 -      limb = (ux ^ up[i]) + uc;
1349 -      uc = limb < uc;
1350 -    }
1351 -  gmp_ctz (cnt, ~limb);
1352 -  return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
1353 +  return mpn_common_scan (limb, i, up, un, ux);
1354  }
1355
1356  \f
1357 @@ -3771,11 +4025,15 @@ mpz_sizeinbase (const mpz_t u, int base)
1358    mpn_copyi (tp, up, un);
1359    mpn_div_qr_1_invert (&bi, base);
1360
1361 -  for (ndigits = 0; un > 0; ndigits++)
1362 +  ndigits = 0;
1363 +  do
1364      {
1365 +      ndigits++;
1366        mpn_div_qr_1_preinv (tp, tp, un, &bi);
1367        un -= (tp[un-1] == 0);
1368      }
1369 +  while (un > 0);
1370 +
1371    gmp_free (tp);
1372    return ndigits;
1373  }
1374 @@ -3852,7 +4110,6 @@ mpz_set_str (mpz_t r, const char *sp, int base)
1375    mp_size_t rn, alloc;
1376    mp_ptr rp;
1377    size_t sn;
1378 -  size_t dn;
1379    int sign;
1380    unsigned char *dp;
1381
1382 @@ -3861,13 +4118,8 @@ mpz_set_str (mpz_t r, const char *sp, int base)
1383    while (isspace( (unsigned char) *sp))
1384      sp++;
1385
1386 -  if (*sp == '-')
1387 -    {
1388 -      sign = 1;
1389 -      sp++;
1390 -    }
1391 -  else
1392 -    sign = 0;
1393 +  sign = (*sp == '-');
1394 +  sp += sign;
1395
1396    if (base == 0)
1397      {
1398 @@ -3894,7 +4146,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
1399    sn = strlen (sp);
1400    dp = gmp_xalloc (sn + (sn == 0));
1401
1402 -  for (dn = 0; *sp; sp++)
1403 +  for (sn = 0; *sp; sp++)
1404      {
1405        unsigned digit;
1406
1407 @@ -3916,7 +4168,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
1408           return -1;
1409         }
1410
1411 -      dp[dn++] = digit;
1412 +      dp[sn++] = digit;
1413      }
1414
1415    bits = mpn_base_power_of_two_p (base);
1416 @@ -3925,7 +4177,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
1417      {
1418        alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
1419        rp = MPZ_REALLOC (r, alloc);
1420 -      rn = mpn_set_str_bits (rp, dp, dn, bits);
1421 +      rn = mpn_set_str_bits (rp, dp, sn, bits);
1422      }
1423    else
1424      {
1425 @@ -3933,7 +4185,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
1426        mpn_get_base_info (&info, base);
1427        alloc = (sn + info.exp - 1) / info.exp;
1428        rp = MPZ_REALLOC (r, alloc);
1429 -      rn = mpn_set_str_other (rp, dp, dn, base, &info);
1430 +      rn = mpn_set_str_other (rp, dp, sn, base, &info);
1431      }
1432    assert (rn <= alloc);
1433    gmp_free (dp);
1434 @@ -3967,14 +4219,9 @@ mpz_out_str (FILE *stream, int base, const mpz_t x)
1435  static int
1436  gmp_detect_endian (void)
1437  {
1438 -  static const int i = 1;
1439 +  static const int i = 2;
1440    const unsigned char *p = (const unsigned char *) &i;
1441 -  if (*p == 1)
1442 -    /* Little endian */
1443 -    return -1;
1444 -  else
1445 -    /* Big endian */
1446 -    return 1;
1447 +  return 1 - *p;
1448  }
1449
1450  /* Import and export. Does not support nails. */
1451 @@ -4037,29 +4284,22 @@ mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
1452             }
1453         }
1454      }
1455 -  if (bytes > 0)
1456 +  assert (i + (bytes > 0) == rn);
1457 +  if (limb != 0)
1458      rp[i++] = limb;
1459 -  assert (i == rn);
1460 +  else
1461 +    i = mpn_normalized_size (rp, i);
1462
1463 -  r->_mp_size = mpn_normalized_size (rp, i);
1464 +  r->_mp_size = i;
1465  }
1466
1467  void *
1468  mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
1469             size_t nails, const mpz_t u)
1470  {
1471 -  unsigned char *p;
1472 -  ptrdiff_t word_step;
1473 -  size_t count, k;
1474 +  size_t count;
1475    mp_size_t un;
1476
1477 -  /* The current (partial) limb. */
1478 -  mp_limb_t limb;
1479 -  /* The number of bytes left to to in this limb. */
1480 -  size_t bytes;
1481 -  /* The index where the limb was read. */
1482 -  mp_size_t i;
1483 -
1484    if (nails != 0)
1485      gmp_die ("mpz_import: Nails not supported.");
1486
1487 @@ -4067,62 +4307,74 @@ mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
1488    assert (endian >= -1 && endian <= 1);
1489    assert (size > 0 || u->_mp_size == 0);
1490
1491 -  un = GMP_ABS (u->_mp_size);
1492 -  if (un == 0)
1493 -    {
1494 -      if (countp)
1495 -       *countp = 0;
1496 -      return r;
1497 -    }
1498 +  un = u->_mp_size;
1499 +  count = 0;
1500 +  if (un != 0)
1501 +    {
1502 +      size_t k;
1503 +      unsigned char *p;
1504 +      ptrdiff_t word_step;
1505 +      /* The current (partial) limb. */
1506 +      mp_limb_t limb;
1507 +      /* The number of bytes left to to in this limb. */
1508 +      size_t bytes;
1509 +      /* The index where the limb was read. */
1510 +      mp_size_t i;
1511
1512 -  /* Count bytes in top limb. */
1513 -  for (limb = u->_mp_d[un-1], k = 0; limb > 0; k++, limb >>= CHAR_BIT)
1514 -    ;
1515 +      un = GMP_ABS (un);
1516
1517 -  assert (k > 0);
1518 +      /* Count bytes in top limb. */
1519 +      limb = u->_mp_d[un-1];
1520 +      assert (limb != 0);
1521
1522 -  count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
1523 +      k = 0;
1524 +      do {
1525 +       k++; limb >>= CHAR_BIT;
1526 +      } while (limb != 0);
1527
1528 -  if (!r)
1529 -    r = gmp_xalloc (count * size);
1530 +      count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
1531
1532 -  if (endian == 0)
1533 -    endian = gmp_detect_endian ();
1534 +      if (!r)
1535 +       r = gmp_xalloc (count * size);
1536
1537 -  p = (unsigned char *) r;
1538 +      if (endian == 0)
1539 +       endian = gmp_detect_endian ();
1540
1541 -  word_step = (order != endian) ? 2 * size : 0;
1542 +      p = (unsigned char *) r;
1543
1544 -  /* Process bytes from the least significant end, so point p at the
1545 -     least significant word. */
1546 -  if (order == 1)
1547 -    {
1548 -      p += size * (count - 1);
1549 -      word_step = - word_step;
1550 -    }
1551 +      word_step = (order != endian) ? 2 * size : 0;
1552
1553 -  /* And at least significant byte of that word. */
1554 -  if (endian == 1)
1555 -    p += (size - 1);
1556 +      /* Process bytes from the least significant end, so point p at the
1557 +        least significant word. */
1558 +      if (order == 1)
1559 +       {
1560 +         p += size * (count - 1);
1561 +         word_step = - word_step;
1562 +       }
1563
1564 -  for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
1565 -      {
1566 -       size_t j;
1567 -       for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
1568 -         {
1569 -           if (bytes == 0)
1570 -             {
1571 -               if (i < un)
1572 -                 limb = u->_mp_d[i++];
1573 -               bytes = sizeof (mp_limb_t);
1574 -             }
1575 -           *p = limb;
1576 -           limb >>= CHAR_BIT;
1577 -           bytes--;
1578 -         }
1579 -      }
1580 -  assert (i == un);
1581 -  assert (k == count);
1582 +      /* And at least significant byte of that word. */
1583 +      if (endian == 1)
1584 +       p += (size - 1);
1585 +
1586 +      for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
1587 +       {
1588 +         size_t j;
1589 +         for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
1590 +           {
1591 +             if (bytes == 0)
1592 +               {
1593 +                 if (i < un)
1594 +                   limb = u->_mp_d[i++];
1595 +                 bytes = sizeof (mp_limb_t);
1596 +               }
1597 +             *p = limb;
1598 +             limb >>= CHAR_BIT;
1599 +             bytes--;
1600 +           }
1601 +       }
1602 +      assert (i == un);
1603 +      assert (k == count);
1604 +    }
1605
1606    if (countp)
1607      *countp = count;
1608 diff --git a/mini-gmp.h b/mini-gmp.h
1609 index 8c94ca2..d8f691f 100644
1610 --- a/mini-gmp.h
1611 +++ b/mini-gmp.h
1612 @@ -1,21 +1,32 @@
1613  /* mini-gmp, a minimalistic implementation of a GNU GMP subset.
1614
1615 -Copyright 2011, 2012, 2013 Free Software Foundation, Inc.
1616 +Copyright 2011-2014 Free Software Foundation, Inc.
1617
1618  This file is part of the GNU MP Library.
1619
1620  The GNU MP Library is free software; you can redistribute it and/or modify
1621 -it under the terms of the GNU Lesser General Public License as published by
1622 -the Free Software Foundation; either version 3 of the License, or (at your
1623 -option) any later version.
1624 +it under the terms of either:
1625 +
1626 +  * the GNU Lesser General Public License as published by the Free
1627 +    Software Foundation; either version 3 of the License, or (at your
1628 +    option) any later version.
1629 +
1630 +or
1631 +
1632 +  * the GNU General Public License as published by the Free Software
1633 +    Foundation; either version 2 of the License, or (at your option) any
1634 +    later version.
1635 +
1636 +or both in parallel, as here.
1637
1638  The GNU MP Library is distributed in the hope that it will be useful, but
1639  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1640 -or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1641 -License for more details.
1642 +or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
1643 +for more details.
1644
1645 -You should have received a copy of the GNU Lesser General Public License
1646 -along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
1647 +You should have received copies of the GNU General Public License and the
1648 +GNU Lesser General Public License along with the GNU MP Library.  If not,
1649 +see https://www.gnu.org/licenses/.  */
1650
1651  /* About mini-gmp: This is a minimal implementation of a subset of the
1652     GMP interface. It is intended for inclusion into applications which
1653 @@ -64,8 +75,11 @@ typedef __mpz_struct mpz_t[1];
1654  typedef __mpz_struct *mpz_ptr;
1655  typedef const __mpz_struct *mpz_srcptr;
1656
1657 +extern const int mp_bits_per_limb;
1658 +
1659  void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t);
1660  void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t);
1661 +void mpn_zero (mp_ptr, mp_size_t);
1662
1663  int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t);
1664
1665 @@ -84,10 +98,17 @@ mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t);
1666  mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
1667  void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
1668  void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t);
1669 +int mpn_perfect_square_p (mp_srcptr, mp_size_t);
1670 +mp_size_t mpn_sqrtrem (mp_ptr, mp_ptr, mp_srcptr, mp_size_t);
1671
1672  mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
1673  mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
1674
1675 +mp_bitcnt_t mpn_scan0 (mp_srcptr, mp_bitcnt_t);
1676 +mp_bitcnt_t mpn_scan1 (mp_srcptr, mp_bitcnt_t);
1677 +
1678 +mp_bitcnt_t mpn_popcount (mp_srcptr, mp_size_t);
1679 +
1680  mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t);
1681  #define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0)
1682
1683 @@ -124,6 +145,10 @@ void mpz_mul_si (mpz_t, const mpz_t, long int);
1684  void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int);
1685  void mpz_mul (mpz_t, const mpz_t, const mpz_t);
1686  void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
1687 +void mpz_addmul_ui (mpz_t, const mpz_t, unsigned long int);
1688 +void mpz_addmul (mpz_t, const mpz_t, const mpz_t);
1689 +void mpz_submul_ui (mpz_t, const mpz_t, unsigned long int);
1690 +void mpz_submul (mpz_t, const mpz_t, const mpz_t);
1691
1692  void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
1693  void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
1694 @@ -147,6 +172,7 @@ void mpz_mod (mpz_t, const mpz_t, const mpz_t);
1695  void mpz_divexact (mpz_t, const mpz_t, const mpz_t);
1696
1697  int mpz_divisible_p (const mpz_t, const mpz_t);
1698 +int mpz_congruent_p (const mpz_t, const mpz_t, const mpz_t);
1699
1700  unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
1701  unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
1702 @@ -176,6 +202,7 @@ int mpz_invert (mpz_t, const mpz_t, const mpz_t);
1703
1704  void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t);
1705  void mpz_sqrt (mpz_t, const mpz_t);
1706 +int mpz_perfect_square_p (const mpz_t);
1707
1708  void mpz_pow_ui (mpz_t, const mpz_t, unsigned long);
1709  void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long);
1710 @@ -188,6 +215,9 @@ int mpz_root (mpz_t, const mpz_t, unsigned long);
1711  void mpz_fac_ui (mpz_t, unsigned long);
1712  void mpz_bin_uiui (mpz_t, unsigned long, unsigned long);
1713
1714 +int
1715 +mpz_probab_prime_p (const mpz_t, int);
1716 +
1717  int mpz_tstbit (const mpz_t, mp_bitcnt_t);
1718  void mpz_setbit (mpz_t, mp_bitcnt_t);
1719  void mpz_clrbit (mpz_t, mp_bitcnt_t);
1720 @@ -211,6 +241,15 @@ double mpz_get_d (const mpz_t);
1721  size_t mpz_size (const mpz_t);
1722  mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t);
1723
1724 +void mpz_realloc2 (mpz_t, mp_bitcnt_t);
1725 +mp_srcptr mpz_limbs_read (mpz_srcptr);
1726 +mp_ptr mpz_limbs_modify (mpz_t, mp_size_t);
1727 +mp_ptr mpz_limbs_write (mpz_t, mp_size_t);
1728 +void mpz_limbs_finish (mpz_t, mp_size_t);
1729 +mpz_srcptr mpz_roinit_n (mpz_t, mp_srcptr, mp_size_t);
1730 +
1731 +#define MPZ_ROINIT_N(xp, xs) {{0, (xs),(xp) }}
1732 +
1733  void mpz_set_si (mpz_t, signed long int);
1734  void mpz_set_ui (mpz_t, unsigned long int);
1735  void mpz_set (mpz_t, const mpz_t);
1736 --
1737 1.9.2
1738