git.iko.soy/Idris-dev
on master help sign in
Fix looped temporary allocations in mini-GMP (#4324)
3be1766yanok
rts/mini-gmp.c+459−167
compare with master
unifiedsplit
@@ -2,24 +2,33 @@
2
2
3 Contributed to the GNU project by Niels Möller
3 Contributed to the GNU project by Niels Möller
4
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.
5+Copyright 1991-1997, 1999-2016 Free Software Foundation, Inc.
8
6
9This file is part of the GNU MP Library.
7This file is part of the GNU MP Library.
10
8
11The GNU MP Library is free software; you can redistribute it and/or modify
9The GNU MP Library is free software; you can redistribute it and/or modify
10+it under the terms of either:
11+
12−it under the terms of the GNU Lesser General Public License as published by
12+ * the GNU Lesser General Public License as published by the Free
13−the Free Software Foundation; either version 3 of the License, or (at your
13+ Software Foundation; either version 3 of the License, or (at your
14−option) any later version.
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.
15
23
16The GNU MP Library is distributed in the hope that it will be useful, but
24The GNU MP Library is distributed in the hope that it will be useful, but
17WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18−or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
26+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19−License for more details.
27+for more details.
20
28
21−You should have received a copy of the GNU Lesser General Public License
29+You should have received copies of the GNU General Public License and the
22−along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
30+GNU Lesser General Public License along with the GNU MP Library. If not,
31+see https://www.gnu.org/licenses/. */
23
32
24/* NOTE: All functions in this file which are not declared in
33/* NOTE: All functions in this file which are not declared in
25 mini-gmp.h are internal, and are not intended to be compatible
34 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/. */
60#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
69#define GMP_MIN(a, b) ((a) < (b) ? (a) : (b))
61#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
70#define GMP_MAX(a, b) ((a) > (b) ? (a) : (b))
62
71
72+#define GMP_CMP(a,b) (((a) > (b)) - ((a) < (b)))
73+
63#define gmp_assert_nocarry(x) do { \
74#define gmp_assert_nocarry(x) do { \
64− mp_limb_t __cy = x; \
75+ mp_limb_t __cy = (x); \
65 assert (__cy == 0); \
76 assert (__cy == 0); \
66 } while (0)
77 } while (0)
67
78
@@ -227,6 +238,8 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
227 (y) = __mpz_srcptr_swap__tmp; \
238 (y) = __mpz_srcptr_swap__tmp; \
228 } while (0)
239 } while (0)
229
240
241+const int mp_bits_per_limb = GMP_LIMB_BITS;
242+
230
243
231/* Memory allocation and other helper functions. */
244/* Memory allocation and other helper functions. */
232static void
245static void
@@ -253,12 +266,12 @@ gmp_default_alloc (size_t size)
253static void *
266static void *
254gmp_default_realloc (void *old, size_t old_size, size_t new_size)
267gmp_default_realloc (void *old, size_t old_size, size_t new_size)
255{
268{
256− mp_ptr p;
269+ void * p;
257
270
258 p = realloc (old, new_size);
271 p = realloc (old, new_size);
259
272
260 if (!p)
273 if (!p)
261− gmp_die("gmp_default_realoc: Virtual memory exhausted.");
274+ gmp_die("gmp_default_realloc: Virtual memory exhausted.");
262
275
263 return p;
276 return p;
264}
277}
@@ -311,14 +324,14 @@ mp_set_memory_functions (void *(*alloc_func) (size_t),
311static mp_ptr
324static mp_ptr
312gmp_xalloc_limbs (mp_size_t size)
325gmp_xalloc_limbs (mp_size_t size)
313{
326{
314− return gmp_xalloc (size * sizeof (mp_limb_t));
327+ return (mp_ptr) gmp_xalloc (size * sizeof (mp_limb_t));
315}
328}
316
329
317static mp_ptr
330static mp_ptr
318gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
331gmp_xrealloc_limbs (mp_ptr old, mp_size_t size)
319{
332{
320 assert (size > 0);
333 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));
322}
335}
323
336
324
337
@@ -335,7 +348,7 @@ mpn_copyi (mp_ptr d, mp_srcptr s, mp_size_t n)
335void
348void
336mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
349mpn_copyd (mp_ptr d, mp_srcptr s, mp_size_t n)
337{
350{
338− while (n-- > 0)
351+ while (--n >= 0)
339 d[n] = s[n];
352 d[n] = s[n];
340}
353}
341
354
@@ -362,12 +375,23 @@ mpn_cmp4 (mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn)
362static mp_size_t
375static mp_size_t
363mpn_normalized_size (mp_srcptr xp, mp_size_t n)
376mpn_normalized_size (mp_srcptr xp, mp_size_t n)
364{
377{
365− for (; n > 0 && xp[n-1] == 0; n--)
378+ while (n > 0 && xp[n-1] == 0)
366− ;
379+ --n;
367 return n;
380 return n;
368}
381}
369
382
383+int
384+mpn_zero_p(mp_srcptr rp, mp_size_t n)
385+{
370−#define mpn_zero_p(xp, n) (mpn_normalized_size ((xp), (n)) == 0)
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+}
371
395
372mp_limb_t
396mp_limb_t
373mpn_add_1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t b)
397mpn_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)
432 {
456 {
433 mp_limb_t a = ap[i];
457 mp_limb_t a = ap[i];
434 /* Carry out */
458 /* Carry out */
435− mp_limb_t cy = a < b;;
459+ mp_limb_t cy = a < b;
436 rp[i] = a - b;
460 rp[i] = a - b;
437 b = cy;
461 b = cy;
438 }
462 }
@@ -558,17 +582,16 @@ mpn_mul (mp_ptr rp, mp_srcptr up, mp_size_t un, mp_srcptr vp, mp_size_t vn)
558 way. */
582 way. */
559
583
560 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
584 rp[un] = mpn_mul_1 (rp, up, un, vp[0]);
561− rp += 1, vp += 1, vn -= 1;
562
585
563 /* Now accumulate the product of up[] and the next higher limb from
586 /* Now accumulate the product of up[] and the next higher limb from
564 vp[]. */
587 vp[]. */
565
588
566− while (vn >= 1)
589+ while (--vn >= 1)
567 {
590 {
591+ rp += 1, vp += 1;
568 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
592 rp[un] = mpn_addmul_1 (rp, up, un, vp[0]);
569− rp += 1, vp += 1, vn -= 1;
570 }
593 }
571− return rp[un - 1];
594+ return rp[un];
572}
595}
573
596
574void
597void
@@ -588,7 +611,6 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
588{
611{
589 mp_limb_t high_limb, low_limb;
612 mp_limb_t high_limb, low_limb;
590 unsigned int tnc;
613 unsigned int tnc;
591− mp_size_t i;
592 mp_limb_t retval;
614 mp_limb_t retval;
593
615
594 assert (n >= 1);
616 assert (n >= 1);
@@ -603,7 +625,7 @@ mpn_lshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
603 retval = low_limb >> tnc;
625 retval = low_limb >> tnc;
604 high_limb = (low_limb << cnt);
626 high_limb = (low_limb << cnt);
605
627
606− for (i = n; --i != 0;)
628+ while (--n != 0)
607 {
629 {
608 low_limb = *--up;
630 low_limb = *--up;
609 *--rp = high_limb | (low_limb >> tnc);
631 *--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)
619{
641{
620 mp_limb_t high_limb, low_limb;
642 mp_limb_t high_limb, low_limb;
621 unsigned int tnc;
643 unsigned int tnc;
622− mp_size_t i;
623 mp_limb_t retval;
644 mp_limb_t retval;
624
645
625 assert (n >= 1);
646 assert (n >= 1);
@@ -631,7 +652,7 @@ mpn_rshift (mp_ptr rp, mp_srcptr up, mp_size_t n, unsigned int cnt)
631 retval = (high_limb << tnc);
652 retval = (high_limb << tnc);
632 low_limb = high_limb >> cnt;
653 low_limb = high_limb >> cnt;
633
654
634− for (i = n; --i != 0;)
655+ while (--n != 0)
635 {
656 {
636 high_limb = *up++;
657 high_limb = *up++;
637 *rp++ = low_limb | (high_limb << tnc);
658 *rp++ = low_limb | (high_limb << tnc);
@@ -682,24 +703,68 @@ mpn_scan0 (mp_srcptr ptr, mp_bitcnt_t bit)
682 i, ptr, i, GMP_LIMB_MAX);
703 i, ptr, i, GMP_LIMB_MAX);
683}
704}
684
705
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+
685
728
686/* MPN division interface. */
729/* MPN division interface. */
730+
731+/* The 3/2 inverse is defined as
732+
733+ m = floor( (B^3-1) / (B u1 + u0)) - B
734+*/
687mp_limb_t
735mp_limb_t
688mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
736mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
689{
737{
690− mp_limb_t r, p, m;
738+ mp_limb_t r, p, m, ql;
691− unsigned ul, uh;
739+ unsigned ul, uh, qh;
692− unsigned ql, qh;
693
740
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 */
697 assert (u1 >= GMP_LIMB_HIGHBIT);
741 assert (u1 >= GMP_LIMB_HIGHBIT);
698
742
743+ /* For notation, let b denote the half-limb base, so that B = b^2.
744+ Split u1 = b uh + ul. */
699 ul = u1 & GMP_LLIMB_MASK;
745 ul = u1 & GMP_LLIMB_MASK;
700 uh = u1 >> (GMP_LIMB_BITS / 2);
746 uh = u1 >> (GMP_LIMB_BITS / 2);
701
747
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. */
702 qh = ~u1 / uh;
751 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+
703 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
768 r = ((~u1 - (mp_limb_t) qh * uh) << (GMP_LIMB_BITS / 2)) | GMP_LLIMB_MASK;
704
769
705 p = (mp_limb_t) qh * ul;
770 p = (mp_limb_t) qh * ul;
@@ -717,11 +782,19 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
717 }
782 }
718 r -= p;
783 r -= p;
719
784
785+ /* Low half of the quotient is
786+
787+ ql = floor ( (b r + b-1) / u1).
788+
720− /* Do a 3/2 division (with half limb size) */
789+ This is a 3/2 division (on half-limbs), for which qh is a
790+ suitable inverse. */
791+
721 p = (r >> (GMP_LIMB_BITS / 2)) * qh + r;
792 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. */
722 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
795 ql = (p >> (GMP_LIMB_BITS / 2)) + 1;
723
796
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. */
725 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
798 r = (r << (GMP_LIMB_BITS / 2)) + GMP_LLIMB_MASK - ql * u1;
726
799
727 if (r >= (p << (GMP_LIMB_BITS / 2)))
800 if (r >= (p << (GMP_LIMB_BITS / 2)))
@@ -736,6 +809,8 @@ mpn_invert_3by2 (mp_limb_t u1, mp_limb_t u0)
736 r -= u1;
809 r -= u1;
737 }
810 }
738
811
812+ /* Now m is the 2/1 invers of u1. If u0 > 0, adjust it to become a
813+ 3/2 inverse. */
739 if (u0 > 0)
814 if (u0 > 0)
740 {
815 {
741 mp_limb_t th, tl;
816 mp_limb_t th, tl;
@@ -847,7 +922,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
847
922
848 if (inv->shift > 0)
923 if (inv->shift > 0)
849 {
924 {
850− tp = gmp_xalloc_limbs (nn);
925+ tp = qp ? qp : gmp_xalloc_limbs (nn);
851 r = mpn_lshift (tp, np, nn, inv->shift);
926 r = mpn_lshift (tp, np, nn, inv->shift);
852 np = tp;
927 np = tp;
853 }
928 }
@@ -856,7 +931,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
856
931
857 d = inv->d1;
932 d = inv->d1;
858 di = inv->di;
933 di = inv->di;
859− while (nn-- > 0)
934+ while (--nn >= 0)
860 {
935 {
861 mp_limb_t q;
936 mp_limb_t q;
862
937
@@ -864,7 +939,7 @@ mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
864 if (qp)
939 if (qp)
865 qp[nn] = q;
940 qp[nn] = q;
866 }
941 }
867− if (inv->shift > 0)
942+ if (inv->shift > 0 && !qp)
868 gmp_free (tp);
943 gmp_free (tp);
869
944
870 return r >> inv->shift;
945 return r >> inv->shift;
@@ -1140,7 +1215,7 @@ mpn_get_str_bits (unsigned char *sp, unsigned bits, mp_srcptr up, mp_size_t un)
1140 unsigned char mask;
1215 unsigned char mask;
1141 size_t sn, j;
1216 size_t sn, j;
1142 mp_size_t i;
1217 mp_size_t i;
1143− int shift;
1218+ unsigned shift;
1144
1219
1145 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1220 sn = ((un - 1) * GMP_LIMB_BITS + mpn_limb_size_in_base_2 (up[un-1])
1146 + bits - 1) / bits;
1221 + bits - 1) / bits;
@@ -1281,6 +1356,8 @@ mpn_set_str_bits (mp_ptr rp, const unsigned char *sp, size_t sn,
1281 return rn;
1356 return rn;
1282}
1357}
1283
1358
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. */
1284static mp_size_t
1361static mp_size_t
1285mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1362mpn_set_str_other (mp_ptr rp, const unsigned char *sp, size_t sn,
1286 mp_limb_t b, const struct mpn_base_info *info)
1363 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,
1290 unsigned k;
1367 unsigned k;
1291 size_t j;
1368 size_t j;
1292
1369
1370+ assert (sn > 0);
1371+
1293 k = 1 + (sn - 1) % info->exp;
1372 k = 1 + (sn - 1) % info->exp;
1294
1373
1295 j = 0;
1374 j = 0;
1296 w = sp[j++];
1375 w = sp[j++];
1297− for (; --k > 0; )
1376+ while (--k != 0)
1298 w = w * b + sp[j++];
1377 w = w * b + sp[j++];
1299
1378
1300 rp[0] = w;
1379 rp[0] = w;
1301
1380
1302− for (rn = (w > 0); j < sn;)
1381+ for (rn = 1; j < sn;)
1303 {
1382 {
1304 mp_limb_t cy;
1383 mp_limb_t cy;
1305
1384
@@ -1342,9 +1421,11 @@ mpn_set_str (mp_ptr rp, const unsigned char *sp, size_t sn, int base)
1342void
1421void
1343mpz_init (mpz_t r)
1422mpz_init (mpz_t r)
1344{
1423{
1424+ static const mp_limb_t dummy_limb = 0xc1a0;
1425+
1345− r->_mp_alloc = 1;
1426+ r->_mp_alloc = 0;
1346 r->_mp_size = 0;
1427 r->_mp_size = 0;
1347− r->_mp_d = gmp_xalloc_limbs (1);
1428+ r->_mp_d = (mp_ptr) &dummy_limb;
1348}
1429}
1349
1430
1350/* The utility of this function is a bit limited, since many functions
1431/* 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)
1365void
1446void
1366mpz_clear (mpz_t r)
1447mpz_clear (mpz_t r)
1367{
1448{
1449+ if (r->_mp_alloc)
1368− gmp_free (r->_mp_d);
1450+ gmp_free (r->_mp_d);
1369}
1451}
1370
1452
1371−static void *
1453+static mp_ptr
1372mpz_realloc (mpz_t r, mp_size_t size)
1454mpz_realloc (mpz_t r, mp_size_t size)
1373{
1455{
1374 size = GMP_MAX (size, 1);
1456 size = GMP_MAX (size, 1);
1375
1457
1458+ if (r->_mp_alloc)
1376− r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1459+ r->_mp_d = gmp_xrealloc_limbs (r->_mp_d, size);
1460+ else
1461+ r->_mp_d = gmp_xalloc_limbs (size);
1377 r->_mp_alloc = size;
1462 r->_mp_alloc = size;
1378
1463
1379 if (GMP_ABS (r->_mp_size) > size)
1464 if (GMP_ABS (r->_mp_size) > size)
@@ -1396,7 +1481,7 @@ mpz_set_si (mpz_t r, signed long int x)
1396 else /* (x < 0) */
1481 else /* (x < 0) */
1397 {
1482 {
1398 r->_mp_size = -1;
1483 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);
1400 }
1485 }
1401}
1486}
1402
1487
@@ -1406,7 +1491,7 @@ mpz_set_ui (mpz_t r, unsigned long int x)
1406 if (x > 0)
1491 if (x > 0)
1407 {
1492 {
1408 r->_mp_size = 1;
1493 r->_mp_size = 1;
1409− r->_mp_d[0] = x;
1494+ MPZ_REALLOC (r, 1)[0] = x;
1410 }
1495 }
1411 else
1496 else
1412 r->_mp_size = 0;
1497 r->_mp_size = 0;
@@ -1455,14 +1540,12 @@ mpz_fits_slong_p (const mpz_t u)
1455{
1540{
1456 mp_size_t us = u->_mp_size;
1541 mp_size_t us = u->_mp_size;
1457
1542
1458− if (us == 0)
1543+ if (us == 1)
1459− return 1;
1460− else if (us == 1)
1461 return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1544 return u->_mp_d[0] < GMP_LIMB_HIGHBIT;
1462 else if (us == -1)
1545 else if (us == -1)
1463 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1546 return u->_mp_d[0] <= GMP_LIMB_HIGHBIT;
1464 else
1547 else
1465− return 0;
1548+ return (us == 0);
1466}
1549}
1467
1550
1468int
1551int
@@ -1476,14 +1559,11 @@ mpz_fits_ulong_p (const mpz_t u)
1476long int
1559long int
1477mpz_get_si (const mpz_t u)
1560mpz_get_si (const mpz_t u)
1478{
1561{
1479− mp_size_t us = u->_mp_size;
1562+ if (u->_mp_size < 0)
1480−
1563+ /* This expression is necessary to properly handle 0x80000000 */
1481− if (us > 0)
1482− return (long) (u->_mp_d[0] & ~GMP_LIMB_HIGHBIT);
1564+ return -1 - (long) ((u->_mp_d[0] - 1) & ~GMP_LIMB_HIGHBIT);
1483− else if (us < 0)
1484− return (long) (- u->_mp_d[0] | GMP_LIMB_HIGHBIT);
1485 else
1565 else
1486− return 0;
1566+ return (long) (mpz_get_ui (u) & ~GMP_LIMB_HIGHBIT);
1487}
1567}
1488
1568
1489unsigned long int
1569unsigned long int
@@ -1507,6 +1587,48 @@ mpz_getlimbn (const mpz_t u, mp_size_t n)
1507 return 0;
1587 return 0;
1508}
1588}
1509
1589
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+
1510
1632
1511/* Conversions and comparison to double. */
1633/* Conversions and comparison to double. */
1512void
1634void
@@ -1654,9 +1776,7 @@ mpz_cmp_d (const mpz_t x, double d)
1654int
1776int
1655mpz_sgn (const mpz_t u)
1777mpz_sgn (const mpz_t u)
1656{
1778{
1657− mp_size_t usize = u->_mp_size;
1779+ return GMP_CMP (u->_mp_size, 0);
1658−
1659− return (usize > 0) - (usize < 0);
1660}
1780}
1661
1781
1662int
1782int
@@ -1671,13 +1791,7 @@ mpz_cmp_si (const mpz_t u, long v)
1671 else if (usize >= 0)
1791 else if (usize >= 0)
1672 return 1;
1792 return 1;
1673 else /* usize == -1 */
1793 else /* usize == -1 */
1674− {
1675− mp_limb_t ul = u->_mp_d[0];
1794+ return GMP_CMP (GMP_NEG_CAST (mp_limb_t, v), 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− }
1681}
1795}
1682
1796
1683int
1797int
@@ -1690,10 +1804,7 @@ mpz_cmp_ui (const mpz_t u, unsigned long v)
1690 else if (usize < 0)
1804 else if (usize < 0)
1691 return -1;
1805 return -1;
1692 else
1806 else
1693− {
1807+ return GMP_CMP (mpz_get_ui (u), v);
1694− mp_limb_t ul = (usize > 0) ? u->_mp_d[0] : 0;
1695− return (ul > v) - (ul < v);
1696− }
1697}
1808}
1698
1809
1699int
1810int
@@ -1713,15 +1824,10 @@ mpz_cmp (const mpz_t a, const mpz_t b)
1713int
1824int
1714mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1825mpz_cmpabs_ui (const mpz_t u, unsigned long v)
1715{
1826{
1716− mp_size_t un = GMP_ABS (u->_mp_size);
1827+ if (GMP_ABS (u->_mp_size) > 1)
1717− mp_limb_t ul;
1718−
1719− if (un > 1)
1720 return 1;
1828 return 1;
1721−
1829+ else
1722− ul = (un == 1) ? u->_mp_d[0] : 0;
1830+ return GMP_CMP (mpz_get_ui (u), v);
1723−
1724− return (ul > v) - (ul < v);
1725}
1831}
1726
1832
1727int
1833int
@@ -1734,18 +1840,14 @@ mpz_cmpabs (const mpz_t u, const mpz_t v)
1734void
1840void
1735mpz_abs (mpz_t r, const mpz_t u)
1841mpz_abs (mpz_t r, const mpz_t u)
1736{
1842{
1737− if (r != u)
1738− mpz_set (r, u);
1843+ mpz_set (r, u);
1739−
1740 r->_mp_size = GMP_ABS (r->_mp_size);
1844 r->_mp_size = GMP_ABS (r->_mp_size);
1741}
1845}
1742
1846
1743void
1847void
1744mpz_neg (mpz_t r, const mpz_t u)
1848mpz_neg (mpz_t r, const mpz_t u)
1745{
1849{
1746− if (r != u)
1747− mpz_set (r, u);
1850+ mpz_set (r, u);
1748−
1749 r->_mp_size = -r->_mp_size;
1851 r->_mp_size = -r->_mp_size;
1750}
1852}
1751
1853
@@ -1771,7 +1873,7 @@ mpz_abs_add_ui (mpz_t r, const mpz_t a, unsigned long b)
1771 an = GMP_ABS (a->_mp_size);
1873 an = GMP_ABS (a->_mp_size);
1772 if (an == 0)
1874 if (an == 0)
1773 {
1875 {
1774− r->_mp_d[0] = b;
1876+ MPZ_REALLOC (r, 1)[0] = b;
1775 return b > 0;
1877 return b > 0;
1776 }
1878 }
1777
1879
@@ -1790,14 +1892,15 @@ static mp_size_t
1790mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1892mpz_abs_sub_ui (mpz_t r, const mpz_t a, unsigned long b)
1791{
1893{
1792 mp_size_t an = GMP_ABS (a->_mp_size);
1894 mp_size_t an = GMP_ABS (a->_mp_size);
1793− mp_ptr rp = MPZ_REALLOC (r, an);
1895+ mp_ptr rp;
1794
1896
1795 if (an == 0)
1897 if (an == 0)
1796 {
1898 {
1797− rp[0] = b;
1899+ MPZ_REALLOC (r, 1)[0] = b;
1798 return -(b > 0);
1900 return -(b > 0);
1799 }
1901 }
1902+ rp = MPZ_REALLOC (r, an);
1800− else if (an == 1 && a->_mp_d[0] < b)
1903+ if (an == 1 && a->_mp_d[0] < b)
1801 {
1904 {
1802 rp[0] = b - a->_mp_d[0];
1905 rp[0] = b - a->_mp_d[0];
1803 return -1;
1906 return -1;
@@ -2015,12 +2118,51 @@ mpz_mul_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bits)
2015 else
2118 else
2016 mpn_copyd (rp + limbs, u->_mp_d, un);
2119 mpn_copyd (rp + limbs, u->_mp_d, un);
2017
2120
2018− while (limbs > 0)
2019− rp[--limbs] = 0;
2121+ mpn_zero (rp, limbs);
2020
2122
2021 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2123 r->_mp_size = (u->_mp_size < 0) ? - rn : rn;
2022}
2124}
2023
2125
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+
2024
2166
2025/* MPZ division */
2167/* MPZ division */
2026enum mpz_div_round_mode { GMP_DIV_FLOOR, GMP_DIV_CEIL, GMP_DIV_TRUNC };
2168enum 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,
2085 mp_size_t qn, rn;
2227 mp_size_t qn, rn;
2086 mpz_t tq, tr;
2228 mpz_t tq, tr;
2087
2229
2088− mpz_init (tr);
2230+ mpz_init_set (tr, n);
2089− mpz_set (tr, n);
2090 np = tr->_mp_d;
2231 np = tr->_mp_d;
2091
2232
2092 qn = nn - dn + 1;
2233 qn = nn - dn + 1;
@@ -2230,7 +2371,6 @@ mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
2230
2371
2231 if (qn <= 0)
2372 if (qn <= 0)
2232 qn = 0;
2373 qn = 0;
2233−
2234 else
2374 else
2235 {
2375 {
2236 qp = MPZ_REALLOC (q, qn);
2376 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,
2284 {
2424 {
2285 /* Have to negate and sign extend. */
2425 /* Have to negate and sign extend. */
2286 mp_size_t i;
2426 mp_size_t i;
2287− mp_limb_t cy;
2288
2427
2289− for (cy = 1, i = 0; i < un; i++)
2428+ gmp_assert_nocarry (! mpn_neg (rp, u->_mp_d, un));
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++)
2429+ for (i = un; i < rn - 1; i++)
2297 rp[i] = GMP_LIMB_MAX;
2430 rp[i] = GMP_LIMB_MAX;
2298
2431
2299 rp[rn-1] = mask;
2432 rp[rn-1] = mask;
@@ -2318,23 +2451,13 @@ mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
2318 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2451 if (mode == ((us > 0) ? GMP_DIV_CEIL : GMP_DIV_FLOOR)) /* us != 0 here. */
2319 {
2452 {
2320 /* If r != 0, compute 2^{bit_count} - r. */
2453 /* If r != 0, compute 2^{bit_count} - r. */
2321− mp_size_t i;
2454+ mpn_neg (rp, rp, rn);
2322
2455
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;
2333
2457
2334− /* us is not used for anything else, so we can modify it
2458+ /* us is not used for anything else, so we can modify it
2335− here to indicate flipped sign. */
2459+ here to indicate flipped sign. */
2336− us = -us;
2460+ us = -us;
2337− }
2338 }
2461 }
2339 }
2462 }
2340 rn = mpn_normalized_size (rp, rn);
2463 rn = mpn_normalized_size (rp, rn);
@@ -2389,6 +2512,24 @@ mpz_divisible_p (const mpz_t n, const mpz_t d)
2389 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2512 return mpz_div_qr (NULL, NULL, n, d, GMP_DIV_TRUNC) == 0;
2390}
2513}
2391
2514
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+
2392static unsigned long
2533static unsigned long
2393mpz_div_qr_ui (mpz_t q, mpz_t r,
2534mpz_div_qr_ui (mpz_t q, mpz_t r,
2394 const mpz_t n, unsigned long d, enum mpz_div_round_mode mode)
2535 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,
2431
2572
2432 if (r)
2573 if (r)
2433 {
2574 {
2434− r->_mp_d[0] = rl;
2575+ MPZ_REALLOC (r, 1)[0] = rl;
2435 r->_mp_size = rs;
2576 r->_mp_size = rs;
2436 }
2577 }
2437 if (q)
2578 if (q)
@@ -2955,9 +3096,7 @@ void
2955mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
3096mpz_ui_pow_ui (mpz_t r, unsigned long blimb, unsigned long e)
2956{
3097{
2957 mpz_t b;
3098 mpz_t b;
2958− mpz_init_set_ui (b, blimb);
2959− mpz_pow_ui (r, b, e);
3099+ mpz_pow_ui (r, mpz_roinit_n (b, &blimb, 1), e);
2960− mpz_clear (b);
2961}
3100}
2962
3101
2963void
3102void
@@ -3029,7 +3168,7 @@ mpz_powm (mpz_t r, const mpz_t b, const mpz_t e, const mpz_t m)
3029 }
3168 }
3030 mpz_init_set_ui (tr, 1);
3169 mpz_init_set_ui (tr, 1);
3031
3170
3032− while (en-- > 0)
3171+ while (--en >= 0)
3033 {
3172 {
3034 mp_limb_t w = e->_mp_d[en];
3173 mp_limb_t w = e->_mp_d[en];
3035 mp_limb_t bit;
3174 mp_limb_t bit;
@@ -3069,9 +3208,7 @@ void
3069mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3208mpz_powm_ui (mpz_t r, const mpz_t b, unsigned long elimb, const mpz_t m)
3070{
3209{
3071 mpz_t e;
3210 mpz_t e;
3072− mpz_init_set_ui (e, elimb);
3073− mpz_powm (r, b, e, m);
3211+ mpz_powm (r, b, mpz_roinit_n (e, &elimb, 1), m);
3074− mpz_clear (e);
3075}
3212}
3076
3213
3077/* x=trunc(y^(1/z)), r=y-x^z */
3214/* 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)
3088 gmp_die ("mpz_rootrem: Zeroth root.");
3225 gmp_die ("mpz_rootrem: Zeroth root.");
3089
3226
3090 if (mpz_cmpabs_ui (y, 1) <= 0) {
3227 if (mpz_cmpabs_ui (y, 1) <= 0) {
3228+ if (x)
3091− mpz_set (x, y);
3229+ mpz_set (x, y);
3092 if (r)
3230 if (r)
3093 r->_mp_size = 0;
3231 r->_mp_size = 0;
3094 return;
3232 return;
3095 }
3233 }
3096
3234
3097 mpz_init (u);
3235 mpz_init (u);
3098− {
3236+ mpz_init (t);
3099− mp_bitcnt_t tb;
3100− tb = mpz_sizeinbase (y, 2) / z + 1;
3237+ mpz_setbit (t, mpz_sizeinbase (y, 2) / z + 1);
3101− mpz_init2 (t, tb);
3102− mpz_setbit (t, tb);
3103− }
3104
3238
3105 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3239 if (z == 2) /* simplify sqrt loop: z-1 == 1 */
3106 do {
3240 do {
@@ -3132,7 +3266,8 @@ mpz_rootrem (mpz_t x, mpz_t r, const mpz_t y, unsigned long z)
3132 mpz_pow_ui (t, u, z);
3266 mpz_pow_ui (t, u, z);
3133 mpz_sub (r, y, t);
3267 mpz_sub (r, y, t);
3134 }
3268 }
3269+ if (x)
3135− mpz_swap (x, u);
3270+ mpz_swap (x, u);
3136 mpz_clear (u);
3271 mpz_clear (u);
3137 mpz_clear (t);
3272 mpz_clear (t);
3138}
3273}
@@ -3164,6 +3299,47 @@ mpz_sqrt (mpz_t s, const mpz_t u)
3164 mpz_rootrem (s, NULL, u, 2);
3299 mpz_rootrem (s, NULL, u, 2);
3165}
3300}
3166
3301
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+}
3167
3343
3168/* Combinatorics */
3344/* Combinatorics */
3169
3345
@@ -3171,7 +3347,7 @@ void
3171mpz_fac_ui (mpz_t x, unsigned long n)
3347mpz_fac_ui (mpz_t x, unsigned long n)
3172{
3348{
3173 mpz_set_ui (x, n + (n == 0));
3349 mpz_set_ui (x, n + (n == 0));
3174− for (;n > 2;)
3350+ while (n > 2)
3175 mpz_mul_ui (x, x, --n);
3351 mpz_mul_ui (x, x, --n);
3176}
3352}
3177
3353
@@ -3195,6 +3371,102 @@ mpz_bin_uiui (mpz_t r, unsigned long n, unsigned long k)
3195 mpz_clear (t);
3371 mpz_clear (t);
3196}
3372}
3197
3373
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+
3198
3470
3199/* Logical operations and bit manipulation. */
3471/* Logical operations and bit manipulation. */
3200
3472
@@ -3246,7 +3518,7 @@ mpz_tstbit (const mpz_t d, mp_bitcnt_t bit_index)
3246 must be complemented. */
3518 must be complemented. */
3247 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3519 if (shift > 0 && (w << (GMP_LIMB_BITS - shift)) > 0)
3248 return bit ^ 1;
3520 return bit ^ 1;
3249− while (limb_index-- > 0)
3521+ while (--limb_index >= 0)
3250 if (d->_mp_d[limb_index] > 0)
3522 if (d->_mp_d[limb_index] > 0)
3251 return bit ^ 1;
3523 return bit ^ 1;
3252 }
3524 }
@@ -3311,7 +3583,7 @@ mpz_abs_sub_bit (mpz_t d, mp_bitcnt_t bit_index)
3311
3583
3312 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3584 gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
3313 dn - limb_index, bit));
3585 dn - limb_index, bit));
3314− dn -= (dp[dn-1] == 0);
3586+ dn = mpn_normalized_size (dp, dn);
3315 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3587 d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
3316}
3588}
3317
3589
@@ -3389,7 +3661,7 @@ mpz_and (mpz_t r, const mpz_t u, const mpz_t v)
3389 /* If the smaller input is positive, higher limbs don't matter. */
3661 /* If the smaller input is positive, higher limbs don't matter. */
3390 rn = vx ? un : vn;
3662 rn = vx ? un : vn;
3391
3663
3392− rp = MPZ_REALLOC (r, rn + rc);
3664+ rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3393
3665
3394 up = u->_mp_d;
3666 up = u->_mp_d;
3395 vp = v->_mp_d;
3667 vp = v->_mp_d;
@@ -3462,7 +3734,7 @@ mpz_ior (mpz_t r, const mpz_t u, const mpz_t v)
3462 don't matter. */
3734 don't matter. */
3463 rn = vx ? vn : un;
3735 rn = vx ? vn : un;
3464
3736
3465− rp = MPZ_REALLOC (r, rn + rc);
3737+ rp = MPZ_REALLOC (r, rn + (mp_size_t) rc);
3466
3738
3467 up = u->_mp_d;
3739 up = u->_mp_d;
3468 vp = v->_mp_d;
3740 vp = v->_mp_d;
@@ -3531,7 +3803,7 @@ mpz_xor (mpz_t r, const mpz_t u, const mpz_t v)
3531 vx = -vc;
3803 vx = -vc;
3532 rx = -rc;
3804 rx = -rc;
3533
3805
3534− rp = MPZ_REALLOC (r, un + rc);
3806+ rp = MPZ_REALLOC (r, un + (mp_size_t) rc);
3535
3807
3536 up = u->_mp_d;
3808 up = u->_mp_d;
3537 vp = v->_mp_d;
3809 vp = v->_mp_d;
@@ -3587,20 +3859,28 @@ gmp_popcount_limb (mp_limb_t x)
3587}
3859}
3588
3860
3589mp_bitcnt_t
3861mp_bitcnt_t
3590−mpz_popcount (const mpz_t u)
3862+mpn_popcount (mp_srcptr p, mp_size_t n)
3591{
3863{
3592− mp_size_t un, i;
3864+ mp_size_t i;
3593 mp_bitcnt_t c;
3865 mp_bitcnt_t c;
3594
3866
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+
3595 un = u->_mp_size;
3878 un = u->_mp_size;
3596
3879
3597 if (un < 0)
3880 if (un < 0)
3598 return ~(mp_bitcnt_t) 0;
3881 return ~(mp_bitcnt_t) 0;
3599
3882
3600− for (c = 0, i = 0; i < un; i++)
3883+ return mpn_popcount (u->_mp_d, un);
3601− c += gmp_popcount_limb (u->_mp_d[i]);
3602−
3603− return c;
3604}
3884}
3605
3885
3606mp_bitcnt_t
3886mp_bitcnt_t
@@ -3799,7 +4079,7 @@ mpz_get_str (char *sp, int base, const mpz_t u)
3799
4079
3800 sn = 1 + mpz_sizeinbase (u, base);
4080 sn = 1 + mpz_sizeinbase (u, base);
3801 if (!sp)
4081 if (!sp)
3802− sp = gmp_xalloc (1 + sn);
4082+ sp = (char *) gmp_xalloc (1 + sn);
3803
4083
3804 un = GMP_ABS (u->_mp_size);
4084 un = GMP_ABS (u->_mp_size);
3805
4085
@@ -3846,7 +4126,7 @@ mpz_set_str (mpz_t r, const char *sp, int base)
3846 unsigned bits;
4126 unsigned bits;
3847 mp_size_t rn, alloc;
4127 mp_size_t rn, alloc;
3848 mp_ptr rp;
4128 mp_ptr rp;
3849− size_t sn;
4129+ size_t dn;
3850 int sign;
4130 int sign;
3851 unsigned char *dp;
4131 unsigned char *dp;
3852
4132
@@ -3860,18 +4140,17 @@ mpz_set_str (mpz_t r, const char *sp, int base)
3860
4140
3861 if (base == 0)
4141 if (base == 0)
3862 {
4142 {
3863− if (*sp == '0')
4143+ if (sp[0] == '0')
3864 {
4144 {
3865− sp++;
3866− if (*sp == 'x' || *sp == 'X')
4145+ if (sp[1] == 'x' || sp[1] == 'X')
3867 {
4146 {
3868 base = 16;
4147 base = 16;
3869− sp++;
4148+ sp += 2;
3870 }
4149 }
3871− else if (*sp == 'b' || *sp == 'B')
4150+ else if (sp[1] == 'b' || sp[1] == 'B')
3872 {
4151 {
3873 base = 2;
4152 base = 2;
3874− sp++;
4153+ sp += 2;
3875 }
4154 }
3876 else
4155 else
3877 base = 8;
4156 base = 8;
@@ -3880,16 +4159,20 @@ mpz_set_str (mpz_t r, const char *sp, int base)
3880 base = 10;
4159 base = 10;
3881 }
4160 }
3882
4161
3883− sn = strlen (sp);
4162+ if (!*sp)
4163+ {
4164+ r->_mp_size = 0;
4165+ return -1;
4166+ }
3884− dp = gmp_xalloc (sn + (sn == 0));
4167+ dp = (unsigned char *) gmp_xalloc (strlen (sp));
3885
4168
3886− for (sn = 0; *sp; sp++)
4169+ for (dn = 0; *sp; sp++)
3887 {
4170 {
3888 unsigned digit;
4171 unsigned digit;
3889
4172
3890 if (isspace ((unsigned char) *sp))
4173 if (isspace ((unsigned char) *sp))
3891 continue;
4174 continue;
3892− if (*sp >= '0' && *sp <= '9')
4175+ else if (*sp >= '0' && *sp <= '9')
3893 digit = *sp - '0';
4176 digit = *sp - '0';
3894 else if (*sp >= 'a' && *sp <= 'z')
4177 else if (*sp >= 'a' && *sp <= 'z')
3895 digit = *sp - 'a' + 10;
4178 digit = *sp - 'a' + 10;
@@ -3898,31 +4181,40 @@ mpz_set_str (mpz_t r, const char *sp, int base)
3898 else
4181 else
3899 digit = base; /* fail */
4182 digit = base; /* fail */
3900
4183
3901− if (digit >= base)
4184+ if (digit >= (unsigned) base)
3902 {
4185 {
3903 gmp_free (dp);
4186 gmp_free (dp);
3904 r->_mp_size = 0;
4187 r->_mp_size = 0;
3905 return -1;
4188 return -1;
3906 }
4189 }
3907
4190
3908− dp[sn++] = digit;
4191+ dp[dn++] = digit;
3909 }
4192 }
3910
4193
4194+ if (!dn)
4195+ {
4196+ gmp_free (dp);
4197+ r->_mp_size = 0;
4198+ return -1;
4199+ }
3911 bits = mpn_base_power_of_two_p (base);
4200 bits = mpn_base_power_of_two_p (base);
3912
4201
3913 if (bits > 0)
4202 if (bits > 0)
3914 {
4203 {
3915− alloc = (sn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
4204+ alloc = (dn * bits + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
3916 rp = MPZ_REALLOC (r, alloc);
4205 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);
3918 }
4207 }
3919 else
4208 else
3920 {
4209 {
3921 struct mpn_base_info info;
4210 struct mpn_base_info info;
3922 mpn_get_base_info (&info, base);
4211 mpn_get_base_info (&info, base);
3923− alloc = (sn + info.exp - 1) / info.exp;
4212+ alloc = (dn + info.exp - 1) / info.exp;
3924 rp = MPZ_REALLOC (r, alloc);
4213 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;
3926 }
4218 }
3927 assert (rn <= alloc);
4219 assert (rn <= alloc);
3928 gmp_free (dp);
4220 gmp_free (dp);
rts/mini-gmp.h+47−8
@@ -1,21 +1,32 @@
1/* mini-gmp, a minimalistic implementation of a GNU GMP subset.
1/* mini-gmp, a minimalistic implementation of a GNU GMP subset.
2
2
3−Copyright 2011, 2012, 2013 Free Software Foundation, Inc.
3+Copyright 2011-2015 Free Software Foundation, Inc.
4
4
5This file is part of the GNU MP Library.
5This file is part of the GNU MP Library.
6
6
7The GNU MP Library is free software; you can redistribute it and/or modify
7The GNU MP Library is free software; you can redistribute it and/or modify
8+it under the terms of either:
9+
8−it under the terms of the GNU Lesser General Public License as published by
10+ * the GNU Lesser General Public License as published by the Free
9−the Free Software Foundation; either version 3 of the License, or (at your
11+ Software Foundation; either version 3 of the License, or (at your
10−option) any later version.
12+ option) any later version.
13+
14+or
15+
16+ * the GNU General Public License as published by the Free Software
17+ Foundation; either version 2 of the License, or (at your option) any
18+ later version.
19+
20+or both in parallel, as here.
11
21
12The GNU MP Library is distributed in the hope that it will be useful, but
22The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
23WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14−or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
24+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15−License for more details.
25+for more details.
16
26
17−You should have received a copy of the GNU Lesser General Public License
27+You should have received copies of the GNU General Public License and the
18−along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
28+GNU Lesser General Public License along with the GNU MP Library. If not,
29+see https://www.gnu.org/licenses/. */
19
30
20/* About mini-gmp: This is a minimal implementation of a subset of the
31/* About mini-gmp: This is a minimal implementation of a subset of the
21 GMP interface. It is intended for inclusion into applications which
32 GMP interface. It is intended for inclusion into applications which
@@ -64,10 +75,14 @@ typedef __mpz_struct mpz_t[1];
64typedef __mpz_struct *mpz_ptr;
75typedef __mpz_struct *mpz_ptr;
65typedef const __mpz_struct *mpz_srcptr;
76typedef const __mpz_struct *mpz_srcptr;
66
77
78+extern const int mp_bits_per_limb;
79+
67void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t);
80void mpn_copyi (mp_ptr, mp_srcptr, mp_size_t);
68void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t);
81void mpn_copyd (mp_ptr, mp_srcptr, mp_size_t);
82+void mpn_zero (mp_ptr, mp_size_t);
69
83
70int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t);
84int mpn_cmp (mp_srcptr, mp_srcptr, mp_size_t);
85+int mpn_zero_p (mp_srcptr, mp_size_t);
71
86
72mp_limb_t mpn_add_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t);
87mp_limb_t mpn_add_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t);
73mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
88mp_limb_t mpn_add_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
@@ -84,6 +99,8 @@ mp_limb_t mpn_submul_1 (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t);
84mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
99mp_limb_t mpn_mul (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
85void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
100void mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t);
86void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t);
101void mpn_sqr (mp_ptr, mp_srcptr, mp_size_t);
102+int mpn_perfect_square_p (mp_srcptr, mp_size_t);
103+mp_size_t mpn_sqrtrem (mp_ptr, mp_ptr, mp_srcptr, mp_size_t);
87
104
88mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
105mp_limb_t mpn_lshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
89mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
106mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
@@ -91,6 +108,11 @@ mp_limb_t mpn_rshift (mp_ptr, mp_srcptr, mp_size_t, unsigned int);
91mp_bitcnt_t mpn_scan0 (mp_srcptr, mp_bitcnt_t);
108mp_bitcnt_t mpn_scan0 (mp_srcptr, mp_bitcnt_t);
92mp_bitcnt_t mpn_scan1 (mp_srcptr, mp_bitcnt_t);
109mp_bitcnt_t mpn_scan1 (mp_srcptr, mp_bitcnt_t);
93
110
111+void mpn_com (mp_ptr, mp_srcptr, mp_size_t);
112+mp_limb_t mpn_neg (mp_ptr, mp_srcptr, mp_size_t);
113+
114+mp_bitcnt_t mpn_popcount (mp_srcptr, mp_size_t);
115+
94mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t);
116mp_limb_t mpn_invert_3by2 (mp_limb_t, mp_limb_t);
95#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0)
117#define mpn_invert_limb(x) mpn_invert_3by2 ((x), 0)
96
118
@@ -127,6 +149,10 @@ void mpz_mul_si (mpz_t, const mpz_t, long int);
127void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int);
149void mpz_mul_ui (mpz_t, const mpz_t, unsigned long int);
128void mpz_mul (mpz_t, const mpz_t, const mpz_t);
150void mpz_mul (mpz_t, const mpz_t, const mpz_t);
129void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
151void mpz_mul_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
152+void mpz_addmul_ui (mpz_t, const mpz_t, unsigned long int);
153+void mpz_addmul (mpz_t, const mpz_t, const mpz_t);
154+void mpz_submul_ui (mpz_t, const mpz_t, unsigned long int);
155+void mpz_submul (mpz_t, const mpz_t, const mpz_t);
130
156
131void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
157void mpz_cdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
132void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
158void mpz_fdiv_qr (mpz_t, mpz_t, const mpz_t, const mpz_t);
@@ -150,6 +176,7 @@ void mpz_mod (mpz_t, const mpz_t, const mpz_t);
150void mpz_divexact (mpz_t, const mpz_t, const mpz_t);
176void mpz_divexact (mpz_t, const mpz_t, const mpz_t);
151
177
152int mpz_divisible_p (const mpz_t, const mpz_t);
178int mpz_divisible_p (const mpz_t, const mpz_t);
179+int mpz_congruent_p (const mpz_t, const mpz_t, const mpz_t);
153
180
154unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
181unsigned long mpz_cdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
155unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
182unsigned long mpz_fdiv_qr_ui (mpz_t, mpz_t, const mpz_t, unsigned long);
@@ -179,6 +206,7 @@ int mpz_invert (mpz_t, const mpz_t, const mpz_t);
179
206
180void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t);
207void mpz_sqrtrem (mpz_t, mpz_t, const mpz_t);
181void mpz_sqrt (mpz_t, const mpz_t);
208void mpz_sqrt (mpz_t, const mpz_t);
209+int mpz_perfect_square_p (const mpz_t);
182
210
183void mpz_pow_ui (mpz_t, const mpz_t, unsigned long);
211void mpz_pow_ui (mpz_t, const mpz_t, unsigned long);
184void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long);
212void mpz_ui_pow_ui (mpz_t, unsigned long, unsigned long);
@@ -191,6 +219,8 @@ int mpz_root (mpz_t, const mpz_t, unsigned long);
191void mpz_fac_ui (mpz_t, unsigned long);
219void mpz_fac_ui (mpz_t, unsigned long);
192void mpz_bin_uiui (mpz_t, unsigned long, unsigned long);
220void mpz_bin_uiui (mpz_t, unsigned long, unsigned long);
193
221
222+int mpz_probab_prime_p (const mpz_t, int);
223+
194int mpz_tstbit (const mpz_t, mp_bitcnt_t);
224int mpz_tstbit (const mpz_t, mp_bitcnt_t);
195void mpz_setbit (mpz_t, mp_bitcnt_t);
225void mpz_setbit (mpz_t, mp_bitcnt_t);
196void mpz_clrbit (mpz_t, mp_bitcnt_t);
226void mpz_clrbit (mpz_t, mp_bitcnt_t);
@@ -214,6 +244,15 @@ double mpz_get_d (const mpz_t);
214size_t mpz_size (const mpz_t);
244size_t mpz_size (const mpz_t);
215mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t);
245mp_limb_t mpz_getlimbn (const mpz_t, mp_size_t);
216
246
247+void mpz_realloc2 (mpz_t, mp_bitcnt_t);
248+mp_srcptr mpz_limbs_read (mpz_srcptr);
249+mp_ptr mpz_limbs_modify (mpz_t, mp_size_t);
250+mp_ptr mpz_limbs_write (mpz_t, mp_size_t);
251+void mpz_limbs_finish (mpz_t, mp_size_t);
252+mpz_srcptr mpz_roinit_n (mpz_t, mp_srcptr, mp_size_t);
253+
254+#define MPZ_ROINIT_N(xp, xs) {{0, (xs),(xp) }}
255+
217void mpz_set_si (mpz_t, signed long int);
256void mpz_set_si (mpz_t, signed long int);
218void mpz_set_ui (mpz_t, unsigned long int);
257void mpz_set_ui (mpz_t, unsigned long int);
219void mpz_set (mpz_t, const mpz_t);
258void mpz_set (mpz_t, const mpz_t);
test/TestData.hs+2−1
@@ -80,7 +80,8 @@ testFamiliesData = [
80 ( 26, ANY )]),
80 ( 26, ANY )]),
81 ("bignum", "Bignum",
81 ("bignum", "Bignum",
82 [ ( 1, ANY ),
82 [ ( 1, ANY ),
83− ( 2, ANY )]),
83+ ( 2, ANY ),
84+ ( 3, C_CG )]),
84 ("bounded", "Bounded",
85 ("bounded", "Bounded",
85 [ ( 1, ANY )]),
86 [ ( 1, ANY )]),
86 ("buffer", "Buffer",
87 ("buffer", "Buffer",
test/bignum003/bignum003.idrnew file+23
@@ -0,0 +1,23 @@
1+allZeroesFromTo : Int -> Int -> String -> Bool
2+allZeroesFromTo n m str
3+ = if n < m
4+ then if strIndex str n == '0'
5+ then allZeroesFromTo (n+1) m str
6+ else False
7+ else True
8+
9+isNthPowerOfTen : Int -> String -> Bool
10+isNthPowerOfTen n str
11+ = let m = cast (length str) in
12+ if m == n+1
13+ then strHead str == '1' && allZeroesFromTo 1 m str
14+ else False
15+
16+main : IO ()
17+main = do
18+ let x = power 10 12468
19+ let sx = show x
20+ printLn (isNthPowerOfTen 12468 sx)
21+ let y = power 10 12470
22+ let sy = show y
23+ printLn (isNthPowerOfTen 12470 sy)
test/bignum003/expectednew file+2
@@ -0,0 +1,2 @@
1+True
2+True
test/bignum003/runnew file+4
@@ -0,0 +1,4 @@
1+#!/usr/bin/env bash
2+${IDRIS:-idris} $@ bignum003.idr -o bignum003
3+./bignum003
4+rm -f bignum003 *.ibc