src/gmp/mini-gmp.c

   1 /* mini-gmp, a minimalistic implementation of a GNU GMP subset.
   2
   3    Contributed to the GNU project by Niels Möller
   4
   5 Copyright 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1999, 2000, 2001,
   6 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013
   7 Free Software Foundation, Inc.
   8
   9 This file is part of the GNU MP Library.
  10
  11 The GNU MP Library is free software; you can redistribute it and/or modify
  12 it under the terms of the GNU Lesser General Public License as published by
  13 the Free Software Foundation; either version 3 of the License, or (at your
  14 option) any later version.
  15
  16 The GNU MP Library is distributed in the hope that it will be useful, but
  17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19 License for more details.
  20
  21 You should have received a copy of the GNU Lesser General Public License
  22 along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  23
  24 /* NOTE: All functions in this file which are not declared in
  25    mini-gmp.h are internal, and are not intended to be compatible
  26    neither with GMP nor with future versions of mini-gmp. */
  27
  28 /* Much of the material copied from GMP files, including: gmp-impl.h,
  29    longlong.h, mpn/generic/add_n.c, mpn/generic/addmul_1.c,
  30    mpn/generic/lshift.c, mpn/generic/mul_1.c,
  31    mpn/generic/mul_basecase.c, mpn/generic/rshift.c,
  32    mpn/generic/sbpi1_div_qr.c, mpn/generic/sub_n.c,
  33    mpn/generic/submul_1.c. */
  34
  35 #include <assert.h>
  36 #include <ctype.h>
  37 #include <limits.h>
  38 #include <stdio.h>
  39 #include <stdlib.h>
  40 #include <string.h>
  41
  42 #include "mini-gmp.h"
  43
  44 \f
  45 /* Macros */
  46 #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
  47
  48 #define GMP_LIMB_MAX (~ (mp_limb_t) 0)
  49 #define GMP_LIMB_HIGHBIT ((mp_limb_t) 1 << (GMP_LIMB_BITS - 1))
  50
  51 #define GMP_HLIMB_BIT ((mp_limb_t) 1 << (GMP_LIMB_BITS / 2))
  52 #define GMP_LLIMB_MASK (GMP_HLIMB_BIT - 1)
  53
  54 #define GMP_ULONG_BITS (sizeof(unsigned long) * CHAR_BIT)
  55 #define GMP_ULONG_HIGHBIT ((unsigned long) 1 << (GMP_ULONG_BITS - 1))
  56
  57 #define GMP_ABS(x) ((x) >= 0 ? (x) : -(x))
  58 #define GMP_NEG_CAST(T,x) (-((T)((x) + 1) - 1))
  59
  60 #define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
  61 #define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
  62
  63 #define gmp_assert_nocarry(x) do { \
  64     mp_limb_t __cy = x;            \
  65     assert (__cy == 0);            \
  66   } while (0)
  67
  68 #define gmp_clz(count, x) do {                                          \
  69     mp_limb_t __clz_x = (x);                                            \
  70     unsigned __clz_c;                                                   \
  71     for (__clz_c = 0;                                                   \
  72          (__clz_x & ((mp_limb_t) 0xff << (GMP_LIMB_BITS - 8))) == 0;    \
  73          __clz_c += 8)                                                  \
  74       __clz_x <<= 8;                                                    \
  75     for (; (__clz_x & GMP_LIMB_HIGHBIT) == 0; __clz_c++)                \
  76       __clz_x <<= 1;                                                    \
  77     (count) = __clz_c;                                                  \
  78   } while (0)
  79
  80 #define gmp_ctz(count, x) do {                                          \
  81     mp_limb_t __ctz_x = (x);                                            \
  82     unsigned __ctz_c = 0;                                               \
  83     gmp_clz (__ctz_c, __ctz_x & - __ctz_x);                             \
  84     (count) = GMP_LIMB_BITS - 1 - __ctz_c;                              \
  85   } while (0)
  86
  87 #define gmp_add_ssaaaa(sh, sl, ah, al, bh, bl) \
  88   do {                                                                  \
  89     mp_limb_t __x;                                                      \
  90     __x = (al) + (bl);                                                  \
  91     (sh) = (ah) + (bh) + (__x < (al));                                  \
  92     (sl) = __x;                                                         \
  93   } while (0)
  94
  95 #define gmp_sub_ddmmss(sh, sl, ah, al, bh, bl) \
  96   do {                                                                  \
  97     mp_limb_t __x;                                                      \
  98     __x = (al) - (bl);                                                  \
  99     (sh) = (ah) - (bh) - ((al) < (bl));                                 \
 100     (sl) = __x;                                                         \
 101   } while (0)
 102
 103 #define gmp_umul_ppmm(w1, w0, u, v)                                     \
 104   do {                                                                  \
 105     mp_limb_t __x0, __x1, __x2, __x3;                                   \
 106     unsigned __ul, __vl, __uh, __vh;                                    \
 107     mp_limb_t __u = (u), __v = (v);                                     \
 108                                                                         \
 109     __ul = __u & GMP_LLIMB_MASK;                                        \
 110     __uh = __u >> (GMP_LIMB_BITS / 2);                                  \
 111     __vl = __v & GMP_LLIMB_MASK;                                        \
 112     __vh = __v >> (GMP_LIMB_BITS / 2);                                  \
 113                                                                         \
 114     __x0 = (mp_limb_t) __ul * __vl;                                     \
 115     __x1 = (mp_limb_t) __ul * __vh;                                     \
 116     __x2 = (mp_limb_t) __uh * __vl;                                     \
 117     __x3 = (mp_limb_t) __uh * __vh;                                     \
 118                                                                         \
 119     __x1 += __x0 >> (GMP_LIMB_BITS / 2);/* this can't give carry */     \
 120     __x1 += __x2;               /* but this indeed can */               \
 121     if (__x1 < __x2)            /* did we get it? */                    \
 122       __x3 += GMP_HLIMB_BIT;    /* yes, add it in the proper pos. */    \
 123                                                                         \
 124     (w1) = __x3 + (__x1 >> (GMP_LIMB_BITS / 2));                        \
 125     (w0) = (__x1 << (GMP_LIMB_BITS / 2)) + (__x0 & GMP_LLIMB_MASK);     \
 126   } while (0)
 127
 128 #define gmp_udiv_qrnnd_preinv(q, r, nh, nl, d, di)                      \
 129   do {                                                                  \
 130     mp_limb_t _qh, _ql, _r, _mask;                                      \
 131     gmp_umul_ppmm (_qh, _ql, (nh), (di));                               \
 132     gmp_add_ssaaaa (_qh, _ql, _qh, _ql, (nh) + 1, (nl));                \
 133     _r = (nl) - _qh * (d);                                              \
 134     _mask = -(mp_limb_t) (_r > _ql); /* both > and >= are OK */         \
 135     _qh += _mask;                                                       \
 136     _r += _mask & (d);                                                  \
 137     if (_r >= (d))                                                      \
 138       {                                                                 \
 139         _r -= (d);                                                      \
 140         _qh++;                                                          \
 141       }                                                                 \
 142                                                                         \
 143     (r) = _r;                                                           \
 144     (q) = _qh;                                                          \
 145   } while (0)
 146
 147 #define gmp_udiv_qr_3by2(q, r1, r0, n2, n1, n0, d1, d0, dinv)           \
 148   do {                                                                  \
 149     mp_limb_t _q0, _t1, _t0, _mask;                                     \
 150     gmp_umul_ppmm ((q), _q0, (n2), (dinv));                             \
 151     gmp_add_ssaaaa ((q), _q0, (q), _q0, (n2), (n1));                    \
 152                                                                         \
 153     /* Compute the two most significant limbs of n - q'd */             \
 154     (r1) = (n1) - (d1) * (q);                                           \
 155     gmp_sub_ddmmss ((r1), (r0), (r1), (n0), (d1), (d0));                \
 156     gmp_umul_ppmm (_t1, _t0, (d0), (q));                                \
 157     gmp_sub_ddmmss ((r1), (r0), (r1), (r0), _t1, _t0);                  \
 158     (q)++;                                                              \
 159                                                                         \
 160     /* Conditionally adjust q and the remainders */                     \
 161     _mask = - (mp_limb_t) ((r1) >= _q0);                                \
 162     (q) += _mask;                                                       \
 163     gmp_add_ssaaaa ((r1), (r0), (r1), (r0), _mask & (d1), _mask & (d0)); \
 164     if ((r1) >= (d1))                                                   \
 165       {                                                                 \
 166         if ((r1) > (d1) || (r0) >= (d0))                                \
 167           {                                                             \
 168             (q)++;                                                      \
 169             gmp_sub_ddmmss ((r1), (r0), (r1), (r0), (d1), (d0));        \
 170           }                                                             \
 171       }                                                                 \
 172   } while (0)
 173
 174 /* Swap macros. */
 175 #define MP_LIMB_T_SWAP(x, y)                                            \
 176   do {                                                                  \
 177     mp_limb_t __mp_limb_t_swap__tmp = (x);                              \
 178     (x) = (y);                                                          \
 179     (y) = __mp_limb_t_swap__tmp;                                        \
 180   } while (0)
 181 #define MP_SIZE_T_SWAP(x, y)                                            \
 182   do {                                                                  \
 183     mp_size_t __mp_size_t_swap__tmp = (x);                              \
 184     (x) = (y);                                                          \
 185     (y) = __mp_size_t_swap__tmp;                                        \
 186   } while (0)
 187 #define MP_BITCNT_T_SWAP(x,y)                   \
 188   do {                                          \
 189     mp_bitcnt_t __mp_bitcnt_t_swap__tmp = (x);  \
 190     (x) = (y);                                  \
 191     (y) = __mp_bitcnt_t_swap__tmp;              \
 192   } while (0)
 193 #define MP_PTR_SWAP(x, y)                                               \
 194   do {                                                                  \
 195     mp_ptr __mp_ptr_swap__tmp = (x);                                    \
 196     (x) = (y);                                                          \
 197     (y) = __mp_ptr_swap__tmp;                                           \
 198   } while (0)
 199 #define MP_SRCPTR_SWAP(x, y)                                            \
 200   do {                                                                  \
 201     mp_srcptr __mp_srcptr_swap__tmp = (x);                              \
 202     (x) = (y);                                                          \
 203     (y) = __mp_srcptr_swap__tmp;                                        \
 204   } while (0)
 205
 206 #define MPN_PTR_SWAP(xp,xs, yp,ys)                                      \
 207   do {                                                                  \
 208     MP_PTR_SWAP (xp, yp);                                               \
 209     MP_SIZE_T_SWAP (xs, ys);                                            \
 210   } while(0)
 211 #define MPN_SRCPTR_SWAP(xp,xs, yp,ys)                                   \
 212   do {                                                                  \
 213     MP_SRCPTR_SWAP (xp, yp);                                            \
 214     MP_SIZE_T_SWAP (xs, ys);                                            \
 215   } while(0)
 216
 217 #define MPZ_PTR_SWAP(x, y)                                              \
 218   do {                                                                  \
 219     mpz_ptr __mpz_ptr_swap__tmp = (x);                                  \
 220     (x) = (y);                                                          \
 221     (y) = __mpz_ptr_swap__tmp;                                          \
 222   } while (0)
 223 #define MPZ_SRCPTR_SWAP(x, y)                                           \
 224   do {                                                                  \
 225     mpz_srcptr __mpz_srcptr_swap__tmp = (x);                    \
 226     (x) = (y);                                                          \
 227     (y) = __mpz_srcptr_swap__tmp;                                       \
 228   } while (0)
 229
 230 \f
 231 /* Memory allocation and other helper functions. */
 232 static void
 233 gmp_die (const char *msg)
 234 {
 235   fprintf (stderr, "%s\n", msg);
 236   abort();
 237 }
 238
 239 static void *
 240 gmp_default_alloc (size_t size)
 241 {
 242   void *p;
 243
 244   assert (size > 0);
 245
 246   p = malloc (size);
 247   if (!p)
 248     gmp_die("gmp_default_alloc: Virtual memory exhausted.");
 249
 250   return p;
 251 }
 252
 253 static void *
 254 gmp_default_realloc (void *old, size_t old_size, size_t new_size)
 255 {
 256   mp_ptr p;
 257
 258   p = realloc (old, new_size);
 259
 260   if (!p)
 261     gmp_die("gmp_default_realoc: Virtual memory exhausted.");
 262
 263   return p;
 264 }
 265
 266 static void
 267 gmp_default_free (void *p, size_t size)
 268 {
 269   free (p);
 270 }
 271
 272 static void * (*gmp_allocate_func) (size_t) = gmp_default_alloc;
 273 static void * (*gmp_reallocate_func) (void *, size_t, size_t) = gmp_default_realloc;
 274 static void (*gmp_free_func) (void *, size_t) = gmp_default_free;
 275
 276 void
 277 mp_get_memory_functions (void *(**alloc_func) (size_t),
 278                          void *(**realloc_func) (void *, size_t, size_t),
 279                          void (**free_func) (void *, size_t))
 280 {
 281   if (alloc_func)
 282     *alloc_func = gmp_allocate_func;
 283
 284   if (realloc_func)
 285     *realloc_func = gmp_reallocate_func;
 286
 287   if (free_func)
 288     *free_func = gmp_free_func;
 289 }
 290
 291 void
 292 mp_set_memory_functions (void *(*alloc_func) (size_t),
 293                          void *(*realloc_func) (void *, size_t, size_t),
 294                          void (*free_func) (void *, size_t))
 295 {
 296   if (!alloc_func)
 297     alloc_func = gmp_default_alloc;
 298   if (!realloc_func)
 299     realloc_func = gmp_default_realloc;
 300   if (!free_func)
 301     free_func = gmp_default_free;
 302
 303   gmp_allocate_func = alloc_func;
 304   gmp_reallocate_func = realloc_func;
 305   gmp_free_func = free_func;
 306 }
 307
 308 #define gmp_xalloc(size) ((*gmp_allocate_func)((size)))
 309 #define gmp_free(p) ((*gmp_free_func) ((p), 0))
 310
 311 static mp_ptr
 312 gmp_xalloc_limbs (mp_size_t size)
 313 {
 314   return gmp_xalloc (size * sizeof (mp_limb_t));
 315 }
 316
 317 static mp_ptr
 318 gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
 319 {
 320   assert (size > 0);
 321   return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
 322 }
 323
 324 \f
 325 /* MPN interface */
 326
 327 void
 328 mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
 329 {
 330   mp_size_t i;
 331   for (i = 0; i < n; i++)
 332     d[i] = s[i];
 333 }
 334
 335 void
 336 mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
 337 {
 338   while (n-- > 0)
 339     d[n] = s[n];
 340 }
 341
 342 int
 343 mpn_cmp (mp_srcptr ap, mp_srcptr bp, mp_size_t n)
 344 {
 345   for (; n > 0; n--)
 346     {
 347       if (ap[n-1] < bp[n-1])
 348         return -1;
 349       else if (ap[n-1] > bp[n-1])
 350         return 1;
 351     }
 352   return 0;
 353 }
 354
 355 static int
 356 mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
 357 {
 358   if (an > bn)
 359     return 1;
 360   else if (an < bn)
 361     return -1;
 362   else
 363     return mpn_cmp (ap, bp, an);
 364 }
 365
 366 static mp_size_t
 367 mpn_normalized_size (mp_srcptr xp, mp_size_t n)
 368 {
 369   for (; n > 0 && xp[n-1] == 0; n--)
 370     ;
 371   return n;
 372 }
 373
 374 #define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0)
 375
 376 mp_limb_t
 377 mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 378 {
 379   mp_size_t i;
 380
 381   assert (n > 0);
 382
 383   for (i = 0; i < n; i++)
 384     {
 385       mp_limb_t r = ap[i] + b;
 386       /* Carry out */
 387       b = (r < b);
 388       rp[i] = r;
 389     }
 390   return b;
 391 }
 392
 393 mp_limb_t
 394 mpn_add_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
 395 {
 396   mp_size_t i;
 397   mp_limb_t cy;
 398
 399   for (i = 0, cy = 0; i < n; i++)
 400     {
 401       mp_limb_t a, b, r;
 402       a = ap[i]; b = bp[i];
 403       r = a + cy;
 404       cy = (r < cy);
 405       r += b;
 406       cy += (r < b);
 407       rp[i] = r;
 408     }
 409   return cy;
 410 }
 411
 412 mp_limb_t
 413 mpn_add (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
 414 {
 415   mp_limb_t cy;
 416
 417   assert (an >= bn);
 418
 419   cy = mpn_add_n (rp, ap, bp, bn);
 420   if (an > bn)
 421     cy = mpn_add_1 (rp + bn, ap + bn, an - bn, cy);
 422   return cy;
 423 }
 424
 425 mp_limb_t
 426 mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
 427 {
 428   mp_size_t i;
 429
 430   assert (n > 0);
 431
 432   for (i = 0; i < n; i++)
 433     {
 434       mp_limb_t a = ap[i];
 435       /* Carry out */
 436       mp_limb_t cy = a < b;;
 437       rp[i] = a - b;
 438       b = cy;
 439     }
 440   return b;
 441 }
 442
 443 mp_limb_t
 444 mpn_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
 445 {
 446   mp_size_t i;
 447   mp_limb_t cy;
 448
 449   for (i = 0, cy = 0; i < n; i++)
 450     {
 451       mp_limb_t a, b;
 452       a = ap[i]; b = bp[i];
 453       b += cy;
 454       cy = (b < cy);
 455       cy += (a < b);
 456       rp[i] = a - b;
 457     }
 458   return cy;
 459 }
 460
 461 mp_limb_t
 462 mpn_sub (mp_ptr rp, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
 463 {
 464   mp_limb_t cy;
 465
 466   assert (an >= bn);
 467
 468   cy = mpn_sub_n (rp, ap, bp, bn);
 469   if (an > bn)
 470     cy = mpn_sub_1 (rp + bn, ap + bn, an - bn, cy);
 471   return cy;
 472 }
 473
 474 mp_limb_t
 475 mpn_mul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
 476 {
 477   mp_limb_t ul, cl, hpl, lpl;
 478
 479   assert (n >= 1);
 480
 481   cl = 0;
 482   do
 483     {
 484       ul = *up++;
 485       gmp_umul_ppmm (hpl, lpl, ul, vl);
 486
 487       lpl += cl;
 488       cl = (lpl < cl) + hpl;
 489
 490       *rp++ = lpl;
 491     }
 492   while (--n != 0);
 493
 494   return cl;
 495 }
 496
 497 mp_limb_t
 498 mpn_addmul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
 499 {
 500   mp_limb_t ul, cl, hpl, lpl, rl;
 501
 502   assert (n >= 1);
 503
 504   cl = 0;
 505   do
 506     {
 507       ul = *up++;
 508       gmp_umul_ppmm (hpl, lpl, ul, vl);
 509
 510       lpl += cl;
 511       cl = (lpl < cl) + hpl;
 512
 513       rl = *rp;
 514       lpl = rl + lpl;
 515       cl += lpl < rl;
 516       *rp++ = lpl;
 517     }
 518   while (--n != 0);
 519
 520   return cl;
 521 }
 522
 523 mp_limb_t
 524 mpn_submul_1 (mp_ptr rp, mp_srcptr up, mp_size_t n, mp_limb_t vl)
 525 {
 526   mp_limb_t ul, cl, hpl, lpl, rl;
 527
 528   assert (n >= 1);
 529
 530   cl = 0;
 531   do
 532     {
 533       ul = *up++;
 534       gmp_umul_ppmm (hpl, lpl, ul, vl);
 535
 536       lpl += cl;
 537       cl = (lpl < cl) + hpl;
 538
 539       rl = *rp;
 540       lpl = rl - lpl;
 541       cl += lpl > rl;
 542       *rp++ = lpl;
 543     }
 544   while (--n != 0);
 545
 546   return cl;
 547 }
 548
 549 mp_limb_t
 550 mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
 551 {
 552   assert (un >= vn);
 553   assert (vn >= 1);
 554
 555   /* We first multiply by the low order limb. This result can be
 556      stored, not added, to rp. We also avoid a loop for zeroing this
 557      way. */
 558
 559   rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
 560   rp += 1, vp += 1, vn -= 1;
 561
 562   /* Now accumulate the product of up[] and the next higher limb from
 563      vp[]. */
 564
 565   while (vn >= 1)
 566     {
 567       rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
 568       rp += 1, vp += 1, vn -= 1;
 569     }
 570   return rp[un - 1];
 571 }
 572
 573 void
 574 mpn_mul_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
 575 {
 576   mpn_mul (rp, ap, n, bp, n);
 577 }
 578
 579 void
 580 mpn_sqr (mp_ptr rp, mp_srcptr ap, mp_size_t n)
 581 {
 582   mpn_mul (rp, ap, n, ap, n);
 583 }
 584
 585 mp_limb_t
 586 mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
 587 {
 588   mp_limb_t high_limb, low_limb;
 589   unsigned int tnc;
 590   mp_size_t i;
 591   mp_limb_t retval;
 592
 593   assert (n >= 1);
 594   assert (cnt >= 1);
 595   assert (cnt < GMP_LIMB_BITS);
 596
 597   up += n;
 598   rp += n;
 599
 600   tnc = GMP_LIMB_BITS - cnt;
 601   low_limb = *--up;
 602   retval = low_limb >> tnc;
 603   high_limb = (low_limb << cnt);
 604
 605   for (i = n - 1; i != 0; i--)
 606     {
 607       low_limb = *--up;
 608       *--rp = high_limb | (low_limb >> tnc);
 609       high_limb = (low_limb << cnt);
 610     }
 611   *--rp = high_limb;
 612
 613   return retval;
 614 }
 615
 616 mp_limb_t
 617 mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
 618 {
 619   mp_limb_t high_limb, low_limb;
 620   unsigned int tnc;
 621   mp_size_t i;
 622   mp_limb_t retval;
 623
 624   assert (n >= 1);
 625   assert (cnt >= 1);
 626   assert (cnt < GMP_LIMB_BITS);
 627
 628   tnc = GMP_LIMB_BITS - cnt;
 629   high_limb = *up++;
 630   retval = (high_limb << tnc);
 631   low_limb = high_limb >> cnt;
 632
 633   for (i = n - 1; i != 0; i--)
 634     {
 635       high_limb = *up++;
 636       *rp++ = low_limb | (high_limb << tnc);
 637       low_limb = high_limb >> cnt;
 638     }
 639   *rp = low_limb;
 640
 641   return retval;
 642 }
 643
 644 \f
 645 /* MPN division interface. */
 646 mp_limb_t
 647 mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
 648 {
 649   mp_limb_t r, p, m;
 650   unsigned ul, uh;
 651   unsigned ql, qh;
 652
 653   /* First, do a 2/1 inverse. */
 654   /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 <
 655    * B^2 - (B + m) u1 <= u1 */
 656   assert (u1 >= GMP_LIMB_HIGHBIT);
 657
 658   ul = u1 & GMP_LLIMB_MASK;
 659   uh = u1 >> (GMP_LIMB_BITS / 2);
 660
 661   qh = ~u1 / uh;
 662   r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
 663
 664   p = (mp_limb_t) qh * ul;
 665   /* Adjustment steps taken from udiv_qrnnd_c */
 666   if (r < p)
 667     {
 668       qh--;
 669       r += u1;
 670       if (r >= u1) /* i.e. we didn't get carry when adding to r */
 671         if (r < p)
 672           {
 673             qh--;
 674             r += u1;
 675           }
 676     }
 677   r -= p;
 678
 679   /* Do a 3/2 division (with half limb size) */
 680   p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
 681   ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
 682
 683   /* By the 3/2 method, we don't need the high half limb. */
 684   r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
 685
 686   if (r >= (p << (GMP_LIMB_BITS / 2)))
 687     {
 688       ql--;
 689       r += u1;
 690     }
 691   m = ((mp_limb_t) qh << (GMP_LIMB_BITS / 2)) + ql;
 692   if (r >= u1)
 693     {
 694       m++;
 695       r -= u1;
 696     }
 697
 698   if (u0 > 0)
 699     {
 700       mp_limb_t th, tl;
 701       r = ~r;
 702       r += u0;
 703       if (r < u0)
 704         {
 705           m--;
 706           if (r >= u1)
 707             {
 708               m--;
 709               r -= u1;
 710             }
 711           r -= u1;
 712         }
 713       gmp_umul_ppmm (th, tl, u0, m);
 714       r += th;
 715       if (r < th)
 716         {
 717           m--;
 718           if (r > u1 || (r == u1 && tl > u0))
 719             m--;
 720         }
 721     }
 722
 723   return m;
 724 }
 725
 726 struct gmp_div_inverse
 727 {
 728   /* Normalization shift count. */
 729   unsigned shift;
 730   /* Normalized divisor (d0 unused for mpn_div_qr_1) */
 731   mp_limb_t d1, d0;
 732   /* Inverse, for 2/1 or 3/2. */
 733   mp_limb_t di;
 734 };
 735
 736 static void
 737 mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
 738 {
 739   unsigned shift;
 740
 741   assert (d > 0);
 742   gmp_clz (shift, d);
 743   inv->shift = shift;
 744   inv->d1 = d << shift;
 745   inv->di = mpn_invert_limb (inv->d1);
 746 }
 747
 748 static void
 749 mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
 750                      mp_limb_t d1, mp_limb_t d0)
 751 {
 752   unsigned shift;
 753
 754   assert (d1 > 0);
 755   gmp_clz (shift, d1);
 756   inv->shift = shift;
 757   if (shift > 0)
 758     {
 759       d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
 760       d0 <<= shift;
 761     }
 762   inv->d1 = d1;
 763   inv->d0 = d0;
 764   inv->di = mpn_invert_3by2 (d1, d0);
 765 }
 766
 767 static void
 768 mpn_div_qr_invert (struct gmp_div_inverse *inv,
 769                    mp_srcptr dp, mp_size_t dn)
 770 {
 771   assert (dn > 0);
 772
 773   if (dn == 1)
 774     mpn_div_qr_1_invert (inv, dp[0]);
 775   else if (dn == 2)
 776     mpn_div_qr_2_invert (inv, dp[1], dp[0]);
 777   else
 778     {
 779       unsigned shift;
 780       mp_limb_t d1, d0;
 781
 782       d1 = dp[dn-1];
 783       d0 = dp[dn-2];
 784       assert (d1 > 0);
 785       gmp_clz (shift, d1);
 786       inv->shift = shift;
 787       if (shift > 0)
 788         {
 789           d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
 790           d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
 791         }
 792       inv->d1 = d1;
 793       inv->d0 = d0;
 794       inv->di = mpn_invert_3by2 (d1, d0);
 795     }
 796 }
 797
 798 /* Not matching current public gmp interface, rather corresponding to
 799    the sbpi1_div_* functions. */
 800 static mp_limb_t
 801 mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
 802                      const struct gmp_div_inverse *inv)
 803 {
 804   mp_limb_t d, di;
 805   mp_limb_t r;
 806   mp_ptr tp = NULL;
 807
 808   if (inv->shift > 0)
 809     {
 810       tp = gmp_xalloc_limbs (nn);
 811       r = mpn_lshift (tp, np, nn, inv->shift);
 812       np = tp;
 813     }
 814   else
 815     r = 0;
 816
 817   d = inv->d1;
 818   di = inv->di;
 819   while (nn-- > 0)
 820     {
 821       mp_limb_t q;
 822
 823       gmp_udiv_qrnnd_preinv (q, r, r, np[nn], d, di);
 824       if (qp)
 825         qp[nn] = q;
 826     }
 827   if (inv->shift > 0)
 828     gmp_free (tp);
 829
 830   return r >> inv->shift;
 831 }
 832
 833 static mp_limb_t
 834 mpn_div_qr_1 (mp_ptr qp, mp_srcptr np, mp_size_t nn, mp_limb_t d)
 835 {
 836   assert (d > 0);
 837
 838   /* Special case for powers of two. */
 839   if (d > 1 && (d & (d-1)) == 0)
 840     {
 841       unsigned shift;
 842       mp_limb_t r = np[0] & (d-1);
 843       gmp_ctz (shift, d);
 844       if (qp)
 845         mpn_rshift (qp, np, nn, shift);
 846
 847       return r;
 848     }
 849   else
 850     {
 851       struct gmp_div_inverse inv;
 852       mpn_div_qr_1_invert (&inv, d);
 853       return mpn_div_qr_1_preinv (qp, np, nn, &inv);
 854     }
 855 }
 856
 857 static void
 858 mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
 859                      const struct gmp_div_inverse *inv)
 860 {
 861   unsigned shift;
 862   mp_size_t i;
 863   mp_limb_t d1, d0, di, r1, r0;
 864   mp_ptr tp;
 865
 866   assert (nn >= 2);
 867   shift = inv->shift;
 868   d1 = inv->d1;
 869   d0 = inv->d0;
 870   di = inv->di;
 871
 872   if (shift > 0)
 873     {
 874       tp = gmp_xalloc_limbs (nn);
 875       r1 = mpn_lshift (tp, np, nn, shift);
 876       np = tp;
 877     }
 878   else
 879     r1 = 0;
 880
 881   r0 = np[nn - 1];
 882
 883   for (i = nn - 2; i >= 0; i--)
 884     {
 885       mp_limb_t n0, q;
 886       n0 = np[i];
 887       gmp_udiv_qr_3by2 (q, r1, r0, r1, r0, n0, d1, d0, di);
 888
 889       if (qp)
 890         qp[i] = q;
 891     }
 892
 893   if (shift > 0)
 894     {
 895       assert ((r0 << (GMP_LIMB_BITS - shift)) == 0);
 896       r0 = (r0 >> shift) | (r1 << (GMP_LIMB_BITS - shift));
 897       r1 >>= shift;
 898
 899       gmp_free (tp);
 900     }
 901
 902   rp[1] = r1;
 903   rp[0] = r0;
 904 }
 905
 906 #if 0
 907 static void
 908 mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
 909               mp_limb_t d1, mp_limb_t d0)
 910 {
 911   struct gmp_div_inverse inv;
 912   assert (nn >= 2);
 913
 914   mpn_div_qr_2_invert (&inv, d1, d0);
 915   mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
 916 }
 917 #endif
 918
 919 static void
 920 mpn_div_qr_pi1 (mp_ptr qp,
 921                 mp_ptr np, mp_size_t nn, mp_limb_t n1,
 922                 mp_srcptr dp, mp_size_t dn,
 923                 mp_limb_t dinv)
 924 {
 925   mp_size_t i;
 926
 927   mp_limb_t d1, d0;
 928   mp_limb_t cy, cy1;
 929   mp_limb_t q;
 930
 931   assert (dn > 2);
 932   assert (nn >= dn);
 933
 934   d1 = dp[dn - 1];
 935   d0 = dp[dn - 2];
 936
 937   assert ((d1 & GMP_LIMB_HIGHBIT) != 0);
 938   /* Iteration variable is the index of the q limb.
 939    *
 940    * We divide <n1, np[dn-1+i], np[dn-2+i], np[dn-3+i],..., np[i]>
 941    * by            <d1,          d0,        dp[dn-3],  ..., dp[0] >
 942    */
 943
 944   for (i = nn - dn; i >= 0; i--)
 945     {
 946       mp_limb_t n0 = np[dn-1+i];
 947
 948       if (n1 == d1 && n0 == d0)
 949         {
 950           q = GMP_LIMB_MAX;
 951           mpn_submul_1 (np+i, dp, dn, q);
 952           n1 = np[dn-1+i];      /* update n1, last loop's value will now be invalid */
 953         }
 954       else
 955         {
 956           gmp_udiv_qr_3by2 (q, n1, n0, n1, n0, np[dn-2+i], d1, d0, dinv);
 957
 958           cy = mpn_submul_1 (np + i, dp, dn-2, q);
 959
 960           cy1 = n0 < cy;
 961           n0 = n0 - cy;
 962           cy = n1 < cy1;
 963           n1 = n1 - cy1;
 964           np[dn-2+i] = n0;
 965
 966           if (cy != 0)
 967             {
 968               n1 += d1 + mpn_add_n (np + i, np + i, dp, dn - 1);
 969               q--;
 970             }
 971         }
 972
 973       if (qp)
 974         qp[i] = q;
 975     }
 976
 977   np[dn - 1] = n1;
 978 }
 979
 980 static void
 981 mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
 982                    mp_srcptr dp, mp_size_t dn,
 983                    const struct gmp_div_inverse *inv)
 984 {
 985   assert (dn > 0);
 986   assert (nn >= dn);
 987
 988   if (dn == 1)
 989     np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
 990   else if (dn == 2)
 991     mpn_div_qr_2_preinv (qp, np, np, nn, inv);
 992   else
 993     {
 994       mp_limb_t nh;
 995       unsigned shift;
 996
 997       assert (inv->d1 == dp[dn-1]);
 998       assert (inv->d0 == dp[dn-2]);
 999       assert ((inv->d1 & GMP_LIMB_HIGHBIT) != 0);
1000
1001       shift = inv->shift;
1002       if (shift > 0)
1003         nh = mpn_lshift (np, np, nn, shift);
1004       else
1005         nh = 0;
1006
1007       mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1008
1009       if (shift > 0)
1010         gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1011     }
1012 }
1013
1014 static void
1015 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1016 {
1017   struct gmp_div_inverse inv;
1018   mp_ptr tp = NULL;
1019
1020   assert (dn > 0);
1021   assert (nn >= dn);
1022
1023   mpn_div_qr_invert (&inv, dp, dn);
1024   if (dn > 2 && inv.shift > 0)
1025     {
1026       tp = gmp_xalloc_limbs (dn);
1027       gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1028       dp = tp;
1029     }
1030   mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1031   if (tp)
1032     gmp_free (tp);
1033 }
1034
1035 \f
1036 /* MPN base conversion. */
1037 static unsigned
1038 mpn_base_power_of_two_p (unsigned b)
1039 {
1040   switch (b)
1041     {
1042     case 2: return 1;
1043     case 4: return 2;
1044     case 8: return 3;
1045     case 16: return 4;
1046     case 32: return 5;
1047     case 64: return 6;
1048     case 128: return 7;
1049     case 256: return 8;
1050     default: return 0;
1051     }
1052 }
1053
1054 struct mpn_base_info
1055 {
1056   /* bb is the largest power of the base which fits in one limb, and
1057      exp is the corresponding exponent. */
1058   unsigned exp;
1059   mp_limb_t bb;
1060 };
1061
1062 static void
1063 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1064 {
1065   mp_limb_t m;
1066   mp_limb_t p;
1067   unsigned exp;
1068
1069   m = GMP_LIMB_MAX / b;
1070   for (exp = 1, p = b; p <= m; exp++)
1071     p *= b;
1072
1073   info->exp = exp;
1074   info->bb = p;
1075 }
1076
1077 static mp_bitcnt_t
1078 mpn_limb_size_in_base_2 (mp_limb_t u)
1079 {
1080   unsigned shift;
1081
1082   assert (u > 0);
1083   gmp_clz (shift, u);
1084   return GMP_LIMB_BITS - shift;
1085 }
1086
1087 static size_t
1088 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1089 {
1090   unsigned char mask;
1091   size_t sn, j;
1092   mp_size_t i;
1093   int shift;
1094
1095   sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1096         + bits - 1) / bits;
1097
1098   mask = (1U << bits) - 1;
1099
1100   for (i = 0, j = sn, shift = 0; j-- > 0;)
1101     {
1102       unsigned char digit = up[i] >> shift;
1103
1104       shift += bits;
1105
1106       if (shift >= GMP_LIMB_BITS && ++i < un)
1107         {
1108           shift -= GMP_LIMB_BITS;
1109           digit |= up[i] << (bits - shift);
1110         }
1111       sp[j] = digit & mask;
1112     }
1113   return sn;
1114 }
1115
1116 /* We generate digits from the least significant end, and reverse at
1117    the end. */
1118 static size_t
1119 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1120                   const struct gmp_div_inverse *binv)
1121 {
1122   mp_size_t i;
1123   for (i = 0; w > 0; i++)
1124     {
1125       mp_limb_t h, l, r;
1126
1127       h = w >> (GMP_LIMB_BITS - binv->shift);
1128       l = w << binv->shift;
1129
1130       gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1131       assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
1132       r >>= binv->shift;
1133
1134       sp[i] = r;
1135     }
1136   return i;
1137 }
1138
1139 static size_t
1140 mpn_get_str_other (unsigned char *sp,
1141                    int base, const struct mpn_base_info *info,
1142                    mp_ptr up, mp_size_t un)
1143 {
1144   struct gmp_div_inverse binv;
1145   size_t sn;
1146   size_t i;
1147
1148   mpn_div_qr_1_invert (&binv, base);
1149
1150   sn = 0;
1151
1152   if (un > 1)
1153     {
1154       struct gmp_div_inverse bbinv;
1155       mpn_div_qr_1_invert (&bbinv, info->bb);
1156
1157       do
1158         {
1159           mp_limb_t w;
1160           size_t done;
1161           w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1162           un -= (up[un-1] == 0);
1163           done = mpn_limb_get_str (sp + sn, w, &binv);
1164
1165           for (sn += done; done < info->exp; done++)
1166             sp[sn++] = 0;
1167         }
1168       while (un > 1);
1169     }
1170   sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1171
1172   /* Reverse order */
1173   for (i = 0; 2*i + 1 < sn; i++)
1174     {
1175       unsigned char t = sp[i];
1176       sp[i] = sp[sn - i - 1];
1177       sp[sn - i - 1] = t;
1178     }
1179
1180   return sn;
1181 }
1182
1183 size_t
1184 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1185 {
1186   unsigned bits;
1187
1188   assert (un > 0);
1189   assert (up[un-1] > 0);
1190
1191   bits = mpn_base_power_of_two_p (base);
1192   if (bits)
1193     return mpn_get_str_bits (sp, bits, up, un);
1194   else
1195     {
1196       struct mpn_base_info info;
1197
1198       mpn_get_base_info (&info, base);
1199       return mpn_get_str_other (sp, base, &info, up, un);
1200     }
1201 }
1202
1203 static mp_size_t
1204 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1205                   unsigned bits)
1206 {
1207   mp_size_t rn;
1208   size_t j;
1209   unsigned shift;
1210
1211   for (j = sn, rn = 0, shift = 0; j-- > 0; )
1212     {
1213       if (shift == 0)
1214         {
1215           rp[rn++] = sp[j];
1216           shift += bits;
1217         }
1218       else
1219         {
1220           rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1221           shift += bits;
1222           if (shift >= GMP_LIMB_BITS)
1223             {
1224               shift -= GMP_LIMB_BITS;
1225               if (shift > 0)
1226                 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1227             }
1228         }
1229     }
1230   rn = mpn_normalized_size (rp, rn);
1231   return rn;
1232 }
1233
1234 static mp_size_t
1235 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1236                    mp_limb_t b, const struct mpn_base_info *info)
1237 {
1238   mp_size_t rn;
1239   mp_limb_t w;
1240   unsigned first;
1241   unsigned k;
1242   size_t j;
1243
1244   first = 1 + (sn - 1) % info->exp;
1245
1246   j = 0;
1247   w = sp[j++];
1248   for (k = 1; k < first; k++)
1249     w = w * b + sp[j++];
1250
1251   rp[0] = w;
1252
1253   for (rn = (w > 0); j < sn;)
1254     {
1255       mp_limb_t cy;
1256
1257       w = sp[j++];
1258       for (k = 1; k < info->exp; k++)
1259         w = w * b + sp[j++];
1260
1261       cy = mpn_mul_1 (rp, rp, rn, info->bb);
1262       cy += mpn_add_1 (rp, rp, rn, w);
1263       if (cy > 0)
1264         rp[rn++] = cy;
1265     }
1266   assert (j == sn);
1267
1268   return rn;
1269 }
1270
1271 mp_size_t
1272 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1273 {
1274   unsigned bits;
1275
1276   if (sn == 0)
1277     return 0;
1278
1279   bits = mpn_base_power_of_two_p (base);
1280   if (bits)
1281     return mpn_set_str_bits (rp, sp, sn, bits);
1282   else
1283     {
1284       struct mpn_base_info info;
1285
1286       mpn_get_base_info (&info, base);
1287       return mpn_set_str_other (rp, sp, sn, base, &info);
1288     }
1289 }
1290
1291 \f
1292 /* MPZ interface */
1293 void
1294 mpz_init (mpz_t r)
1295 {
1296   r->_mp_alloc = 1;
1297   r->_mp_size = 0;
1298   r->_mp_d = gmp_xalloc_limbs (1);
1299 }
1300
1301 /* The utility of this function is a bit limited, since many functions
1302    assings the result variable using mpz_swap. */
1303 void
1304 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1305 {
1306   mp_size_t rn;
1307
1308   bits -= (bits != 0);          /* Round down, except if 0 */
1309   rn = 1 + bits / GMP_LIMB_BITS;
1310
1311   r->_mp_alloc = rn;
1312   r->_mp_size = 0;
1313   r->_mp_d = gmp_xalloc_limbs (rn);
1314 }
1315
1316 void
1317 mpz_clear (mpz_t r)
1318 {
1319   gmp_free (r->_mp_d);
1320 }
1321
1322 static void *
1323 mpz_realloc (mpz_t r, mp_size_t size)
1324 {
1325   size = GMP_MAX (size, 1);
1326
1327   r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1328   r->_mp_alloc = size;
1329
1330   if (GMP_ABS (r->_mp_size) > size)
1331     r->_mp_size = 0;
1332
1333   return r->_mp_d;
1334 }
1335
1336 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs.  */
1337 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc                  \
1338                           ? mpz_realloc(z,n)                    \
1339                           : (z)->_mp_d)
1340 \f
1341 /* MPZ assignment and basic conversions. */
1342 void
1343 mpz_set_si (mpz_t r, signed long int x)
1344 {
1345   if (x >= 0)
1346     mpz_set_ui (r, x);
1347   else /* (x < 0) */
1348     {
1349       r->_mp_size = -1;
1350       r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x);
1351     }
1352 }
1353
1354 void
1355 mpz_set_ui (mpz_t r, unsigned long int x)
1356 {
1357   if (x > 0)
1358     {
1359       r->_mp_size = 1;
1360       r->_mp_d[0] = x;
1361     }
1362   else
1363     r->_mp_size = 0;
1364 }
1365
1366 void
1367 mpz_set (mpz_t r, const mpz_t x)
1368 {
1369   /* Allow the NOP r == x */
1370   if (r != x)
1371     {
1372       mp_size_t n;
1373       mp_ptr rp;
1374
1375       n = GMP_ABS (x->_mp_size);
1376       rp = MPZ_REALLOC (r, n);
1377
1378       mpn_copyi (rp, x->_mp_d, n);
1379       r->_mp_size = x->_mp_size;
1380     }
1381 }
1382
1383 void
1384 mpz_init_set_si (mpz_t r, signed long int x)
1385 {
1386   mpz_init (r);
1387   mpz_set_si (r, x);
1388 }
1389
1390 void
1391 mpz_init_set_ui (mpz_t r, unsigned long int x)
1392 {
1393   mpz_init (r);
1394   mpz_set_ui (r, x);
1395 }
1396
1397 void
1398 mpz_init_set (mpz_t r, const mpz_t x)
1399 {
1400   mpz_init (r);
1401   mpz_set (r, x);
1402 }
1403
1404 int
1405 mpz_fits_slong_p (const mpz_t u)
1406 {
1407   mp_size_t us = u->_mp_size;
1408
1409   if (us == 0)
1410     return 1;
1411   else if (us == 1)
1412     return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1413   else if (us == -1)
1414     return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1415   else
1416     return 0;
1417 }
1418
1419 int
1420 mpz_fits_ulong_p (const mpz_t u)
1421 {
1422   mp_size_t us = u->_mp_size;
1423
1424   return us == 0 || us == 1;
1425 }
1426
1427 long int
1428 mpz_get_si (const mpz_t u)
1429 {
1430   mp_size_t us = u->_mp_size;
1431
1432   if (us > 0)
1433     return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT);
1434   else if (us < 0)
1435     return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT);
1436   else
1437     return 0;
1438 }
1439
1440 unsigned long int
1441 mpz_get_ui (const mpz_t u)
1442 {
1443   return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1444 }
1445
1446 size_t
1447 mpz_size (const mpz_t u)
1448 {
1449   return GMP_ABS (u->_mp_size);
1450 }
1451
1452 mp_limb_t
1453 mpz_getlimbn (const mpz_t u, mp_size_t n)
1454 {
1455   if (n >= 0 && n < GMP_ABS (u->_mp_size))
1456     return u->_mp_d[n];
1457   else
1458     return 0;
1459 }
1460
1461 \f
1462 /* Conversions and comparison to double. */
1463 void
1464 mpz_set_d (mpz_t r, double x)
1465 {
1466   int sign;
1467   mp_ptr rp;
1468   mp_size_t rn, i;
1469   double B;
1470   double Bi;
1471   mp_limb_t f;
1472
1473   /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1474      zero or infinity. */
1475   if (x == 0.0 || x != x || x == x * 0.5)
1476     {
1477       r->_mp_size = 0;
1478       return;
1479     }
1480
1481   if (x < 0.0)
1482     {
1483       x = - x;
1484       sign = 1;
1485     }
1486   else
1487     sign = 0;
1488
1489   if (x < 1.0)
1490     {
1491       r->_mp_size = 0;
1492       return;
1493     }
1494   B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1495   Bi = 1.0 / B;
1496   for (rn = 1; x >= B; rn++)
1497     x *= Bi;
1498
1499   rp = MPZ_REALLOC (r, rn);
1500
1501   f = (mp_limb_t) x;
1502   x -= f;
1503   assert (x < 1.0);
1504   rp[rn-1] = f;
1505   for (i = rn-1; i-- > 0; )
1506     {
1507       x = B * x;
1508       f = (mp_limb_t) x;
1509       x -= f;
1510       assert (x < 1.0);
1511       rp[i] = f;
1512     }
1513
1514   r->_mp_size = sign ? - rn : rn;
1515 }
1516
1517 void
1518 mpz_init_set_d (mpz_t r, double x)
1519 {
1520   mpz_init (r);
1521   mpz_set_d (r, x);
1522 }
1523
1524 double
1525 mpz_get_d (const mpz_t u)
1526 {
1527   mp_size_t un;
1528   double x;
1529   double B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1530
1531   un = GMP_ABS (u->_mp_size);
1532
1533   if (un == 0)
1534     return 0.0;
1535
1536   x = u->_mp_d[--un];
1537   while (un > 0)
1538     x = B*x + u->_mp_d[--un];
1539
1540   if (u->_mp_size < 0)
1541     x = -x;
1542
1543   return x;
1544 }
1545
1546 int
1547 mpz_cmpabs_d (const mpz_t x, double d)
1548 {
1549   mp_size_t xn;
1550   double B, Bi;
1551   mp_size_t i;
1552
1553   xn = x->_mp_size;
1554   d = GMP_ABS (d);
1555
1556   if (xn != 0)
1557     {
1558       xn = GMP_ABS (xn);
1559
1560       B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1561       Bi = 1.0 / B;
1562
1563       /* Scale d so it can be compared with the top limb. */
1564       for (i = 1; i < xn; i++)
1565         d *= Bi;
1566
1567       if (d >= B)
1568         return -1;
1569
1570       /* Compare floor(d) to top limb, subtract and cancel when equal. */
1571       for (i = xn; i-- > 0;)
1572         {
1573           mp_limb_t f, xl;
1574
1575           f = (mp_limb_t) d;
1576           xl = x->_mp_d[i];
1577           if (xl > f)
1578             return 1;
1579           else if (xl < f)
1580             return -1;
1581           d = B * (d - f);
1582         }
1583     }
1584   return - (d > 0.0);
1585 }
1586
1587 int
1588 mpz_cmp_d (const mpz_t x, double d)
1589 {
1590   if (x->_mp_size < 0)
1591     {
1592       if (d >= 0.0)
1593         return -1;
1594       else
1595         return -mpz_cmpabs_d (x, d);
1596     }
1597   else
1598     {
1599       if (d < 0.0)
1600         return 1;
1601       else
1602         return mpz_cmpabs_d (x, d);
1603     }
1604 }
1605
1606 \f
1607 /* MPZ comparisons and the like. */
1608 int
1609 mpz_sgn (const mpz_t u)
1610 {
1611   mp_size_t usize = u->_mp_size;
1612
1613   if (usize > 0)
1614     return 1;
1615   else if (usize < 0)
1616     return -1;
1617   else
1618     return 0;
1619 }
1620
1621 int
1622 mpz_cmp_si (const mpz_t u, long v)
1623 {
1624   mp_size_t usize = u->_mp_size;
1625
1626   if (usize < -1)
1627     return -1;
1628   else if (v >= 0)
1629     return mpz_cmp_ui (u, v);
1630   else if (usize >= 0)
1631     return 1;
1632   else /* usize == -1 */
1633     {
1634       mp_limb_t ul = u->_mp_d[0];
1635       if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
1636         return -1;
1637       else if ( (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul)
1638         return 1;
1639     }
1640   return 0;
1641 }
1642
1643 int
1644 mpz_cmp_ui (const mpz_t u, unsigned long v)
1645 {
1646   mp_size_t usize = u->_mp_size;
1647
1648   if (usize > 1)
1649     return 1;
1650   else if (usize < 0)
1651     return -1;
1652   else
1653     {
1654       mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
1655       if (ul > v)
1656         return 1;
1657       else if (ul < v)
1658         return -1;
1659     }
1660   return 0;
1661 }
1662
1663 int
1664 mpz_cmp (const mpz_t a, const mpz_t b)
1665 {
1666   mp_size_t asize = a->_mp_size;
1667   mp_size_t bsize = b->_mp_size;
1668
1669   if (asize > bsize)
1670     return 1;
1671   else if (asize < bsize)
1672     return -1;
1673   else if (asize > 0)
1674     return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1675   else if (asize < 0)
1676     return -mpn_cmp (a->_mp_d, b->_mp_d, -asize);
1677   else
1678     return 0;
1679 }
1680
1681 int
1682 mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1683 {
1684   mp_size_t un = GMP_ABS (u->_mp_size);
1685   mp_limb_t ul;
1686
1687   if (un > 1)
1688     return 1;
1689
1690   ul = (un == 1) ? u->_mp_d[0] : 0;
1691
1692   if (ul > v)
1693     return 1;
1694   else if (ul < v)
1695     return -1;
1696   else
1697     return 0;
1698 }
1699
1700 int
1701 mpz_cmpabs (const mpz_t u, const mpz_t v)
1702 {
1703   return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1704                    v->_mp_d, GMP_ABS (v->_mp_size));
1705 }
1706
1707 void
1708 mpz_abs (mpz_t r, const mpz_t u)
1709 {
1710   if (r != u)
1711     mpz_set (r, u);
1712
1713   r->_mp_size = GMP_ABS (r->_mp_size);
1714 }
1715
1716 void
1717 mpz_neg (mpz_t r, const mpz_t u)
1718 {
1719   if (r != u)
1720     mpz_set (r, u);
1721
1722   r->_mp_size = -r->_mp_size;
1723 }
1724
1725 void
1726 mpz_swap (mpz_t u, mpz_t v)
1727 {
1728   MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1729   MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1730   MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1731 }
1732
1733 \f
1734 /* MPZ addition and subtraction */
1735
1736 /* Adds to the absolute value. Returns new size, but doesn't store it. */
1737 static mp_size_t
1738 mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1739 {
1740   mp_size_t an;
1741   mp_ptr rp;
1742   mp_limb_t cy;
1743
1744   an = GMP_ABS (a->_mp_size);
1745   if (an == 0)
1746     {
1747       r->_mp_d[0] = b;
1748       return b > 0;
1749     }
1750
1751   rp = MPZ_REALLOC (r, an + 1);
1752
1753   cy = mpn_add_1 (rp, a->_mp_d, an, b);
1754   rp[an] = cy;
1755   an += (cy > 0);
1756
1757   return an;
1758 }
1759
1760 /* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1761    but doesn't store it. */
1762 static mp_size_t
1763 mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1764 {
1765   mp_size_t an = GMP_ABS (a->_mp_size);
1766   mp_ptr rp = MPZ_REALLOC (r, an);
1767
1768   if (an == 0)
1769     {
1770       rp[0] = b;
1771       return -(b > 0);
1772     }
1773   else if (an == 1 && a->_mp_d[0] < b)
1774     {
1775       rp[0] = b - a->_mp_d[0];
1776       return -1;
1777     }
1778   else
1779     {
1780       gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b));
1781       return mpn_normalized_size (rp, an);
1782     }
1783 }
1784
1785 void
1786 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1787 {
1788   if (a->_mp_size >= 0)
1789     r->_mp_size = mpz_abs_add_ui (r, a, b);
1790   else
1791     r->_mp_size = -mpz_abs_sub_ui (r, a, b);
1792 }
1793
1794 void
1795 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1796 {
1797   if (a->_mp_size < 0)
1798     r->_mp_size = -mpz_abs_add_ui (r, a, b);
1799   else
1800     r->_mp_size = mpz_abs_sub_ui (r, a, b);
1801 }
1802
1803 void
1804 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1805 {
1806   if (b->_mp_size < 0)
1807     r->_mp_size = mpz_abs_add_ui (r, b, a);
1808   else
1809     r->_mp_size = -mpz_abs_sub_ui (r, b, a);
1810 }
1811
1812 static mp_size_t
1813 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1814 {
1815   mp_size_t an = GMP_ABS (a->_mp_size);
1816   mp_size_t bn = GMP_ABS (b->_mp_size);
1817   mp_size_t rn;
1818   mp_ptr rp;
1819   mp_limb_t cy;
1820
1821   rn = GMP_MAX (an, bn);
1822   rp = MPZ_REALLOC (r, rn + 1);
1823   if (an >= bn)
1824     cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1825   else
1826     cy = mpn_add (rp, b->_mp_d, bn, a->_mp_d, an);
1827
1828   rp[rn] = cy;
1829
1830   return rn + (cy > 0);
1831 }
1832
1833 static mp_size_t
1834 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1835 {
1836   mp_size_t an = GMP_ABS (a->_mp_size);
1837   mp_size_t bn = GMP_ABS (b->_mp_size);
1838   int cmp;
1839   mp_ptr rp;
1840
1841   cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1842   if (cmp > 0)
1843     {
1844       rp = MPZ_REALLOC (r, an);
1845       gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1846       return mpn_normalized_size (rp, an);
1847     }
1848   else if (cmp < 0)
1849     {
1850       rp = MPZ_REALLOC (r, bn);
1851       gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1852       return -mpn_normalized_size (rp, bn);
1853     }
1854   else
1855     return 0;
1856 }
1857
1858 void
1859 mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
1860 {
1861   mp_size_t rn;
1862
1863   if ( (a->_mp_size ^ b->_mp_size) >= 0)
1864     rn = mpz_abs_add (r, a, b);
1865   else
1866     rn = mpz_abs_sub (r, a, b);
1867
1868   r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1869 }
1870
1871 void
1872 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
1873 {
1874   mp_size_t rn;
1875
1876   if ( (a->_mp_size ^ b->_mp_size) >= 0)
1877     rn = mpz_abs_sub (r, a, b);
1878   else
1879     rn = mpz_abs_add (r, a, b);
1880
1881   r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1882 }
1883
1884 \f
1885 /* MPZ multiplication */
1886 void
1887 mpz_mul_si (mpz_t r, const mpz_t u, long int v)
1888 {
1889   if (v < 0)
1890     {
1891       mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
1892       mpz_neg (r, r);
1893     }
1894   else
1895     mpz_mul_ui (r, u, (unsigned long int) v);
1896 }
1897
1898 void
1899 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
1900 {
1901   mp_size_t un;
1902   mpz_t t;
1903   mp_ptr tp;
1904   mp_limb_t cy;
1905
1906   un = GMP_ABS (u->_mp_size);
1907
1908   if (un == 0 || v == 0)
1909     {
1910       r->_mp_size = 0;
1911       return;
1912     }
1913
1914   mpz_init2 (t, (un + 1) * GMP_LIMB_BITS);
1915
1916   tp = t->_mp_d;
1917   cy = mpn_mul_1 (tp, u->_mp_d, un, v);
1918   tp[un] = cy;
1919
1920   t->_mp_size = un + (cy > 0);
1921   if (u->_mp_size < 0)
1922     t->_mp_size = - t->_mp_size;
1923
1924   mpz_swap (r, t);
1925   mpz_clear (t);
1926 }
1927
1928 void
1929 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
1930 {
1931   int sign;
1932   mp_size_t un, vn, rn;
1933   mpz_t t;
1934   mp_ptr tp;
1935
1936   un = GMP_ABS (u->_mp_size);
1937   vn = GMP_ABS (v->_mp_size);
1938
1939   if (un == 0 || vn == 0)
1940     {
1941       r->_mp_size = 0;
1942       return;
1943     }
1944
1945   sign = (u->_mp_size ^ v->_mp_size) < 0;
1946
1947   mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
1948
1949   tp = t->_mp_d;
1950   if (un >= vn)
1951     mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
1952   else
1953     mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
1954
1955   rn = un + vn;
1956   rn -= tp[rn-1] == 0;
1957
1958   t->_mp_size = sign ? - rn : rn;
1959   mpz_swap (r, t);
1960   mpz_clear (t);
1961 }
1962
1963 void
1964 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
1965 {
1966   mp_size_t un, rn;
1967   mp_size_t limbs;
1968   unsigned shift;
1969   mp_ptr rp;
1970
1971   un = GMP_ABS (u->_mp_size);
1972   if (un == 0)
1973     {
1974       r->_mp_size = 0;
1975       return;
1976     }
1977
1978   limbs = bits / GMP_LIMB_BITS;
1979   shift = bits % GMP_LIMB_BITS;
1980
1981   rn = un + limbs + (shift > 0);
1982   rp = MPZ_REALLOC (r, rn);
1983   if (shift > 0)
1984     {
1985       mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
1986       rp[rn-1] = cy;
1987       rn -= (cy == 0);
1988     }
1989   else
1990     mpn_copyd (rp + limbs, u->_mp_d, un);
1991
1992   while (limbs > 0)
1993     rp[--limbs] = 0;
1994
1995   r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
1996 }
1997
1998 \f
1999 /* MPZ division */
2000 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2001
2002 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2003 static int
2004 mpz_div_qr (mpz_t q, mpz_t r,
2005             const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2006 {
2007   mp_size_t ns, ds, nn, dn, qs;
2008   ns = n->_mp_size;
2009   ds = d->_mp_size;
2010
2011   if (ds == 0)
2012     gmp_die("mpz_div_qr: Divide by zero.");
2013
2014   if (ns == 0)
2015     {
2016       if (q)
2017         q->_mp_size = 0;
2018       if (r)
2019         r->_mp_size = 0;
2020       return 0;
2021     }
2022
2023   nn = GMP_ABS (ns);
2024   dn = GMP_ABS (ds);
2025
2026   qs = ds ^ ns;
2027
2028   if (nn < dn)
2029     {
2030       if (mode == GMP_DIV_CEIL && qs >= 0)
2031         {
2032           /* q = 1, r = n - d */
2033           if (r)
2034             mpz_sub (r, n, d);
2035           if (q)
2036             mpz_set_ui (q, 1);
2037         }
2038       else if (mode == GMP_DIV_FLOOR && qs < 0)
2039         {
2040           /* q = -1, r = n + d */
2041           if (r)
2042             mpz_add (r, n, d);
2043           if (q)
2044             mpz_set_si (q, -1);
2045         }
2046       else
2047         {
2048           /* q = 0, r = d */
2049           if (r)
2050             mpz_set (r, n);
2051           if (q)
2052             q->_mp_size = 0;
2053         }
2054       return 1;
2055     }
2056   else
2057     {
2058       mp_ptr np, qp;
2059       mp_size_t qn, rn;
2060       mpz_t tq, tr;
2061
2062       mpz_init (tr);
2063       mpz_set (tr, n);
2064       np = tr->_mp_d;
2065
2066       qn = nn - dn + 1;
2067
2068       if (q)
2069         {
2070           mpz_init2 (tq, qn * GMP_LIMB_BITS);
2071           qp = tq->_mp_d;
2072         }
2073       else
2074         qp = NULL;
2075
2076       mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2077
2078       if (qp)
2079         {
2080           qn -= (qp[qn-1] == 0);
2081
2082           tq->_mp_size = qs < 0 ? -qn : qn;
2083         }
2084       rn = mpn_normalized_size (np, dn);
2085       tr->_mp_size = ns < 0 ? - rn : rn;
2086
2087       if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2088         {
2089           if (q)
2090             mpz_sub_ui (tq, tq, 1);
2091           if (r)
2092             mpz_add (tr, tr, d);
2093         }
2094       else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2095         {
2096           if (q)
2097             mpz_add_ui (tq, tq, 1);
2098           if (r)
2099             mpz_sub (tr, tr, d);
2100         }
2101
2102       if (q)
2103         {
2104           mpz_swap (tq, q);
2105           mpz_clear (tq);
2106         }
2107       if (r)
2108         mpz_swap (tr, r);
2109
2110       mpz_clear (tr);
2111
2112       return rn != 0;
2113     }
2114 }
2115
2116 void
2117 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2118 {
2119   mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2120 }
2121
2122 void
2123 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2124 {
2125   mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2126 }
2127
2128 void
2129 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2130 {
2131   mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2132 }
2133
2134 void
2135 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2136 {
2137   mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2138 }
2139
2140 void
2141 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2142 {
2143   mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2144 }
2145
2146 void
2147 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2148 {
2149   mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2150 }
2151
2152 void
2153 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2154 {
2155   mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2156 }
2157
2158 void
2159 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2160 {
2161   mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2162 }
2163
2164 void
2165 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2166 {
2167   mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2168 }
2169
2170 void
2171 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2172 {
2173   if (d->_mp_size >= 0)
2174     mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2175   else
2176     mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2177 }
2178
2179 static void
2180 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2181                 enum mpz_div_round_mode mode)
2182 {
2183   mp_size_t un, qn;
2184   mp_size_t limb_cnt;
2185   mp_ptr qp;
2186   mp_limb_t adjust;
2187
2188   un = u->_mp_size;
2189   if (un == 0)
2190     {
2191       q->_mp_size = 0;
2192       return;
2193     }
2194   limb_cnt = bit_index / GMP_LIMB_BITS;
2195   qn = GMP_ABS (un) - limb_cnt;
2196   bit_index %= GMP_LIMB_BITS;
2197
2198   if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2199     /* Note: Below, the final indexing at limb_cnt is valid because at
2200        that point we have qn > 0. */
2201     adjust = (qn <= 0
2202               || !mpn_zero_p (u->_mp_d, limb_cnt)
2203               || (u->_mp_d[limb_cnt]
2204                   & (((mp_limb_t) 1 << bit_index) - 1)));
2205   else
2206     adjust = 0;
2207
2208   if (qn <= 0)
2209     qn = 0;
2210
2211   else
2212     {
2213       qp = MPZ_REALLOC (q, qn);
2214
2215       if (bit_index != 0)
2216         {
2217           mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2218           qn -= qp[qn - 1] == 0;
2219         }
2220       else
2221         {
2222           mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2223         }
2224     }
2225
2226   q->_mp_size = qn;
2227
2228   mpz_add_ui (q, q, adjust);
2229   if (un < 0)
2230     mpz_neg (q, q);
2231 }
2232
2233 static void
2234 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2235                 enum mpz_div_round_mode mode)
2236 {
2237   mp_size_t us, un, rn;
2238   mp_ptr rp;
2239   mp_limb_t mask;
2240
2241   us = u->_mp_size;
2242   if (us == 0 || bit_index == 0)
2243     {
2244       r->_mp_size = 0;
2245       return;
2246     }
2247   rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2248   assert (rn > 0);
2249
2250   rp = MPZ_REALLOC (r, rn);
2251   un = GMP_ABS (us);
2252
2253   mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2254
2255   if (rn > un)
2256     {
2257       /* Quotient (with truncation) is zero, and remainder is
2258          non-zero */
2259       if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2260         {
2261           /* Have to negate and sign extend. */
2262           mp_size_t i;
2263           mp_limb_t cy;
2264
2265           for (cy = 1, i = 0; i < un; i++)
2266             {
2267               mp_limb_t s = ~u->_mp_d[i] + cy;
2268               cy = s < cy;
2269               rp[i] = s;
2270             }
2271           assert (cy == 0);
2272           for (; i < rn - 1; i++)
2273             rp[i] = GMP_LIMB_MAX;
2274
2275           rp[rn-1] = mask;
2276           us = -us;
2277         }
2278       else
2279         {
2280           /* Just copy */
2281           if (r != u)
2282             mpn_copyi (rp, u->_mp_d, un);
2283
2284           rn = un;
2285         }
2286     }
2287   else
2288     {
2289       if (r != u)
2290         mpn_copyi (rp, u->_mp_d, rn - 1);
2291
2292       rp[rn-1] = u->_mp_d[rn-1] & mask;
2293
2294       if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2295         {
2296           /* If r != 0, compute 2^{bit_count} - r. */
2297           mp_size_t i;
2298
2299           for (i = 0; i < rn && rp[i] == 0; i++)
2300             ;
2301           if (i < rn)
2302             {
2303               /* r > 0, need to flip sign. */
2304               rp[i] = ~rp[i] + 1;
2305               for (i++; i < rn; i++)
2306                 rp[i] = ~rp[i];
2307
2308               rp[rn-1] &= mask;
2309
2310               /* us is not used for anything else, so we can modify it
2311                  here to indicate flipped sign. */
2312               us = -us;
2313             }
2314         }
2315     }
2316   rn = mpn_normalized_size (rp, rn);
2317   r->_mp_size = us < 0 ? -rn : rn;
2318 }
2319
2320 void
2321 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2322 {
2323   mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2324 }
2325
2326 void
2327 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2328 {
2329   mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2330 }
2331
2332 void
2333 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2334 {
2335   mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2336 }
2337
2338 void
2339 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2340 {
2341   mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2342 }
2343
2344 void
2345 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2346 {
2347   mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2348 }
2349
2350 void
2351 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2352 {
2353   mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2354 }
2355
2356 void
2357 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2358 {
2359   gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2360 }
2361
2362 int
2363 mpz_divisible_p (const mpz_t n, const mpz_t d)
2364 {
2365   return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2366 }
2367
2368 static unsigned long
2369 mpz_div_qr_ui (mpz_t q, mpz_t r,
2370                const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2371 {
2372   mp_size_t ns, qn;
2373   mp_ptr qp;
2374   mp_limb_t rl;
2375   mp_size_t rs;
2376
2377   ns = n->_mp_size;
2378   if (ns == 0)
2379     {
2380       if (q)
2381         q->_mp_size = 0;
2382       if (r)
2383         r->_mp_size = 0;
2384       return 0;
2385     }
2386
2387   qn = GMP_ABS (ns);
2388   if (q)
2389     qp = MPZ_REALLOC (q, qn);
2390   else
2391     qp = NULL;
2392
2393   rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d);
2394   assert (rl < d);
2395
2396   rs = rl > 0;
2397   rs = (ns < 0) ? -rs : rs;
2398
2399   if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0)
2400                   || (mode == GMP_DIV_CEIL && ns >= 0)))
2401     {
2402       if (q)
2403         gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1));
2404       rl = d - rl;
2405       rs = -rs;
2406     }
2407
2408   if (r)
2409     {
2410       r->_mp_d[0] = rl;
2411       r->_mp_size = rs;
2412     }
2413   if (q)
2414     {
2415       qn -= (qp[qn-1] == 0);
2416       assert (qn == 0 || qp[qn-1] > 0);
2417
2418       q->_mp_size = (ns < 0) ? - qn : qn;
2419     }
2420
2421   return rl;
2422 }
2423
2424 unsigned long
2425 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2426 {
2427   return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2428 }
2429
2430 unsigned long
2431 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2432 {
2433   return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2434 }
2435
2436 unsigned long
2437 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2438 {
2439   return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2440 }
2441
2442 unsigned long
2443 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2444 {
2445   return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2446 }
2447
2448 unsigned long
2449 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2450 {
2451   return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2452 }
2453
2454 unsigned long
2455 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2456 {
2457   return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2458 }
2459
2460 unsigned long
2461 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2462 {
2463   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2464 }
2465 unsigned long
2466 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2467 {
2468   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2469 }
2470 unsigned long
2471 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2472 {
2473   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2474 }
2475
2476 unsigned long
2477 mpz_cdiv_ui (const mpz_t n, unsigned long d)
2478 {
2479   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2480 }
2481
2482 unsigned long
2483 mpz_fdiv_ui (const mpz_t n, unsigned long d)
2484 {
2485   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2486 }
2487
2488 unsigned long
2489 mpz_tdiv_ui (const mpz_t n, unsigned long d)
2490 {
2491   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2492 }
2493
2494 unsigned long
2495 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2496 {
2497   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2498 }
2499
2500 void
2501 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2502 {
2503   gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2504 }
2505
2506 int
2507 mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2508 {
2509   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2510 }
2511
2512 \f
2513 /* GCD */
2514 static mp_limb_t
2515 mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2516 {
2517   unsigned shift;
2518
2519   assert ( (u | v) > 0);
2520
2521   if (u == 0)
2522     return v;
2523   else if (v == 0)
2524     return u;
2525
2526   gmp_ctz (shift, u | v);
2527
2528   u >>= shift;
2529   v >>= shift;
2530
2531   if ( (u & 1) == 0)
2532     MP_LIMB_T_SWAP (u, v);
2533
2534   while ( (v & 1) == 0)
2535     v >>= 1;
2536
2537   while (u != v)
2538     {
2539       if (u > v)
2540         {
2541           u -= v;
2542           do
2543             u >>= 1;
2544           while ( (u & 1) == 0);
2545         }
2546       else
2547         {
2548           v -= u;
2549           do
2550             v >>= 1;
2551           while ( (v & 1) == 0);
2552         }
2553     }
2554   return u << shift;
2555 }
2556
2557 unsigned long
2558 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2559 {
2560   mp_size_t un;
2561
2562   if (v == 0)
2563     {
2564       if (g)
2565         mpz_abs (g, u);
2566     }
2567   else
2568     {
2569       un = GMP_ABS (u->_mp_size);
2570       if (un != 0)
2571         v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v);
2572
2573       if (g)
2574         mpz_set_ui (g, v);
2575     }
2576
2577   return v;
2578 }
2579
2580 static mp_bitcnt_t
2581 mpz_make_odd (mpz_t r, const mpz_t u)
2582 {
2583   mp_size_t un, rn, i;
2584   mp_ptr rp;
2585   unsigned shift;
2586
2587   un = GMP_ABS (u->_mp_size);
2588   assert (un > 0);
2589
2590   for (i = 0; u->_mp_d[i] == 0; i++)
2591     ;
2592
2593   gmp_ctz (shift, u->_mp_d[i]);
2594
2595   rn = un - i;
2596   rp = MPZ_REALLOC (r, rn);
2597   if (shift > 0)
2598     {
2599       mpn_rshift (rp, u->_mp_d + i, rn, shift);
2600       rn -= (rp[rn-1] == 0);
2601     }
2602   else
2603     mpn_copyi (rp, u->_mp_d + i, rn);
2604
2605   r->_mp_size = rn;
2606   return i * GMP_LIMB_BITS + shift;
2607 }
2608
2609 void
2610 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2611 {
2612   mpz_t tu, tv;
2613   mp_bitcnt_t uz, vz, gz;
2614
2615   if (u->_mp_size == 0)
2616     {
2617       mpz_abs (g, v);
2618       return;
2619     }
2620   if (v->_mp_size == 0)
2621     {
2622       mpz_abs (g, u);
2623       return;
2624     }
2625
2626   mpz_init (tu);
2627   mpz_init (tv);
2628
2629   uz = mpz_make_odd (tu, u);
2630   vz = mpz_make_odd (tv, v);
2631   gz = GMP_MIN (uz, vz);
2632
2633   if (tu->_mp_size < tv->_mp_size)
2634     mpz_swap (tu, tv);
2635
2636   mpz_tdiv_r (tu, tu, tv);
2637   if (tu->_mp_size == 0)
2638     {
2639       mpz_swap (g, tv);
2640     }
2641   else
2642     for (;;)
2643       {
2644         int c;
2645
2646         mpz_make_odd (tu, tu);
2647         c = mpz_cmp (tu, tv);
2648         if (c == 0)
2649           {
2650             mpz_swap (g, tu);
2651             break;
2652           }
2653         if (c < 0)
2654           mpz_swap (tu, tv);
2655
2656         if (tv->_mp_size == 1)
2657           {
2658             mp_limb_t vl = tv->_mp_d[0];
2659             mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2660             mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2661             break;
2662           }
2663         mpz_sub (tu, tu, tv);
2664       }
2665   mpz_clear (tu);
2666   mpz_clear (tv);
2667   mpz_mul_2exp (g, g, gz);
2668 }
2669
2670 void
2671 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2672 {
2673   mpz_t tu, tv, s0, s1, t0, t1;
2674   mp_bitcnt_t uz, vz, gz;
2675   mp_bitcnt_t power;
2676
2677   if (u->_mp_size == 0)
2678     {
2679       /* g = 0 u + sgn(v) v */
2680       signed long sign = mpz_sgn (v);
2681       mpz_abs (g, v);
2682       if (s)
2683         mpz_set_ui (s, 0);
2684       if (t)
2685         mpz_set_si (t, sign);
2686       return;
2687     }
2688
2689   if (v->_mp_size == 0)
2690     {
2691       /* g = sgn(u) u + 0 v */
2692       signed long sign = mpz_sgn (u);
2693       mpz_abs (g, u);
2694       if (s)
2695         mpz_set_si (s, sign);
2696       if (t)
2697         mpz_set_ui (t, 0);
2698       return;
2699     }
2700
2701   mpz_init (tu);
2702   mpz_init (tv);
2703   mpz_init (s0);
2704   mpz_init (s1);
2705   mpz_init (t0);
2706   mpz_init (t1);
2707
2708   uz = mpz_make_odd (tu, u);
2709   vz = mpz_make_odd (tv, v);
2710   gz = GMP_MIN (uz, vz);
2711
2712   uz -= gz;
2713   vz -= gz;
2714
2715   /* Cofactors corresponding to odd gcd. gz handled later. */
2716   if (tu->_mp_size < tv->_mp_size)
2717     {
2718       mpz_swap (tu, tv);
2719       MPZ_SRCPTR_SWAP (u, v);
2720       MPZ_PTR_SWAP (s, t);
2721       MP_BITCNT_T_SWAP (uz, vz);
2722     }
2723
2724   /* Maintain
2725    *
2726    * u = t0 tu + t1 tv
2727    * v = s0 tu + s1 tv
2728    *
2729    * where u and v denote the inputs with common factors of two
2730    * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2731    *
2732    * 2^p tu =  s1 u - t1 v
2733    * 2^p tv = -s0 u + t0 v
2734    */
2735
2736   /* After initial division, tu = q tv + tu', we have
2737    *
2738    * u = 2^uz (tu' + q tv)
2739    * v = 2^vz tv
2740    *
2741    * or
2742    *
2743    * t0 = 2^uz, t1 = 2^uz q
2744    * s0 = 0,    s1 = 2^vz
2745    */
2746
2747   mpz_setbit (t0, uz);
2748   mpz_tdiv_qr (t1, tu, tu, tv);
2749   mpz_mul_2exp (t1, t1, uz);
2750
2751   mpz_setbit (s1, vz);
2752   power = uz + vz;
2753
2754   if (tu->_mp_size > 0)
2755     {
2756       mp_bitcnt_t shift;
2757       shift = mpz_make_odd (tu, tu);
2758       mpz_mul_2exp (t0, t0, shift);
2759       mpz_mul_2exp (s0, s0, shift);
2760       power += shift;
2761
2762       for (;;)
2763         {
2764           int c;
2765           c = mpz_cmp (tu, tv);
2766           if (c == 0)
2767             break;
2768
2769           if (c < 0)
2770             {
2771               /* tv = tv' + tu
2772                *
2773                * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2774                * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2775
2776               mpz_sub (tv, tv, tu);
2777               mpz_add (t0, t0, t1);
2778               mpz_add (s0, s0, s1);
2779
2780               shift = mpz_make_odd (tv, tv);
2781               mpz_mul_2exp (t1, t1, shift);
2782               mpz_mul_2exp (s1, s1, shift);
2783             }
2784           else
2785             {
2786               mpz_sub (tu, tu, tv);
2787               mpz_add (t1, t0, t1);
2788               mpz_add (s1, s0, s1);
2789
2790               shift = mpz_make_odd (tu, tu);
2791               mpz_mul_2exp (t0, t0, shift);
2792               mpz_mul_2exp (s0, s0, shift);
2793             }
2794           power += shift;
2795         }
2796     }
2797
2798   /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2799      cofactors. */
2800
2801   mpz_mul_2exp (tv, tv, gz);
2802   mpz_neg (s0, s0);
2803
2804   /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2805      adjust cofactors, we need u / g and v / g */
2806
2807   mpz_divexact (s1, v, tv);
2808   mpz_abs (s1, s1);
2809   mpz_divexact (t1, u, tv);
2810   mpz_abs (t1, t1);
2811
2812   while (power-- > 0)
2813     {
2814       /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2815       if (mpz_odd_p (s0) || mpz_odd_p (t0))
2816         {
2817           mpz_sub (s0, s0, s1);
2818           mpz_add (t0, t0, t1);
2819         }
2820       mpz_divexact_ui (s0, s0, 2);
2821       mpz_divexact_ui (t0, t0, 2);
2822     }
2823
2824   /* Arrange so that |s| < |u| / 2g */
2825   mpz_add (s1, s0, s1);
2826   if (mpz_cmpabs (s0, s1) > 0)
2827     {
2828       mpz_swap (s0, s1);
2829       mpz_sub (t0, t0, t1);
2830     }
2831   if (u->_mp_size < 0)
2832     mpz_neg (s0, s0);
2833   if (v->_mp_size < 0)
2834     mpz_neg (t0, t0);
2835
2836   mpz_swap (g, tv);
2837   if (s)
2838     mpz_swap (s, s0);
2839   if (t)
2840     mpz_swap (t, t0);
2841
2842   mpz_clear (tu);
2843   mpz_clear (tv);
2844   mpz_clear (s0);
2845   mpz_clear (s1);
2846   mpz_clear (t0);
2847   mpz_clear (t1);
2848 }
2849
2850 void
2851 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
2852 {
2853   mpz_t g;
2854
2855   if (u->_mp_size == 0 || v->_mp_size == 0)
2856     {
2857       r->_mp_size = 0;
2858       return;
2859     }
2860
2861   mpz_init (g);
2862
2863   mpz_gcd (g, u, v);
2864   mpz_divexact (g, u, g);
2865   mpz_mul (r, g, v);
2866
2867   mpz_clear (g);
2868   mpz_abs (r, r);
2869 }
2870
2871 void
2872 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
2873 {
2874   if (v == 0 || u->_mp_size == 0)
2875     {
2876       r->_mp_size = 0;
2877       return;
2878     }
2879
2880   v /= mpz_gcd_ui (NULL, u, v);
2881   mpz_mul_ui (r, u, v);
2882
2883   mpz_abs (r, r);
2884 }
2885
2886 int
2887 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
2888 {
2889   mpz_t g, tr;
2890   int invertible;
2891
2892   if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
2893     return 0;
2894
2895   mpz_init (g);
2896   mpz_init (tr);
2897
2898   mpz_gcdext (g, tr, NULL, u, m);
2899   invertible = (mpz_cmp_ui (g, 1) == 0);
2900
2901   if (invertible)
2902     {
2903       if (tr->_mp_size < 0)
2904         {
2905           if (m->_mp_size >= 0)
2906             mpz_add (tr, tr, m);
2907           else
2908             mpz_sub (tr, tr, m);
2909         }
2910       mpz_swap (r, tr);
2911     }
2912
2913   mpz_clear (g);
2914   mpz_clear (tr);
2915   return invertible;
2916 }
2917
2918 \f
2919 /* Higher level operations (sqrt, pow and root) */
2920
2921 void
2922 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
2923 {
2924   unsigned long bit;
2925   mpz_t tr;
2926   mpz_init_set_ui (tr, 1);
2927
2928   for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1)
2929     {
2930       mpz_mul (tr, tr, tr);
2931       if (e & bit)
2932         mpz_mul (tr, tr, b);
2933     }
2934   mpz_swap (r, tr);
2935   mpz_clear (tr);
2936 }
2937
2938 void
2939 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
2940 {
2941   mpz_t b;
2942   mpz_init_set_ui (b, blimb);
2943   mpz_pow_ui (r, b, e);
2944   mpz_clear (b);
2945 }
2946
2947 void
2948 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
2949 {
2950   mpz_t tr;
2951   mpz_t base;
2952   mp_size_t en, mn;
2953   mp_srcptr mp;
2954   struct gmp_div_inverse minv;
2955   unsigned shift;
2956   mp_ptr tp = NULL;
2957
2958   en = GMP_ABS (e->_mp_size);
2959   mn = GMP_ABS (m->_mp_size);
2960   if (mn == 0)
2961     gmp_die ("mpz_powm: Zero modulo.");
2962
2963   if (en == 0)
2964     {
2965       mpz_set_ui (r, 1);
2966       return;
2967     }
2968
2969   mp = m->_mp_d;
2970   mpn_div_qr_invert (&minv, mp, mn);
2971   shift = minv.shift;
2972
2973   if (shift > 0)
2974     {
2975       /* To avoid shifts, we do all our reductions, except the final
2976          one, using a *normalized* m. */
2977       minv.shift = 0;
2978
2979       tp = gmp_xalloc_limbs (mn);
2980       gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
2981       mp = tp;
2982     }
2983
2984   mpz_init (base);
2985
2986   if (e->_mp_size < 0)
2987     {
2988       if (!mpz_invert (base, b, m))
2989         gmp_die ("mpz_powm: Negative exponent and non-invertibe base.");
2990     }
2991   else
2992     {
2993       mp_size_t bn;
2994       mpz_abs (base, b);
2995
2996       bn = base->_mp_size;
2997       if (bn >= mn)
2998         {
2999           mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3000           bn = mn;
3001         }
3002
3003       /* We have reduced the absolute value. Now take care of the
3004          sign. Note that we get zero represented non-canonically as
3005          m. */
3006       if (b->_mp_size < 0)
3007         {
3008           mp_ptr bp = MPZ_REALLOC (base, mn);
3009           gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3010           bn = mn;
3011         }
3012       base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3013     }
3014   mpz_init_set_ui (tr, 1);
3015
3016   while (en-- > 0)
3017     {
3018       mp_limb_t w = e->_mp_d[en];
3019       mp_limb_t bit;
3020
3021       for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1)
3022         {
3023           mpz_mul (tr, tr, tr);
3024           if (w & bit)
3025             mpz_mul (tr, tr, base);
3026           if (tr->_mp_size > mn)
3027             {
3028               mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3029               tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3030             }
3031         }
3032     }
3033
3034   /* Final reduction */
3035   if (tr->_mp_size >= mn)
3036     {
3037       minv.shift = shift;
3038       mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3039       tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3040     }
3041   if (tp)
3042     gmp_free (tp);
3043
3044   mpz_swap (r, tr);
3045   mpz_clear (tr);
3046   mpz_clear (base);
3047 }
3048
3049 void
3050 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3051 {
3052   mpz_t e;
3053   mpz_init_set_ui (e, elimb);
3054   mpz_powm (r, b, e, m);
3055   mpz_clear (e);
3056 }
3057
3058 /* x=trunc(y^(1/z)), r=y-x^z */
3059 void
3060 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3061 {
3062   int sgn;
3063   mpz_t t, u;
3064
3065   sgn = y->_mp_size < 0;
3066   if (sgn && (z & 1) == 0)
3067     gmp_die ("mpz_rootrem: Negative argument, with even root.");
3068   if (z == 0)
3069     gmp_die ("mpz_rootrem: Zeroth root.");
3070
3071   if (mpz_cmpabs_ui (y, 1) <= 0) {
3072     mpz_set (x, y);
3073     if (r)
3074       r->_mp_size = 0;
3075     return;
3076   }
3077
3078   mpz_init (t);
3079   mpz_init (u);
3080   mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3081
3082   if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3083     do {
3084       mpz_swap (u, t);                  /* u = x */
3085       mpz_tdiv_q (t, y, u);             /* t = y/x */
3086       mpz_add (t, t, u);                /* t = y/x + x */
3087       mpz_tdiv_q_2exp (t, t, 1);        /* x'= (y/x + x)/2 */
3088     } while (mpz_cmpabs (t, u) < 0);    /* |x'| < |x| */
3089   else /* z != 2 */ {
3090     mpz_t v;
3091
3092     mpz_init (v);
3093     if (sgn)
3094       mpz_neg (t, t);
3095
3096     do {
3097       mpz_swap (u, t);                  /* u = x */
3098       mpz_pow_ui (t, u, z - 1);         /* t = x^(z-1) */
3099       mpz_tdiv_q (t, y, t);             /* t = y/x^(z-1) */
3100       mpz_mul_ui (v, u, z - 1);         /* v = x*(z-1) */
3101       mpz_add (t, t, v);                /* t = y/x^(z-1) + x*(z-1) */
3102       mpz_tdiv_q_ui (t, t, z);          /* x'=(y/x^(z-1) + x*(z-1))/z */
3103     } while (mpz_cmpabs (t, u) < 0);    /* |x'| < |x| */
3104
3105     mpz_clear (v);
3106   }
3107
3108   if (r) {
3109     mpz_pow_ui (t, u, z);
3110     mpz_sub (r, y, t);
3111   }
3112   mpz_swap (x, u);
3113   mpz_clear (u);
3114   mpz_clear (t);
3115 }
3116
3117 int
3118 mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3119 {
3120   int res;
3121   mpz_t r;
3122
3123   mpz_init (r);
3124   mpz_rootrem (x, r, y, z);
3125   res = r->_mp_size == 0;
3126   mpz_clear (r);
3127
3128   return res;
3129 }
3130
3131 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3132 void
3133 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3134 {
3135   mpz_rootrem (s, r, u, 2);
3136 }
3137
3138 void
3139 mpz_sqrt (mpz_t s, const mpz_t u)
3140 {
3141   mpz_rootrem (s, NULL, u, 2);
3142 }
3143
3144 \f
3145 /* Combinatorics */
3146
3147 void
3148 mpz_fac_ui (mpz_t x, unsigned long n)
3149 {
3150   if (n < 2) {
3151     mpz_set_ui (x, 1);
3152     return;
3153   }
3154   mpz_set_ui (x, n);
3155   for (;--n > 1;)
3156     mpz_mul_ui (x, x, n);
3157 }
3158
3159 void
3160 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3161 {
3162   mpz_t t;
3163
3164   if (k > n) {
3165     r->_mp_size = 0;
3166     return;
3167   }
3168   mpz_fac_ui (r, n);
3169   mpz_init (t);
3170   mpz_fac_ui (t, k);
3171   mpz_divexact (r, r, t);
3172   mpz_fac_ui (t, n - k);
3173   mpz_divexact (r, r, t);
3174   mpz_clear (t);
3175 }
3176
3177 \f
3178 /* Logical operations and bit manipulation. */
3179
3180 /* Numbers are treated as if represented in two's complement (and
3181    infinitely sign extended). For a negative values we get the two's
3182    complement from -x = ~x + 1, where ~ is bitwise complementt.
3183    Negation transforms
3184
3185      xxxx10...0
3186
3187    into
3188
3189      yyyy10...0
3190
3191    where yyyy is the bitwise complement of xxxx. So least significant
3192    bits, up to and including the first one bit, are unchanged, and
3193    the more significant bits are all complemented.
3194
3195    To change a bit from zero to one in a negative number, subtract the
3196    corresponding power of two from the absolute value. This can never
3197    underflow. To change a bit from one to zero, add the corresponding
3198    power of two, and this might overflow. E.g., if x = -001111, the
3199    two's complement is 110001. Clearing the least significant bit, we
3200    get two's complement 110000, and -010000. */
3201
3202 int
3203 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3204 {
3205   mp_size_t limb_index;
3206   unsigned shift;
3207   mp_size_t ds;
3208   mp_size_t dn;
3209   mp_limb_t w;
3210   int bit;
3211
3212   ds = d->_mp_size;
3213   dn = GMP_ABS (ds);
3214   limb_index = bit_index / GMP_LIMB_BITS;
3215   if (limb_index >= dn)
3216     return ds < 0;
3217
3218   shift = bit_index % GMP_LIMB_BITS;
3219   w = d->_mp_d[limb_index];
3220   bit = (w >> shift) & 1;
3221
3222   if (ds < 0)
3223     {
3224       /* d < 0. Check if any of the bits below is set: If so, our bit
3225          must be complemented. */
3226       if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3227         return bit ^ 1;
3228       while (limb_index-- > 0)
3229         if (d->_mp_d[limb_index] > 0)
3230           return bit ^ 1;
3231     }
3232   return bit;
3233 }
3234
3235 static void
3236 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3237 {
3238   mp_size_t dn, limb_index;
3239   mp_limb_t bit;
3240   mp_ptr dp;
3241
3242   dn = GMP_ABS (d->_mp_size);
3243
3244   limb_index = bit_index / GMP_LIMB_BITS;
3245   bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3246
3247   if (limb_index >= dn)
3248     {
3249       mp_size_t i;
3250       /* The bit should be set outside of the end of the number.
3251          We have to increase the size of the number. */
3252       dp = MPZ_REALLOC (d, limb_index + 1);
3253
3254       dp[limb_index] = bit;
3255       for (i = dn; i < limb_index; i++)
3256         dp[i] = 0;
3257       dn = limb_index + 1;
3258     }
3259   else
3260     {
3261       mp_limb_t cy;
3262
3263       dp = d->_mp_d;
3264
3265       cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3266       if (cy > 0)
3267         {
3268           dp = MPZ_REALLOC (d, dn + 1);
3269           dp[dn++] = cy;
3270         }
3271     }
3272
3273   d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3274 }
3275
3276 static void
3277 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3278 {
3279   mp_size_t dn, limb_index;
3280   mp_ptr dp;
3281   mp_limb_t bit;
3282
3283   dn = GMP_ABS (d->_mp_size);
3284   dp = d->_mp_d;
3285
3286   limb_index = bit_index / GMP_LIMB_BITS;
3287   bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3288
3289   assert (limb_index < dn);
3290
3291   gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3292                                  dn - limb_index, bit));
3293   dn -= (dp[dn-1] == 0);
3294   d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3295 }
3296
3297 void
3298 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3299 {
3300   if (!mpz_tstbit (d, bit_index))
3301     {
3302       if (d->_mp_size >= 0)
3303         mpz_abs_add_bit (d, bit_index);
3304       else
3305         mpz_abs_sub_bit (d, bit_index);
3306     }
3307 }
3308
3309 void
3310 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3311 {
3312   if (mpz_tstbit (d, bit_index))
3313     {
3314       if (d->_mp_size >= 0)
3315         mpz_abs_sub_bit (d, bit_index);
3316       else
3317         mpz_abs_add_bit (d, bit_index);
3318     }
3319 }
3320
3321 void
3322 mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3323 {
3324   if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3325     mpz_abs_sub_bit (d, bit_index);
3326   else
3327     mpz_abs_add_bit (d, bit_index);
3328 }
3329
3330 void
3331 mpz_com (mpz_t r, const mpz_t u)
3332 {
3333   mpz_neg (r, u);
3334   mpz_sub_ui (r, r, 1);
3335 }
3336
3337 void
3338 mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3339 {
3340   mp_size_t un, vn, rn, i;
3341   mp_ptr up, vp, rp;
3342
3343   mp_limb_t ux, vx, rx;
3344   mp_limb_t uc, vc, rc;
3345   mp_limb_t ul, vl, rl;
3346
3347   un = GMP_ABS (u->_mp_size);
3348   vn = GMP_ABS (v->_mp_size);
3349   if (un < vn)
3350     {
3351       MPZ_SRCPTR_SWAP (u, v);
3352       MP_SIZE_T_SWAP (un, vn);
3353     }
3354   if (vn == 0)
3355     {
3356       r->_mp_size = 0;
3357       return;
3358     }
3359
3360   uc = u->_mp_size < 0;
3361   vc = v->_mp_size < 0;
3362   rc = uc & vc;
3363
3364   ux = -uc;
3365   vx = -vc;
3366   rx = -rc;
3367
3368   /* If the smaller input is positive, higher limbs don't matter. */
3369   rn = vx ? un : vn;
3370
3371   rp = MPZ_REALLOC (r, rn + rc);
3372
3373   up = u->_mp_d;
3374   vp = v->_mp_d;
3375
3376   for (i = 0; i < vn; i++)
3377     {
3378       ul = (up[i] ^ ux) + uc;
3379       uc = ul < uc;
3380
3381       vl = (vp[i] ^ vx) + vc;
3382       vc = vl < vc;
3383
3384       rl = ( (ul & vl) ^ rx) + rc;
3385       rc = rl < rc;
3386       rp[i] = rl;
3387     }
3388   assert (vc == 0);
3389
3390   for (; i < rn; i++)
3391     {
3392       ul = (up[i] ^ ux) + uc;
3393       uc = ul < uc;
3394
3395       rl = ( (ul & vx) ^ rx) + rc;
3396       rc = rl < rc;
3397       rp[i] = rl;
3398     }
3399   if (rc)
3400     rp[rn++] = rc;
3401   else
3402     rn = mpn_normalized_size (rp, rn);
3403
3404   r->_mp_size = rx ? -rn : rn;
3405 }
3406
3407 void
3408 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3409 {
3410   mp_size_t un, vn, rn, i;
3411   mp_ptr up, vp, rp;
3412
3413   mp_limb_t ux, vx, rx;
3414   mp_limb_t uc, vc, rc;
3415   mp_limb_t ul, vl, rl;
3416
3417   un = GMP_ABS (u->_mp_size);
3418   vn = GMP_ABS (v->_mp_size);
3419   if (un < vn)
3420     {
3421       MPZ_SRCPTR_SWAP (u, v);
3422       MP_SIZE_T_SWAP (un, vn);
3423     }
3424   if (vn == 0)
3425     {
3426       mpz_set (r, u);
3427       return;
3428     }
3429
3430   uc = u->_mp_size < 0;
3431   vc = v->_mp_size < 0;
3432   rc = uc | vc;
3433
3434   ux = -uc;
3435   vx = -vc;
3436   rx = -rc;
3437
3438   /* If the smaller input is negative, by sign extension higher limbs
3439      don't matter. */
3440   rn = vx ? vn : un;
3441
3442   rp = MPZ_REALLOC (r, rn + rc);
3443
3444   up = u->_mp_d;
3445   vp = v->_mp_d;
3446
3447   for (i = 0; i < vn; i++)
3448     {
3449       ul = (up[i] ^ ux) + uc;
3450       uc = ul < uc;
3451
3452       vl = (vp[i] ^ vx) + vc;
3453       vc = vl < vc;
3454
3455       rl = ( (ul | vl) ^ rx) + rc;
3456       rc = rl < rc;
3457       rp[i] = rl;
3458     }
3459   assert (vc == 0);
3460
3461   for (; i < rn; i++)
3462     {
3463       ul = (up[i] ^ ux) + uc;
3464       uc = ul < uc;
3465
3466       rl = ( (ul | vx) ^ rx) + rc;
3467       rc = rl < rc;
3468       rp[i] = rl;
3469     }
3470   if (rc)
3471     rp[rn++] = rc;
3472   else
3473     rn = mpn_normalized_size (rp, rn);
3474
3475   r->_mp_size = rx ? -rn : rn;
3476 }
3477
3478 void
3479 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3480 {
3481   mp_size_t un, vn, i;
3482   mp_ptr up, vp, rp;
3483
3484   mp_limb_t ux, vx, rx;
3485   mp_limb_t uc, vc, rc;
3486   mp_limb_t ul, vl, rl;
3487
3488   un = GMP_ABS (u->_mp_size);
3489   vn = GMP_ABS (v->_mp_size);
3490   if (un < vn)
3491     {
3492       MPZ_SRCPTR_SWAP (u, v);
3493       MP_SIZE_T_SWAP (un, vn);
3494     }
3495   if (vn == 0)
3496     {
3497       mpz_set (r, u);
3498       return;
3499     }
3500
3501   uc = u->_mp_size < 0;
3502   vc = v->_mp_size < 0;
3503   rc = uc ^ vc;
3504
3505   ux = -uc;
3506   vx = -vc;
3507   rx = -rc;
3508
3509   rp = MPZ_REALLOC (r, un + rc);
3510
3511   up = u->_mp_d;
3512   vp = v->_mp_d;
3513
3514   for (i = 0; i < vn; i++)
3515     {
3516       ul = (up[i] ^ ux) + uc;
3517       uc = ul < uc;
3518
3519       vl = (vp[i] ^ vx) + vc;
3520       vc = vl < vc;
3521
3522       rl = (ul ^ vl ^ rx) + rc;
3523       rc = rl < rc;
3524       rp[i] = rl;
3525     }
3526   assert (vc == 0);
3527
3528   for (; i < un; i++)
3529     {
3530       ul = (up[i] ^ ux) + uc;
3531       uc = ul < uc;
3532
3533       rl = (ul ^ ux) + rc;
3534       rc = rl < rc;
3535       rp[i] = rl;
3536     }
3537   if (rc)
3538     rp[un++] = rc;
3539   else
3540     un = mpn_normalized_size (rp, un);
3541
3542   r->_mp_size = rx ? -un : un;
3543 }
3544
3545 static unsigned
3546 gmp_popcount_limb (mp_limb_t x)
3547 {
3548   unsigned c;
3549
3550   /* Do 16 bits at a time, to avoid limb-sized constants. */
3551   for (c = 0; x > 0; x >>= 16)
3552     {
3553       unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555);
3554       w = ((w >> 2) & 0x3333) + (w & 0x3333);
3555       w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f);
3556       w = (w >> 8) + (w & 0x00ff);
3557       c += w;
3558     }
3559   return c;
3560 }
3561
3562 mp_bitcnt_t
3563 mpz_popcount (const mpz_t u)
3564 {
3565   mp_size_t un, i;
3566   mp_bitcnt_t c;
3567
3568   un = u->_mp_size;
3569
3570   if (un < 0)
3571     return ~(mp_bitcnt_t) 0;
3572
3573   for (c = 0, i = 0; i < un; i++)
3574     c += gmp_popcount_limb (u->_mp_d[i]);
3575
3576   return c;
3577 }
3578
3579 mp_bitcnt_t
3580 mpz_hamdist (const mpz_t u, const mpz_t v)
3581 {
3582   mp_size_t un, vn, i;
3583   mp_limb_t uc, vc, ul, vl, comp;
3584   mp_srcptr up, vp;
3585   mp_bitcnt_t c;
3586
3587   un = u->_mp_size;
3588   vn = v->_mp_size;
3589
3590   if ( (un ^ vn) < 0)
3591     return ~(mp_bitcnt_t) 0;
3592
3593   if (un < 0)
3594     {
3595       assert (vn < 0);
3596       un = -un;
3597       vn = -vn;
3598       uc = vc = 1;
3599       comp = - (mp_limb_t) 1;
3600     }
3601   else
3602     uc = vc = comp = 0;
3603
3604   up = u->_mp_d;
3605   vp = v->_mp_d;
3606
3607   if (un < vn)
3608     MPN_SRCPTR_SWAP (up, un, vp, vn);
3609
3610   for (i = 0, c = 0; i < vn; i++)
3611     {
3612       ul = (up[i] ^ comp) + uc;
3613       uc = ul < uc;
3614
3615       vl = (vp[i] ^ comp) + vc;
3616       vc = vl < vc;
3617
3618       c += gmp_popcount_limb (ul ^ vl);
3619     }
3620   assert (vc == 0);
3621
3622   for (; i < un; i++)
3623     {
3624       ul = (up[i] ^ comp) + uc;
3625       uc = ul < uc;
3626
3627       c += gmp_popcount_limb (ul ^ comp);
3628     }
3629
3630   return c;
3631 }
3632
3633 mp_bitcnt_t
3634 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
3635 {
3636   mp_ptr up;
3637   mp_size_t us, un, i;
3638   mp_limb_t limb, ux, uc;
3639   unsigned cnt;
3640
3641   up = u->_mp_d;
3642   us = u->_mp_size;
3643   un = GMP_ABS (us);
3644   i = starting_bit / GMP_LIMB_BITS;
3645
3646   /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3647      for u<0. Notice this test picks up any u==0 too. */
3648   if (i >= un)
3649     return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
3650
3651   if (us < 0)
3652     {
3653       ux = GMP_LIMB_MAX;
3654       uc = mpn_zero_p (up, i);
3655     }
3656   else
3657     ux = uc = 0;
3658
3659   limb = (ux ^ up[i]) + uc;
3660   uc = limb < uc;
3661
3662   /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3663   limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3664
3665   while (limb == 0)
3666     {
3667       i++;
3668       if (i == un)
3669         {
3670           assert (uc == 0);
3671           /* For the u > 0 case, this can happen only for the first
3672              masked limb. For the u < 0 case, it happens when the
3673              highest limbs of the absolute value are all ones. */
3674           return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
3675         }
3676       limb = (ux ^ up[i]) + uc;
3677       uc = limb < uc;
3678     }
3679   gmp_ctz (cnt, limb);
3680   return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
3681 }
3682
3683 mp_bitcnt_t
3684 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
3685 {
3686   mp_ptr up;
3687   mp_size_t us, un, i;
3688   mp_limb_t limb, ux, uc;
3689   unsigned cnt;
3690
3691   up = u->_mp_d;
3692   us = u->_mp_size;
3693   un = GMP_ABS (us);
3694   i = starting_bit / GMP_LIMB_BITS;
3695
3696   /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3697      u<0.  Notice this test picks up all cases of u==0 too. */
3698   if (i >= un)
3699     return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0);
3700
3701   if (us < 0)
3702     {
3703       ux = GMP_LIMB_MAX;
3704       uc = mpn_zero_p (up, i);
3705     }
3706   else
3707     ux = uc = 0;
3708
3709   limb = (ux ^ up[i]) + uc;
3710   uc = limb < uc;
3711
3712   /* Mask to 1 all bits before starting_bit, thus ignoring them. */
3713   limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1;
3714
3715   while (limb == GMP_LIMB_MAX)
3716     {
3717       i++;
3718       if (i == un)
3719         {
3720           assert (uc == 0);
3721           return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0);
3722         }
3723       limb = (ux ^ up[i]) + uc;
3724       uc = limb < uc;
3725     }
3726   gmp_ctz (cnt, ~limb);
3727   return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
3728 }
3729
3730 \f
3731 /* MPZ base conversion. */
3732
3733 size_t
3734 mpz_sizeinbase (const mpz_t u, int base)
3735 {
3736   mp_size_t un;
3737   mp_srcptr up;
3738   mp_ptr tp;
3739   mp_bitcnt_t bits;
3740   struct gmp_div_inverse bi;
3741   size_t ndigits;
3742
3743   assert (base >= 2);
3744   assert (base <= 36);
3745
3746   un = GMP_ABS (u->_mp_size);
3747   if (un == 0)
3748     return 1;
3749
3750   up = u->_mp_d;
3751
3752   bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
3753   switch (base)
3754     {
3755     case 2:
3756       return bits;
3757     case 4:
3758       return (bits + 1) / 2;
3759     case 8:
3760       return (bits + 2) / 3;
3761     case 16:
3762       return (bits + 3) / 4;
3763     case 32:
3764       return (bits + 4) / 5;
3765       /* FIXME: Do something more clever for the common case of base
3766          10. */
3767     }
3768
3769   tp = gmp_xalloc_limbs (un);
3770   mpn_copyi (tp, up, un);
3771   mpn_div_qr_1_invert (&bi, base);
3772
3773   for (ndigits = 0; un > 0; ndigits++)
3774     {
3775       mpn_div_qr_1_preinv (tp, tp, un, &bi);
3776       un -= (tp[un-1] == 0);
3777     }
3778   gmp_free (tp);
3779   return ndigits;
3780 }
3781
3782 char *
3783 mpz_get_str (char *sp, int base, const mpz_t u)
3784 {
3785   unsigned bits;
3786   const char *digits;
3787   mp_size_t un;
3788   size_t i, sn;
3789
3790   if (base >= 0)
3791     {
3792       digits = "0123456789abcdefghijklmnopqrstuvwxyz";
3793     }
3794   else
3795     {
3796       base = -base;
3797       digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
3798     }
3799   if (base <= 1)
3800     base = 10;
3801   if (base > 36)
3802     return NULL;
3803
3804   sn = 1 + mpz_sizeinbase (u, base);
3805   if (!sp)
3806     sp = gmp_xalloc (1 + sn);
3807
3808   un = GMP_ABS (u->_mp_size);
3809
3810   if (un == 0)
3811     {
3812       sp[0] = '0';
3813       sp[1] = '\0';
3814       return sp;
3815     }
3816
3817   i = 0;
3818
3819   if (u->_mp_size < 0)
3820     sp[i++] = '-';
3821
3822   bits = mpn_base_power_of_two_p (base);
3823
3824   if (bits)
3825     /* Not modified in this case. */
3826     sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
3827   else
3828     {
3829       struct mpn_base_info info;
3830       mp_ptr tp;
3831
3832       mpn_get_base_info (&info, base);
3833       tp = gmp_xalloc_limbs (un);
3834       mpn_copyi (tp, u->_mp_d, un);
3835
3836       sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
3837       gmp_free (tp);
3838     }
3839
3840   for (; i < sn; i++)
3841     sp[i] = digits[(unsigned char) sp[i]];
3842
3843   sp[sn] = '\0';
3844   return sp;
3845 }
3846
3847 int
3848 mpz_set_str (mpz_t r, const char *sp, int base)
3849 {
3850   unsigned bits;
3851   mp_size_t rn, alloc;
3852   mp_ptr rp;
3853   size_t sn;
3854   size_t dn;
3855   int sign;
3856   unsigned char *dp;
3857
3858   assert (base == 0 || (base >= 2 && base <= 36));
3859
3860   while (isspace( (unsigned char) *sp))
3861     sp++;
3862
3863   if (*sp == '-')
3864     {
3865       sign = 1;
3866       sp++;
3867     }
3868   else
3869     sign = 0;
3870
3871   if (base == 0)
3872     {
3873       if (*sp == '0')
3874         {
3875           sp++;
3876           if (*sp == 'x' || *sp == 'X')
3877             {
3878               base = 16;
3879               sp++;
3880             }
3881           else if (*sp == 'b' || *sp == 'B')
3882             {
3883               base = 2;
3884               sp++;
3885             }
3886           else
3887             base = 8;
3888         }
3889       else
3890         base = 10;
3891     }
3892
3893   sn = strlen (sp);
3894   dp = gmp_xalloc (sn + (sn == 0));
3895
3896   for (dn = 0; *sp; sp++)
3897     {
3898       unsigned digit;
3899
3900       if (isspace ((unsigned char) *sp))
3901         continue;
3902       if (*sp >= '0' && *sp <= '9')
3903         digit = *sp - '0';
3904       else if (*sp >= 'a' && *sp <= 'z')
3905         digit = *sp - 'a' + 10;
3906       else if (*sp >= 'A' && *sp <= 'Z')
3907         digit = *sp - 'A' + 10;
3908       else
3909         digit = base; /* fail */
3910
3911       if (digit >= base)
3912         {
3913           gmp_free (dp);
3914           r->_mp_size = 0;
3915           return -1;
3916         }
3917
3918       dp[dn++] = digit;
3919     }
3920
3921   bits = mpn_base_power_of_two_p (base);
3922
3923   if (bits > 0)
3924     {
3925       alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
3926       rp = MPZ_REALLOC (r, alloc);
3927       rn = mpn_set_str_bits (rp, dp, dn, bits);
3928     }
3929   else
3930     {
3931       struct mpn_base_info info;
3932       mpn_get_base_info (&info, base);
3933       alloc = (sn + info.exp - 1) / info.exp;
3934       rp = MPZ_REALLOC (r, alloc);
3935       rn = mpn_set_str_other (rp, dp, dn, base, &info);
3936     }
3937   assert (rn <= alloc);
3938   gmp_free (dp);
3939
3940   r->_mp_size = sign ? - rn : rn;
3941
3942   return 0;
3943 }
3944
3945 int
3946 mpz_init_set_str (mpz_t r, const char *sp, int base)
3947 {
3948   mpz_init (r);
3949   return mpz_set_str (r, sp, base);
3950 }
3951
3952 size_t
3953 mpz_out_str (FILE *stream, int base, const mpz_t x)
3954 {
3955   char *str;
3956   size_t len;
3957
3958   str = mpz_get_str (NULL, base, x);
3959   len = strlen (str);
3960   len = fwrite (str, 1, len, stream);
3961   gmp_free (str);
3962   return len;
3963 }
3964
3965 \f
3966 static int
3967 gmp_detect_endian (void)
3968 {
3969   static const int i = 1;
3970   const unsigned char *p = (const unsigned char *) &i;
3971   if (*p == 1)
3972     /* Little endian */
3973     return -1;
3974   else
3975     /* Big endian */
3976     return 1;
3977 }
3978
3979 /* Import and export. Does not support nails. */
3980 void
3981 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
3982             size_t nails, const void *src)
3983 {
3984   const unsigned char *p;
3985   ptrdiff_t word_step;
3986   mp_ptr rp;
3987   mp_size_t rn;
3988
3989   /* The current (partial) limb. */
3990   mp_limb_t limb;
3991   /* The number of bytes already copied to this limb (starting from
3992      the low end). */
3993   size_t bytes;
3994   /* The index where the limb should be stored, when completed. */
3995   mp_size_t i;
3996
3997   if (nails != 0)
3998     gmp_die ("mpz_import: Nails not supported.");
3999
4000   assert (order == 1 || order == -1);
4001   assert (endian >= -1 && endian <= 1);
4002
4003   if (endian == 0)
4004     endian = gmp_detect_endian ();
4005
4006   p = (unsigned char *) src;
4007
4008   word_step = (order != endian) ? 2 * size : 0;
4009
4010   /* Process bytes from the least significant end, so point p at the
4011      least significant word. */
4012   if (order == 1)
4013     {
4014       p += size * (count - 1);
4015       word_step = - word_step;
4016     }
4017
4018   /* And at least significant byte of that word. */
4019   if (endian == 1)
4020     p += (size - 1);
4021
4022   rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4023   rp = MPZ_REALLOC (r, rn);
4024
4025   for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4026     {
4027       size_t j;
4028       for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4029         {
4030           limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4031           if (bytes == sizeof(mp_limb_t))
4032             {
4033               rp[i++] = limb;
4034               bytes = 0;
4035               limb = 0;
4036             }
4037         }
4038     }
4039   if (bytes > 0)
4040     rp[i++] = limb;
4041   assert (i == rn);
4042
4043   r->_mp_size = mpn_normalized_size (rp, i);
4044 }
4045
4046 void *
4047 mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4048             size_t nails, const mpz_t u)
4049 {
4050   unsigned char *p;
4051   ptrdiff_t word_step;
4052   size_t count, k;
4053   mp_size_t un;
4054
4055   /* The current (partial) limb. */
4056   mp_limb_t limb;
4057   /* The number of bytes left to to in this limb. */
4058   size_t bytes;
4059   /* The index where the limb was read. */
4060   mp_size_t i;
4061
4062   if (nails != 0)
4063     gmp_die ("mpz_import: Nails not supported.");
4064
4065   assert (order == 1 || order == -1);
4066   assert (endian >= -1 && endian <= 1);
4067   assert (size > 0 || u->_mp_size == 0);
4068
4069   un = GMP_ABS (u->_mp_size);
4070   if (un == 0)
4071     {
4072       if (countp)
4073         *countp = 0;
4074       return r;
4075     }
4076
4077   /* Count bytes in top limb. */
4078   for (limb = u->_mp_d[un-1], k = 0; limb > 0; k++, limb >>= CHAR_BIT)
4079     ;
4080
4081   assert (k > 0);
4082
4083   count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4084
4085   if (!r)
4086     r = gmp_xalloc (count * size);
4087
4088   if (endian == 0)
4089     endian = gmp_detect_endian ();
4090
4091   p = (unsigned char *) r;
4092
4093   word_step = (order != endian) ? 2 * size : 0;
4094
4095   /* Process bytes from the least significant end, so point p at the
4096      least significant word. */
4097   if (order == 1)
4098     {
4099       p += size * (count - 1);
4100       word_step = - word_step;
4101     }
4102
4103   /* And at least significant byte of that word. */
4104   if (endian == 1)
4105     p += (size - 1);
4106
4107   for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4108       {
4109         size_t j;
4110         for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4111           {
4112             if (bytes == 0)
4113               {
4114                 if (i < un)
4115                   limb = u->_mp_d[i++];
4116                 bytes = sizeof (mp_limb_t);
4117               }
4118             *p = limb;
4119             limb >>= CHAR_BIT;
4120             bytes--;
4121           }
4122       }
4123   assert (i == un);
4124   assert (k == count);
4125
4126   if (countp)
4127     *countp = count;
4128
4129   return r;
4130 }