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

mercurial at gmplib.org mercurial at gmplib.org
Tue Aug 13 17:48:34 UTC 2019


details:   /var/hg/gmp/rev/01835ae182d1
changeset: 17798:01835ae182d1
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:32:39 2019 +0200
description:
Update mpz_millerrabin documentation, by Seth Troisi

details:   /var/hg/gmp/rev/35602cbb22db
changeset: 17799:35602cbb22db
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:35:18 2019 +0200
description:
mpn/x86_64/bd2/gcd_11.asm: Micro-optimisation

details:   /var/hg/gmp/rev/d8d558260918
changeset: 17800:d8d558260918
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:37:50 2019 +0200
description:
doc/gmp.texi: Further update in mpz_millerrabin documentation

details:   /var/hg/gmp/rev/75ec05102858
changeset: 17801:75ec05102858
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:38:32 2019 +0200
description:
mini-gmp/mini-gmp.c: Silence a couple of warnings

details:   /var/hg/gmp/rev/8e5a2741ed59
changeset: 17802:8e5a2741ed59
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:40:10 2019 +0200
description:
tests/misc.c: silence a warning

details:   /var/hg/gmp/rev/08ec43bdeb2d
changeset: 17803:08ec43bdeb2d
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:42:55 2019 +0200
description:
tests/mpz/t-pprime_p.c: Test one more prime

details:   /var/hg/gmp/rev/d23f58a30336
changeset: 17804:d23f58a30336
user:      Marco Bodrato <bodrato at mail.dm.unipi.it>
date:      Tue Aug 13 19:43:43 2019 +0200
description:
mpz/millerrabin.c: Return 2 for surely primes (BPSW checked)

diffstat:

 doc/gmp.texi              |  16 ++++++---
 mini-gmp/mini-gmp.c       |   8 ++--
 mpn/x86_64/bd2/gcd_11.asm |  21 ++++++-------
 mpz/millerrabin.c         |  72 +++++++++++++++++++++++++++++++++++++++++++++-
 tests/misc.c              |   2 +-
 tests/mpz/t-pprime_p.c    |   2 +-
 6 files changed, 96 insertions(+), 25 deletions(-)

diffs (236 lines):

diff -r b252c7e4f9b6 -r d23f58a30336 doc/gmp.texi
--- a/doc/gmp.texi	Thu Aug 08 16:29:36 2019 +0200
+++ b/doc/gmp.texi	Tue Aug 13 19:43:43 2019 +0200
@@ -3548,12 +3548,16 @@
 return 1 if @var{n} is probably prime (without being certain), or return 0 if
 @var{n} is definitely non-prime.
 
-This function performs some trial divisions, then @var{reps} Miller-Rabin
-probabilistic primality tests.  A higher @var{reps} value will reduce the
-chances of a non-prime being identified as ``probably prime''.  A composite
-number will be identified as a prime with a probability of less than
- at m{4^{-reps},4^(- at var{reps})}.  Reasonable values of @var{reps} are between 15
-and 50.
+This function performs some trial divisions, a Baillie-PSW probable prime
+test, then @var{reps-24} Miller-Rabin probabilistic primality tests.  A
+higher @var{reps} value will reduce the chances of a non-prime being
+identified as ``probably prime''.  A composite number will be identified as a
+prime with an asymptotic probability of less than @m{4^{-reps},4^(- at var{reps})}.
+Reasonable values of @var{reps} are between 15 and 50.
+
+GMP versions up to and including 6.1.2 did not use the Baillie-PSW
+primality test. In those older versions of GMP, this function performed
+ at var{reps} Miller-Rabin tests.
 @end deftypefun
 
 @deftypefun void mpz_nextprime (mpz_t @var{rop}, const mpz_t @var{op})
diff -r b252c7e4f9b6 -r d23f58a30336 mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c	Thu Aug 08 16:29:36 2019 +0200
+++ b/mini-gmp/mini-gmp.c	Tue Aug 13 19:43:43 2019 +0200
@@ -295,7 +295,7 @@
 }
 
 static void *
-gmp_default_realloc (void *old, size_t old_size, size_t new_size)
+gmp_default_realloc (void *old, size_t unused_old_size, size_t new_size)
 {
   void * p;
 
@@ -308,7 +308,7 @@
 }
 
 static void
