@@ -2,24 +2,33 @@
22
33 Contributed to the GNU project by Niels Möller
44
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.
5+Copyright 1991-1997, 1999-2016 Free Software Foundation, Inc.
86
97This file is part of the GNU MP Library.
108
119The 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.
10+it under the terms of either:
11+
12+ * the GNU Lesser General Public License as published by the Free
13+ Software Foundation; either version 3 of the License, or (at your
14+ option) any later version.
15+
16+or
17+
18+ * the GNU General Public License as published by the Free Software
19+ Foundation; either version 2 of the License, or (at your option) any
20+ later version.
21+
22+or both in parallel, as here.
1523
1624The GNU MP Library is distributed in the hope that it will be useful, but
1725WITHOUT 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.
26+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
27+for more details.
2028
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/. */
29+You should have received copies of the GNU General Public License and the
30+GNU Lesser General Public License along with the GNU MP Library. If not,
31+see https://www.gnu.org/licenses/. */
2332
2433/* NOTE: All functions in this file which are not declared in
2534 mini-gmp.h are internal, and are not intended to be compatible
@@ -60,8 +69,10 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
6069#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
6170#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
6271
72+#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
73+
6374#define gmp_assert_nocarry(x) do { \
64− mp_limb_t __cy = x; \
75+ mp_limb_t __cy = (x); \
6576 assert (__cy == 0); \
6677 } while (0)
6778
@@ -227,6 +238,8 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
227238 (y) = __mpz_srcptr_swap__tmp; \
228239 } while (0)
229240
241+const int mp_bits_per_limb = GMP_LIMB_BITS;
242+
230243
231244/* Memory allocation and other helper functions. */
232245static void
@@ -253,12 +266,12 @@ gmp_default_alloc (size_t size)
253266static void *
254267gmp_default_realloc (void *old, size_t old_size, size_t new_size)
255268{
256− mp_ptr p;
269+ void * p;
257270
258271 p = realloc (old, new_size);
259272
260273 if (!p)
261− gmp_die("gmp_default_realoc: Virtual memory exhausted.");
274+ gmp_die("gmp_default_realloc: Virtual memory exhausted.");
262275
263276 return p;
264277}
@@ -311,14 +324,14 @@ mp_set_memory_functions (void *(*alloc_func) (size_t),
311324static mp_ptr
312325gmp_xalloc_limbs (mp_size_t size)
313326{
314− return gmp_xalloc (size * sizeof (mp_limb_t));
327+ return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
315328}
316329
317330static mp_ptr
318331gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
319332{
320333 assert (size > 0);
321− return (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
334+ return (mp_ptr) (*gmp_reallocate_func) (old, 0, size * sizeof (mp_limb_t));
322335}
323336
324337
@@ -335,7 +348,7 @@ mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
335348void
336349mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
337350{
338− while (n-- > 0)
351+ while (--n >= 0)
339352 d[n] = s[n];
340353}
341354
@@ -362,12 +375,23 @@ mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
362375static mp_size_t
363376mpn_normalized_size (mp_srcptr xp, mp_size_t n)
364377{
365− for (; n > 0 && xp[n-1] == 0; n--)
366− ;
378+ while (n > 0 && xp[n-1] == 0)
379+ --n;
367380 return n;
368381}
369382
370−#define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0)
383+int
384+mpn_zero_p(mp_srcptr rp, mp_size_t n)
385+{
386+ return mpn_normalized_size (rp, n) == 0;
387+}
388+
389+void
390+mpn_zero (mp_ptr rp, mp_size_t n)
391+{
392+ while (--n >= 0)
393+ rp[n] = 0;
394+}
371395
372396mp_limb_t
373397mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
@@ -432,7 +456,7 @@ mpn_sub_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
432456 {
433457 mp_limb_t a = ap[i];
434458 /* Carry out */
435− mp_limb_t cy = a < b;;
459+ mp_limb_t cy = a < b;
436460 rp[i] = a - b;
437461 b = cy;
438462 }
@@ -558,17 +582,16 @@ mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
558582 way. */
559583
560584 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
561− rp += 1, vp += 1, vn -= 1;
562585
563586 /* Now accumulate the product of up[] and the next higher limb from
564587 vp[]. */
565588
566− while (vn >= 1)
589+ while (--vn >= 1)
567590 {
591+ rp += 1, vp += 1;
568592 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
569− rp += 1, vp += 1, vn -= 1;
570593 }
571− return rp[un - 1];
594+ return rp[un];
572595}
573596
574597void
@@ -588,7 +611,6 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
588611{
589612 mp_limb_t high_limb, low_limb;
590613 unsigned int tnc;
591− mp_size_t i;
592614 mp_limb_t retval;
593615
594616 assert (n >= 1);
@@ -603,7 +625,7 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
603625 retval = low_limb >> tnc;
604626 high_limb = (low_limb << cnt);
605627
606− for (i = n; --i != 0;)
628+ while (--n != 0)
607629 {
608630 low_limb = *--up;
609631 *--rp = high_limb | (low_limb >> tnc);
@@ -619,7 +641,6 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
619641{
620642 mp_limb_t high_limb, low_limb;
621643 unsigned int tnc;
622− mp_size_t i;
623644 mp_limb_t retval;
624645
625646 assert (n >= 1);
@@ -631,7 +652,7 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
631652 retval = (high_limb << tnc);
632653 low_limb = high_limb >> cnt;
633654
634− for (i = n; --i != 0;)
655+ while (--n != 0)
635656 {
636657 high_limb = *up++;
637658 *rp++ = low_limb | (high_limb << tnc);
@@ -682,24 +703,68 @@ mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
682703 i, ptr, i, GMP_LIMB_MAX);
683704}
684705
706+void
707+mpn_com (mp_ptr rp, mp_srcptr up, mp_size_t n)
708+{
709+ while (--n >= 0)
710+ *rp++ = ~ *up++;
711+}
712+
713+mp_limb_t
714+mpn_neg (mp_ptr rp, mp_srcptr up, mp_size_t n)
715+{
716+ while (*up == 0)
717+ {
718+ *rp = 0;
719+ if (!--n)
720+ return 0;
721+ ++up; ++rp;
722+ }
723+ *rp = - *up;
724+ mpn_com (++rp, ++up, --n);
725+ return 1;
726+}
727+
685728
686729/* MPN division interface. */
730+
731+/* The 3/2 inverse is defined as
732+
733+ m = floor( (B^3-1) / (B u1 + u0)) - B
734+*/
687735mp_limb_t
688736mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
689737{
690− mp_limb_t r, p, m;
691− unsigned ul, uh;
692− unsigned ql, qh;
738+ mp_limb_t r, p, m, ql;
739+ unsigned ul, uh, qh;
693740
694− /* First, do a 2/1 inverse. */
695− /* The inverse m is defined as floor( (B^2 - 1 - u1)/u1 ), so that 0 <
696− * B^2 - (B + m) u1 <= u1 */
697741 assert (u1 >= GMP_LIMB_HIGHBIT);
698742
743+ /* For notation, let b denote the half-limb base, so that B = b^2.
744+ Split u1 = b uh + ul. */
699745 ul = u1 & GMP_LLIMB_MASK;
700746 uh = u1 >> (GMP_LIMB_BITS / 2);
701747
748+ /* Approximation of the high half of quotient. Differs from the 2/1
749+ inverse of the half limb uh, since we have already subtracted
750+ u0. */
702751 qh = ~u1 / uh;
752+
753+ /* Adjust to get a half-limb 3/2 inverse, i.e., we want
754+
755+ qh' = floor( (b^3 - 1) / u) - b = floor ((b^3 - b u - 1) / u
756+ = floor( (b (~u) + b-1) / u),
757+
758+ and the remainder
759+
760+ r = b (~u) + b-1 - qh (b uh + ul)
761+ = b (~u - qh uh) + b-1 - qh ul
762+
763+ Subtraction of qh ul may underflow, which implies adjustments.
764+ But by normalization, 2 u >= B > qh ul, so we need to adjust by
765+ at most 2.
766+ */
767+
703768 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
704769
705770 p = (mp_limb_t) qh * ul;
@@ -717,11 +782,19 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
717782 }
718783 r -= p;
719784
720− /* Do a 3/2 division (with half limb size) */
785+ /* Low half of the quotient is
786+
787+ ql = floor ( (b r + b-1) / u1).
788+
789+ This is a 3/2 division (on half-limbs), for which qh is a
790+ suitable inverse. */
791+
721792 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
793+ /* Unlike full-limb 3/2, we can add 1 without overflow. For this to
794+ work, it is essential that ql is a full mp_limb_t. */
722795 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
723796
724− /* By the 3/2 method, we don't need the high half limb. */
797+ /* By the 3/2 trick, we don't need the high half limb. */
725798 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
726799
727800 if (r >= (p << (GMP_LIMB_BITS / 2)))
@@ -736,6 +809,8 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
736809 r -= u1;
737810 }
738811
812+ /* Now m is the 2/1 invers of u1. If u0 > 0, adjust it to become a
813+ 3/2 inverse. */
739814 if (u0 > 0)
740815 {
741816 mp_limb_t th, tl;
@@ -847,7 +922,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
847922
848923 if (inv->shift > 0)
849924 {
850− tp = gmp_xalloc_limbs (nn);
925+ tp = qp ? qp : gmp_xalloc_limbs (nn);
851926 r = mpn_lshift (tp, np, nn, inv->shift);
852927 np = tp;
853928 }
@@ -856,7 +931,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
856931
857932 d = inv->d1;
858933 di = inv->di;
859− while (nn-- > 0)
934+ while (--nn >= 0)
860935 {
861936 mp_limb_t q;
862937
@@ -864,7 +939,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
864939 if (qp)
865940 qp[nn] = q;
866941 }
867− if (inv->shift > 0)
942+ if (inv->shift > 0 && !qp)
868943 gmp_free (tp);
869944
870945 return r >> inv->shift;
@@ -1140,7 +1215,7 @@ mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
11401215 unsigned char mask;
11411216 size_t sn, j;
11421217 mp_size_t i;
1143− int shift;
1218+ unsigned shift;
11441219
11451220 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
11461221 + bits - 1) / bits;
@@ -1281,6 +1356,8 @@ mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
12811356 return rn;
12821357}
12831358
1359+/* Result is usually normalized, except for all-zero input, in which
1360+ case a single zero limb is written at *RP, and 1 is returned. */
12841361static mp_size_t
12851362mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
12861363 mp_limb_t b, const struct mpn_base_info *info)
@@ -1290,16 +1367,18 @@ mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
12901367 unsigned k;
12911368 size_t j;
12921369
1370+ assert (sn > 0);
1371+
12931372 k = 1 + (sn - 1) % info->exp;
12941373
12951374 j = 0;
12961375 w = sp[j++];
1297− for (; --k > 0; )
1376+ while (--k != 0)
12981377 w = w * b + sp[j++];
12991378
13001379 rp[0] = w;
13011380
1302− for (rn = (w > 0); j < sn;)
1381+ for (rn = 1; j < sn;)
13031382 {
13041383 mp_limb_t cy;
13051384
@@ -1342,9 +1421,11 @@ mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
13421421void
13431422mpz_init (mpz_t r)
13441423{
1345− r->_mp_alloc = 1;
1424+ static const mp_limb_t dummy_limb = 0xc1a0;
1425+
1426+ r->_mp_alloc = 0;
13461427 r->_mp_size = 0;
1347− r->_mp_d = gmp_xalloc_limbs (1);
1428+ r->_mp_d = (mp_ptr) &dummy_limb;
13481429}
13491430
13501431/* The utility of this function is a bit limited, since many functions
@@ -1365,15 +1446,19 @@ mpz_init2 (mpz_t r, mp_bitcnt_t bits)
13651446void
13661447mpz_clear (mpz_t r)
13671448{
1368− gmp_free (r->_mp_d);
1449+ if (r->_mp_alloc)
1450+ gmp_free (r->_mp_d);
13691451}
13701452
1371−static void *
1453+static mp_ptr
13721454mpz_realloc (mpz_t r, mp_size_t size)
13731455{
13741456 size = GMP_MAX (size, 1);
13751457
1376− r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1458+ if (r->_mp_alloc)
1459+ r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1460+ else
1461+ r->_mp_d = gmp_xalloc_limbs (size);
13771462 r->_mp_alloc = size;
13781463
13791464 if (GMP_ABS (r->_mp_size) > size)
@@ -1396,7 +1481,7 @@ mpz_set_si (mpz_t r, signed long int x)
13961481 else /* (x < 0) */
13971482 {
13981483 r->_mp_size = -1;
1399− r->_mp_d[0] = GMP_NEG_CAST (unsigned long int, x);
1484+ MPZ_REALLOC (r, 1)[0] = GMP_NEG_CAST (unsigned long int, x);
14001485 }
14011486}
14021487
@@ -1406,7 +1491,7 @@ mpz_set_ui (mpz_t r, unsigned long int x)
14061491 if (x > 0)
14071492 {
14081493 r->_mp_size = 1;
1409− r->_mp_d[0] = x;
1494+ MPZ_REALLOC (r, 1)[0] = x;
14101495 }
14111496 else
14121497 r->_mp_size = 0;
@@ -1455,14 +1540,12 @@ mpz_fits_slong_p (const mpz_t u)
14551540{
14561541 mp_size_t us = u->_mp_size;
14571542
1458− if (us == 0)
1459− return 1;
1460− else if (us == 1)
1543+ if (us == 1)
14611544 return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
14621545 else if (us == -1)
14631546 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
14641547 else
1465− return 0;
1548+ return (us == 0);
14661549}
14671550
14681551int
@@ -1476,14 +1559,11 @@ mpz_fits_ulong_p (const mpz_t u)
14761559long int
14771560mpz_get_si (const mpz_t u)
14781561{
1479− mp_size_t us = u->_mp_size;
1480−
1481− if (us > 0)
1482− return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT);
1483− else if (us < 0)
1484− return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT);
1562+ if (u->_mp_size < 0)
1563+ /* This expression is necessary to properly handle 0x80000000 */
1564+ return -1 - (long) ((u->_mp_d[0] - 1) & ~GMP_LIMB_HIGHBIT);
14851565 else
1486− return 0;
1566+ return (long) (mpz_get_ui (u) & ~GMP_LIMB_HIGHBIT);
14871567}
14881568
14891569unsigned long int
@@ -1507,6 +1587,48 @@ mpz_getlimbn (const mpz_t u, mp_size_t n)
15071587 return 0;
15081588}
15091589
1590+void
1591+mpz_realloc2 (mpz_t x, mp_bitcnt_t n)
1592+{
1593+ mpz_realloc (x, 1 + (n - (n != 0)) / GMP_LIMB_BITS);
1594+}
1595+
1596+mp_srcptr
1597+mpz_limbs_read (mpz_srcptr x)
1598+{
1599+ return x->_mp_d;
1600+}
1601+
1602+mp_ptr
1603+mpz_limbs_modify (mpz_t x, mp_size_t n)
1604+{
1605+ assert (n > 0);
1606+ return MPZ_REALLOC (x, n);
1607+}
1608+
1609+mp_ptr
1610+mpz_limbs_write (mpz_t x, mp_size_t n)
1611+{
1612+ return mpz_limbs_modify (x, n);
1613+}
1614+
1615+void
1616+mpz_limbs_finish (mpz_t x, mp_size_t xs)
1617+{
1618+ mp_size_t xn;
1619+ xn = mpn_normalized_size (x->_mp_d, GMP_ABS (xs));
1620+ x->_mp_size = xs < 0 ? -xn : xn;
1621+}
1622+
1623+mpz_srcptr
1624+mpz_roinit_n (mpz_t x, mp_srcptr xp, mp_size_t xs)
1625+{
1626+ x->_mp_alloc = 0;
1627+ x->_mp_d = (mp_ptr) xp;
1628+ mpz_limbs_finish (x, xs);
1629+ return x;
1630+}
1631+
15101632
15111633/* Conversions and comparison to double. */
15121634void
@@ -1654,9 +1776,7 @@ mpz_cmp_d (const mpz_t x, double d)
16541776int
16551777mpz_sgn (const mpz_t u)
16561778{
1657− mp_size_t usize = u->_mp_size;
1658−
1659− return (usize > 0) - (usize < 0);
1779+ return GMP_CMP (u->_mp_size, 0);
16601780}
16611781
16621782int
@@ -1671,13 +1791,7 @@ mpz_cmp_si (const mpz_t u, long v)
16711791 else if (usize >= 0)
16721792 return 1;
16731793 else /* usize == -1 */
1674− {
1675− mp_limb_t ul = u->_mp_d[0];
1676− if ((mp_limb_t)GMP_NEG_CAST (unsigned long int, v) < ul)
1677− return -1;
1678− else
1679− return (mp_limb_t)GMP_NEG_CAST (unsigned long int, v) > ul;
1680− }
1794+ return GMP_CMP (GMP_NEG_CAST (mp_limb_t, v), u->_mp_d[0]);
16811795}
16821796
16831797int
@@ -1690,10 +1804,7 @@ mpz_cmp_ui (const mpz_t u, unsigned long v)
16901804 else if (usize < 0)
16911805 return -1;
16921806 else
1693− {
1694− mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
1695− return (ul > v) - (ul < v);
1696− }
1807+ return GMP_CMP (mpz_get_ui (u), v);
16971808}
16981809
16991810int
@@ -1713,15 +1824,10 @@ mpz_cmp (const mpz_t a, const mpz_t b)
17131824int
17141825mpz_cmpabs_ui (const mpz_t u, unsigned long v)
17151826{
1716− mp_size_t un = GMP_ABS (u->_mp_size);
1717− mp_limb_t ul;
1718−
1719− if (un > 1)
1827+ if (GMP_ABS (u->_mp_size) > 1)
17201828 return 1;
1721−
1722− ul = (un == 1) ? u->_mp_d[0] : 0;
1723−
1724− return (ul > v) - (ul < v);
1829+ else
1830+ return GMP_CMP (mpz_get_ui (u), v);
17251831}
17261832
17271833int
@@ -1734,18 +1840,14 @@ mpz_cmpabs (const mpz_t u, const mpz_t v)
17341840void
17351841mpz_abs (mpz_t r, const mpz_t u)
17361842{
1737− if (r != u)
1738− mpz_set (r, u);
1739−
1843+ mpz_set (r, u);
17401844 r->_mp_size = GMP_ABS (r->_mp_size);
17411845}
17421846
17431847void
17441848mpz_neg (mpz_t r, const mpz_t u)
17451849{
1746− if (r != u)
1747− mpz_set (r, u);
1748−
1850+ mpz_set (r, u);
17491851 r->_mp_size = -r->_mp_size;
17501852}
17511853
@@ -1771,7 +1873,7 @@ mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
17711873 an = GMP_ABS (a->_mp_size);
17721874 if (an == 0)
17731875 {
1774− r->_mp_d[0] = b;
1876+ MPZ_REALLOC (r, 1)[0] = b;
17751877 return b > 0;
17761878 }
17771879
@@ -1790,14 +1892,15 @@ static mp_size_t
17901892mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
17911893{
17921894 mp_size_t an = GMP_ABS (a->_mp_size);
1793− mp_ptr rp = MPZ_REALLOC (r, an);
1895+ mp_ptr rp;
17941896
17951897 if (an == 0)
17961898 {
1797− rp[0] = b;
1899+ MPZ_REALLOC (r, 1)[0] = b;
17981900 return -(b > 0);
17991901 }
1800− else if (an == 1 && a->_mp_d[0] < b)
1902+ rp = MPZ_REALLOC (r, an);
1903+ if (an == 1 && a->_mp_d[0] < b)
18011904 {
18021905 rp[0] = b - a->_mp_d[0];
18031906 return -1;
@@ -2015,12 +2118,51 @@ mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
20152118 else
20162119 mpn_copyd (rp + limbs, u->_mp_d, un);
20172120
2018− while (limbs > 0)
2019− rp[--limbs] = 0;
2121+ mpn_zero (rp, limbs);
20202122
20212123 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
20222124}
20232125
2126+void
2127+mpz_addmul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2128+{
2129+ mpz_t t;
2130+ mpz_init (t);
2131+ mpz_mul_ui (t, u, v);
2132+ mpz_add (r, r, t);
2133+ mpz_clear (t);
2134+}
2135+
2136+void
2137+mpz_submul_ui (mpz_t r, const mpz_t u, unsigned long int v)
2138+{
2139+ mpz_t t;
2140+ mpz_init (t);
2141+ mpz_mul_ui (t, u, v);
2142+ mpz_sub (r, r, t);
2143+ mpz_clear (t);
2144+}
2145+
2146+void
2147+mpz_addmul (mpz_t r, const mpz_t u, const mpz_t v)
2148+{
2149+ mpz_t t;
2150+ mpz_init (t);
2151+ mpz_mul (t, u, v);
2152+ mpz_add (r, r, t);
2153+ mpz_clear (t);
2154+}
2155+
2156+void
2157+mpz_submul (mpz_t r, const mpz_t u, const mpz_t v)
2158+{
2159+ mpz_t t;
2160+ mpz_init (t);
2161+ mpz_mul (t, u, v);
2162+ mpz_sub (r, r, t);
2163+ mpz_clear (t);
2164+}
2165+
20242166
20252167/* MPZ division */
20262168enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
@@ -2085,8 +2227,7 @@ mpz_div_qr (mpz_t q, mpz_t r,
20852227 mp_size_t qn, rn;
20862228 mpz_t tq, tr;
20872229
2088− mpz_init (tr);
2089− mpz_set (tr, n);
2230+ mpz_init_set (tr, n);
20902231 np = tr->_mp_d;
20912232
20922233 qn = nn - dn + 1;
@@ -2230,7 +2371,6 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
22302371
22312372 if (qn <= 0)
22322373 qn = 0;
2233−
22342374 else
22352375 {
22362376 qp = MPZ_REALLOC (q, qn);
@@ -2284,16 +2424,9 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
22842424 {
22852425 /* Have to negate and sign extend. */
22862426 mp_size_t i;
2287− mp_limb_t cy;
22882427
2289− for (cy = 1, i = 0; i < un; i++)
2290− {
2291− mp_limb_t s = ~u->_mp_d[i] + cy;
2292− cy = s < cy;
2293− rp[i] = s;
2294− }
2295− assert (cy == 0);
2296− for (; i < rn - 1; i++)
2428+ gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
2429+ for (i = un; i < rn - 1; i++)
22972430 rp[i] = GMP_LIMB_MAX;
22982431
22992432 rp[rn-1] = mask;
@@ -2318,23 +2451,13 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
23182451 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
23192452 {
23202453 /* If r != 0, compute 2^{bit_count} - r. */
2321− mp_size_t i;
2454+ mpn_neg (rp, rp, rn);
23222455
2323− for (i = 0; i < rn && rp[i] == 0; i++)
2324− ;
2325− if (i < rn)
2326− {
2327− /* r > 0, need to flip sign. */
2328− rp[i] = ~rp[i] + 1;
2329− while (++i < rn)
2330− rp[i] = ~rp[i];
2331−
2332− rp[rn-1] &= mask;
2456+ rp[rn-1] &= mask;
23332457
2334− /* us is not used for anything else, so we can modify it
2335− here to indicate flipped sign. */
2336− us = -us;
2337− }
2458+ /* us is not used for anything else, so we can modify it
2459+ here to indicate flipped sign. */
2460+ us = -us;
23382461 }
23392462 }
23402463 rn = mpn_normalized_size (rp, rn);
@@ -2389,6 +2512,24 @@ mpz_divisible_p (const mpz_t n, const mpz_t d)
23892512 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
23902513}
23912514
2515+int
2516+mpz_congruent_p (const mpz_t a, const mpz_t b, const mpz_t m)
2517+{
2518+ mpz_t t;
2519+ int res;
2520+
2521+ /* a == b (mod 0) iff a == b */
2522+ if (mpz_sgn (m) == 0)
2523+ return (mpz_cmp (a, b) == 0);
2524+
2525+ mpz_init (t);
2526+ mpz_sub (t, a, b);
2527+ res = mpz_divisible_p (t, m);
2528+ mpz_clear (t);
2529+
2530+ return res;
2531+}
2532+
23922533static unsigned long
23932534mpz_div_qr_ui (mpz_t q, mpz_t r,
23942535 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
@@ -2431,7 +2572,7 @@ mpz_div_qr_ui (mpz_t q, mpz_t r,
24312572
24322573 if (r)
24332574 {
2434− r->_mp_d[0] = rl;
2575+ MPZ_REALLOC (r, 1)[0] = rl;
24352576 r->_mp_size = rs;
24362577 }
24372578 if (q)
@@ -2955,9 +3096,7 @@ void
29553096mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
29563097{
29573098 mpz_t b;
2958− mpz_init_set_ui (b, blimb);
2959− mpz_pow_ui (r, b, e);
2960− mpz_clear (b);
3099+ mpz_pow_ui (r, mpz_roinit_n (b, &blimb, 1), e);
29613100}
29623101
29633102void
@@ -3029,7 +3168,7 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
30293168 }
30303169 mpz_init_set_ui (tr, 1);
30313170
3032− while (en-- > 0)
3171+ while (--en >= 0)
30333172 {
30343173 mp_limb_t w = e->_mp_d[en];
30353174 mp_limb_t bit;
@@ -3069,9 +3208,7 @@ void
30693208mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
30703209{
30713210 mpz_t e;
3072− mpz_init_set_ui (e, elimb);
3073− mpz_powm (r, b, e, m);
3074− mpz_clear (e);
3211+ mpz_powm (r, b, mpz_roinit_n (e, &elimb, 1), m);
30753212}
30763213
30773214/* x=trunc(y^(1/z)), r=y-x^z */
@@ -3088,19 +3225,16 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
30883225 gmp_die ("mpz_rootrem: Zeroth root.");
30893226
30903227 if (mpz_cmpabs_ui (y, 1) <= 0) {
3091− mpz_set (x, y);
3228+ if (x)
3229+ mpz_set (x, y);
30923230 if (r)
30933231 r->_mp_size = 0;
30943232 return;
30953233 }
30963234
30973235 mpz_init (u);
3098− {
3099− mp_bitcnt_t tb;
3100− tb = mpz_sizeinbase (y, 2) / z + 1;
3101− mpz_init2 (t, tb);
3102− mpz_setbit (t, tb);
3103− }
3236+ mpz_init (t);
3237+ mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
31043238
31053239 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
31063240 do {
@@ -3132,7 +3266,8 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
31323266 mpz_pow_ui (t, u, z);
31333267 mpz_sub (r, y, t);
31343268 }
3135− mpz_swap (x, u);
3269+ if (x)
3270+ mpz_swap (x, u);
31363271 mpz_clear (u);
31373272 mpz_clear (t);
31383273}
@@ -3164,6 +3299,47 @@ mpz_sqrt (mpz_t s, const mpz_t u)
31643299 mpz_rootrem (s, NULL, u, 2);
31653300}
31663301
3302+int
3303+mpz_perfect_square_p (const mpz_t u)
3304+{
3305+ if (u->_mp_size <= 0)
3306+ return (u->_mp_size == 0);
3307+ else
3308+ return mpz_root (NULL, u, 2);
3309+}
3310+
3311+int
3312+mpn_perfect_square_p (mp_srcptr p, mp_size_t n)
3313+{
3314+ mpz_t t;
3315+
3316+ assert (n > 0);
3317+ assert (p [n-1] != 0);
3318+ return mpz_root (NULL, mpz_roinit_n (t, p, n), 2);
3319+}
3320+
3321+mp_size_t
3322+mpn_sqrtrem (mp_ptr sp, mp_ptr rp, mp_srcptr p, mp_size_t n)
3323+{
3324+ mpz_t s, r, u;
3325+ mp_size_t res;
3326+
3327+ assert (n > 0);
3328+ assert (p [n-1] != 0);
3329+
3330+ mpz_init (r);
3331+ mpz_init (s);
3332+ mpz_rootrem (s, r, mpz_roinit_n (u, p, n), 2);
3333+
3334+ assert (s->_mp_size == (n+1)/2);
3335+ mpn_copyd (sp, s->_mp_d, s->_mp_size);
3336+ mpz_clear (s);
3337+ res = r->_mp_size;
3338+ if (rp)
3339+ mpn_copyd (rp, r->_mp_d, res);
3340+ mpz_clear (r);
3341+ return res;
3342+}
31673343
31683344/* Combinatorics */
31693345
@@ -3171,7 +3347,7 @@ void
31713347mpz_fac_ui (mpz_t x, unsigned long n)
31723348{
31733349 mpz_set_ui (x, n + (n == 0));
3174− for (;n > 2;)
3350+ while (n > 2)
31753351 mpz_mul_ui (x, x, --n);
31763352}
31773353
@@ -3195,6 +3371,102 @@ mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
31953371 mpz_clear (t);
31963372}
31973373
3374+
3375+/* Primality testing */
3376+static int
3377+gmp_millerrabin (const mpz_t n, const mpz_t nm1, mpz_t y,
3378+ const mpz_t q, mp_bitcnt_t k)
3379+{
3380+ assert (k > 0);
3381+
3382+ /* Caller must initialize y to the base. */
3383+ mpz_powm (y, y, q, n);
3384+
3385+ if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
3386+ return 1;
3387+
3388+ while (--k > 0)
3389+ {
3390+ mpz_powm_ui (y, y, 2, n);
3391+ if (mpz_cmp (y, nm1) == 0)
3392+ return 1;
3393+ /* y == 1 means that the previous y was a non-trivial square root
3394+ of 1 (mod n). y == 0 means that n is a power of the base.
3395+ In either case, n is not prime. */
3396+ if (mpz_cmp_ui (y, 1) <= 0)
3397+ return 0;
3398+ }
3399+ return 0;
3400+}
3401+
3402+/* This product is 0xc0cfd797, and fits in 32 bits. */
3403+#define GMP_PRIME_PRODUCT \
3404+ (3UL*5UL*7UL*11UL*13UL*17UL*19UL*23UL*29UL)
3405+
3406+/* Bit (p+1)/2 is set, for each odd prime <= 61 */
3407+#define GMP_PRIME_MASK 0xc96996dcUL
3408+
3409+int
3410+mpz_probab_prime_p (const mpz_t n, int reps)
3411+{
3412+ mpz_t nm1;
3413+ mpz_t q;
3414+ mpz_t y;
3415+ mp_bitcnt_t k;
3416+ int is_prime;
3417+ int j;
3418+
3419+ /* Note that we use the absolute value of n only, for compatibility
3420+ with the real GMP. */
3421+ if (mpz_even_p (n))
3422+ return (mpz_cmpabs_ui (n, 2) == 0) ? 2 : 0;
3423+
3424+ /* Above test excludes n == 0 */
3425+ assert (n->_mp_size != 0);
3426+
3427+ if (mpz_cmpabs_ui (n, 64) < 0)
3428+ return (GMP_PRIME_MASK >> (n->_mp_d[0] >> 1)) & 2;
3429+
3430+ if (mpz_gcd_ui (NULL, n, GMP_PRIME_PRODUCT) != 1)
3431+ return 0;
3432+
3433+ /* All prime factors are >= 31. */
3434+ if (mpz_cmpabs_ui (n, 31*31) < 0)
3435+ return 2;
3436+
3437+ /* Use Miller-Rabin, with a deterministic sequence of bases, a[j] =
3438+ j^2 + j + 41 using Euler's polynomial. We potentially stop early,
3439+ if a[j] >= n - 1. Since n >= 31*31, this can happen only if reps >
3440+ 30 (a[30] == 971 > 31*31 == 961). */
3441+
3442+ mpz_init (nm1);
3443+ mpz_init (q);
3444+ mpz_init (y);
3445+
3446+ /* Find q and k, where q is odd and n = 1 + 2**k * q. */
3447+ nm1->_mp_size = mpz_abs_sub_ui (nm1, n, 1);
3448+ k = mpz_scan1 (nm1, 0);
3449+ mpz_tdiv_q_2exp (q, nm1, k);
3450+
3451+ for (j = 0, is_prime = 1; is_prime & (j < reps); j++)
3452+ {
3453+ mpz_set_ui (y, (unsigned long) j*j+j+41);
3454+ if (mpz_cmp (y, nm1) >= 0)
3455+ {
3456+ /* Don't try any further bases. This "early" break does not affect
3457+ the result for any reasonable reps value (<=5000 was tested) */
3458+ assert (j >= 30);
3459+ break;
3460+ }
3461+ is_prime = gmp_millerrabin (n, nm1, y, q, k);
3462+ }
3463+ mpz_clear (nm1);
3464+ mpz_clear (q);
3465+ mpz_clear (y);
3466+
3467+ return is_prime;
3468+}
3469+
31983470
31993471/* Logical operations and bit manipulation. */
32003472
@@ -3246,7 +3518,7 @@ mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
32463518 must be complemented. */
32473519 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
32483520 return bit ^ 1;
3249− while (limb_index-- > 0)
3521+ while (--limb_index >= 0)
32503522 if (d->_mp_d[limb_index] > 0)
32513523 return bit ^ 1;
32523524 }
@@ -3311,7 +3583,7 @@ mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
33113583
33123584 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
33133585 dn - limb_index, bit));
3314− dn -= (dp[dn-1] == 0);
3586+ dn = mpn_normalized_size (dp, dn);
33153587 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
33163588}
33173589
@@ -3389,7 +3661,7 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
33893661 /* If the smaller input is positive, higher limbs don't matter. */
33903662 rn = vx ? un : vn;
33913663
3392− rp = MPZ_REALLOC (r, rn + rc);
3664+ rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
33933665
33943666 up = u->_mp_d;
33953667 vp = v->_mp_d;
@@ -3462,7 +3734,7 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
34623734 don't matter. */
34633735 rn = vx ? vn : un;
34643736
3465− rp = MPZ_REALLOC (r, rn + rc);
3737+ rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
34663738
34673739 up = u->_mp_d;
34683740 vp = v->_mp_d;
@@ -3531,7 +3803,7 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
35313803 vx = -vc;
35323804 rx = -rc;
35333805
3534− rp = MPZ_REALLOC (r, un + rc);
3806+ rp = MPZ_REALLOC (r, un + (mp_size_t) rc);
35353807
35363808 up = u->_mp_d;
35373809 vp = v->_mp_d;
@@ -3587,20 +3859,28 @@ gmp_popcount_limb (mp_limb_t x)
35873859}
35883860
35893861mp_bitcnt_t
3590−mpz_popcount (const mpz_t u)
3862+mpn_popcount (mp_srcptr p, mp_size_t n)
35913863{
3592− mp_size_t un, i;
3864+ mp_size_t i;
35933865 mp_bitcnt_t c;
35943866
3867+ for (c = 0, i = 0; i < n; i++)
3868+ c += gmp_popcount_limb (p[i]);
3869+
3870+ return c;
3871+}
3872+
3873+mp_bitcnt_t
3874+mpz_popcount (const mpz_t u)
3875+{
3876+ mp_size_t un;
3877+
35953878 un = u->_mp_size;
35963879
35973880 if (un < 0)
35983881 return ~(mp_bitcnt_t) 0;
35993882
3600− for (c = 0, i = 0; i < un; i++)
3601− c += gmp_popcount_limb (u->_mp_d[i]);
3602−
3603− return c;
3883+ return mpn_popcount (u->_mp_d, un);
36043884}
36053885
36063886mp_bitcnt_t
@@ -3799,7 +4079,7 @@ mpz_get_str (char *sp, int base, const mpz_t u)
37994079
38004080 sn = 1 + mpz_sizeinbase (u, base);
38014081 if (!sp)
3802− sp = gmp_xalloc (1 + sn);
4082+ sp = (char *) gmp_xalloc (1 + sn);
38034083
38044084 un = GMP_ABS (u->_mp_size);
38054085
@@ -3846,7 +4126,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
38464126 unsigned bits;
38474127 mp_size_t rn, alloc;
38484128 mp_ptr rp;
3849− size_t sn;
4129+ size_t dn;
38504130 int sign;
38514131 unsigned char *dp;
38524132
@@ -3860,18 +4140,17 @@ mpz_set_str (mpz_t r, const char *sp, int base)
38604140
38614141 if (base == 0)
38624142 {
3863− if (*sp == '0')
4143+ if (sp[0] == '0')
38644144 {
3865− sp++;
3866− if (*sp == 'x' || *sp == 'X')
4145+ if (sp[1] == 'x' || sp[1] == 'X')
38674146 {
38684147 base = 16;
3869− sp++;
4148+ sp += 2;
38704149 }
3871− else if (*sp == 'b' || *sp == 'B')
4150+ else if (sp[1] == 'b' || sp[1] == 'B')
38724151 {
38734152 base = 2;
3874− sp++;
4153+ sp += 2;
38754154 }
38764155 else
38774156 base = 8;
@@ -3880,16 +4159,20 @@ mpz_set_str (mpz_t r, const char *sp, int base)
38804159 base = 10;
38814160 }
38824161
3883− sn = strlen (sp);
3884− dp = gmp_xalloc (sn + (sn == 0));
4162+ if (!*sp)
4163+ {
4164+ r->_mp_size = 0;
4165+ return -1;
4166+ }
4167+ dp = (unsigned char *) gmp_xalloc (strlen (sp));
38854168
3886− for (sn = 0; *sp; sp++)
4169+ for (dn = 0; *sp; sp++)
38874170 {
38884171 unsigned digit;
38894172
38904173 if (isspace ((unsigned char) *sp))
38914174 continue;
3892− if (*sp >= '0' && *sp <= '9')
4175+ else if (*sp >= '0' && *sp <= '9')
38934176 digit = *sp - '0';
38944177 else if (*sp >= 'a' && *sp <= 'z')
38954178 digit = *sp - 'a' + 10;
@@ -3898,31 +4181,40 @@ mpz_set_str (mpz_t r, const char *sp, int base)
38984181 else
38994182 digit = base; /* fail */
39004183
3901− if (digit >= base)
4184+ if (digit >= (unsigned) base)
39024185 {
39034186 gmp_free (dp);
39044187 r->_mp_size = 0;
39054188 return -1;
39064189 }
39074190
3908− dp[sn++] = digit;
4191+ dp[dn++] = digit;
39094192 }
39104193
4194+ if (!dn)
4195+ {
4196+ gmp_free (dp);
4197+ r->_mp_size = 0;
4198+ return -1;
4199+ }
39114200 bits = mpn_base_power_of_two_p (base);
39124201
39134202 if (bits > 0)
39144203 {
3915− alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
4204+ alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
39164205 rp = MPZ_REALLOC (r, alloc);
3917− rn = mpn_set_str_bits (rp, dp, sn, bits);
4206+ rn = mpn_set_str_bits (rp, dp, dn, bits);
39184207 }
39194208 else
39204209 {
39214210 struct mpn_base_info info;
39224211 mpn_get_base_info (&info, base);
3923− alloc = (sn + info.exp - 1) / info.exp;
4212+ alloc = (dn + info.exp - 1) / info.exp;
39244213 rp = MPZ_REALLOC (r, alloc);
3925− rn = mpn_set_str_other (rp, dp, sn, base, &info);
4214+ rn = mpn_set_str_other (rp, dp, dn, base, &info);
4215+ /* Normalization, needed for all-zero input. */
4216+ assert (rn > 0);
4217+ rn -= rp[rn-1] == 0;
39264218 }
39274219 assert (rn <= alloc);
39284220 gmp_free (dp);