[Gmp-commit] /var/hg/gmp: 9 new changesets

mercurial at gmplib.org mercurial at gmplib.org
Mon Aug 17 15:09:54 CEST 2026


details:   /var/hg/gmp/rev/fb00206f046f
changeset: 18575:fb00206f046f
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 12:41:00 2026 +0200
description:
mpz/millerrabin.c (millerrabin): Slightly smaller mempry usage.

details:   /var/hg/gmp/rev/b33c2330c419
changeset: 18576:b33c2330c419
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 12:51:09 2026 +0200
description:
mpz/millerrabin.c: Special (simple) test for Proth numbers.

details:   /var/hg/gmp/rev/801ea8481774
changeset: 18577:801ea8481774
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 14:13:29 2026 +0200
description:
tests/mpz/t-pprime_p.c (check_one): One more argument,
to know in which section a test failed.

details:   /var/hg/gmp/rev/fc59cf70d6fa
changeset: 18578:fc59cf70d6fa
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 14:24:15 2026 +0200
description:
tests/mpz/t-pprime_p.c (check_fermat_mersenne): Check squared numbers,
they are (non-prime) Proth numbers.

details:   /var/hg/gmp/rev/0de55c998d8a
changeset: 18579:0de55c998d8a
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 14:30:04 2026 +0200
description:
tests/mpz/t-pprime_p.c (check_proth_fixed): New tests.

details:   /var/hg/gmp/rev/0090896cc3cf
changeset: 18580:0090896cc3cf
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 14:37:11 2026 +0200
description:
tests/mpz/t-pprime_p.c: Declare rands in main(), to pass it to functions.

details:   /var/hg/gmp/rev/810cc6ebcb92
changeset: 18581:810cc6ebcb92
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 14:45:33 2026 +0200
description:
tests/mpz/t-pprime_p.c (check_proth): New randomized test.

details:   /var/hg/gmp/rev/e66521cb5d68
changeset: 18582:e66521cb5d68
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 14:50:40 2026 +0200
description:
Copyright years

details:   /var/hg/gmp/rev/b74bd2160c17
changeset: 18583:b74bd2160c17
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Mon Aug 17 15:09:43 2026 +0200
description:
ChangeLog

diffstat:

 ChangeLog              |    5 +
 mpz/millerrabin.c      |   68 +++++++++++++++++--
 tests/mpz/t-pprime_p.c |  164 +++++++++++++++++++++++++++++++++++++++++++++---
 3 files changed, 215 insertions(+), 22 deletions(-)

diffs (truncated from 401 to 300 lines):

diff -r b4d0ee7d9773 -r b74bd2160c17 ChangeLog
--- a/ChangeLog	Sun Jul 26 11:54:09 2026 +0200
+++ b/ChangeLog	Mon Aug 17 15:09:43 2026 +0200
@@ -1,3 +1,8 @@
+2026-08-17  Marco Bodrato <bodrato at anjara.org>
+
+	* mpz/millerrabin.c: Special (simple) test for Proth numbers.
+	* tests/mpz/t-pprime_p.c: Exercise the new test.
+
 2026-07-20  Marco Bodrato <marco.bodrato at protonmail.com>
 
 	* mpn/generic/bsqrt.c: Use a last iteration to increase precision
diff -r b4d0ee7d9773 -r b74bd2160c17 mpz/millerrabin.c
--- a/mpz/millerrabin.c	Sun Jul 26 11:54:09 2026 +0200
+++ b/mpz/millerrabin.c	Mon Aug 17 15:09:43 2026 +0200
@@ -19,8 +19,8 @@
    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
    FUTURE GNU MP RELEASES.
 
-Copyright 1991, 1993, 1994, 1996-2002, 2005, 2014, 2018-2022, 2024 Free
-Software Foundation, Inc.
+Copyright 1991, 1993, 1994, 1996-2002, 2005, 2014, 2018-2022, 2024,
+2026 Free Software Foundation, Inc.
 
 Contributed by John Amanatides.
 Changed to "BPSW, then Miller Rabin if required" by Marco Bodrato.
@@ -57,6 +57,10 @@
 #define GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS 0
 #endif
 
