[Gmp-commit] /var/hg/gmp: 4 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Mon Aug 31 20:22:11 CEST 2026
details: /var/hg/gmp/rev/c3eb728d8467
changeset: 18584:c3eb728d8467
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Wed Aug 26 18:06:21 2026 +0200
description:
primesieve.c (gmp_primesieve): Use MPN_COPY.
details: /var/hg/gmp/rev/914a0905ad30
changeset: 18585:914a0905ad30
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Wed Aug 26 18:12:23 2026 +0200
description:
mpz/millerrabin.c (mpz_millerrabin): Save higest limb and use it for checking sizes.
details: /var/hg/gmp/rev/e066bc2f906b
changeset: 18586:e066bc2f906b
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Aug 31 20:20:31 2026 +0200
description:
mpn/generic/bsqrtinv.c: Smaller table and simpler code for initialization.
details: /var/hg/gmp/rev/ab4a7408e2ce
changeset: 18587:ab4a7408e2ce
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Aug 31 20:21:46 2026 +0200
description:
mpn/generic/bsqrt.c: Add a comment.
diffstat:
ChangeLog | 2 +
mpn/generic/bsqrt.c | 5 +-
mpn/generic/bsqrtinv.c | 109 ++++++++++++++++++++----------------------------
mpz/millerrabin.c | 29 ++++++++----
primesieve.c | 4 +-
5 files changed, 72 insertions(+), 77 deletions(-)
diffs (truncated from 307 to 300 lines):
diff -r b74bd2160c17 -r ab4a7408e2ce ChangeLog
--- a/ChangeLog Mon Aug 17 15:09:43 2026 +0200
+++ b/ChangeLog Mon Aug 31 20:21:46 2026 +0200
@@ -3,6 +3,8 @@
* mpz/millerrabin.c: Special (simple) test for Proth numbers.
* tests/mpz/t-pprime_p.c: Exercise the new test.
+ * primesieve.c (gmp_primesieve): Use MPN_COPY.
+
2026-07-20 Marco Bodrato <marco.bodrato at protonmail.com>
* mpn/generic/bsqrt.c: Use a last iteration to increase precision
diff -r b74bd2160c17 -r ab4a7408e2ce mpn/generic/bsqrt.c
--- a/mpn/generic/bsqrt.c Mon Aug 17 15:09:43 2026 +0200
+++ b/mpn/generic/bsqrt.c Mon Aug 31 20:21:46 2026 +0200
@@ -34,7 +34,8 @@
used by mpn_bsqrtinv, i.e 3*(1 + nb / GMP_NUMB_BITS)
For large enough sizes, it calls mpn_bsqrtinv for just half the
- needed precision, the last half is computed here.
+ needed precision, the last half is computed here, using a
+ Karp–Markstein step.
If T is the result of mpn_bsqrtinv, and A the input,
R = A * T;
@@ -64,7 +65,7 @@
{
mp_ptr sp = tp + nn;
- MPN_FILL (tp, nn, CNST_LIMB(0));
+ MPN_FILL (tp + 1, nn - 1, CNST_LIMB(0));
if (! mpn_bsqrtinv (tp, ap, nb, sp))
return 0;
diff -r b74bd2160c17 -r ab4a7408e2ce mpn/generic/bsqrtinv.c
--- a/mpn/generic/bsqrtinv.c Mon Aug 17 15:09:43 2026 +0200
+++ b/mpn/generic/bsqrtinv.c Mon Aug 31 20:21:46 2026 +0200
@@ -57,27 +57,19 @@
#ifndef BSQRTINV_DONT_USE_TABLE
/* Generated with GP-Pari:
- b=8;v=Vecsmall(-binary(2^2^b-1));
- forstep(i=1,2^(b+1),2,v[lift(Mod(i,2^(b+3))^-2)>>3+1]=i\2);
+ b=7;v=Vecsmall(-binary(2^2^b-1));
+ forstep(i=1,2^(b+1),2,v[lift(Mod(i,2^(b+3))^-2)>>3+1]=i);
for(i=1,2^b,print1(v[i],",");if(i%16==0,print(),print1(" ")))
*/
-static const unsigned char binvsqrttab[256] = /* The least significant 1 was removed */
- { 0, 170, 172, 102, 184, 253, 219, 129, 240, 218, 227, 22, 168, 50, 148, 46,
- 31, 245, 115, 57, 152, 157, 4, 222, 208, 197, 3, 137, 136, 146, 139, 113,
- 63, 149, 108, 217, 120, 61, 228, 62, 176, 101, 220, 214, 104, 242, 84, 238,
- 95, 53, 179, 134, 88, 34, 59, 97, 144, 5, 67, 54, 72, 173, 203, 78,
- 127, 42, 44, 25, 56, 130, 164, 254, 112, 90, 156, 105, 40, 77, 20, 81,
- 159, 138, 243, 185, 24, 226, 123, 94, 80, 186, 131, 246, 8, 18, 244, 241,
- 191, 234, 19, 166, 7, 189, 100, 65, 48, 229, 92, 86, 23, 114, 43, 110,
- 223, 181, 204, 6, 39, 93, 187, 225, 16, 133, 195, 73, 55, 210, 180, 49,
- 255, 85, 83, 153, 71, 2, 36, 126, 15, 37, 28, 233, 87, 205, 107, 209,
- 224, 10, 140, 198, 103, 98, 251, 33, 47, 58, 252, 118, 119, 109, 116, 142,
- 192, 106, 147, 38, 135, 194, 27, 193, 79, 154, 35, 41, 151, 13, 171, 17,
- 160, 202, 76, 121, 167, 221, 196, 158, 111, 250, 188, 201, 183, 82, 52, 177,
- 128, 213, 211, 230, 199, 125, 91, 1, 143, 165, 99, 150, 215, 178, 235, 174,
- 96, 117, 12, 70, 231, 29, 132, 161, 175, 69, 124, 9, 247, 237, 11, 14,
- 64, 21, 236, 89, 248, 66, 155, 190, 207, 26, 163, 169, 232, 141, 212, 145,
- 32, 74, 51, 249, 216, 162, 68, 30, 239, 122, 60, 182, 200, 45, 75, 206};
+static const unsigned char binvsqrttab[128] =
+ { 1, 171, 167, 205, 143, 5, 73, 253, 31, 75, 57, 45, 175, 101, 215, 93,
+ 63, 21, 231, 115, 207, 197, 9, 67, 95, 117, 7, 237, 239, 219, 233, 227,
+ 127, 213, 217, 77, 241, 123, 55, 125, 159, 203, 71, 83, 209, 27, 169, 35,
+ 191, 107, 153, 243, 177, 69, 119, 195, 223, 11, 135, 109, 145, 165, 105, 157,
+ 255, 85, 89, 51, 113, 251, 183, 3, 225, 181, 199, 211, 81, 155, 41, 163,
+ 193, 235, 25, 141, 49, 59, 247, 189, 161, 139, 249, 19, 17, 37, 23, 29,
+ 129, 43, 39, 179, 15, 133, 201, 131, 97, 53, 185, 173, 47, 229, 87, 221,
+ 65, 149, 103, 13, 79, 187, 137, 61, 33, 245, 121, 147, 111, 91, 151, 99};
#endif
/* tp needs 2*(1 + bnb / GMP_NUMB_BITS) limbs of space */
@@ -88,7 +80,7 @@
ASSERT (bnb > 0);
#ifndef BSQRTINV_RP_NOT_ZEROED
- ASSERT ((bnb <= GMP_NUMB_BITS) || mpn_zero_p (rp, 1 + bnb / GMP_NUMB_BITS));
+ ASSERT ((bnb <= GMP_NUMB_BITS) || mpn_zero_p (rp + 1, bnb / GMP_NUMB_BITS));
#endif
if (UNLIKELY (bnb == 1))
{
@@ -99,63 +91,49 @@
{
mp_ptr tp2 = tp + 1 + bnb / GMP_NUMB_BITS;
mp_size_t bn, order[GMP_LIMB_BITS + 1];
- mp_limb_t t0, r0;
+ mp_limb_t r0;
int i;
if ((y0 & 7) != 1)
return 0;
- /* Rekte : 4 */
- /* 4 -4x-> 11 -4x-> 32 (ne 33, 8x) */
- /* 4 -4x-> 11 -3x-> 21 -3x-> preter 33 (10x) [4,33]*/
-
- /* 4 -4x-> 11 -4x-> 32 -3x-> 63 (ne 64, 11x) [3,4]*/
- /* 4 -4x-> 11 -4x-> 32 -4x-> preter 65 (12x) */
-
- /* 4 -3x-> 7 -3x-> 13 -3x-> 25 -4x-> preter 65 (13x) */
- /* 4 -3x-> 7 -3x-> 13 -3x-> 25 -3x-> preter 33 (12x) */
- /* 4 -3x-> 7 -3x-> 13 -4x-> preter 33 (10x) */
-
- /* Tabulo: 10 */
- /* 10 -4x-> 29 -4x-> preter 65, 8x [4,4]*/
- /* 10 -3x-> 19 -3x-> preter 33, 6x [3,3]*/
-
#ifdef BSQRTINV_DONT_USE_TABLE
/* 16-bits computations are enough */
- unsigned ru = y0 + ((y0 & 8) >> 2) + ((y0 & 16) >> 1);
+ unsigned ru = 1 + ((y0 & 8) >> 2) + ((y0 & 16) >> 1);
unsigned tu = ru * ru * (unsigned) y0 >> 1;
ASSERT ((tu & (GMP_NUMB_MAX >> (GMP_NUMB_BITS - 4))) == 0);
ru += ru * tu * ((tu >> 1) + tu - 1); /* Halley -> 11 */
+ // ru += ru * (tu >> 1) * (3 * tu - 2); /* Halley -> 11 */
/* Better sequences are possible from size 11, but this
code is currently not used. */
+ r0 = ru;
#else /* ! defined(BSQRTINV_DONT_USE_TABLE) */
- unsigned ru = binvsqrttab[(y0 >> 3) & 0xff];
- ru = (ru << 1) + 1;
+ r0 = binvsqrttab[(y0 >> 3) & 0x7f];
#endif
- r0 = ru;
-#if GMP_NUMB_BITS < 10 * 2 - 2
- const mp_bitcnt_t precomputed_bits = 10;
-#else /* GMP_NUMB_BITS >= 10 * 2 - 2 */
- t0 = r0 * r0 * y0 >> 1;
- ASSERT ((t0 & (GMP_NUMB_MAX >> (GMP_NUMB_BITS - 10))) == 0);
-#if GMP_NUMB_BITS < 19 * 2 - 2
- r0 -= r0 * t0;
- const mp_bitcnt_t precomputed_bits = 19;
-#else /* GMP_NUMB_BITS >= 19 * 2 - 2 */
- mp_limb_t t2 = t0 * t0;
- r0 += r0 * t0 * ((t0 >> 1) + t0 - 1 - (t2 >> 1) - (t2 << 1)); /* Halley -> 37*/
-#if GMP_NUMB_BITS < 37 * 2 - 2
- const mp_bitcnt_t precomputed_bits = 37;
-#else /* GMP_NUMB_BITS >= 37 * 2 - 2 */
-#error Not implemented, yet.
+#if GMP_NUMB_BITS < 9 * 2 - 1
+ const mp_bitcnt_t precomputed_bits = 9;
+#else /* GMP_NUMB_BITS > 9 * 2 - 2 */
+ r0 -= r0 * (y0 * r0 * r0 >> 1); /* 9 -> 17 */
+#if GMP_NUMB_BITS < 17 * 2 - 1
+ const mp_bitcnt_t precomputed_bits = 17;
+#else /* GMP_NUMB_BITS > 17 * 2 - 2 */
+ r0 -= r0 * (y0 * r0 * r0 >> 1); /* 17 -> 33 */
+#if GMP_NUMB_BITS < 33 * 2 - 1
+ const mp_bitcnt_t precomputed_bits = 33;
+#else /* GMP_NUMB_BITS > 33 * 2 - 2 */
+ mp_bitcnt_t precomputed_bits = 33;
+ do {
+ r0 -= r0 * (y0 * r0 * r0 >> 1); /* n -> 2*n-1 */
+ precomputed_bits = precomputed_bits * 2 - 1;
+ } while ((GMP_NUMB_BITS + 3) / 2 > precomputed_bits);
#endif
#endif
#endif
i = 0;
- for (; bnb > GMP_NUMB_BITS + 1; bnb = (bnb + 2) >> 1)
+ for (; bnb > GMP_NUMB_BITS + 1; bnb = (bnb >> 1) + 1)
order[i++] = bnb;
if (bnb > precomputed_bits) {
if (bnb >= GMP_NUMB_BITS) {
@@ -172,15 +150,12 @@
rp[1] = 0;
#endif
} else {
- t0 = r0 * r0 * y0 >> 1;
- r0 -= r0 * t0;
- /* r0 += r0 * t0 * ((t0 >> 1) + t0 - 1); /\* Halley, x3 - 1 *\/ */
- ASSERT ((t0 & (GMP_NUMB_MAX >> (GMP_NUMB_BITS - precomputed_bits))) == 0);
+ r0 -= r0 * (y0 * r0 * r0 >> 1); /* -> GMP_NUMB_BITS - 1 */
}
}
if (i) {
- mp_limb_t t4, t3, t2, t1, r1;
+ mp_limb_t t4, t3, t2, t1, t0, r1;
umul_ppmm (t1, t0, r0, r0); /* [t1,t0] <- r^2 */
if (bnb <= GMP_NUMB_BITS) {
@@ -196,7 +171,7 @@
/* r (r^2 y - 1) / 2 - r */
sub_ddmmss(rp[1], rp[0], r1, t4, 0, r0);
} else {
- t0 = (t0 >> 2) | (t1 << (GMP_NUMB_BITS -2));
+ t0 = (t0 >> 2) | (t1 << (GMP_NUMB_BITS -2)) & GMP_NUMB_MAX;
t1 = (t1 >> 2); /* [t1,t0] = r0*r0 >> 2 */
umul_ppmm (t3, t2, y0, t0);
t3 += y0 * t1 + yp[1] * t0;
@@ -222,7 +197,15 @@
for (bn = 2 + (bnb > GMP_NUMB_BITS); --i >= 0;)
{
mp_size_t pbn = bn;
- /* sqr may partially overwrite tp2, but that part is unused here */
+ /* The portion of the result of sqr that overlaps with
+ tp2, is not relevant anyway. */
+
+ /* FIXME: Could maybe be updated:
+ - use wraparound, low quarter known from previous
+ iteration; or
+ - the current r = prev_r + 2^n*d,
+ r^2 = prev_r^2 + 2^{n+1}*prev_r*d + 2^{2n}*d^2 .
+ */
mpn_sqr (tp, rp, bn); /* tp <- r^2 */
bnb = order[i];
diff -r b74bd2160c17 -r ab4a7408e2ce mpz/millerrabin.c
--- a/mpz/millerrabin.c Mon Aug 17 15:09:43 2026 +0200
+++ b/mpz/millerrabin.c Mon Aug 31 20:21:46 2026 +0200
@@ -121,24 +121,34 @@
mpz_millerrabin (mpz_srcptr n, int reps)
{
mpz_t nm, x, y, q;
- mp_bitcnt_t k, l;
+ mp_bitcnt_t k;
+ mp_limb_t h;
int is_prime;
TMP_DECL;
ASSERT (SIZ (n) > 0);
- ASSERT ((SIZ (n) > 1) || (*PTR(n) > 3));
+ h = PTR (n) [SIZ (n) - 1];
+ ASSERT ((SIZ (n) > 1) || (h > 9));
/* Find q and k, where q is odd and n = 1 + 2**k * q. */
k = mpn_scan1 (PTR (n), 1);
- if (GMP_ENABLE_PROTH_TEST && ((l = mpz_sizeinbase (n, 2) - k) <= k)) {
+ if (GMP_ENABLE_PROTH_TEST &&
+ (SIZ (n) <= (k - 1) / (GMP_NUMB_BITS >> 1) + 1) &&
+ ((SIZ (n) < (k - 1) / (GMP_NUMB_BITS >> 1) + 1) ||
+ (k % (GMP_NUMB_BITS >> 1) == 0) ||
+ !(h >> k % (GMP_NUMB_BITS >> 1) * 2))) {
+ /* if (GMP_ENABLE_PROTH_TEST && ((l = mpz_sizeinbase (n, 2) - k) <= k)) { */
ASSERT (k > 1);
/* The number n is a Proth number: 2**k > q. */
+
/* The next search with _kronecker_ would fail with a square n,
- we detect possible squares (of the given form) here. */
- if (((l == k - 1) && /* (2^(k-1)+1)^2 */
+ we detect possible squares (of the given form, > 9) here. */
+ if (((SIZ (n) == (k - 1) / (GMP_NUMB_BITS >> 1) + 1) && /* (2^(k-1)+1)^2 */
+ ((h >> (k - 1) % (GMP_NUMB_BITS >> 1) * 2) == CNST_LIMB (1)) &&
(mpn_scan1 (PTR (n), k + 1) == k - 1 << 1)) ||
- ((l == k - 2) && /* (2^(k-1)-1)^2 */
+ ((SIZ (n) == (k - 2) / (GMP_NUMB_BITS >> 1) + 1) && /* (2^(k-1)-1)^2 */
+ ((h >> (k - 2) % (GMP_NUMB_BITS >> 1) * 2) == CNST_LIMB (3)) &&
(mpz_scan0 (n, k + 1) == k - 1 << 1)))
return 0; /* n is a square => it is a composite */
@@ -172,6 +182,7 @@
/* k == 1, continue */
b += 2; /* FIXME: Loop on primes only. */
} while (b < MIN (ULONG_MAX, GMP_NUMB_MAX));
+ /* FIXME: We should impose a smaller limit. */
}
TMP_MARK;
@@ -211,7 +222,7 @@
0
#endif
#if 64 % GMP_NUMB_BITS != 0
- || SIZ (n) - 64 / GMP_NUMB_BITS == (PTR (n) [64 / GMP_NUMB_BITS] < CNST_LIMB(1) << 64 % GMP_NUMB_BITS)
+ || SIZ (n) - 64 / GMP_NUMB_BITS == (h < CNST_LIMB(1) << 64 % GMP_NUMB_BITS)
#endif
#else
/* Consider numbers that pass the BPSW test as primes, if
@@ -221,14 +232,14 @@
SIZ (n) <= GMP_BPSW_LIMBS_LIMIT ||
#endif
#if GMP_BPSW_BITS_MOD != 0
- SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] <
+ SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (h <
#if GMP_BPSW_BITS_MOD >= GMP_BPSW_BITS_CONST
GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST))
#else
GMP_BPSW_LIMB_CONST >> (GMP_BPSW_BITS_CONST - GMP_BPSW_BITS_MOD))
#endif
#else /* GMP_BPSW_BITS_MOD == 0 */
- SIZ (nm) - GMP_BPSW_LIMBS_LIMIT + 1 == (PTR (nm) [GMP_BPSW_LIMBS_LIMIT - 1] <
+ SIZ (nm) - GMP_BPSW_LIMBS_LIMIT + 1 == (h <
#if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST
GMP_BPSW_LIMB_CONST << (GMP_NUMB_BITS - 1 - GMP_BPSW_BITS_CONST))
#else
diff -r b74bd2160c17 -r ab4a7408e2ce primesieve.c
--- a/primesieve.c Mon Aug 17 15:09:43 2026 +0200
+++ b/primesieve.c Mon Aug 31 20:21:46 2026 +0200
@@ -286,9 +286,7 @@
bits = n_fto_bit(n);
size = bits / GMP_LIMB_BITS + 1;
More information about the gmp-commit
mailing list