-gmp_default_free (void *p, size_t size)
+gmp_default_free (void *p, size_t unused_size)
 {
   free (p);
 }
@@ -1595,7 +1595,7 @@
       int LOCAL_GMP_LIMB_BITS = GMP_LIMB_BITS;
       unsigned long r = 0;
       mp_size_t n = GMP_ABS (u->_mp_size);
-      n = GMP_MIN (n, 1 + (GMP_ULONG_BITS - 1) / GMP_LIMB_BITS);
+      n = GMP_MIN (n, 1 + (mp_size_t) (GMP_ULONG_BITS - 1) / GMP_LIMB_BITS);
       while (--n >= 0)
 	r = (r << LOCAL_GMP_LIMB_BITS) + u->_mp_d[n];
       return r;
@@ -3499,7 +3499,7 @@
   b0 = mpz_scan0 (n, 0);
 
   /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
-  Q = (D & 2) ? (D >> 2) + 1 : -(long) (D >> 2);
+  Q = (D & 2) ? (long) (D >> 2) + 1 : -(long) (D >> 2);
 
   if (! gmp_lucas_mod (V, Qk, Q, b0, n))	/* If Ud != 0 */
     while (V->_mp_size != 0 && --b0 != 0)	/* while Vk != 0 */
diff -r b252c7e4f9b6 -r d23f58a30336 mpn/x86_64/bd2/gcd_11.asm
--- a/mpn/x86_64/bd2/gcd_11.asm	Thu Aug 08 16:29:36 2019 +0200
+++ b/mpn/x86_64/bd2/gcd_11.asm	Tue Aug 13 19:43:43 2019 +0200
@@ -72,22 +72,21 @@
 	ALIGN(16)
 PROLOGUE(mpn_gcd_11)
 	FUNC_ENTRY(2)
-	mov	v0, %r10	C
-	sub	u0, %r10	C
+	mov	v0, %rax	C
+	sub	u0, v0		C
 	jz	L(end)		C
 
 	ALIGN(16)		C              K10 BD1 BD2 ZEN CNR NHM SBR
-L(top):	rep;bsf	%r10, %rcx	C tzcnt!       3   3   3   2   6   5   5
+L(top):	rep;bsf	v0, %rcx	C tzcnt!       3   3   3   2   6   5   5	(TODO)
 	mov	u0, %r9		C              2   2   2   2   3   3   4
-	sub	v0, u0		C              2   2   2   2   4   3   4
-	cmovc	%r10, u0	C if x-y < 0   0,3 0,3 0,3 0,3 0,6 0,5 0,5
-	cmovc	%r9, v0		C use x,y-x    0,3 0,3 0,3 0,3 2,8 1,7 1,7
+	sub	%rax, u0	C              2   2   2   2   4   3   4	(TODO)
+	cmovc	v0, u0		C if x-y < 0   0,3 0,3 0,3 0,3 0,6 0,5 0,5	(TODO)
+	cmovc	%r9, %rax	C use x,y-x    0,3 0,3 0,3 0,3 2,8 1,7 1,7	(TODO)
 	shr	R8(%rcx), u0	C              1,7 1,6 1,5 1,4 2,8 2,8 2,8
-	mov	v0, %r10	C              1   1   1   1   4   3   3
-	sub	u0, %r10	C              2   2   2   1   5   4   4
-	jnz	L(top)		C
+	mov	%rax, v0	C              1   1   1   1   4   3   3	(TODO)
+	sub	u0, v0		C              2   2   2   1   5   4   4	(TODO)
+	jnz	L(top)		C	(TODO: registers was changed, cycles not updated)
 
-L(end):	mov	v0, %rax
-	FUNC_EXIT()
+L(end):	FUNC_EXIT()
 	ret
 EPILOGUE()
diff -r b252c7e4f9b6 -r d23f58a30336 mpz/millerrabin.c
--- a/mpz/millerrabin.c	Thu Aug 08 16:29:36 2019 +0200
+++ b/mpz/millerrabin.c	Tue Aug 13 19:43:43 2019 +0200
@@ -5,6 +5,16 @@
    reps is the number of internal passes of the probabilistic algorithm.  Knuth
    indicates that 25 passes are reasonable.
 