+#ifndef GMP_ENABLE_PROTH_TEST
+#define GMP_ENABLE_PROTH_TEST 1
+#endif
+
 static int
 mod_eq_m1 (mpz_srcptr x, mpz_srcptr m)
 {
@@ -101,7 +105,7 @@
 {
   mpz_powm (y, x, q, n);
 
-  if (mpz_cmp_ui (y, 1L) == 0 || mod_eq_m1 (y, n))
+  if (((SIZ (y) == 1) & (*PTR (y) == 1)) || mod_eq_m1 (y, n))
     return 1;
 
   for (mp_bitcnt_t i = 1; i < k; ++i)
@@ -117,19 +121,65 @@
 mpz_millerrabin (mpz_srcptr n, int reps)
 {
   mpz_t nm, x, y, q;
-  mp_bitcnt_t k;
+  mp_bitcnt_t k, l;
   int is_prime;
   TMP_DECL;
-  TMP_MARK;
 
   ASSERT (SIZ (n) > 0);
+  ASSERT ((SIZ (n) > 1) || (*PTR(n) > 3));
+
+  /* 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)) {
+    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 */
+	 (mpn_scan1 (PTR (n), k + 1) == k - 1 << 1)) ||
+	((l == k - 2) && /* (2^(k-1)-1)^2 */
+	 (mpz_scan0 (n, k + 1) == k - 1 << 1)))
+      return 0; /* n is a square => it is a composite */
+
+    unsigned long b = 3;
+    do {
+      int knb = mpz_kronecker_ui (n, b);
+      if (knb <= 0) {
+	if (knb == 0) /* knb == 0, gcd(b,n) != 1. */
+	  is_prime = 0; /* Composite. */
+	else {
+	  mp_limb_t const xp[1] = {b};
+	  TMP_MARK;
+
+	  MPZ_TMP_INIT (y, SIZ (n));
+	  MPZ_TMP_INIT (q, SIZ (n));
+
+	  mpz_tdiv_q_2exp (q, n, 1);
+
+	  mpz_roinit_n (x, xp, 1);
+	  mpz_powm (y, x, q, n);
+	  /* n is prime if and only if x^{(n-1)/2} (mod n) = n - 1.
+	     In case, n is surely prime, we can return 2. */
+	  is_prime = mod_eq_m1 (y, n) << 1;
+
+	  TMP_FREE;
+	}
+
+	return is_prime;
+      }
+
+      /* k == 1, continue */
+      b += 2; /* FIXME: Loop on primes only. */
+    } while (b < MIN (ULONG_MAX, GMP_NUMB_MAX));
+  }
+
+  TMP_MARK;
 
   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 (n), 1);
   mpz_tdiv_q_2exp (q, n, k);
 
   /* BPSW test */
