[Gmp-commit] /var/hg/gmp: 5 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Thu Nov 15 17:49:55 UTC 2018
details: /var/hg/gmp/rev/0c123def7bf0
changeset: 17702:0c123def7bf0
user: "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date: Mon Nov 12 18:54:02 2018 +0100
description:
mpn/generic/strongfibo.c: New file.
details: /var/hg/gmp/rev/04b7365f7a7b
changeset: 17703:04b7365f7a7b
user: "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date: Mon Nov 12 18:56:37 2018 +0100
description:
mpz/stronglucas.c: New file
details: /var/hg/gmp/rev/d417bcef18e6
changeset: 17704:d417bcef18e6
user: "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date: Mon Nov 12 23:04:51 2018 +0100
description:
mpz/millerrabin.c: Implement BPSW test for primality.
details: /var/hg/gmp/rev/873b9552435f
changeset: 17705:873b9552435f
user: "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date: Mon Nov 12 19:24:14 2018 +0100
description:
ChangeLog
details: /var/hg/gmp/rev/204e6427d536
changeset: 17706:204e6427d536
user: "Marco Bodrato <bodrato at mail.dm.unipi.it>"
date: Mon Nov 12 23:13:39 2018 +0100
description:
AUTHORS and NEWS
diffstat:
AUTHORS | 6 +-
ChangeLog | 11 ++
Makefile.am | 3 +-
NEWS | 3 +
configure.ac | 2 +-
gmp-impl.h | 6 +
mpn/generic/strongfibo.c | 217 +++++++++++++++++++++++++++++++++++++++++++++++
mpz/Makefile.am | 3 +-
mpz/millerrabin.c | 53 +++++------
mpz/stronglucas.c | 178 ++++++++++++++++++++++++++++++++++++++
10 files changed, 448 insertions(+), 34 deletions(-)
diffs (truncated from 614 to 300 lines):
diff -r 3ead655ddc63 -r 204e6427d536 AUTHORS
--- a/AUTHORS Mon Nov 12 23:11:26 2018 +0100
+++ b/AUTHORS Mon Nov 12 23:13:39 2018 +0100
@@ -36,7 +36,7 @@
gcd.c, gcdext.c, matrix22_mul.c,
gcdext_1.c, gcd_subdiv_step.c, gcd_lehmer.c,
gcdext_subdiv_step.c, gcdext_lehmer.c,
- jacobi_2.c, jacbase.c, hgcd_jacobi.c, hgcd2_jacobi.c
+ jacobi_2.c, jacbase.c, hgcd_jacobi.c, hgcd2_jacobi.c,
matrix22_mul1_inverse_vector.c,
toom_interpolate_7pts, mulmod_bnm1.c, dcpi1_bdiv_qr.c,
dcpi1_bdiv_q.c, sbpi1_bdiv_qr.c, sbpi1_bdiv_q.c,
@@ -62,12 +62,14 @@
toom8h_mul.c, toom8_sqr.c, toom_interpolate_16pts.c,
mulmod_bnm1.c, sqrmod_bnm1.c, nussbaumer_mul.c,
toom_eval_pm2.c, toom_eval_pm2rexp.c,
+ fib2m.c, strongfibo.c,
mullo_n.c, sqrlo.c, invert.c, invertappr.c;
mpn/x86/atom/aors_n.asm, aorslshC_n.asm,
aorrlsh{1,2,C}_n.asm, aorsmul_1.asm, logops_n.asm,
sublsh2_n.asm, rshift.asm; primesieve.c;
mpz/fac_ui.c, 2fac_ui.c, mfac_uiui.c, oddfac_1.c,
- primorial_ui.c, prodlimbs.c, bin_ui.c,
+ primorial_ui.c, prodlimbs.c, bin_ui.c,
+ lucmod.c, stronglucas.c,
goetgheluck_bin_uiui.c; mini-gmp/mini-mpq.c.
David Harvey mpn/generic/add_err1_n.c, add_err2_n.c,
diff -r 3ead655ddc63 -r 204e6427d536 ChangeLog
--- a/ChangeLog Mon Nov 12 23:11:26 2018 +0100
+++ b/ChangeLog Mon Nov 12 23:13:39 2018 +0100
@@ -42,6 +42,17 @@
* tests/mpz/t-pprime_p.c (check_primes): Two more primes.
* tests/mp?: Use TESTS_REPS in many files.
+ * mpn/generic/strongfibo.c: New file, Fibonacci primality test.
+ * configure.ac (gmp_mpn_functions): Add it.
+ * gmp-impl.h: Declare mpn_strongfibo.
+
+ * mpz/stronglucas.c: New file, strong Lucas primality test.
+ * Makefile.am (MPZ_OBJECTS): Add it.
+ * mpz/Makefile.am (libmpz_la_SOURCES): Add it.
+ * gmp-impl.h: Declare mpz_stronglucas.
+
+ * mpz/millerrabin.c: Implement BPSW test for primality.
+
2018-11-07 Torbjörn Granlund <tg at gmplib.org>
* configure.ac (arm): Support a12 and a17.
diff -r 3ead655ddc63 -r 204e6427d536 Makefile.am
--- a/Makefile.am Mon Nov 12 23:11:26 2018 +0100
+++ b/Makefile.am Mon Nov 12 23:13:39 2018 +0100
@@ -210,7 +210,8 @@
mpz/set_q$U.lo mpz/set_si$U.lo mpz/set_str$U.lo mpz/set_ui$U.lo \
mpz/setbit$U.lo \
mpz/size$U.lo mpz/sizeinbase$U.lo mpz/sqrt$U.lo \
- mpz/sqrtrem$U.lo mpz/sub$U.lo mpz/sub_ui$U.lo mpz/swap$U.lo \
+ mpz/sqrtrem$U.lo mpz/stronglucas$U.lo mpz/sub$U.lo \
+ mpz/sub_ui$U.lo mpz/swap$U.lo \
mpz/tdiv_ui$U.lo mpz/tdiv_q$U.lo mpz/tdiv_q_2exp$U.lo \
mpz/tdiv_q_ui$U.lo mpz/tdiv_qr$U.lo mpz/tdiv_qr_ui$U.lo \
mpz/tdiv_r$U.lo mpz/tdiv_r_2exp$U.lo mpz/tdiv_r_ui$U.lo \
diff -r 3ead655ddc63 -r 204e6427d536 NEWS
--- a/NEWS Mon Nov 12 23:11:26 2018 +0100
+++ b/NEWS Mon Nov 12 23:13:39 2018 +0100
@@ -23,6 +23,9 @@
* Mini-GMP: added support for the mpq_t layer.
+ * Functions to detect primality now substitute 24 Miller-Rabin
+ iterations with the BPSW test.
+
SPEEDUPS
* The n-over-k function have been reimplemented for great
speedups for large operands handled by mpz_bin_ui .
diff -r 3ead655ddc63 -r 204e6427d536 configure.ac
--- a/configure.ac Mon Nov 12 23:11:26 2018 +0100
+++ b/configure.ac Mon Nov 12 23:13:39 2018 +0100
@@ -2953,7 +2953,7 @@
random random2 pow_1 \
rootrem sqrtrem sizeinbase get_str set_str compute_powtab \
scan0 scan1 popcount hamdist cmp zero_p \
- perfsqr perfpow \
+ perfsqr perfpow strongfibo \
gcd_1 gcd gcdext_1 gcdext gcd_subdiv_step \
gcdext_lehmer \
div_q tdiv_qr jacbase jacobi_2 jacobi get_d \
diff -r 3ead655ddc63 -r 204e6427d536 gmp-impl.h
--- a/gmp-impl.h Mon Nov 12 23:11:26 2018 +0100
+++ b/gmp-impl.h Mon Nov 12 23:13:39 2018 +0100
@@ -1136,6 +1136,9 @@
#define mpn_fib2m __MPN(fib2m)
__GMP_DECLSPEC int mpn_fib2m (mp_ptr, mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
+#define mpn_strongfibo __MPN(strongfibo)
+__GMP_DECLSPEC int mpn_strongfibo (mp_srcptr, mp_size_t, mp_ptr);
+
/* Remap names of internal mpn functions. */
#define __clz_tab __MPN(clz_tab)
#define mpn_udiv_w_sdiv __MPN(udiv_w_sdiv)
@@ -1704,6 +1707,9 @@
#define mpz_oddfac_1 __gmpz_oddfac_1
__GMP_DECLSPEC void mpz_oddfac_1 (mpz_ptr, mp_limb_t, unsigned);
+#define mpz_stronglucas __gmpz_stronglucas
+__GMP_DECLSPEC int mpz_stronglucas (mpz_srcptr, mpz_ptr, mpz_ptr);
+
#define mpz_lucas_mod __gmpz_lucas_mod
__GMP_DECLSPEC int mpz_lucas_mod (mpz_ptr, mpz_ptr, long, mp_bitcnt_t, mpz_srcptr, mpz_ptr, mpz_ptr);
diff -r 3ead655ddc63 -r 204e6427d536 mpn/generic/strongfibo.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/strongfibo.c Mon Nov 12 23:13:39 2018 +0100
@@ -0,0 +1,217 @@
+/* mpn_fib2m -- calculate Fibonacci numbers, modulo m.
+
+Contributed to the GNU project by Marco Bodrato.
+
+ 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.
+
+Copyright 2001, 2002, 2005, 2009, 2018 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of either:
+
+ * the GNU Lesser General Public License as published by the Free
+ Software Foundation; either version 3 of the License, or (at your
+ option) any later version.
+
+or
+
+ * the GNU General Public License as published by the Free Software
+ Foundation; either version 2 of the License, or (at your option) any
+ later version.
+
+or both in parallel, as here.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+for more details.
+
+You should have received copies of the GNU General Public License and the
+GNU Lesser General Public License along with the GNU MP Library. If not,
+see https://www.gnu.org/licenses/. */
+
+#include <stdio.h>
+#include "gmp-impl.h"
+
+/* Stores |{ap,n}-{bp,n}| in {rp,n},
+ returns the sign of {ap,n}-{bp,n}. */
+static int
+abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
+{
+ mp_limb_t x, y;
+ while (--n >= 0)
+ {
+ x = ap[n];
+ y = bp[n];
+ if (x != y)
+ {
+ ++n;
+ if (x > y)
+ {
+ mpn_sub_n (rp, ap, bp, n);
+ return 1;
+ }
+ else
+ {
+ mpn_sub_n (rp, bp, ap, n);
+ return -1;
+ }
+ }
+ rp[n] = 0;
+ }
+ return 0;
+}
+
+/* Computes at most count terms of the sequence needed by the
+ Lucas-Lehmer-Riesel test, indexing backward:
+ L_i = L_{i+1}^2 - 2
+
+ The sequence is computed modulo M = {mp, mn}.
+ The starting point is given in L_{count+1} = {lp, mn}.
+ The scratch pointed by sp, needs a space of at least 3 * mn + 1 limbs.
+
+ Returns the index i>0 if L_i = 0 (mod M) is found within the
+ computed count terms of the sequence. Otherwise it returns zero.
+
+ Note: (+/-2)^2-2=2, (+/-1)^2-2=-1, 0^2-2=-2
+ */
+
+static mp_bitcnt_t
+mpn_llriter (mp_ptr lp, mp_srcptr mp, mp_size_t mn, mp_bitcnt_t count, mp_ptr sp)
+{
+ do
+ {
+ mpn_sqr (sp, lp, mn);
+ mpn_tdiv_qr (sp + 2 * mn, lp, 0, sp, 2 * mn, mp, mn);
+ if (lp[0] < 5)
+ {
+ /* If L^2 % M < 5, |L^2 % M - 2| <= 2 */
+ if (mn == 1 || mpn_zero_p (lp + 1, mn - 1))
+ return (lp[0] == 2) ? count : 0;
+ else
+ MPN_DECR_U (lp, mn, 2);
+ }
+ else
+ lp[0] -= 2;
+ } while (--count != 0);
+ return 0;
+}
+
+/* Store the Lucas' number L[n] at lp (maybe), computed modulo m. lp
+ and scratch should have room for mn*2+1 limbs.
+
+ Returns the size of L[n] normally.
+
+ If F[n] is zero modulo m, or L[n] is, returns 0 and lp is
+ undefined.
+*/
+
+static mp_size_t
+mpn_lucm (mp_ptr lp, mp_srcptr np, mp_size_t nn, mp_srcptr mp, mp_size_t mn, mp_ptr scratch)
+{
+ int neg;
+ mp_limb_t cy;
+
+ ASSERT (! MPN_OVERLAP_P (lp, MAX(2*mn+1,5), scratch, MAX(2*mn+1,5)));
+ ASSERT (nn > 0);
+
+ neg = mpn_fib2m (lp, scratch, np, nn, mp, mn);
+
+ /* F[n] = +/-{lp, mn}, F[n-1] = +/-{scratch, mn} */
+ if (mpn_zero_p (lp, mn))
+ return 0;
+
+ if (neg) /* One sign is opposite, use sub instead of add. */
+ {
+#if HAVE_NATIVE_mpn_rsblsh1_n || HAVE_NATIVE_mpn_sublsh1_n
+#if HAVE_NATIVE_mpn_rsblsh1_n
+ cy = mpn_rsblsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */
+#else
+ cy = mpn_sublsh1_n (lp, lp, scratch, mn); /* L[n] = -/+(F[n]-(-2F[n-1])) */
+ if (cy != 0)
+ cy = mpn_add_n (lp, lp, mp, mn) - cy;
+#endif
+ if (cy > 1)
+ cy += mpn_add_n (lp, lp, mp, mn);
+#else
+ cy = mpn_lshift (scratch, scratch, mn, 1); /* 2F[n-1] */
+ if (cy)
+ cy -= mpn_sub_n (lp, scratch, lp, mn); /* L[n] = +/-(2F[n-1]-(-F[n])) */
+ else
+ abs_sub_n (lp, lp, scratch, mn);
+#endif
+ ASSERT (cy <= 1);
+ }
+ else
+ {
+#if HAVE_NATIVE_mpn_addlsh1_n
+ cy = mpn_addlsh1_n (lp, lp, scratch, mn); /* L[n] = +/-(2F[n-1]+F[n])) */
+#else
+ cy = mpn_lshift (scratch, scratch, mn, 1);
+ cy+= mpn_add_n (lp, lp, scratch, mn);
+#endif
+ ASSERT (cy <= 2);
+ }
+ while (cy || mpn_cmp (lp, mp, mn) >= 0)
+ cy -= mpn_sub_n (lp, lp, mp, mn);
+ MPN_NORMALIZE (lp, mn);
+ return mn;
+}
+
+int
+mpn_strongfibo (mp_srcptr mp, mp_size_t mn, mp_ptr scratch)
+{
+ mp_ptr lp, sp;
+ mp_size_t en;
+ mp_bitcnt_t b0;
+ int neg;
+ TMP_DECL;
+
+#if GMP_NUMB_BITS % 4 == 0
+ b0 = mpn_scan0 (mp, 0);
+#else
+ {
+ mpz_t m = MPZ_ROINIT_N(mp, mn);
+ b0 = mpz_scan0 (m, 0);
+ }
+ if (b0 == mn * GMP_NUMB_BITS)
+ {
+ en = 1;
+ scratch [0] = 1;
+ }
+ else
More information about the gmp-commit
mailing list