+   With the current implementation, the first 24 MR-tests are substituted by a
+   Baillie-PSW probable prime test.
+
+   This implementation the Baillie-PSW test was checked up to 37*2^45,
+   for smaller values no MR-test is performed, regardless of reps, and
+   2 ("surely prime") is returned if the number was not proved composite.
+
+   If GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS is defined as non-zero,
+   the code assumes that the Baillie-PSW test was checked up to 2^64.
+
    THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
    CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
    FUTURE GNU MP RELEASES.
@@ -42,6 +52,10 @@
 
 #include "gmp-impl.h"
 
+#ifndef GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
+#define GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS 0
+#endif
+
 static int millerrabin (mpz_srcptr,
 			mpz_ptr, mpz_ptr,
 			mpz_srcptr, unsigned long int);
@@ -56,6 +70,7 @@
   TMP_DECL;
   TMP_MARK;
 
+  ASSERT (SIZ (n) > 0);
   MPZ_TMP_INIT (nm, SIZ (n) + 1);
   mpz_tdiv_q_2exp (nm, n, 1);
 
@@ -72,8 +87,60 @@
   mpz_set_ui (x, 2);
   is_prime = millerrabin (n, x, y, q, k) && mpz_stronglucas (n, x, y);
 
+  if (is_prime)
+    {
+      if (
+#if GMP_BPSW_NOFALSEPOSITIVES_UPTO_64BITS
+	  /* Consider numbers up to 2^64 that pass the BPSW test as primes. */
+#if GMP_NUMB_BITS <= 64
+	  SIZ (n) <= 64 / GMP_NUMB_BITS
+#else
+	  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)
+#endif
+#else
+	  /* Consider numbers up to 37*2^45 that pass the BPSW test as primes.
+	     This implementation was tested up to 37*2^45 = 2^50+2^47+2^45 */
+	  /* 2^5 < 37 = 0b100101 < 2^6 */
+#define GMP_BPSW_LIMB_CONST CNST_LIMB(37)
+#define GMP_BPSW_BITS_CONST (LOG2C(37) - 1)
+#define GMP_BPSW_BITS_LIMIT (45 + 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 GMP_NUMB_BITS <=  GMP_BPSW_BITS_LIMIT
+	  SIZ (n) <= GMP_BPSW_LIMBS_LIMIT
+#else
+	  0
+#endif
+#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))
+#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))
+#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))
+#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
+	{
   reps -= 24;
-  if (reps > 0 && is_prime)
+  if (reps > 0)
     {
       /* (n-5)/2 */
       mpz_sub_ui (nm, nm, 2L);
@@ -92,7 +159,8 @@
 
       gmp_randclear (rstate);
     }
-
+	}
+    }
   TMP_FREE;
   return is_prime;
 }
diff -r b252c7e4f9b6 -r d23f58a30336 tests/misc.c
--- a/tests/misc.c	Thu Aug 08 16:29:36 2019 +0200
+++ b/tests/misc.c	Tue Aug 13 19:43:43 2019 +0200
@@ -209,7 +209,7 @@
 {
   char  *s;
   for (s = s_orig; *s != '\0'; s++)
-    if (isascii (*s))
+    if (islower (*s))
       *s = toupper (*s);
   return s_orig;
 }
diff -r b252c7e4f9b6 -r d23f58a30336 tests/mpz/t-pprime_p.c
--- a/tests/mpz/t-pprime_p.c	Thu Aug 08 16:29:36 2019 +0200
+++ b/tests/mpz/t-pprime_p.c	Tue Aug 13 19:43:43 2019 +0200
@@ -147,7 +147,7 @@
 {
   static const char * const primes[] = {
     "2", "53", "1234567891",
-    "2055693949", "16412292043871650369",
+    "2055693949", "1125899906842597", "16412292043871650369",
     /* diffie-hellman-group1-sha1, also "Well known group 2" in RFC
        2412, 2^1024 - 2^960 - 1 + 2^64 * { [2^894 pi] + 129093 } */
     "0xFFFFFFFFFFFFFFFFC90FDAA22168C234C4C6628B80DC1CD1"


More information about the gmp-commit mailing list