@@ -149,7 +199,7 @@
   if (is_prime)
     {
 #if !GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS && GMP_BPSW_BITS_MOD == 0
-      MPZ_TMP_INIT (nm, SIZ (n) + 1);
+      MPZ_TMP_INIT (nm, SIZ (n));
       mpz_tdiv_q_2exp (nm, n, 1);
 #endif
       if (
@@ -191,7 +241,7 @@
       else
 	{
 #if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS || GMP_BPSW_BITS_MOD != 0
-	  MPZ_TMP_INIT (nm, SIZ (n) + 1);
+	  MPZ_TMP_INIT (nm, SIZ (n));
 	  mpz_tdiv_q_2exp (nm, n, 1);
 #endif
 	  reps -= 24;
diff -r b4d0ee7d9773 -r b74bd2160c17 tests/mpz/t-pprime_p.c
--- a/tests/mpz/t-pprime_p.c	Sun Jul 26 11:54:09 2026 +0200
+++ b/tests/mpz/t-pprime_p.c	Mon Aug 17 15:09:43 2026 +0200
@@ -1,6 +1,6 @@
 /* Exercise mpz_probab_prime_p.
 
-Copyright 2002, 2018-2019, 2022 Free Software Foundation, Inc.
+Copyright 2002, 2018-2019, 2022, 2026 Free Software Foundation, Inc.
 
 This file is part of the GNU MP Library test suite.
 
@@ -22,6 +22,9 @@
 #include "gmp-impl.h"
 #include "tests.h"
 
+#ifndef GMP_ENABLE_PROTH_TEST
+#define GMP_ENABLE_PROTH_TEST 1
+#endif
 
 /* Enhancements:
 
@@ -48,7 +51,7 @@
 }
 
 void
-check_one (mpz_srcptr n, int want)
+check_one (mpz_srcptr n, int want, char a)
 {
   int  got;
 
@@ -57,7 +60,7 @@
   /* "definitely prime" (2) is fine if we only wanted "probably prime" (1) */
   if ((got != want) && (got != want * 2))
     {
-      printf ("mpz_probab_prime_p\n");
+      printf ("mpz_probab_prime_p (%c)\n", a);
       mpz_trace ("  n    ", n);
       printf    ("  got =%d", got);
       printf    ("  want=%d", want);
@@ -68,9 +71,9 @@
 void
 check_pn (mpz_ptr n, int want)
 {
-  check_one (n, want);
+  check_one (n, want, '+');
   mpz_neg (n, n);
-  check_one (n, want);
+  check_one (n, want, '-');
 }
 
 /* expect certainty for small n */
@@ -92,12 +95,133 @@
 }
 
 void
-check_composites (int count)
+check_proth_fixed ()
+{
+  static const struct {
+    char*        k;
+    mp_bitcnt_t  n;
+    int       want;
+  } data [] = {
+    {"18975", 16, 1}, /* prime, b=61 */
+    {"642497427", 34, 1}, /* prime, b=107 */
+    {"430341165", 32, 0}, /* composite, b=113 */
+    {"3432101253", 32, 1}, /* prime, b=113 */
+    {"67067655", 64, 1}, /* prime, b=103 */
+    {"117377265", 64, 0}, /* composite, b=113 */
+    {"4654135305", 63, 1}, /* prime, b=107 */
+    {"3836987091", 63, 0}, /* composite, b=131 */
+    {"11838137235", 63, 1}, /* prime, b=109 */
+    {"6269603373", 64, 0}, /* composite, b=137 */
+    {"1300166691", 67, 1}, /* prime, b=113 */
+    {"11257223805", 65, 0}, /* composite, b=151 */
+    {"71459624811", 63, 1}, /* prime, b=137 */
+    {"59242853985", 64, 0}, /* composite, b=157 */
+    {"25739539989", 65, 1}, /* prime, b=139 */
+    {"59242853985", 65, 0}, /* composite, b=179 */
+    {"6148401", 103, 0}, /* composite, b=101 */
+    {"69646689", 101, 1}, /* prime, b=97 */
+  };
+  mpz_t n;
+  mpz_init (n);
+
+  for (int i = 0; i < numberof (data); ++i) {
+    mpz_set_str (n, data[i].k, 10);
+    mpz_mul_2exp (n, n, data[i].n);
+    int want = data[i].want << GMP_ENABLE_PROTH_TEST;
+    mpz_add_ui (n, n, 1);
+
+    check_one (n, want, 'h');
+  }
+
+  mpz_clear (n);
+}
+
+void
+check_proth (gmp_randstate_ptr rands, int count)
+{
+  /* Exponents such that k*2^n + 1 is prime, with
+     odd 2 < k < 18, and k < 2^n < 2^(2^14 + 10).
+     Sequences where recomputed, but checked with OEIS. */
+#define DATA_END 32767
+  /* 3 * 2^n + 1, OEIS A002253 */
+  static const unsigned data3[] = {
+    2, 5, 6, 8, 12, 18, 30, 36, 41, 66, 189, 201, 209, 276, 353,
+    408, 438, 534, 2208, 2816, 3168, 3189, 3912, DATA_END};
+  /* 5 * 2^n + 1, OEIS A002254 */
+  static const unsigned data5[] = {
+    3, 7, 13, 15, 25, 39, 55, 75, 85, 127, 1947,
+    3313, 4687, 5947, 13165, DATA_END};
+  /* 7 * 2^n + 1, OEIS A002255 */
+  static const unsigned data7[] = {
+    4, 6, 14, 20, 26, 50, 52, 92, 120, 174, 180, 190, 290, 320,
+    390, 432, 616, 830, 1804, 2256, 6614, 13496, 15494, DATA_END};
+  /* 9 * 2^n + 1, OEIS A002256 */
+  static const unsigned data9[] = {
+    6, 7, 11, 14, 17, 33, 42, 43, 63, 65, 67, 81, 134, 162, 206,
+    211, 366, 663, 782, 1305, 1411, 1494, 2297, 2826, 3230, 3354,
+    3417, 3690, 4842, 5802, 6937, 7967, 9431, 13903, DATA_END};
+  /* 11 * 2^n + 1, OEIS A002261 */
+  static const unsigned data11[] = {
+    5, 7, 19, 21, 43, 81, 125, 127, 209, 211, 3225,
+    4543, 10179, 15329, DATA_END};
+  /* 13 * 2^n + 1, OEIS A002257 */
+  static const unsigned data13[] = {
+    8, 10, 20, 28, 82, 188, 308, 316, 1000, DATA_END};
+  /* 15 * 2^n + 1, OEIS A002258 */
+  static const unsigned data15[] = {
+    4, 9, 10, 12, 27, 37, 38, 44, 48, 78, 112, 168, 229, 297,
+    339, 517, 522, 654, 900, 1518, 2808, 2875, 3128, 3888,
+    4410, 6804, 7050, 7392, DATA_END};
+  /* 17 * 2^n + 1, OEIS A002259 */
+  static const unsigned data17[] = {
+    15, 27, 51, 147, 243, 267, 347, 471, 747, 2163,
+    3087, 5355, 6539, 7311, DATA_END};
+  static const unsigned *data[8] = {
+    data3, data5, data7, data9, data11, data13, data15, data17
+  };
+
+  mpz_t n;
+  mpz_init (n);
+
+  unsigned bits = 10;
+  while ((count >> bits) > 2)
+    {
+      count >>= 1;
+      ++bits;
+      if (bits > 13)
+	break;
+    }
+
+  for (int i = count; i != 0; --i) {
+    unsigned long k = gmp_urandomb_ui (rands, 3); /* 0..7 */
+    const unsigned *exponents = data[k];
+    mpz_set_ui (n, 2 * k + 3); /* 3..17, odd */
+
+    unsigned shift = gmp_urandomb_ui (rands, bits) + mpz_sizeinbase (n, 2);
+    mpz_mul_2exp (n, n, shift);
+    mpz_add_ui (n, n, 1);
+
+    int want;
+    for (unsigned j = 0; ; ++j)
+      if (shift <= exponents[j])
+	{


More information about the gmp-commit mailing list