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   assert ((dp[dn-1] & GMP_LIMB_HIGHBIT) != 0);
 934
 935   d1 = dp[dn - 1];
 936   d0 = dp[dn - 2];
 937
 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 (dp[dn-1] & GMP_LIMB_HIGHBIT);
 998
 999       shift = inv->shift;
1000       if (shift > 0)
1001         nh = mpn_lshift (np, np, nn, shift);
1002       else
1003         nh = 0;
1004
1005       assert (inv->d1 == dp[dn-1]);
1006       assert (inv->d0 == dp[dn-2]);
1007
1008       mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
1009
1010       if (shift > 0)
1011         gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
1012     }
1013 }
1014
1015 static void
1016 mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
1017 {
1018   struct gmp_div_inverse inv;
1019   mp_ptr tp = NULL;
1020
1021   assert (dn > 0);
1022   assert (nn >= dn);
1023
1024   mpn_div_qr_invert (&inv, dp, dn);
1025   if (dn > 2 && inv.shift > 0)
1026     {
1027       tp = gmp_xalloc_limbs (dn);
1028       gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
1029       dp = tp;
1030     }
1031   mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
1032   if (tp)
1033     gmp_free (tp);
1034 }
1035
1036 \f
1037 /* MPN base conversion. */
1038 static unsigned
1039 mpn_base_power_of_two_p (unsigned b)
1040 {
1041   switch (b)
1042     {
1043     case 2: return 1;
1044     case 4: return 2;
1045     case 8: return 3;
1046     case 16: return 4;
1047     case 32: return 5;
1048     case 64: return 6;
1049     case 128: return 7;
1050     case 256: return 8;
1051     default: return 0;
1052     }
1053 }
1054
1055 struct mpn_base_info
1056 {
1057   /* bb is the largest power of the base which fits in one limb, and
1058      exp is the corresponding exponent. */
1059   unsigned exp;
1060   mp_limb_t bb;
1061 };
1062
1063 static void
1064 mpn_get_base_info (struct mpn_base_info *info, mp_limb_t b)
1065 {
1066   mp_limb_t m;
1067   mp_limb_t p;
1068   unsigned exp;
1069
1070   m = GMP_LIMB_MAX / b;
1071   for (exp = 1, p = b; p <= m; exp++)
1072     p *= b;
1073
1074   info->exp = exp;
1075   info->bb = p;
1076 }
1077
1078 static mp_bitcnt_t
1079 mpn_limb_size_in_base_2 (mp_limb_t u)
1080 {
1081   unsigned shift;
1082
1083   assert (u > 0);
1084   gmp_clz (shift, u);
1085   return GMP_LIMB_BITS - shift;
1086 }
1087
1088 static size_t
1089 mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1090 {
1091   unsigned char mask;
1092   size_t sn, j;
1093   mp_size_t i;
1094   int shift;
1095
1096   sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1097         + bits - 1) / bits;
1098
1099   mask = (1U << bits) - 1;
1100
1101   for (i = 0, j = sn, shift = 0; j-- > 0;)
1102     {
1103       unsigned char digit = up[i] >> shift;
1104
1105       shift += bits;
1106
1107       if (shift >= GMP_LIMB_BITS && ++i < un)
1108         {
1109           shift -= GMP_LIMB_BITS;
1110           digit |= up[i] << (bits - shift);
1111         }
1112       sp[j] = digit & mask;
1113     }
1114   return sn;
1115 }
1116
1117 /* We generate digits from the least significant end, and reverse at
1118    the end. */
1119 static size_t
1120 mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
1121                   const struct gmp_div_inverse *binv)
1122 {
1123   mp_size_t i;
1124   for (i = 0; w > 0; i++)
1125     {
1126       mp_limb_t h, l, r;
1127
1128       h = w >> (GMP_LIMB_BITS - binv->shift);
1129       l = w << binv->shift;
1130
1131       gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
1132       assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
1133       r >>= binv->shift;
1134
1135       sp[i] = r;
1136     }
1137   return i;
1138 }
1139
1140 static size_t
1141 mpn_get_str_other (unsigned char *sp,
1142                    int base, const struct mpn_base_info *info,
1143                    mp_ptr up, mp_size_t un)
1144 {
1145   struct gmp_div_inverse binv;
1146   size_t sn;
1147   size_t i;
1148
1149   mpn_div_qr_1_invert (&binv, base);
1150
1151   sn = 0;
1152
1153   if (un > 1)
1154     {
1155       struct gmp_div_inverse bbinv;
1156       mpn_div_qr_1_invert (&bbinv, info->bb);
1157
1158       do
1159         {
1160           mp_limb_t w;
1161           size_t done;
1162           w = mpn_div_qr_1_preinv (up, up, un, &bbinv);
1163           un -= (up[un-1] == 0);
1164           done = mpn_limb_get_str (sp + sn, w, &binv);
1165
1166           for (sn += done; done < info->exp; done++)
1167             sp[sn++] = 0;
1168         }
1169       while (un > 1);
1170     }
1171   sn += mpn_limb_get_str (sp + sn, up[0], &binv);
1172
1173   /* Reverse order */
1174   for (i = 0; 2*i + 1 < sn; i++)
1175     {
1176       unsigned char t = sp[i];
1177       sp[i] = sp[sn - i - 1];
1178       sp[sn - i - 1] = t;
1179     }
1180
1181   return sn;
1182 }
1183
1184 size_t
1185 mpn_get_str (unsigned char *sp, int base, mp_ptr up, mp_size_t un)
1186 {
1187   unsigned bits;
1188
1189   assert (un > 0);
1190   assert (up[un-1] > 0);
1191
1192   bits = mpn_base_power_of_two_p (base);
1193   if (bits)
1194     return mpn_get_str_bits (sp, bits, up, un);
1195   else
1196     {
1197       struct mpn_base_info info;
1198
1199       mpn_get_base_info (&info, base);
1200       return mpn_get_str_other (sp, base, &info, up, un);
1201     }
1202 }
1203
1204 static mp_size_t
1205 mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1206                   unsigned bits)
1207 {
1208   mp_size_t rn;
1209   size_t j;
1210   unsigned shift;
1211
1212   for (j = sn, rn = 0, shift = 0; j-- > 0; )
1213     {
1214       if (shift == 0)
1215         {
1216           rp[rn++] = sp[j];
1217           shift += bits;
1218         }
1219       else
1220         {
1221           rp[rn-1] |= (mp_limb_t) sp[j] << shift;
1222           shift += bits;
1223           if (shift >= GMP_LIMB_BITS)
1224             {
1225               shift -= GMP_LIMB_BITS;
1226               if (shift > 0)
1227                 rp[rn++] = (mp_limb_t) sp[j] >> (bits - shift);
1228             }
1229         }
1230     }
1231   rn = mpn_normalized_size (rp, rn);
1232   return rn;
1233 }
1234
1235 static mp_size_t
1236 mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1237                    mp_limb_t b, const struct mpn_base_info *info)
1238 {
1239   mp_size_t rn;
1240   mp_limb_t w;
1241   unsigned first;
1242   unsigned k;
1243   size_t j;
1244
1245   first = 1 + (sn - 1) % info->exp;
1246
1247   j = 0;
1248   w = sp[j++];
1249   for (k = 1; k < first; k++)
1250     w = w * b + sp[j++];
1251
1252   rp[0] = w;
1253
1254   for (rn = (w > 0); j < sn;)
1255     {
1256       mp_limb_t cy;
1257
1258       w = sp[j++];
1259       for (k = 1; k < info->exp; k++)
1260         w = w * b + sp[j++];
1261
1262       cy = mpn_mul_1 (rp, rp, rn, info->bb);
1263       cy += mpn_add_1 (rp, rp, rn, w);
1264       if (cy > 0)
1265         rp[rn++] = cy;
1266     }
1267   assert (j == sn);
1268
1269   return rn;
1270 }
1271
1272 mp_size_t
1273 mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1274 {
1275   unsigned bits;
1276
1277   if (sn == 0)
1278     return 0;
1279
1280   bits = mpn_base_power_of_two_p (base);
1281   if (bits)
1282     return mpn_set_str_bits (rp, sp, sn, bits);
1283   else
1284     {
1285       struct mpn_base_info info;
1286
1287       mpn_get_base_info (&info, base);
1288       return mpn_set_str_other (rp, sp, sn, base, &info);
1289     }
1290 }
1291
1292 \f
1293 /* MPZ interface */
1294 void
1295 mpz_init (mpz_t r)
1296 {
1297   r->_mp_alloc = 1;
1298   r->_mp_size = 0;
1299   r->_mp_d = gmp_xalloc_limbs (1);
1300 }
1301
1302 /* The utility of this function is a bit limited, since many functions
1303    assings the result variable using mpz_swap. */
1304 void
1305 mpz_init2 (mpz_t r, mp_bitcnt_t bits)
1306 {
1307   mp_size_t rn;
1308
1309   bits -= (bits != 0);          /* Round down, except if 0 */
1310   rn = 1 + bits / GMP_LIMB_BITS;
1311
1312   r->_mp_alloc = rn;
1313   r->_mp_size = 0;
1314   r->_mp_d = gmp_xalloc_limbs (rn);
1315 }
1316
1317 void
1318 mpz_clear (mpz_t r)
1319 {
1320   gmp_free (r->_mp_d);
1321 }
1322
1323 static void *
1324 mpz_realloc (mpz_t r, mp_size_t size)
1325 {
1326   size = GMP_MAX (size, 1);
1327
1328   r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1329   r->_mp_alloc = size;
1330
1331   if (GMP_ABS (r->_mp_size) > size)
1332     r->_mp_size = 0;
1333
1334   return r->_mp_d;
1335 }
1336
1337 /* Realloc for an mpz_t WHAT if it has less than NEEDED limbs.  */
1338 #define MPZ_REALLOC(z,n) ((n) > (z)->_mp_alloc                  \
1339                           ? mpz_realloc(z,n)                    \
1340                           : (z)->_mp_d)
1341 \f
1342 /* MPZ assignment and basic conversions. */
1343 void
1344 mpz_set_si (mpz_t r, signed long int x)
1345 {
1346   if (x >= 0)
1347     mpz_set_ui (r, x);
1348   else /* (x < 0) */
1349     {
1350       r->_mp_size = -1;
1351       r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x);
1352     }
1353 }
1354
1355 void
1356 mpz_set_ui (mpz_t r, unsigned long int x)
1357 {
1358   if (x > 0)
1359     {
1360       r->_mp_size = 1;
1361       r->_mp_d[0] = x;
1362     }
1363   else
1364     r->_mp_size = 0;
1365 }
1366
1367 void
1368 mpz_set (mpz_t r, const mpz_t x)
1369 {
1370   /* Allow the NOP r == x */
1371   if (r != x)
1372     {
1373       mp_size_t n;
1374       mp_ptr rp;
1375
1376       n = GMP_ABS (x->_mp_size);
1377       rp = MPZ_REALLOC (r, n);
1378
1379       mpn_copyi (rp, x->_mp_d, n);
1380       r->_mp_size = x->_mp_size;
1381     }
1382 }
1383
1384 void
1385 mpz_init_set_si (mpz_t r, signed long int x)
1386 {
1387   mpz_init (r);
1388   mpz_set_si (r, x);
1389 }
1390
1391 void
1392 mpz_init_set_ui (mpz_t r, unsigned long int x)
1393 {
1394   mpz_init (r);
1395   mpz_set_ui (r, x);
1396 }
1397
1398 void
1399 mpz_init_set (mpz_t r, const mpz_t x)
1400 {
1401   mpz_init (r);
1402   mpz_set (r, x);
1403 }
1404
1405 int
1406 mpz_fits_slong_p (const mpz_t u)
1407 {
1408   mp_size_t us = u->_mp_size;
1409
1410   if (us == 0)
1411     return 1;
1412   else if (us == 1)
1413     return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1414   else if (us == -1)
1415     return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1416   else
1417     return 0;
1418 }
1419
1420 int
1421 mpz_fits_ulong_p (const mpz_t u)
1422 {
1423   mp_size_t us = u->_mp_size;
1424
1425   return us == 0 || us == 1;
1426 }
1427
1428 long int
1429 mpz_get_si (const mpz_t u)
1430 {
1431   mp_size_t us = u->_mp_size;
1432
1433   if (us > 0)
1434     return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT);
1435   else if (us < 0)
1436     return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT);
1437   else
1438     return 0;
1439 }
1440
1441 unsigned long int
1442 mpz_get_ui (const mpz_t u)
1443 {
1444   return u->_mp_size == 0 ? 0 : u->_mp_d[0];
1445 }
1446
1447 size_t
1448 mpz_size (const mpz_t u)
1449 {
1450   return GMP_ABS (u->_mp_size);
1451 }
1452
1453 mp_limb_t
1454 mpz_getlimbn (const mpz_t u, mp_size_t n)
1455 {
1456   if (n >= 0 && n < GMP_ABS (u->_mp_size))
1457     return u->_mp_d[n];
1458   else
1459     return 0;
1460 }
1461
1462 \f
1463 /* Conversions and comparison to double. */
1464 void
1465 mpz_set_d (mpz_t r, double x)
1466 {
1467   int sign;
1468   mp_ptr rp;
1469   mp_size_t rn, i;
1470   double B;
1471   double Bi;
1472   mp_limb_t f;
1473
1474   /* x != x is true when x is a NaN, and x == x * 0.5 is true when x is
1475      zero or infinity. */
1476   if (x == 0.0 || x != x || x == x * 0.5)
1477     {
1478       r->_mp_size = 0;
1479       return;
1480     }
1481
1482   if (x < 0.0)
1483     {
1484       x = - x;
1485       sign = 1;
1486     }
1487   else
1488     sign = 0;
1489
1490   if (x < 1.0)
1491     {
1492       r->_mp_size = 0;
1493       return;
1494     }
1495   B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1496   Bi = 1.0 / B;
1497   for (rn = 1; x >= B; rn++)
1498     x *= Bi;
1499
1500   rp = MPZ_REALLOC (r, rn);
1501
1502   f = (mp_limb_t) x;
1503   x -= f;
1504   assert (x < 1.0);
1505   rp[rn-1] = f;
1506   for (i = rn-1; i-- > 0; )
1507     {
1508       x = B * x;
1509       f = (mp_limb_t) x;
1510       x -= f;
1511       assert (x < 1.0);
1512       rp[i] = f;
1513     }
1514
1515   r->_mp_size = sign ? - rn : rn;
1516 }
1517
1518 void
1519 mpz_init_set_d (mpz_t r, double x)
1520 {
1521   mpz_init (r);
1522   mpz_set_d (r, x);
1523 }
1524
1525 double
1526 mpz_get_d (const mpz_t u)
1527 {
1528   mp_size_t un;
1529   double x;
1530   double B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1531
1532   un = GMP_ABS (u->_mp_size);
1533
1534   if (un == 0)
1535     return 0.0;
1536
1537   x = u->_mp_d[--un];
1538   while (un > 0)
1539     x = B*x + u->_mp_d[--un];
1540
1541   if (u->_mp_size < 0)
1542     x = -x;
1543
1544   return x;
1545 }
1546
1547 int
1548 mpz_cmpabs_d (const mpz_t x, double d)
1549 {
1550   mp_size_t xn;
1551   double B, Bi;
1552   mp_size_t i;
1553
1554   xn = x->_mp_size;
1555   d = GMP_ABS (d);
1556
1557   if (xn != 0)
1558     {
1559       xn = GMP_ABS (xn);
1560
1561       B = 2.0 * (double) GMP_LIMB_HIGHBIT;
1562       Bi = 1.0 / B;
1563
1564       /* Scale d so it can be compared with the top limb. */
1565       for (i = 1; i < xn; i++)
1566         d *= Bi;
1567
1568       if (d >= B)
1569         return -1;
1570
1571       /* Compare floor(d) to top limb, subtract and cancel when equal. */
1572       for (i = xn; i-- > 0;)
1573         {
1574           mp_limb_t f, xl;
1575
1576           f = (mp_limb_t) d;
1577           xl = x->_mp_d[i];
1578           if (xl > f)
1579             return 1;
1580           else if (xl < f)
1581             return -1;
1582           d = B * (d - f);
1583         }
1584     }
1585   return - (d > 0.0);
1586 }
1587
1588 int
1589 mpz_cmp_d (const mpz_t x, double d)
1590 {
1591   if (x->_mp_size < 0)
1592     {
1593       if (d >= 0.0)
1594         return -1;
1595       else
1596         return -mpz_cmpabs_d (x, d);
1597     }
1598   else
1599     {
1600       if (d < 0.0)
1601         return 1;
1602       else
1603         return mpz_cmpabs_d (x, d);
1604     }
1605 }
1606
1607 \f
1608 /* MPZ comparisons and the like. */
1609 int
1610 mpz_sgn (const mpz_t u)
1611 {
1612   mp_size_t usize = u->_mp_size;
1613
1614   if (usize > 0)
1615     return 1;
1616   else if (usize < 0)
1617     return -1;
1618   else
1619     return 0;
1620 }
1621
1622 int
1623 mpz_cmp_si (const mpz_t u, long v)
1624 {
1625   mp_size_t usize = u->_mp_size;
1626
1627   if (usize < -1)
1628     return -1;
1629   else if (v >= 0)
1630     return mpz_cmp_ui (u, v);
1631   else if (usize >= 0)
1632     return 1;
1633   else /* usize == -1 */
1634     {
1635       mp_limb_t ul = u->_mp_d[0];
1636       if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
1637         return -1;
1638       else if ( (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul)
1639         return 1;
1640     }
1641   return 0;
1642 }
1643
1644 int
1645 mpz_cmp_ui (const mpz_t u, unsigned long v)
1646 {
1647   mp_size_t usize = u->_mp_size;
1648
1649   if (usize > 1)
1650     return 1;
1651   else if (usize < 0)
1652     return -1;
1653   else
1654     {
1655       mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
1656       if (ul > v)
1657         return 1;
1658       else if (ul < v)
1659         return -1;
1660     }
1661   return 0;
1662 }
1663
1664 int
1665 mpz_cmp (const mpz_t a, const mpz_t b)
1666 {
1667   mp_size_t asize = a->_mp_size;
1668   mp_size_t bsize = b->_mp_size;
1669
1670   if (asize > bsize)
1671     return 1;
1672   else if (asize < bsize)
1673     return -1;
1674   else if (asize > 0)
1675     return mpn_cmp (a->_mp_d, b->_mp_d, asize);
1676   else if (asize < 0)
1677     return -mpn_cmp (a->_mp_d, b->_mp_d, -asize);
1678   else
1679     return 0;
1680 }
1681
1682 int
1683 mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1684 {
1685   mp_size_t un = GMP_ABS (u->_mp_size);
1686   mp_limb_t ul;
1687
1688   if (un > 1)
1689     return 1;
1690
1691   ul = (un == 1) ? u->_mp_d[0] : 0;
1692
1693   if (ul > v)
1694     return 1;
1695   else if (ul < v)
1696     return -1;
1697   else
1698     return 0;
1699 }
1700
1701 int
1702 mpz_cmpabs (const mpz_t u, const mpz_t v)
1703 {
1704   return mpn_cmp4 (u->_mp_d, GMP_ABS (u->_mp_size),
1705                    v->_mp_d, GMP_ABS (v->_mp_size));
1706 }
1707
1708 void
1709 mpz_abs (mpz_t r, const mpz_t u)
1710 {
1711   if (r != u)
1712     mpz_set (r, u);
1713
1714   r->_mp_size = GMP_ABS (r->_mp_size);
1715 }
1716
1717 void
1718 mpz_neg (mpz_t r, const mpz_t u)
1719 {
1720   if (r != u)
1721     mpz_set (r, u);
1722
1723   r->_mp_size = -r->_mp_size;
1724 }
1725
1726 void
1727 mpz_swap (mpz_t u, mpz_t v)
1728 {
1729   MP_SIZE_T_SWAP (u->_mp_size, v->_mp_size);
1730   MP_SIZE_T_SWAP (u->_mp_alloc, v->_mp_alloc);
1731   MP_PTR_SWAP (u->_mp_d, v->_mp_d);
1732 }
1733
1734 \f
1735 /* MPZ addition and subtraction */
1736
1737 /* Adds to the absolute value. Returns new size, but doesn't store it. */
1738 static mp_size_t
1739 mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1740 {
1741   mp_size_t an;
1742   mp_ptr rp;
1743   mp_limb_t cy;
1744
1745   an = GMP_ABS (a->_mp_size);
1746   if (an == 0)
1747     {
1748       r->_mp_d[0] = b;
1749       return b > 0;
1750     }
1751
1752   rp = MPZ_REALLOC (r, an + 1);
1753
1754   cy = mpn_add_1 (rp, a->_mp_d, an, b);
1755   rp[an] = cy;
1756   an += (cy > 0);
1757
1758   return an;
1759 }
1760
1761 /* Subtract from the absolute value. Returns new size, (or -1 on underflow),
1762    but doesn't store it. */
1763 static mp_size_t
1764 mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1765 {
1766   mp_size_t an = GMP_ABS (a->_mp_size);
1767   mp_ptr rp = MPZ_REALLOC (r, an);
1768
1769   if (an == 0)
1770     {
1771       rp[0] = b;
1772       return -(b > 0);
1773     }
1774   else if (an == 1 && a->_mp_d[0] < b)
1775     {
1776       rp[0] = b - a->_mp_d[0];
1777       return -1;
1778     }
1779   else
1780     {
1781       gmp_assert_nocarry (mpn_sub_1 (rp, a->_mp_d, an, b));
1782       return mpn_normalized_size (rp, an);
1783     }
1784 }
1785
1786 void
1787 mpz_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1788 {
1789   if (a->_mp_size >= 0)
1790     r->_mp_size = mpz_abs_add_ui (r, a, b);
1791   else
1792     r->_mp_size = -mpz_abs_sub_ui (r, a, b);
1793 }
1794
1795 void
1796 mpz_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1797 {
1798   if (a->_mp_size < 0)
1799     r->_mp_size = -mpz_abs_add_ui (r, a, b);
1800   else
1801     r->_mp_size = mpz_abs_sub_ui (r, a, b);
1802 }
1803
1804 void
1805 mpz_ui_sub (mpz_t r, unsigned long a, const mpz_t b)
1806 {
1807   if (b->_mp_size < 0)
1808     r->_mp_size = mpz_abs_add_ui (r, b, a);
1809   else
1810     r->_mp_size = -mpz_abs_sub_ui (r, b, a);
1811 }
1812
1813 static mp_size_t
1814 mpz_abs_add (mpz_t r, const mpz_t a, const mpz_t b)
1815 {
1816   mp_size_t an = GMP_ABS (a->_mp_size);
1817   mp_size_t bn = GMP_ABS (b->_mp_size);
1818   mp_size_t rn;
1819   mp_ptr rp;
1820   mp_limb_t cy;
1821
1822   rn = GMP_MAX (an, bn);
1823   rp = MPZ_REALLOC (r, rn + 1);
1824   if (an >= bn)
1825     cy = mpn_add (rp, a->_mp_d, an, b->_mp_d, bn);
1826   else
1827     cy = mpn_add (rp, b->_mp_d, bn, a->_mp_d, an);
1828
1829   rp[rn] = cy;
1830
1831   return rn + (cy > 0);
1832 }
1833
1834 static mp_size_t
1835 mpz_abs_sub (mpz_t r, const mpz_t a, const mpz_t b)
1836 {
1837   mp_size_t an = GMP_ABS (a->_mp_size);
1838   mp_size_t bn = GMP_ABS (b->_mp_size);
1839   int cmp;
1840   mp_ptr rp;
1841
1842   cmp = mpn_cmp4 (a->_mp_d, an, b->_mp_d, bn);
1843   if (cmp > 0)
1844     {
1845       rp = MPZ_REALLOC (r, an);
1846       gmp_assert_nocarry (mpn_sub (rp, a->_mp_d, an, b->_mp_d, bn));
1847       return mpn_normalized_size (rp, an);
1848     }
1849   else if (cmp < 0)
1850     {
1851       rp = MPZ_REALLOC (r, bn);
1852       gmp_assert_nocarry (mpn_sub (rp, b->_mp_d, bn, a->_mp_d, an));
1853       return -mpn_normalized_size (rp, bn);
1854     }
1855   else
1856     return 0;
1857 }
1858
1859 void
1860 mpz_add (mpz_t r, const mpz_t a, const mpz_t b)
1861 {
1862   mp_size_t rn;
1863
1864   if ( (a->_mp_size ^ b->_mp_size) >= 0)
1865     rn = mpz_abs_add (r, a, b);
1866   else
1867     rn = mpz_abs_sub (r, a, b);
1868
1869   r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1870 }
1871
1872 void
1873 mpz_sub (mpz_t r, const mpz_t a, const mpz_t b)
1874 {
1875   mp_size_t rn;
1876
1877   if ( (a->_mp_size ^ b->_mp_size) >= 0)
1878     rn = mpz_abs_sub (r, a, b);
1879   else
1880     rn = mpz_abs_add (r, a, b);
1881
1882   r->_mp_size = a->_mp_size >= 0 ? rn : - rn;
1883 }
1884
1885 \f
1886 /* MPZ multiplication */
1887 void
1888 mpz_mul_si (mpz_t r, const mpz_t u, long int v)
1889 {
1890   if (v < 0)
1891     {
1892       mpz_mul_ui (r, u, GMP_NEG_CAST (unsigned long int, v));
1893       mpz_neg (r, r);
1894     }
1895   else
1896     mpz_mul_ui (r, u, (unsigned long int) v);
1897 }
1898
1899 void
1900 mpz_mul_ui (mpz_t r, const mpz_t u, unsigned long int v)
1901 {
1902   mp_size_t un;
1903   mpz_t t;
1904   mp_ptr tp;
1905   mp_limb_t cy;
1906
1907   un = GMP_ABS (u->_mp_size);
1908
1909   if (un == 0 || v == 0)
1910     {
1911       r->_mp_size = 0;
1912       return;
1913     }
1914
1915   mpz_init2 (t, (un + 1) * GMP_LIMB_BITS);
1916
1917   tp = t->_mp_d;
1918   cy = mpn_mul_1 (tp, u->_mp_d, un, v);
1919   tp[un] = cy;
1920
1921   t->_mp_size = un + (cy > 0);
1922   if (u->_mp_size < 0)
1923     t->_mp_size = - t->_mp_size;
1924
1925   mpz_swap (r, t);
1926   mpz_clear (t);
1927 }
1928
1929 void
1930 mpz_mul (mpz_t r, const mpz_t u, const mpz_t v)
1931 {
1932   int sign;
1933   mp_size_t un, vn, rn;
1934   mpz_t t;
1935   mp_ptr tp;
1936
1937   un = GMP_ABS (u->_mp_size);
1938   vn = GMP_ABS (v->_mp_size);
1939
1940   if (un == 0 || vn == 0)
1941     {
1942       r->_mp_size = 0;
1943       return;
1944     }
1945
1946   sign = (u->_mp_size ^ v->_mp_size) < 0;
1947
1948   mpz_init2 (t, (un + vn) * GMP_LIMB_BITS);
1949
1950   tp = t->_mp_d;
1951   if (un >= vn)
1952     mpn_mul (tp, u->_mp_d, un, v->_mp_d, vn);
1953   else
1954     mpn_mul (tp, v->_mp_d, vn, u->_mp_d, un);
1955
1956   rn = un + vn;
1957   rn -= tp[rn-1] == 0;
1958
1959   t->_mp_size = sign ? - rn : rn;
1960   mpz_swap (r, t);
1961   mpz_clear (t);
1962 }
1963
1964 void
1965 mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
1966 {
1967   mp_size_t un, rn;
1968   mp_size_t limbs;
1969   unsigned shift;
1970   mp_ptr rp;
1971
1972   un = GMP_ABS (u->_mp_size);
1973   if (un == 0)
1974     {
1975       r->_mp_size = 0;
1976       return;
1977     }
1978
1979   limbs = bits / GMP_LIMB_BITS;
1980   shift = bits % GMP_LIMB_BITS;
1981
1982   rn = un + limbs + (shift > 0);
1983   rp = MPZ_REALLOC (r, rn);
1984   if (shift > 0)
1985     {
1986       mp_limb_t cy = mpn_lshift (rp + limbs, u->_mp_d, un, shift);
1987       rp[rn-1] = cy;
1988       rn -= (cy == 0);
1989     }
1990   else
1991     mpn_copyd (rp + limbs, u->_mp_d, un);
1992
1993   while (limbs > 0)
1994     rp[--limbs] = 0;
1995
1996   r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
1997 }
1998
1999 \f
2000 /* MPZ division */
2001 enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2002
2003 /* Allows q or r to be zero. Returns 1 iff remainder is non-zero. */
2004 static int
2005 mpz_div_qr (mpz_t q, mpz_t r,
2006             const mpz_t n, const mpz_t d, enum mpz_div_round_mode mode)
2007 {
2008   mp_size_t ns, ds, nn, dn, qs;
2009   ns = n->_mp_size;
2010   ds = d->_mp_size;
2011
2012   if (ds == 0)
2013     gmp_die("mpz_div_qr: Divide by zero.");
2014
2015   if (ns == 0)
2016     {
2017       if (q)
2018         q->_mp_size = 0;
2019       if (r)
2020         r->_mp_size = 0;
2021       return 0;
2022     }
2023
2024   nn = GMP_ABS (ns);
2025   dn = GMP_ABS (ds);
2026
2027   qs = ds ^ ns;
2028
2029   if (nn < dn)
2030     {
2031       if (mode == GMP_DIV_CEIL && qs >= 0)
2032         {
2033           /* q = 1, r = n - d */
2034           if (r)
2035             mpz_sub (r, n, d);
2036           if (q)
2037             mpz_set_ui (q, 1);
2038         }
2039       else if (mode == GMP_DIV_FLOOR && qs < 0)
2040         {
2041           /* q = -1, r = n + d */
2042           if (r)
2043             mpz_add (r, n, d);
2044           if (q)
2045             mpz_set_si (q, -1);
2046         }
2047       else
2048         {
2049           /* q = 0, r = d */
2050           if (r)
2051             mpz_set (r, n);
2052           if (q)
2053             q->_mp_size = 0;
2054         }
2055       return 1;
2056     }
2057   else
2058     {
2059       mp_ptr np, qp;
2060       mp_size_t qn, rn;
2061       mpz_t tq, tr;
2062
2063       mpz_init (tr);
2064       mpz_set (tr, n);
2065       np = tr->_mp_d;
2066
2067       qn = nn - dn + 1;
2068
2069       if (q)
2070         {
2071           mpz_init2 (tq, qn * GMP_LIMB_BITS);
2072           qp = tq->_mp_d;
2073         }
2074       else
2075         qp = NULL;
2076
2077       mpn_div_qr (qp, np, nn, d->_mp_d, dn);
2078
2079       if (qp)
2080         {
2081           qn -= (qp[qn-1] == 0);
2082
2083           tq->_mp_size = qs < 0 ? -qn : qn;
2084         }
2085       rn = mpn_normalized_size (np, dn);
2086       tr->_mp_size = ns < 0 ? - rn : rn;
2087
2088       if (mode == GMP_DIV_FLOOR && qs < 0 && rn != 0)
2089         {
2090           if (q)
2091             mpz_sub_ui (tq, tq, 1);
2092           if (r)
2093             mpz_add (tr, tr, d);
2094         }
2095       else if (mode == GMP_DIV_CEIL && qs >= 0 && rn != 0)
2096         {
2097           if (q)
2098             mpz_add_ui (tq, tq, 1);
2099           if (r)
2100             mpz_sub (tr, tr, d);
2101         }
2102
2103       if (q)
2104         {
2105           mpz_swap (tq, q);
2106           mpz_clear (tq);
2107         }
2108       if (r)
2109         mpz_swap (tr, r);
2110
2111       mpz_clear (tr);
2112
2113       return rn != 0;
2114     }
2115 }
2116
2117 void
2118 mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2119 {
2120   mpz_div_qr (q, r, n, d, GMP_DIV_CEIL);
2121 }
2122
2123 void
2124 mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2125 {
2126   mpz_div_qr (q, r, n, d, GMP_DIV_FLOOR);
2127 }
2128
2129 void
2130 mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
2131 {
2132   mpz_div_qr (q, r, n, d, GMP_DIV_TRUNC);
2133 }
2134
2135 void
2136 mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2137 {
2138   mpz_div_qr (q, NULL, n, d, GMP_DIV_CEIL);
2139 }
2140
2141 void
2142 mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2143 {
2144   mpz_div_qr (q, NULL, n, d, GMP_DIV_FLOOR);
2145 }
2146
2147 void
2148 mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
2149 {
2150   mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC);
2151 }
2152
2153 void
2154 mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2155 {
2156   mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2157 }
2158
2159 void
2160 mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2161 {
2162   mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2163 }
2164
2165 void
2166 mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
2167 {
2168   mpz_div_qr (NULL, r, n, d, GMP_DIV_TRUNC);
2169 }
2170
2171 void
2172 mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
2173 {
2174   if (d->_mp_size >= 0)
2175     mpz_div_qr (NULL, r, n, d, GMP_DIV_FLOOR);
2176   else
2177     mpz_div_qr (NULL, r, n, d, GMP_DIV_CEIL);
2178 }
2179
2180 static void
2181 mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2182                 enum mpz_div_round_mode mode)
2183 {
2184   mp_size_t un, qn;
2185   mp_size_t limb_cnt;
2186   mp_ptr qp;
2187   mp_limb_t adjust;
2188
2189   un = u->_mp_size;
2190   if (un == 0)
2191     {
2192       q->_mp_size = 0;
2193       return;
2194     }
2195   limb_cnt = bit_index / GMP_LIMB_BITS;
2196   qn = GMP_ABS (un) - limb_cnt;
2197   bit_index %= GMP_LIMB_BITS;
2198
2199   if (mode == ((un > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* un != 0 here. */
2200     /* Note: Below, the final indexing at limb_cnt is valid because at
2201        that point we have qn > 0. */
2202     adjust = (qn <= 0
2203               || !mpn_zero_p (u->_mp_d, limb_cnt)
2204               || (u->_mp_d[limb_cnt]
2205                   & (((mp_limb_t) 1 << bit_index) - 1)));
2206   else
2207     adjust = 0;
2208
2209   if (qn <= 0)
2210     qn = 0;
2211
2212   else
2213     {
2214       qp = MPZ_REALLOC (q, qn);
2215
2216       if (bit_index != 0)
2217         {
2218           mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
2219           qn -= qp[qn - 1] == 0;
2220         }
2221       else
2222         {
2223           mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
2224         }
2225     }
2226
2227   q->_mp_size = qn;
2228
2229   mpz_add_ui (q, q, adjust);
2230   if (un < 0)
2231     mpz_neg (q, q);
2232 }
2233
2234 static void
2235 mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2236                 enum mpz_div_round_mode mode)
2237 {
2238   mp_size_t us, un, rn;
2239   mp_ptr rp;
2240   mp_limb_t mask;
2241
2242   us = u->_mp_size;
2243   if (us == 0 || bit_index == 0)
2244     {
2245       r->_mp_size = 0;
2246       return;
2247     }
2248   rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
2249   assert (rn > 0);
2250
2251   rp = MPZ_REALLOC (r, rn);
2252   un = GMP_ABS (us);
2253
2254   mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
2255
2256   if (rn > un)
2257     {
2258       /* Quotient (with truncation) is zero, and remainder is
2259          non-zero */
2260       if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2261         {
2262           /* Have to negate and sign extend. */
2263           mp_size_t i;
2264           mp_limb_t cy;
2265
2266           for (cy = 1, i = 0; i < un; i++)
2267             {
2268               mp_limb_t s = ~u->_mp_d[i] + cy;
2269               cy = s < cy;
2270               rp[i] = s;
2271             }
2272           assert (cy == 0);
2273           for (; i < rn - 1; i++)
2274             rp[i] = GMP_LIMB_MAX;
2275
2276           rp[rn-1] = mask;
2277           us = -us;
2278         }
2279       else
2280         {
2281           /* Just copy */
2282           if (r != u)
2283             mpn_copyi (rp, u->_mp_d, un);
2284
2285           rn = un;
2286         }
2287     }
2288   else
2289     {
2290       if (r != u)
2291         mpn_copyi (rp, u->_mp_d, rn - 1);
2292
2293       rp[rn-1] = u->_mp_d[rn-1] & mask;
2294
2295       if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2296         {
2297           /* If r != 0, compute 2^{bit_count} - r. */
2298           mp_size_t i;
2299
2300           for (i = 0; i < rn && rp[i] == 0; i++)
2301             ;
2302           if (i < rn)
2303             {
2304               /* r > 0, need to flip sign. */
2305               rp[i] = ~rp[i] + 1;
2306               for (i++; i < rn; i++)
2307                 rp[i] = ~rp[i];
2308
2309               rp[rn-1] &= mask;
2310
2311               /* us is not used for anything else, so we can modify it
2312                  here to indicate flipped sign. */
2313               us = -us;
2314             }
2315         }
2316     }
2317   rn = mpn_normalized_size (rp, rn);
2318   r->_mp_size = us < 0 ? -rn : rn;
2319 }
2320
2321 void
2322 mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2323 {
2324   mpz_div_q_2exp (r, u, cnt, GMP_DIV_CEIL);
2325 }
2326
2327 void
2328 mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2329 {
2330   mpz_div_q_2exp (r, u, cnt, GMP_DIV_FLOOR);
2331 }
2332
2333 void
2334 mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2335 {
2336   mpz_div_q_2exp (r, u, cnt, GMP_DIV_TRUNC);
2337 }
2338
2339 void
2340 mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2341 {
2342   mpz_div_r_2exp (r, u, cnt, GMP_DIV_CEIL);
2343 }
2344
2345 void
2346 mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2347 {
2348   mpz_div_r_2exp (r, u, cnt, GMP_DIV_FLOOR);
2349 }
2350
2351 void
2352 mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
2353 {
2354   mpz_div_r_2exp (r, u, cnt, GMP_DIV_TRUNC);
2355 }
2356
2357 void
2358 mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d)
2359 {
2360   gmp_assert_nocarry (mpz_div_qr (q, NULL, n, d, GMP_DIV_TRUNC));
2361 }
2362
2363 int
2364 mpz_divisible_p (const mpz_t n, const mpz_t d)
2365 {
2366   return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2367 }
2368
2369 static unsigned long
2370 mpz_div_qr_ui (mpz_t q, mpz_t r,
2371                const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2372 {
2373   mp_size_t ns, qn;
2374   mp_ptr qp;
2375   mp_limb_t rl;
2376   mp_size_t rs;
2377
2378   ns = n->_mp_size;
2379   if (ns == 0)
2380     {
2381       if (q)
2382         q->_mp_size = 0;
2383       if (r)
2384         r->_mp_size = 0;
2385       return 0;
2386     }
2387
2388   qn = GMP_ABS (ns);
2389   if (q)
2390     qp = MPZ_REALLOC (q, qn);
2391   else
2392     qp = NULL;
2393
2394   rl = mpn_div_qr_1 (qp, n->_mp_d, qn, d);
2395   assert (rl < d);
2396
2397   rs = rl > 0;
2398   rs = (ns < 0) ? -rs : rs;
2399
2400   if (rl > 0 && ( (mode == GMP_DIV_FLOOR && ns < 0)
2401                   || (mode == GMP_DIV_CEIL && ns >= 0)))
2402     {
2403       if (q)
2404         gmp_assert_nocarry (mpn_add_1 (qp, qp, qn, 1));
2405       rl = d - rl;
2406       rs = -rs;
2407     }
2408
2409   if (r)
2410     {
2411       r->_mp_d[0] = rl;
2412       r->_mp_size = rs;
2413     }
2414   if (q)
2415     {
2416       qn -= (qp[qn-1] == 0);
2417       assert (qn == 0 || qp[qn-1] > 0);
2418
2419       q->_mp_size = (ns < 0) ? - qn : qn;
2420     }
2421
2422   return rl;
2423 }
2424
2425 unsigned long
2426 mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2427 {
2428   return mpz_div_qr_ui (q, r, n, d, GMP_DIV_CEIL);
2429 }
2430
2431 unsigned long
2432 mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2433 {
2434   return mpz_div_qr_ui (q, r, n, d, GMP_DIV_FLOOR);
2435 }
2436
2437 unsigned long
2438 mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long d)
2439 {
2440   return mpz_div_qr_ui (q, r, n, d, GMP_DIV_TRUNC);
2441 }
2442
2443 unsigned long
2444 mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2445 {
2446   return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_CEIL);
2447 }
2448
2449 unsigned long
2450 mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2451 {
2452   return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_FLOOR);
2453 }
2454
2455 unsigned long
2456 mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long d)
2457 {
2458   return mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC);
2459 }
2460
2461 unsigned long
2462 mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2463 {
2464   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_CEIL);
2465 }
2466 unsigned long
2467 mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2468 {
2469   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2470 }
2471 unsigned long
2472 mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long d)
2473 {
2474   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_TRUNC);
2475 }
2476
2477 unsigned long
2478 mpz_cdiv_ui (const mpz_t n, unsigned long d)
2479 {
2480   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_CEIL);
2481 }
2482
2483 unsigned long
2484 mpz_fdiv_ui (const mpz_t n, unsigned long d)
2485 {
2486   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_FLOOR);
2487 }
2488
2489 unsigned long
2490 mpz_tdiv_ui (const mpz_t n, unsigned long d)
2491 {
2492   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC);
2493 }
2494
2495 unsigned long
2496 mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long d)
2497 {
2498   return mpz_div_qr_ui (NULL, r, n, d, GMP_DIV_FLOOR);
2499 }
2500
2501 void
2502 mpz_divexact_ui (mpz_t q, const mpz_t n, unsigned long d)
2503 {
2504   gmp_assert_nocarry (mpz_div_qr_ui (q, NULL, n, d, GMP_DIV_TRUNC));
2505 }
2506
2507 int
2508 mpz_divisible_ui_p (const mpz_t n, unsigned long d)
2509 {
2510   return mpz_div_qr_ui (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2511 }
2512
2513 \f
2514 /* GCD */
2515 static mp_limb_t
2516 mpn_gcd_11 (mp_limb_t u, mp_limb_t v)
2517 {
2518   unsigned shift;
2519
2520   assert ( (u | v) > 0);
2521
2522   if (u == 0)
2523     return v;
2524   else if (v == 0)
2525     return u;
2526
2527   gmp_ctz (shift, u | v);
2528
2529   u >>= shift;
2530   v >>= shift;
2531
2532   if ( (u & 1) == 0)
2533     MP_LIMB_T_SWAP (u, v);
2534
2535   while ( (v & 1) == 0)
2536     v >>= 1;
2537
2538   while (u != v)
2539     {
2540       if (u > v)
2541         {
2542           u -= v;
2543           do
2544             u >>= 1;
2545           while ( (u & 1) == 0);
2546         }
2547       else
2548         {
2549           v -= u;
2550           do
2551             v >>= 1;
2552           while ( (v & 1) == 0);
2553         }
2554     }
2555   return u << shift;
2556 }
2557
2558 unsigned long
2559 mpz_gcd_ui (mpz_t g, const mpz_t u, unsigned long v)
2560 {
2561   mp_size_t un;
2562
2563   if (v == 0)
2564     {
2565       if (g)
2566         mpz_abs (g, u);
2567     }
2568   else
2569     {
2570       un = GMP_ABS (u->_mp_size);
2571       if (un != 0)
2572         v = mpn_gcd_11 (mpn_div_qr_1 (NULL, u->_mp_d, un, v), v);
2573
2574       if (g)
2575         mpz_set_ui (g, v);
2576     }
2577
2578   return v;
2579 }
2580
2581 static mp_bitcnt_t
2582 mpz_make_odd (mpz_t r, const mpz_t u)
2583 {
2584   mp_size_t un, rn, i;
2585   mp_ptr rp;
2586   unsigned shift;
2587
2588   un = GMP_ABS (u->_mp_size);
2589   assert (un > 0);
2590
2591   for (i = 0; u->_mp_d[i] == 0; i++)
2592     ;
2593
2594   gmp_ctz (shift, u->_mp_d[i]);
2595
2596   rn = un - i;
2597   rp = MPZ_REALLOC (r, rn);
2598   if (shift > 0)
2599     {
2600       mpn_rshift (rp, u->_mp_d + i, rn, shift);
2601       rn -= (rp[rn-1] == 0);
2602     }
2603   else
2604     mpn_copyi (rp, u->_mp_d + i, rn);
2605
2606   r->_mp_size = rn;
2607   return i * GMP_LIMB_BITS + shift;
2608 }
2609
2610 void
2611 mpz_gcd (mpz_t g, const mpz_t u, const mpz_t v)
2612 {
2613   mpz_t tu, tv;
2614   mp_bitcnt_t uz, vz, gz;
2615
2616   if (u->_mp_size == 0)
2617     {
2618       mpz_abs (g, v);
2619       return;
2620     }
2621   if (v->_mp_size == 0)
2622     {
2623       mpz_abs (g, u);
2624       return;
2625     }
2626
2627   mpz_init (tu);
2628   mpz_init (tv);
2629
2630   uz = mpz_make_odd (tu, u);
2631   vz = mpz_make_odd (tv, v);
2632   gz = GMP_MIN (uz, vz);
2633
2634   if (tu->_mp_size < tv->_mp_size)
2635     mpz_swap (tu, tv);
2636
2637   mpz_tdiv_r (tu, tu, tv);
2638   if (tu->_mp_size == 0)
2639     {
2640       mpz_swap (g, tv);
2641     }
2642   else
2643     for (;;)
2644       {
2645         int c;
2646
2647         mpz_make_odd (tu, tu);
2648         c = mpz_cmp (tu, tv);
2649         if (c == 0)
2650           {
2651             mpz_swap (g, tu);
2652             break;
2653           }
2654         if (c < 0)
2655           mpz_swap (tu, tv);
2656
2657         if (tv->_mp_size == 1)
2658           {
2659             mp_limb_t vl = tv->_mp_d[0];
2660             mp_limb_t ul = mpz_tdiv_ui (tu, vl);
2661             mpz_set_ui (g, mpn_gcd_11 (ul, vl));
2662             break;
2663           }
2664         mpz_sub (tu, tu, tv);
2665       }
2666   mpz_clear (tu);
2667   mpz_clear (tv);
2668   mpz_mul_2exp (g, g, gz);
2669 }
2670
2671 void
2672 mpz_gcdext (mpz_t g, mpz_t s, mpz_t t, const mpz_t u, const mpz_t v)
2673 {
2674   mpz_t tu, tv, s0, s1, t0, t1;
2675   mp_bitcnt_t uz, vz, gz;
2676   mp_bitcnt_t power;
2677
2678   if (u->_mp_size == 0)
2679     {
2680       /* g = 0 u + sgn(v) v */
2681       signed long sign = mpz_sgn (v);
2682       mpz_abs (g, v);
2683       if (s)
2684         mpz_set_ui (s, 0);
2685       if (t)
2686         mpz_set_si (t, sign);
2687       return;
2688     }
2689
2690   if (v->_mp_size == 0)
2691     {
2692       /* g = sgn(u) u + 0 v */
2693       signed long sign = mpz_sgn (u);
2694       mpz_abs (g, u);
2695       if (s)
2696         mpz_set_si (s, sign);
2697       if (t)
2698         mpz_set_ui (t, 0);
2699       return;
2700     }
2701
2702   mpz_init (tu);
2703   mpz_init (tv);
2704   mpz_init (s0);
2705   mpz_init (s1);
2706   mpz_init (t0);
2707   mpz_init (t1);
2708
2709   uz = mpz_make_odd (tu, u);
2710   vz = mpz_make_odd (tv, v);
2711   gz = GMP_MIN (uz, vz);
2712
2713   uz -= gz;
2714   vz -= gz;
2715
2716   /* Cofactors corresponding to odd gcd. gz handled later. */
2717   if (tu->_mp_size < tv->_mp_size)
2718     {
2719       mpz_swap (tu, tv);
2720       MPZ_SRCPTR_SWAP (u, v);
2721       MPZ_PTR_SWAP (s, t);
2722       MP_BITCNT_T_SWAP (uz, vz);
2723     }
2724
2725   /* Maintain
2726    *
2727    * u = t0 tu + t1 tv
2728    * v = s0 tu + s1 tv
2729    *
2730    * where u and v denote the inputs with common factors of two
2731    * eliminated, and det (s0, t0; s1, t1) = 2^p. Then
2732    *
2733    * 2^p tu =  s1 u - t1 v
2734    * 2^p tv = -s0 u + t0 v
2735    */
2736
2737   /* After initial division, tu = q tv + tu', we have
2738    *
2739    * u = 2^uz (tu' + q tv)
2740    * v = 2^vz tv
2741    *
2742    * or
2743    *
2744    * t0 = 2^uz, t1 = 2^uz q
2745    * s0 = 0,    s1 = 2^vz
2746    */
2747
2748   mpz_setbit (t0, uz);
2749   mpz_tdiv_qr (t1, tu, tu, tv);
2750   mpz_mul_2exp (t1, t1, uz);
2751
2752   mpz_setbit (s1, vz);
2753   power = uz + vz;
2754
2755   if (tu->_mp_size > 0)
2756     {
2757       mp_bitcnt_t shift;
2758       shift = mpz_make_odd (tu, tu);
2759       mpz_mul_2exp (t0, t0, shift);
2760       mpz_mul_2exp (s0, s0, shift);
2761       power += shift;
2762
2763       for (;;)
2764         {
2765           int c;
2766           c = mpz_cmp (tu, tv);
2767           if (c == 0)
2768             break;
2769
2770           if (c < 0)
2771             {
2772               /* tv = tv' + tu
2773                *
2774                * u = t0 tu + t1 (tv' + tu) = (t0 + t1) tu + t1 tv'
2775                * v = s0 tu + s1 (tv' + tu) = (s0 + s1) tu + s1 tv' */
2776
2777               mpz_sub (tv, tv, tu);
2778               mpz_add (t0, t0, t1);
2779               mpz_add (s0, s0, s1);
2780
2781               shift = mpz_make_odd (tv, tv);
2782               mpz_mul_2exp (t1, t1, shift);
2783               mpz_mul_2exp (s1, s1, shift);
2784             }
2785           else
2786             {
2787               mpz_sub (tu, tu, tv);
2788               mpz_add (t1, t0, t1);
2789               mpz_add (s1, s0, s1);
2790
2791               shift = mpz_make_odd (tu, tu);
2792               mpz_mul_2exp (t0, t0, shift);
2793               mpz_mul_2exp (s0, s0, shift);
2794             }
2795           power += shift;
2796         }
2797     }
2798
2799   /* Now tv = odd part of gcd, and -s0 and t0 are corresponding
2800      cofactors. */
2801
2802   mpz_mul_2exp (tv, tv, gz);
2803   mpz_neg (s0, s0);
2804
2805   /* 2^p g = s0 u + t0 v. Eliminate one factor of two at a time. To
2806      adjust cofactors, we need u / g and v / g */
2807
2808   mpz_divexact (s1, v, tv);
2809   mpz_abs (s1, s1);
2810   mpz_divexact (t1, u, tv);
2811   mpz_abs (t1, t1);
2812
2813   while (power-- > 0)
2814     {
2815       /* s0 u + t0 v = (s0 - v/g) u - (t0 + u/g) v */
2816       if (mpz_odd_p (s0) || mpz_odd_p (t0))
2817         {
2818           mpz_sub (s0, s0, s1);
2819           mpz_add (t0, t0, t1);
2820         }
2821       mpz_divexact_ui (s0, s0, 2);
2822       mpz_divexact_ui (t0, t0, 2);
2823     }
2824
2825   /* Arrange so that |s| < |u| / 2g */
2826   mpz_add (s1, s0, s1);
2827   if (mpz_cmpabs (s0, s1) > 0)
2828     {
2829       mpz_swap (s0, s1);
2830       mpz_sub (t0, t0, t1);
2831     }
2832   if (u->_mp_size < 0)
2833     mpz_neg (s0, s0);
2834   if (v->_mp_size < 0)
2835     mpz_neg (t0, t0);
2836
2837   mpz_swap (g, tv);
2838   if (s)
2839     mpz_swap (s, s0);
2840   if (t)
2841     mpz_swap (t, t0);
2842
2843   mpz_clear (tu);
2844   mpz_clear (tv);
2845   mpz_clear (s0);
2846   mpz_clear (s1);
2847   mpz_clear (t0);
2848   mpz_clear (t1);
2849 }
2850
2851 void
2852 mpz_lcm (mpz_t r, const mpz_t u, const mpz_t v)
2853 {
2854   mpz_t g;
2855
2856   if (u->_mp_size == 0 || v->_mp_size == 0)
2857     {
2858       r->_mp_size = 0;
2859       return;
2860     }
2861
2862   mpz_init (g);
2863
2864   mpz_gcd (g, u, v);
2865   mpz_divexact (g, u, g);
2866   mpz_mul (r, g, v);
2867
2868   mpz_clear (g);
2869   mpz_abs (r, r);
2870 }
2871
2872 void
2873 mpz_lcm_ui (mpz_t r, const mpz_t u, unsigned long v)
2874 {
2875   if (v == 0 || u->_mp_size == 0)
2876     {
2877       r->_mp_size = 0;
2878       return;
2879     }
2880
2881   v /= mpz_gcd_ui (NULL, u, v);
2882   mpz_mul_ui (r, u, v);
2883
2884   mpz_abs (r, r);
2885 }
2886
2887 int
2888 mpz_invert (mpz_t r, const mpz_t u, const mpz_t m)
2889 {
2890   mpz_t g, tr;
2891   int invertible;
2892
2893   if (u->_mp_size == 0 || mpz_cmpabs_ui (m, 1) <= 0)
2894     return 0;
2895
2896   mpz_init (g);
2897   mpz_init (tr);
2898
2899   mpz_gcdext (g, tr, NULL, u, m);
2900   invertible = (mpz_cmp_ui (g, 1) == 0);
2901
2902   if (invertible)
2903     {
2904       if (tr->_mp_size < 0)
2905         {
2906           if (m->_mp_size >= 0)
2907             mpz_add (tr, tr, m);
2908           else
2909             mpz_sub (tr, tr, m);
2910         }
2911       mpz_swap (r, tr);
2912     }
2913
2914   mpz_clear (g);
2915   mpz_clear (tr);
2916   return invertible;
2917 }
2918
2919 \f
2920 /* Higher level operations (sqrt, pow and root) */
2921
2922 void
2923 mpz_pow_ui (mpz_t r, const mpz_t b, unsigned long e)
2924 {
2925   unsigned long bit;
2926   mpz_t tr;
2927   mpz_init_set_ui (tr, 1);
2928
2929   for (bit = GMP_ULONG_HIGHBIT; bit > 0; bit >>= 1)
2930     {
2931       mpz_mul (tr, tr, tr);
2932       if (e & bit)
2933         mpz_mul (tr, tr, b);
2934     }
2935   mpz_swap (r, tr);
2936   mpz_clear (tr);
2937 }
2938
2939 void
2940 mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
2941 {
2942   mpz_t b;
2943   mpz_init_set_ui (b, blimb);
2944   mpz_pow_ui (r, b, e);
2945   mpz_clear (b);
2946 }
2947
2948 void
2949 mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
2950 {
2951   mpz_t tr;
2952   mpz_t base;
2953   mp_size_t en, mn;
2954   mp_srcptr mp;
2955   struct gmp_div_inverse minv;
2956   unsigned shift;
2957   mp_ptr tp = NULL;
2958
2959   en = GMP_ABS (e->_mp_size);
2960   mn = GMP_ABS (m->_mp_size);
2961   if (mn == 0)
2962     gmp_die ("mpz_powm: Zero modulo.");
2963
2964   if (en == 0)
2965     {
2966       mpz_set_ui (r, 1);
2967       return;
2968     }
2969
2970   mp = m->_mp_d;
2971   mpn_div_qr_invert (&minv, mp, mn);
2972   shift = minv.shift;
2973
2974   if (shift > 0)
2975     {
2976       /* To avoid shifts, we do all our reductions, except the final
2977          one, using a *normalized* m. */
2978       minv.shift = 0;
2979
2980       tp = gmp_xalloc_limbs (mn);
2981       gmp_assert_nocarry (mpn_lshift (tp, mp, mn, shift));
2982       mp = tp;
2983     }
2984
2985   mpz_init (base);
2986
2987   if (e->_mp_size < 0)
2988     {
2989       if (!mpz_invert (base, b, m))
2990         gmp_die ("mpz_powm: Negative exponent and non-invertibe base.");
2991     }
2992   else
2993     {
2994       mp_size_t bn;
2995       mpz_abs (base, b);
2996
2997       bn = base->_mp_size;
2998       if (bn >= mn)
2999         {
3000           mpn_div_qr_preinv (NULL, base->_mp_d, base->_mp_size, mp, mn, &minv);
3001           bn = mn;
3002         }
3003
3004       /* We have reduced the absolute value. Now take care of the
3005          sign. Note that we get zero represented non-canonically as
3006          m. */
3007       if (b->_mp_size < 0)
3008         {
3009           mp_ptr bp = MPZ_REALLOC (base, mn);
3010           gmp_assert_nocarry (mpn_sub (bp, mp, mn, bp, bn));
3011           bn = mn;
3012         }
3013       base->_mp_size = mpn_normalized_size (base->_mp_d, bn);
3014     }
3015   mpz_init_set_ui (tr, 1);
3016
3017   while (en-- > 0)
3018     {
3019       mp_limb_t w = e->_mp_d[en];
3020       mp_limb_t bit;
3021
3022       for (bit = GMP_LIMB_HIGHBIT; bit > 0; bit >>= 1)
3023         {
3024           mpz_mul (tr, tr, tr);
3025           if (w & bit)
3026             mpz_mul (tr, tr, base);
3027           if (tr->_mp_size > mn)
3028             {
3029               mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3030               tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3031             }
3032         }
3033     }
3034
3035   /* Final reduction */
3036   if (tr->_mp_size >= mn)
3037     {
3038       minv.shift = shift;
3039       mpn_div_qr_preinv (NULL, tr->_mp_d, tr->_mp_size, mp, mn, &minv);
3040       tr->_mp_size = mpn_normalized_size (tr->_mp_d, mn);
3041     }
3042   if (tp)
3043     gmp_free (tp);
3044
3045   mpz_swap (r, tr);
3046   mpz_clear (tr);
3047   mpz_clear (base);
3048 }
3049
3050 void
3051 mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3052 {
3053   mpz_t e;
3054   mpz_init_set_ui (e, elimb);
3055   mpz_powm (r, b, e, m);
3056   mpz_clear (e);
3057 }
3058
3059 /* x=trunc(y^(1/z)), r=y-x^z */
3060 void
3061 mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3062 {
3063   int sgn;
3064   mpz_t t, u;
3065
3066   sgn = y->_mp_size < 0;
3067   if (sgn && (z & 1) == 0)
3068     gmp_die ("mpz_rootrem: Negative argument, with even root.");
3069   if (z == 0)
3070     gmp_die ("mpz_rootrem: Zeroth root.");
3071
3072   if (mpz_cmpabs_ui (y, 1) <= 0) {
3073     mpz_set (x, y);
3074     if (r)
3075       r->_mp_size = 0;
3076     return;
3077   }
3078
3079   mpz_init (t);
3080   mpz_init (u);
3081   mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3082
3083   if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3084     do {
3085       mpz_swap (u, t);                  /* u = x */
3086       mpz_tdiv_q (t, y, u);             /* t = y/x */
3087       mpz_add (t, t, u);                /* t = y/x + x */
3088       mpz_tdiv_q_2exp (t, t, 1);        /* x'= (y/x + x)/2 */
3089     } while (mpz_cmpabs (t, u) < 0);    /* |x'| < |x| */
3090   else /* z != 2 */ {
3091     mpz_t v;
3092
3093     mpz_init (v);
3094     if (sgn)
3095       mpz_neg (t, t);
3096
3097     do {
3098       mpz_swap (u, t);                  /* u = x */
3099       mpz_pow_ui (t, u, z - 1);         /* t = x^(z-1) */
3100       mpz_tdiv_q (t, y, t);             /* t = y/x^(z-1) */
3101       mpz_mul_ui (v, u, z - 1);         /* v = x*(z-1) */
3102       mpz_add (t, t, v);                /* t = y/x^(z-1) + x*(z-1) */
3103       mpz_tdiv_q_ui (t, t, z);          /* x'=(y/x^(z-1) + x*(z-1))/z */
3104     } while (mpz_cmpabs (t, u) < 0);    /* |x'| < |x| */
3105
3106     mpz_clear (v);
3107   }
3108
3109   if (r) {
3110     mpz_pow_ui (t, u, z);
3111     mpz_sub (r, y, t);
3112   }
3113   mpz_swap (x, u);
3114   mpz_clear (u);
3115   mpz_clear (t);
3116 }
3117
3118 int
3119 mpz_root (mpz_t x, const mpz_t y, unsigned long z)
3120 {
3121   int res;
3122   mpz_t r;
3123
3124   mpz_init (r);
3125   mpz_rootrem (x, r, y, z);
3126   res = r->_mp_size == 0;
3127   mpz_clear (r);
3128
3129   return res;
3130 }
3131
3132 /* Compute s = floor(sqrt(u)) and r = u - s^2. Allows r == NULL */
3133 void
3134 mpz_sqrtrem (mpz_t s, mpz_t r, const mpz_t u)
3135 {
3136   mpz_rootrem (s, r, u, 2);
3137 }
3138
3139 void
3140 mpz_sqrt (mpz_t s, const mpz_t u)
3141 {
3142   mpz_rootrem (s, NULL, u, 2);
3143 }
3144
3145 \f
3146 /* Combinatorics */
3147
3148 void
3149 mpz_fac_ui (mpz_t x, unsigned long n)
3150 {
3151   if (n < 2) {
3152     mpz_set_ui (x, 1);
3153     return;
3154   }
3155   mpz_set_ui (x, n);
3156   for (;--n > 1;)
3157     mpz_mul_ui (x, x, n);
3158 }
3159
3160 void
3161 mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3162 {
3163   mpz_t t;
3164
3165   if (k > n) {
3166     r->_mp_size = 0;
3167     return;
3168   }
3169   mpz_fac_ui (r, n);
3170   mpz_init (t);
3171   mpz_fac_ui (t, k);
3172   mpz_divexact (r, r, t);
3173   mpz_fac_ui (t, n - k);
3174   mpz_divexact (r, r, t);
3175   mpz_clear (t);
3176 }
3177
3178 \f
3179 /* Logical operations and bit manipulation. */
3180
3181 /* Numbers are treated as if represented in two's complement (and
3182    infinitely sign extended). For a negative values we get the two's
3183    complement from -x = ~x + 1, where ~ is bitwise complementt.
3184    Negation transforms
3185
3186      xxxx10...0
3187
3188    into
3189
3190      yyyy10...0
3191
3192    where yyyy is the bitwise complement of xxxx. So least significant
3193    bits, up to and including the first one bit, are unchanged, and
3194    the more significant bits are all complemented.
3195
3196    To change a bit from zero to one in a negative number, subtract the
3197    corresponding power of two from the absolute value. This can never
3198    underflow. To change a bit from one to zero, add the corresponding
3199    power of two, and this might overflow. E.g., if x = -001111, the
3200    two's complement is 110001. Clearing the least significant bit, we
3201    get two's complement 110000, and -010000. */
3202
3203 int
3204 mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3205 {
3206   mp_size_t limb_index;
3207   unsigned shift;
3208   mp_size_t ds;
3209   mp_size_t dn;
3210   mp_limb_t w;
3211   int bit;
3212
3213   ds = d->_mp_size;
3214   dn = GMP_ABS (ds);
3215   limb_index = bit_index / GMP_LIMB_BITS;
3216   if (limb_index >= dn)
3217     return ds < 0;
3218
3219   shift = bit_index % GMP_LIMB_BITS;
3220   w = d->_mp_d[limb_index];
3221   bit = (w >> shift) & 1;
3222
3223   if (ds < 0)
3224     {
3225       /* d < 0. Check if any of the bits below is set: If so, our bit
3226          must be complemented. */
3227       if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3228         return bit ^ 1;
3229       while (limb_index-- > 0)
3230         if (d->_mp_d[limb_index] > 0)
3231           return bit ^ 1;
3232     }
3233   return bit;
3234 }
3235
3236 static void
3237 mpz_abs_add_bit (mpz_t d, mp_bitcnt_t bit_index)
3238 {
3239   mp_size_t dn, limb_index;
3240   mp_limb_t bit;
3241   mp_ptr dp;
3242
3243   dn = GMP_ABS (d->_mp_size);
3244
3245   limb_index = bit_index / GMP_LIMB_BITS;
3246   bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3247
3248   if (limb_index >= dn)
3249     {
3250       mp_size_t i;
3251       /* The bit should be set outside of the end of the number.
3252          We have to increase the size of the number. */
3253       dp = MPZ_REALLOC (d, limb_index + 1);
3254
3255       dp[limb_index] = bit;
3256       for (i = dn; i < limb_index; i++)
3257         dp[i] = 0;
3258       dn = limb_index + 1;
3259     }
3260   else
3261     {
3262       mp_limb_t cy;
3263
3264       dp = d->_mp_d;
3265
3266       cy = mpn_add_1 (dp + limb_index, dp + limb_index, dn - limb_index, bit);
3267       if (cy > 0)
3268         {
3269           dp = MPZ_REALLOC (d, dn + 1);
3270           dp[dn++] = cy;
3271         }
3272     }
3273
3274   d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3275 }
3276
3277 static void
3278 mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3279 {
3280   mp_size_t dn, limb_index;
3281   mp_ptr dp;
3282   mp_limb_t bit;
3283
3284   dn = GMP_ABS (d->_mp_size);
3285   dp = d->_mp_d;
3286
3287   limb_index = bit_index / GMP_LIMB_BITS;
3288   bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
3289
3290   assert (limb_index < dn);
3291
3292   gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3293                                  dn - limb_index, bit));
3294   dn -= (dp[dn-1] == 0);
3295   d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3296 }
3297
3298 void
3299 mpz_setbit (mpz_t d, mp_bitcnt_t bit_index)
3300 {
3301   if (!mpz_tstbit (d, bit_index))
3302     {
3303       if (d->_mp_size >= 0)
3304         mpz_abs_add_bit (d, bit_index);
3305       else
3306         mpz_abs_sub_bit (d, bit_index);
3307     }
3308 }
3309
3310 void
3311 mpz_clrbit (mpz_t d, mp_bitcnt_t bit_index)
3312 {
3313   if (mpz_tstbit (d, bit_index))
3314     {
3315       if (d->_mp_size >= 0)
3316         mpz_abs_sub_bit (d, bit_index);
3317       else
3318         mpz_abs_add_bit (d, bit_index);
3319     }
3320 }
3321
3322 void
3323 mpz_combit (mpz_t d, mp_bitcnt_t bit_index)
3324 {
3325   if (mpz_tstbit (d, bit_index) ^ (d->_mp_size < 0))
3326     mpz_abs_sub_bit (d, bit_index);
3327   else
3328     mpz_abs_add_bit (d, bit_index);
3329 }
3330
3331 void
3332 mpz_com (mpz_t r, const mpz_t u)
3333 {
3334   mpz_neg (r, u);
3335   mpz_sub_ui (r, r, 1);
3336 }
3337
3338 void
3339 mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3340 {
3341   mp_size_t un, vn, rn, i;
3342   mp_ptr up, vp, rp;
3343
3344   mp_limb_t ux, vx, rx;
3345   mp_limb_t uc, vc, rc;
3346   mp_limb_t ul, vl, rl;
3347
3348   un = GMP_ABS (u->_mp_size);
3349   vn = GMP_ABS (v->_mp_size);
3350   if (un < vn)
3351     {
3352       MPZ_SRCPTR_SWAP (u, v);
3353       MP_SIZE_T_SWAP (un, vn);
3354     }
3355   if (vn == 0)
3356     {
3357       r->_mp_size = 0;
3358       return;
3359     }
3360
3361   uc = u->_mp_size < 0;
3362   vc = v->_mp_size < 0;
3363   rc = uc & vc;
3364
3365   ux = -uc;
3366   vx = -vc;
3367   rx = -rc;
3368
3369   /* If the smaller input is positive, higher limbs don't matter. */
3370   rn = vx ? un : vn;
3371
3372   rp = MPZ_REALLOC (r, rn + rc);
3373
3374   up = u->_mp_d;
3375   vp = v->_mp_d;
3376
3377   for (i = 0; i < vn; i++)
3378     {
3379       ul = (up[i] ^ ux) + uc;
3380       uc = ul < uc;
3381
3382       vl = (vp[i] ^ vx) + vc;
3383       vc = vl < vc;
3384
3385       rl = ( (ul & vl) ^ rx) + rc;
3386       rc = rl < rc;
3387       rp[i] = rl;
3388     }
3389   assert (vc == 0);
3390
3391   for (; i < rn; i++)
3392     {
3393       ul = (up[i] ^ ux) + uc;
3394       uc = ul < uc;
3395
3396       rl = ( (ul & vx) ^ rx) + rc;
3397       rc = rl < rc;
3398       rp[i] = rl;
3399     }
3400   if (rc)
3401     rp[rn++] = rc;
3402   else
3403     rn = mpn_normalized_size (rp, rn);
3404
3405   r->_mp_size = rx ? -rn : rn;
3406 }
3407
3408 void
3409 mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3410 {
3411   mp_size_t un, vn, rn, i;
3412   mp_ptr up, vp, rp;
3413
3414   mp_limb_t ux, vx, rx;
3415   mp_limb_t uc, vc, rc;
3416   mp_limb_t ul, vl, rl;
3417
3418   un = GMP_ABS (u->_mp_size);
3419   vn = GMP_ABS (v->_mp_size);
3420   if (un < vn)
3421     {
3422       MPZ_SRCPTR_SWAP (u, v);
3423       MP_SIZE_T_SWAP (un, vn);
3424     }
3425   if (vn == 0)
3426     {
3427       mpz_set (r, u);
3428       return;
3429     }
3430
3431   uc = u->_mp_size < 0;
3432   vc = v->_mp_size < 0;
3433   rc = uc | vc;
3434
3435   ux = -uc;
3436   vx = -vc;
3437   rx = -rc;
3438
3439   /* If the smaller input is negative, by sign extension higher limbs
3440      don't matter. */
3441   rn = vx ? vn : un;
3442
3443   rp = MPZ_REALLOC (r, rn + rc);
3444
3445   up = u->_mp_d;
3446   vp = v->_mp_d;
3447
3448   for (i = 0; i < vn; i++)
3449     {
3450       ul = (up[i] ^ ux) + uc;
3451       uc = ul < uc;
3452
3453       vl = (vp[i] ^ vx) + vc;
3454       vc = vl < vc;
3455
3456       rl = ( (ul | vl) ^ rx) + rc;
3457       rc = rl < rc;
3458       rp[i] = rl;
3459     }
3460   assert (vc == 0);
3461
3462   for (; i < rn; i++)
3463     {
3464       ul = (up[i] ^ ux) + uc;
3465       uc = ul < uc;
3466
3467       rl = ( (ul | vx) ^ rx) + rc;
3468       rc = rl < rc;
3469       rp[i] = rl;
3470     }
3471   if (rc)
3472     rp[rn++] = rc;
3473   else
3474     rn = mpn_normalized_size (rp, rn);
3475
3476   r->_mp_size = rx ? -rn : rn;
3477 }
3478
3479 void
3480 mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3481 {
3482   mp_size_t un, vn, i;
3483   mp_ptr up, vp, rp;
3484
3485   mp_limb_t ux, vx, rx;
3486   mp_limb_t uc, vc, rc;
3487   mp_limb_t ul, vl, rl;
3488
3489   un = GMP_ABS (u->_mp_size);
3490   vn = GMP_ABS (v->_mp_size);
3491   if (un < vn)
3492     {
3493       MPZ_SRCPTR_SWAP (u, v);
3494       MP_SIZE_T_SWAP (un, vn);
3495     }
3496   if (vn == 0)
3497     {
3498       mpz_set (r, u);
3499       return;
3500     }
3501
3502   uc = u->_mp_size < 0;
3503   vc = v->_mp_size < 0;
3504   rc = uc ^ vc;
3505
3506   ux = -uc;
3507   vx = -vc;
3508   rx = -rc;
3509
3510   rp = MPZ_REALLOC (r, un + rc);
3511
3512   up = u->_mp_d;
3513   vp = v->_mp_d;
3514
3515   for (i = 0; i < vn; i++)
3516     {
3517       ul = (up[i] ^ ux) + uc;
3518       uc = ul < uc;
3519
3520       vl = (vp[i] ^ vx) + vc;
3521       vc = vl < vc;
3522
3523       rl = (ul ^ vl ^ rx) + rc;
3524       rc = rl < rc;
3525       rp[i] = rl;
3526     }
3527   assert (vc == 0);
3528
3529   for (; i < un; i++)
3530     {
3531       ul = (up[i] ^ ux) + uc;
3532       uc = ul < uc;
3533
3534       rl = (ul ^ ux) + rc;
3535       rc = rl < rc;
3536       rp[i] = rl;
3537     }
3538   if (rc)
3539     rp[un++] = rc;
3540   else
3541     un = mpn_normalized_size (rp, un);
3542
3543   r->_mp_size = rx ? -un : un;
3544 }
3545
3546 static unsigned
3547 gmp_popcount_limb (mp_limb_t x)
3548 {
3549   unsigned c;
3550
3551   /* Do 16 bits at a time, to avoid limb-sized constants. */
3552   for (c = 0; x > 0; x >>= 16)
3553     {
3554       unsigned w = ((x >> 1) & 0x5555) + (x & 0x5555);
3555       w = ((w >> 2) & 0x3333) + (w & 0x3333);
3556       w = ((w >> 4) & 0x0f0f) + (w & 0x0f0f);
3557       w = (w >> 8) + (w & 0x00ff);
3558       c += w;
3559     }
3560   return c;
3561 }
3562
3563 mp_bitcnt_t
3564 mpz_popcount (const mpz_t u)
3565 {
3566   mp_size_t un, i;
3567   mp_bitcnt_t c;
3568
3569   un = u->_mp_size;
3570
3571   if (un < 0)
3572     return ~(mp_bitcnt_t) 0;
3573
3574   for (c = 0, i = 0; i < un; i++)
3575     c += gmp_popcount_limb (u->_mp_d[i]);
3576
3577   return c;
3578 }
3579
3580 mp_bitcnt_t
3581 mpz_hamdist (const mpz_t u, const mpz_t v)
3582 {
3583   mp_size_t un, vn, i;
3584   mp_limb_t uc, vc, ul, vl, comp;
3585   mp_srcptr up, vp;
3586   mp_bitcnt_t c;
3587
3588   un = u->_mp_size;
3589   vn = v->_mp_size;
3590
3591   if ( (un ^ vn) < 0)
3592     return ~(mp_bitcnt_t) 0;
3593
3594   if (un < 0)
3595     {
3596       assert (vn < 0);
3597       un = -un;
3598       vn = -vn;
3599       uc = vc = 1;
3600       comp = - (mp_limb_t) 1;
3601     }
3602   else
3603     uc = vc = comp = 0;
3604
3605   up = u->_mp_d;
3606   vp = v->_mp_d;
3607
3608   if (un < vn)
3609     MPN_SRCPTR_SWAP (up, un, vp, vn);
3610
3611   for (i = 0, c = 0; i < vn; i++)
3612     {
3613       ul = (up[i] ^ comp) + uc;
3614       uc = ul < uc;
3615
3616       vl = (vp[i] ^ comp) + vc;
3617       vc = vl < vc;
3618
3619       c += gmp_popcount_limb (ul ^ vl);
3620     }
3621   assert (vc == 0);
3622
3623   for (; i < un; i++)
3624     {
3625       ul = (up[i] ^ comp) + uc;
3626       uc = ul < uc;
3627
3628       c += gmp_popcount_limb (ul ^ comp);
3629     }
3630
3631   return c;
3632 }
3633
3634 mp_bitcnt_t
3635 mpz_scan1 (const mpz_t u, mp_bitcnt_t starting_bit)
3636 {
3637   mp_ptr up;
3638   mp_size_t us, un, i;
3639   mp_limb_t limb, ux, uc;
3640   unsigned cnt;
3641
3642   up = u->_mp_d;
3643   us = u->_mp_size;
3644   un = GMP_ABS (us);
3645   i = starting_bit / GMP_LIMB_BITS;
3646
3647   /* Past the end there's no 1 bits for u>=0, or an immediate 1 bit
3648      for u<0. Notice this test picks up any u==0 too. */
3649   if (i >= un)
3650     return (us >= 0 ? ~(mp_bitcnt_t) 0 : starting_bit);
3651
3652   if (us < 0)
3653     {
3654       ux = GMP_LIMB_MAX;
3655       uc = mpn_zero_p (up, i);
3656     }
3657   else
3658     ux = uc = 0;
3659
3660   limb = (ux ^ up[i]) + uc;
3661   uc = limb < uc;
3662
3663   /* Mask to 0 all bits before starting_bit, thus ignoring them. */
3664   limb &= (GMP_LIMB_MAX << (starting_bit % GMP_LIMB_BITS));
3665
3666   while (limb == 0)
3667     {
3668       i++;
3669       if (i == un)
3670         {
3671           assert (uc == 0);
3672           /* For the u > 0 case, this can happen only for the first
3673              masked limb. For the u < 0 case, it happens when the
3674              highest limbs of the absolute value are all ones. */
3675           return (us >= 0 ? ~(mp_bitcnt_t) 0 : un * GMP_LIMB_BITS);
3676         }
3677       limb = (ux ^ up[i]) + uc;
3678       uc = limb < uc;
3679     }
3680   gmp_ctz (cnt, limb);
3681   return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
3682 }
3683
3684 mp_bitcnt_t
3685 mpz_scan0 (const mpz_t u, mp_bitcnt_t starting_bit)
3686 {
3687   mp_ptr up;
3688   mp_size_t us, un, i;
3689   mp_limb_t limb, ux, uc;
3690   unsigned cnt;
3691
3692   up = u->_mp_d;
3693   us = u->_mp_size;
3694   un = GMP_ABS (us);
3695   i = starting_bit / GMP_LIMB_BITS;
3696
3697   /* When past end, there's an immediate 0 bit for u>=0, or no 0 bits for
3698      u<0.  Notice this test picks up all cases of u==0 too. */
3699   if (i >= un)
3700     return (us >= 0 ? starting_bit : ~(mp_bitcnt_t) 0);
3701
3702   if (us < 0)
3703     {
3704       ux = GMP_LIMB_MAX;
3705       uc = mpn_zero_p (up, i);
3706     }
3707   else
3708     ux = uc = 0;
3709
3710   limb = (ux ^ up[i]) + uc;
3711   uc = limb < uc;
3712
3713   /* Mask to 1 all bits before starting_bit, thus ignoring them. */
3714   limb |= ((mp_limb_t) 1 << (starting_bit % GMP_LIMB_BITS)) - 1;
3715
3716   while (limb == GMP_LIMB_MAX)
3717     {
3718       i++;
3719       if (i == un)
3720         {
3721           assert (uc == 0);
3722           return (us >= 0 ? un * GMP_LIMB_BITS : ~(mp_bitcnt_t) 0);
3723         }
3724       limb = (ux ^ up[i]) + uc;
3725       uc = limb < uc;
3726     }
3727   gmp_ctz (cnt, ~limb);
3728   return (mp_bitcnt_t) i * GMP_LIMB_BITS + cnt;
3729 }
3730
3731 \f
3732 /* MPZ base conversion. */
3733
3734 size_t
3735 mpz_sizeinbase (const mpz_t u, int base)
3736 {
3737   mp_size_t un;
3738   mp_srcptr up;
3739   mp_ptr tp;
3740   mp_bitcnt_t bits;
3741   struct gmp_div_inverse bi;
3742   size_t ndigits;
3743
3744   assert (base >= 2);
3745   assert (base <= 36);
3746
3747   un = GMP_ABS (u->_mp_size);
3748   if (un == 0)
3749     return 1;
3750
3751   up = u->_mp_d;
3752
3753   bits = (un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1]);
3754   switch (base)
3755     {
3756     case 2:
3757       return bits;
3758     case 4:
3759       return (bits + 1) / 2;
3760     case 8:
3761       return (bits + 2) / 3;
3762     case 16:
3763       return (bits + 3) / 4;
3764     case 32:
3765       return (bits + 4) / 5;
3766       /* FIXME: Do something more clever for the common case of base
3767          10. */
3768     }
3769
3770   tp = gmp_xalloc_limbs (un);
3771   mpn_copyi (tp, up, un);
3772   mpn_div_qr_1_invert (&bi, base);
3773
3774   for (ndigits = 0; un > 0; ndigits++)
3775     {
3776       mpn_div_qr_1_preinv (tp, tp, un, &bi);
3777       un -= (tp[un-1] == 0);
3778     }
3779   gmp_free (tp);
3780   return ndigits;
3781 }
3782
3783 char *
3784 mpz_get_str (char *sp, int base, const mpz_t u)
3785 {
3786   unsigned bits;
3787   const char *digits;
3788   mp_size_t un;
3789   size_t i, sn;
3790
3791   if (base >= 0)
3792     {
3793       digits = "0123456789abcdefghijklmnopqrstuvwxyz";
3794     }
3795   else
3796     {
3797       base = -base;
3798       digits = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
3799     }
3800   if (base <= 1)
3801     base = 10;
3802   if (base > 36)
3803     return NULL;
3804
3805   sn = 1 + mpz_sizeinbase (u, base);
3806   if (!sp)
3807     sp = gmp_xalloc (1 + sn);
3808
3809   un = GMP_ABS (u->_mp_size);
3810
3811   if (un == 0)
3812     {
3813       sp[0] = '0';
3814       sp[1] = '\0';
3815       return sp;
3816     }
3817
3818   i = 0;
3819
3820   if (u->_mp_size < 0)
3821     sp[i++] = '-';
3822
3823   bits = mpn_base_power_of_two_p (base);
3824
3825   if (bits)
3826     /* Not modified in this case. */
3827     sn = i + mpn_get_str_bits ((unsigned char *) sp + i, bits, u->_mp_d, un);
3828   else
3829     {
3830       struct mpn_base_info info;
3831       mp_ptr tp;
3832
3833       mpn_get_base_info (&info, base);
3834       tp = gmp_xalloc_limbs (un);
3835       mpn_copyi (tp, u->_mp_d, un);
3836
3837       sn = i + mpn_get_str_other ((unsigned char *) sp + i, base, &info, tp, un);
3838       gmp_free (tp);
3839     }
3840
3841   for (; i < sn; i++)
3842     sp[i] = digits[(unsigned char) sp[i]];
3843
3844   sp[sn] = '\0';
3845   return sp;
3846 }
3847
3848 int
3849 mpz_set_str (mpz_t r, const char *sp, int base)
3850 {
3851   unsigned bits;
3852   mp_size_t rn, alloc;
3853   mp_ptr rp;
3854   size_t sn;
3855   size_t dn;
3856   int sign;
3857   unsigned char *dp;
3858
3859   assert (base == 0 || (base >= 2 && base <= 36));
3860
3861   while (isspace( (unsigned char) *sp))
3862     sp++;
3863
3864   if (*sp == '-')
3865     {
3866       sign = 1;
3867       sp++;
3868     }
3869   else
3870     sign = 0;
3871
3872   if (base == 0)
3873     {
3874       if (*sp == '0')
3875         {
3876           sp++;
3877           if (*sp == 'x' || *sp == 'X')
3878             {
3879               base = 16;
3880               sp++;
3881             }
3882           else if (*sp == 'b' || *sp == 'B')
3883             {
3884               base = 2;
3885               sp++;
3886             }
3887           else
3888             base = 8;
3889         }
3890       else
3891         base = 10;
3892     }
3893
3894   sn = strlen (sp);
3895   dp = gmp_xalloc (sn + (sn == 0));
3896
3897   for (dn = 0; *sp; sp++)
3898     {
3899       unsigned digit;
3900
3901       if (isspace ((unsigned char) *sp))
3902         continue;
3903       if (*sp >= '0' && *sp <= '9')
3904         digit = *sp - '0';
3905       else if (*sp >= 'a' && *sp <= 'z')
3906         digit = *sp - 'a' + 10;
3907       else if (*sp >= 'A' && *sp <= 'Z')
3908         digit = *sp - 'A' + 10;
3909       else
3910         digit = base; /* fail */
3911
3912       if (digit >= base)
3913         {
3914           gmp_free (dp);
3915           r->_mp_size = 0;
3916           return -1;
3917         }
3918
3919       dp[dn++] = digit;
3920     }
3921
3922   bits = mpn_base_power_of_two_p (base);
3923
3924   if (bits > 0)
3925     {
3926       alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
3927       rp = MPZ_REALLOC (r, alloc);
3928       rn = mpn_set_str_bits (rp, dp, dn, bits);
3929     }
3930   else
3931     {
3932       struct mpn_base_info info;
3933       mpn_get_base_info (&info, base);
3934       alloc = (sn + info.exp - 1) / info.exp;
3935       rp = MPZ_REALLOC (r, alloc);
3936       rn = mpn_set_str_other (rp, dp, dn, base, &info);
3937     }
3938   assert (rn <= alloc);
3939   gmp_free (dp);
3940
3941   r->_mp_size = sign ? - rn : rn;
3942
3943   return 0;
3944 }
3945
3946 int
3947 mpz_init_set_str (mpz_t r, const char *sp, int base)
3948 {
3949   mpz_init (r);
3950   return mpz_set_str (r, sp, base);
3951 }
3952
3953 size_t
3954 mpz_out_str (FILE *stream, int base, const mpz_t x)
3955 {
3956   char *str;
3957   size_t len;
3958
3959   str = mpz_get_str (NULL, base, x);
3960   len = strlen (str);
3961   len = fwrite (str, 1, len, stream);
3962   gmp_free (str);
3963   return len;
3964 }
3965
3966 \f
3967 static int
3968 gmp_detect_endian (void)
3969 {
3970   static const int i = 1;
3971   const unsigned char *p = (const unsigned char *) &i;
3972   if (*p == 1)
3973     /* Little endian */
3974     return -1;
3975   else
3976     /* Big endian */
3977     return 1;
3978 }
3979
3980 /* Import and export. Does not support nails. */
3981 void
3982 mpz_import (mpz_t r, size_t count, int order, size_t size, int endian,
3983             size_t nails, const void *src)
3984 {
3985   const unsigned char *p;
3986   ptrdiff_t word_step;
3987   mp_ptr rp;
3988   mp_size_t rn;
3989
3990   /* The current (partial) limb. */
3991   mp_limb_t limb;
3992   /* The number of bytes already copied to this limb (starting from
3993      the low end). */
3994   size_t bytes;
3995   /* The index where the limb should be stored, when completed. */
3996   mp_size_t i;
3997
3998   if (nails != 0)
3999     gmp_die ("mpz_import: Nails not supported.");
4000
4001   assert (order == 1 || order == -1);
4002   assert (endian >= -1 && endian <= 1);
4003
4004   if (endian == 0)
4005     endian = gmp_detect_endian ();
4006
4007   p = (unsigned char *) src;
4008
4009   word_step = (order != endian) ? 2 * size : 0;
4010
4011   /* Process bytes from the least significant end, so point p at the
4012      least significant word. */
4013   if (order == 1)
4014     {
4015       p += size * (count - 1);
4016       word_step = - word_step;
4017     }
4018
4019   /* And at least significant byte of that word. */
4020   if (endian == 1)
4021     p += (size - 1);
4022
4023   rn = (size * count + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
4024   rp = MPZ_REALLOC (r, rn);
4025
4026   for (limb = 0, bytes = 0, i = 0; count > 0; count--, p += word_step)
4027     {
4028       size_t j;
4029       for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4030         {
4031           limb |= (mp_limb_t) *p << (bytes++ * CHAR_BIT);
4032           if (bytes == sizeof(mp_limb_t))
4033             {
4034               rp[i++] = limb;
4035               bytes = 0;
4036               limb = 0;
4037             }
4038         }
4039     }
4040   if (bytes > 0)
4041     rp[i++] = limb;
4042   assert (i == rn);
4043
4044   r->_mp_size = mpn_normalized_size (rp, i);
4045 }
4046
4047 void *
4048 mpz_export (void *r, size_t *countp, int order, size_t size, int endian,
4049             size_t nails, const mpz_t u)
4050 {
4051   unsigned char *p;
4052   ptrdiff_t word_step;
4053   size_t count, k;
4054   mp_size_t un;
4055
4056   /* The current (partial) limb. */
4057   mp_limb_t limb;
4058   /* The number of bytes left to to in this limb. */
4059   size_t bytes;
4060   /* The index where the limb was read. */
4061   mp_size_t i;
4062
4063   if (nails != 0)
4064     gmp_die ("mpz_import: Nails not supported.");
4065
4066   assert (order == 1 || order == -1);
4067   assert (endian >= -1 && endian <= 1);
4068   assert (size > 0 || u->_mp_size == 0);
4069
4070   un = GMP_ABS (u->_mp_size);
4071   if (un == 0)
4072     {
4073       if (countp)
4074         *countp = 0;
4075       return r;
4076     }
4077
4078   /* Count bytes in top limb. */
4079   for (limb = u->_mp_d[un-1], k = 0; limb > 0; k++, limb >>= CHAR_BIT)
4080     ;
4081
4082   assert (k > 0);
4083
4084   count = (k + (un-1) * sizeof (mp_limb_t) + size - 1) / size;
4085
4086   if (!r)
4087     r = gmp_xalloc (count * size);
4088
4089   if (endian == 0)
4090     endian = gmp_detect_endian ();
4091
4092   p = (unsigned char *) r;
4093
4094   word_step = (order != endian) ? 2 * size : 0;
4095
4096   /* Process bytes from the least significant end, so point p at the
4097      least significant word. */
4098   if (order == 1)
4099     {
4100       p += size * (count - 1);
4101       word_step = - word_step;
4102     }
4103
4104   /* And at least significant byte of that word. */
4105   if (endian == 1)
4106     p += (size - 1);
4107
4108   for (bytes = 0, i = 0, k = 0; k < count; k++, p += word_step)
4109       {
4110         size_t j;
4111         for (j = 0; j < size; j++, p -= (ptrdiff_t) endian)
4112           {
4113             if (bytes == 0)
4114               {
4115                 if (i < un)
4116                   limb = u->_mp_d[i++];
4117                 bytes = sizeof (mp_limb_t);
4118               }
4119             *p = limb;
4120             limb >>= CHAR_BIT;
4121             bytes--;
4122           }
4123       }
4124   assert (i == un);
4125   assert (k == count);
4126
4127   if (countp)
4128     *countp = count;
4129
4130   return r;
4131 }