[Gmp-commit] /var/hg/gmp: 5 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Mon Jun 29 15:22:08 CEST 2026
details: /var/hg/gmp/rev/a9d5c11bebe6
changeset: 18543:a9d5c11bebe6
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 29 15:16:20 2026 +0200
description:
primesieve.c: Move a branch outside the loop.
details: /var/hg/gmp/rev/9f85154c5fa9
changeset: 18544:9f85154c5fa9
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 29 15:17:12 2026 +0200
description:
mpn/generic/bsqrtinv.c: Higher order iteration step on half-limb.
details: /var/hg/gmp/rev/c8968306a20c
changeset: 18545:c8968306a20c
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 29 15:18:08 2026 +0200
description:
mpz/millerrabin.c: Compute nm only when needed.
details: /var/hg/gmp/rev/0741e4cb376f
changeset: 18546:0741e4cb376f
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 29 15:18:52 2026 +0200
description:
mpz/stronglucas.c: Use mpn_perfect_square_p instead of mpz_
details: /var/hg/gmp/rev/b5f6b754d5bc
changeset: 18547:b5f6b754d5bc
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 29 15:21:15 2026 +0200
description:
ChangeLog
diffstat:
ChangeLog | 5 +++
mpn/generic/bsqrtinv.c | 19 ++++--------
mpz/millerrabin.c | 73 ++++++++++++++++++++++++++++---------------------
mpz/stronglucas.c | 4 +-
primesieve.c | 38 ++++++++++++-------------
5 files changed, 73 insertions(+), 66 deletions(-)
diffs (256 lines):
diff -r 05ed7e126732 -r b5f6b754d5bc ChangeLog
--- a/ChangeLog Fri Jun 26 13:34:23 2026 +0200
+++ b/ChangeLog Mon Jun 29 15:21:15 2026 +0200
@@ -12,6 +12,11 @@
* gen-trialdivtab.c: Likewise.
* bootstrap.c: Likewise.
+ * mpz/stronglucas.c: Use mpn_perfect_square_p instead of mpz_.
+ * mpz/millerrabin.c: Compute nm only when needed.
+ * mpn/generic/bsqrtinv.c: Higher order iteration step on half-limb.
+ * primesieve.c: Move a branch outside the loop.
+
2026-06-19 Marco Bodrato <marco.bodrato at protonmail.com>
* mpn/generic/perfpow.c (is_kth_power): Check only one possible
diff -r 05ed7e126732 -r b5f6b754d5bc mpn/generic/bsqrtinv.c
--- a/mpn/generic/bsqrtinv.c Fri Jun 26 13:34:23 2026 +0200
+++ b/mpn/generic/bsqrtinv.c Mon Jun 29 15:21:15 2026 +0200
@@ -144,10 +144,11 @@
r0 -= r0 * t0;
const mp_bitcnt_t precomputed_bits = 19;
#else /* GMP_NUMB_BITS >= 19 * 2 - 2 */
- r0 += r0 * t0 * ((t0 >> 1) + t0 - 1); /* Halley -> 29*/
-#if GMP_NUMB_BITS < 29 * 3 - 2
- const mp_bitcnt_t precomputed_bits = 29;
-#else /* GMP_NUMB_BITS >= 29 * 3 - 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.
#endif
#endif
@@ -164,22 +165,16 @@
mp_limb_t yh = (y0 >> 2) + (yp[1] << (GMP_NUMB_BITS - 2));
mp_limb_t yrrm1d4 = y0 * r0sqm1 + yh; /* (r0*r0*y0-1) >> 2 */
ASSERT ((yrrm1d4 & (GMP_NUMB_MAX >> (GMP_NUMB_BITS - precomputed_bits + 1))) == 0);
-#if GMP_NUMB_BITS < 19 * 2 - 2
mp_limb_t rt = r0h - r0 * yrrm1d4; /* (r - r*(r0*r0*y0-1)/2) >>1 */
-#else /* GMP_NUMB_BITS >= 19 * 2 - 2 */
- mp_limb_t rt = r0h + r0 * yrrm1d4 * (yrrm1d4 * 3 - 1); /* Halley, x3 - 1 */
-#endif
+ /* mp_limb_t rt = r0h + r0 * yrrm1d4 * (yrrm1d4 * 3 - 1); /\* Halley, x3 - 1 *\/ */
r0 = (rt << 1) ^ ((rt & GMP_LIMB_HIGHBIT) ? GMP_NUMB_MAX : CNST_LIMB(1));
#ifdef BSQRTINV_RP_NOT_ZEROED
rp[1] = 0;
#endif
} else {
t0 = r0 * r0 * y0 >> 1;
-#if GMP_NUMB_BITS < 19 * 2 - 2
r0 -= r0 * t0;
-#else /* GMP_NUMB_BITS >= 19 * 2 - 2 */
- r0 += r0 * t0 * ((t0 >> 1) + t0 - 1); /* Halley, x3 - 1 */
-#endif
+ /* r0 += r0 * t0 * ((t0 >> 1) + t0 - 1); /\* Halley, x3 - 1 *\/ */
ASSERT ((t0 & (GMP_NUMB_MAX >> (GMP_NUMB_BITS - precomputed_bits))) == 0);
}
}
diff -r 05ed7e126732 -r b5f6b754d5bc mpz/millerrabin.c
--- a/mpz/millerrabin.c Fri Jun 26 13:34:23 2026 +0200
+++ b/mpz/millerrabin.c Mon Jun 29 15:21:15 2026 +0200
@@ -91,7 +91,7 @@
/* Performs a Miller-Rabin test, on the number n, with base x.
* The value q is the odd number such that (q<<k) + 1 = n.
*
- * The variable y is overwritten, it's only role is to reuse the same
+ * The variable y is overwritten, its only role is to reuse the same
* temp variable.
*/
@@ -123,24 +123,35 @@
TMP_MARK;
ASSERT (SIZ (n) > 0);
- MPZ_TMP_INIT (nm, SIZ (n) + 1);
- mpz_tdiv_q_2exp (nm, n, 1);
MPZ_TMP_INIT (x, SIZ (n) + 1);
MPZ_TMP_INIT (y, 2 * SIZ (n)); /* mpz_powm_ui needs excessive memory!!! */
MPZ_TMP_INIT (q, SIZ (n));
/* Find q and k, where q is odd and n = 1 + 2**k * q. */
- k = mpn_scan1 (PTR (nm), 0);
- mpz_tdiv_q_2exp (q, nm, k);
- ++k;
+ k = mpn_scan1 (PTR (n), 1);
+ mpz_tdiv_q_2exp (q, n, k);
/* BPSW test */
mpz_set_ui (x, 2);
is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y);
+ /* Consider numbers up to 41*2^46 that pass the BPSW test as primes.
+ This implementation was tested up to 289*10^13 > 2^51+2^49+2^46 */
+ /* 2^5 < 41 = 0b101001 < 2^6 */
+#define GMP_BPSW_LIMB_CONST CNST_LIMB(41)
+#define GMP_BPSW_BITS_CONST (LOG2C(41) - 1)
+#define GMP_BPSW_BITS_LIMIT (46 + GMP_BPSW_BITS_CONST)
+
+#define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS)
+#define GMP_BPSW_BITS_MOD (GMP_BPSW_BITS_LIMIT % GMP_NUMB_BITS)
+
if (is_prime)
{
+#if !GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS && GMP_BPSW_BITS_MOD == 0
+ MPZ_TMP_INIT (nm, SIZ (n) + 1);
+ mpz_tdiv_q_2exp (nm, n, 1);
+#endif
if (
#if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
/* Consider numbers up to 2^64 that pass the BPSW test as primes. */
@@ -153,44 +164,36 @@
|| SIZ (n) - 64 / GMP_NUMB_BITS == (PTR (n) [64 / GMP_NUMB_BITS] < CNST_LIMB(1) << 64 % GMP_NUMB_BITS)
#endif
#else
- /* Consider numbers up to 41*2^46 that pass the BPSW test as primes.
- This implementation was tested up to 289*10^13 > 2^51+2^49+2^46 */
- /* 2^5 < 41 = 0b101001 < 2^6 */
-#define GMP_BPSW_LIMB_CONST CNST_LIMB(41)
-#define GMP_BPSW_BITS_CONST (LOG2C(41) - 1)
-#define GMP_BPSW_BITS_LIMIT (46 + GMP_BPSW_BITS_CONST)
-
-#define GMP_BPSW_LIMBS_LIMIT (GMP_BPSW_BITS_LIMIT / GMP_NUMB_BITS)
-#define GMP_BPSW_BITS_MOD (GMP_BPSW_BITS_LIMIT % GMP_NUMB_BITS)
-
+ /* Consider numbers that pass the BPSW test as primes, if
+ they are in the range where this implementation of the
+ test has been fully tested. */
#if GMP_NUMB_BITS <= GMP_BPSW_BITS_LIMIT
- SIZ (n) <= GMP_BPSW_LIMBS_LIMIT
-#else
- 0
+ 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] <
#if GMP_BPSW_BITS_MOD >= GMP_BPSW_BITS_CONST
- || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST))
+ GMP_BPSW_LIMB_CONST << (GMP_BPSW_BITS_MOD - GMP_BPSW_BITS_CONST))
#else
-#if GMP_BPSW_BITS_MOD != 0
- || SIZ (n) - GMP_BPSW_LIMBS_LIMIT == (PTR (n) [GMP_BPSW_LIMBS_LIMIT] < GMP_BPSW_LIMB_CONST >> (GMP_BPSW_BITS_CONST - GMP_BPSW_BITS_MOD))
+ 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] <
+#if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST
+ GMP_BPSW_LIMB_CONST << (GMP_NUMB_BITS - 1 - GMP_BPSW_BITS_CONST))
#else
-#if GMP_NUMB_BITS > GMP_BPSW_BITS_CONST
- || SIZ (nm) - GMP_BPSW_LIMBS_LIMIT + 1 == (PTR (nm) [GMP_BPSW_LIMBS_LIMIT - 1] < GMP_BPSW_LIMB_CONST << (GMP_NUMB_BITS - 1 - GMP_BPSW_BITS_CONST))
+ GMP_BPSW_LIMB_CONST >> (1 + GMP_BPSW_BITS_CONST - GMP_NUMB_BITS))
#endif
#endif
#endif
-
-#undef GMP_BPSW_BITS_LIMIT
-#undef GMP_BPSW_LIMB_CONST
-#undef GMP_BPSW_BITS_CONST
-#undef GMP_BPSW_LIMBS_LIMIT
-#undef GMP_BPSW_BITS_MOD
-
-#endif
)
is_prime = 2;
else
{
+#if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS || GMP_BPSW_BITS_MOD != 0
+ MPZ_TMP_INIT (nm, SIZ (n) + 1);
+ mpz_tdiv_q_2exp (nm, n, 1);
+#endif
reps -= 24;
if (reps > 0)
{
@@ -217,3 +220,9 @@
TMP_FREE;
return is_prime;
}
+
+#undef GMP_BPSW_BITS_LIMIT
+#undef GMP_BPSW_LIMB_CONST
+#undef GMP_BPSW_BITS_CONST
+#undef GMP_BPSW_LIMBS_LIMIT
+#undef GMP_BPSW_BITS_MOD
diff -r 05ed7e126732 -r b5f6b754d5bc mpz/stronglucas.c
--- a/mpz/stronglucas.c Fri Jun 26 13:34:23 2026 +0200
+++ b/mpz/stronglucas.c Mon Jun 29 15:21:15 2026 +0200
@@ -135,8 +135,8 @@
int jac;
/* n is odd, to possibly be a square, n % 8 = 1 is needed. */
- if (((*PTR (n) & 6) == 0) && UNLIKELY (mpz_perfect_square_p (n)))
- return 0; /* A square is composite. */
+ if (((*PTR (n) & 6) == 0) && UNLIKELY (mpn_perfect_square_p (PTR (n), SIZ (n))))
+ return 0; /* A square is not prime. */
/* Check Ds up to square root (in case, n is prime)
or avoid overflows */
diff -r 05ed7e126732 -r b5f6b754d5bc primesieve.c
--- a/primesieve.c Fri Jun 26 13:34:23 2026 +0200
+++ b/primesieve.c Mon Jun 29 15:21:15 2026 +0200
@@ -154,18 +154,17 @@
offset %= 13 * 7 * 2;
SET_OFF2 (m21, m22, m23, SIEVE_2MSK1, SIEVE_2MSK2, SIEVE_2MSKT, offset, 13 * 7 * 2);
}
- /* THINK: Consider handling odd values of 'limbs' outside the loop,
- to have a single exit condition. */
- do {
+ for (; (limbs -= 2) >= 0; )
+ {
+ bit_array[0] = m11 | m21;
+ ROTATE1 (m11, m12, 11 * 5 * 2);
+ bit_array[1] = m11 | m22;
+ bit_array += 2;
+ ROTATE1 (m11, m12, 11 * 5 * 2);
+ ROTATE2 (m21, m22, m23, 13 * 7 * 2);
+ }
+ if (limbs == -1) /* limbs+2 != 0, limbs&1 != 0 */
bit_array[0] = m11 | m21;
- if (--limbs == 0)
- break;
- ROTATE1 (m11, m12, 11 * 5 * 2);
- bit_array[1] = m11 | m22;
- bit_array += 2;
- ROTATE1 (m11, m12, 11 * 5 * 2);
- ROTATE2 (m21, m22, m23, 13 * 7 * 2);
- } while (--limbs != 0);
return n_cto_bit (13 + 1);
#else
#ifdef SIEVE_MASK2
@@ -175,16 +174,15 @@
offset %= 7 * 5 * 2;
SET_OFF2 (mask, mask2, tail, SIEVE_MASK1, SIEVE_MASK2, SIEVE_MASKT, offset, 7 * 5 * 2);
}
- /* THINK: Consider handling odd values of 'limbs' outside the loop,
- to have a single exit condition. */
- do {
+ for (; (limbs -= 2) >= 0; )
+ {
+ bit_array[0] = mask;
+ bit_array[1] = mask2;
+ bit_array += 2;
+ ROTATE2 (mask, mask2, tail, 7 * 5 * 2);
+ }
+ if (limbs == -1)
bit_array[0] = mask;
- if (--limbs == 0)
- break;
- bit_array[1] = mask2;
- bit_array += 2;
- ROTATE2 (mask, mask2, tail, 7 * 5 * 2);
- } while (--limbs != 0);
return n_cto_bit (7 + 1);
#else
MPN_FILL (bit_array, limbs, CNST_LIMB(0));
More information about the gmp-commit
mailing list