[Gmp-commit] /var/hg/gmp: New function mpn_gcd_22.
mercurial at gmplib.org
mercurial at gmplib.org
Fri Aug 16 06:06:04 UTC 2019
details: /var/hg/gmp/rev/dc279a8fb16c
changeset: 17815:dc279a8fb16c
user: Niels M?ller <nisse at lysator.liu.se>
date: Fri Aug 16 08:00:46 2019 +0200
description:
New function mpn_gcd_22.
* mpn/generic/gcd.c (gcd_2): Moved to gcd_22.c below.
(mpn_gcd): Adapt for calling gcd_22.
* mpn/generic/gcd_22.c (mpn_gcd_22): New file and function.
* gmp-impl.h (mp_double_limb_t): New (typedef) struct.
* configure.ac (gmp_mpn_functions): Added gcd_22.
* tests/mpn/t-gcd_22.c: New test.
* tests/mpn/Makefile.am (check_PROGRAMS): Add t-gcd_22.
* tests/refmpz.c (refmpz_gcd): New function (plain binary gcd).
diffstat:
ChangeLog | 12 +++++++
configure.ac | 2 +-
gmp-impl.h | 8 ++++
mpn/generic/gcd.c | 63 +++-----------------------------------
mpn/generic/gcd_22.c | 81 +++++++++++++++++++++++++++++++++++++++++++++++++
tests/mpn/Makefile.am | 2 +-
tests/mpn/t-gcd_22.c | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++
tests/refmpz.c | 46 ++++++++++++++++++++++++++++
tests/tests.h | 1 +
9 files changed, 239 insertions(+), 59 deletions(-)
diffs (truncated from 380 to 300 lines):
diff -r fad0937a0af8 -r dc279a8fb16c ChangeLog
--- a/ChangeLog Thu Aug 15 22:45:45 2019 +0200
+++ b/ChangeLog Fri Aug 16 08:00:46 2019 +0200
@@ -1,3 +1,15 @@
+2019-08-16 Niels Möller <nisse at lysator.liu.se>
+
+ * mpn/generic/gcd.c (gcd_2): Moved to gcd_22.c below.
+ (mpn_gcd): Adapt for calling gcd_22.
+ * mpn/generic/gcd_22.c (mpn_gcd_22): New file and function.
+ * gmp-impl.h (mp_double_limb_t): New (typedef) struct.
+ * configure.ac (gmp_mpn_functions): Added gcd_22.
+
+ * tests/mpn/t-gcd_22.c: New test.
+ * tests/mpn/Makefile.am (check_PROGRAMS): Add t-gcd_22.
+ * tests/refmpz.c (refmpz_gcd): New function (plain binary gcd).
+
2019-08-13 Torbjörn Granlund <tg at gmplib.org>
* mpn/x86_64: Add more gcd_11 variants of of x86_64 gcd_11.asm and
diff -r fad0937a0af8 -r dc279a8fb16c configure.ac
--- a/configure.ac Thu Aug 15 22:45:45 2019 +0200
+++ b/configure.ac Fri Aug 16 08:00:46 2019 +0200
@@ -2970,7 +2970,7 @@
rootrem sqrtrem sizeinbase get_str set_str compute_powtab \
scan0 scan1 popcount hamdist cmp zero_p \
perfsqr perfpow strongfibo \
- gcd_11 gcd_1 gcd gcdext_1 gcdext gcd_subdiv_step \
+ gcd_11 gcd_22 gcd_1 gcd gcdext_1 gcdext gcd_subdiv_step \
gcdext_lehmer \
div_q tdiv_qr jacbase jacobi_2 jacobi get_d \
matrix22_mul matrix22_mul1_inverse_vector \
diff -r fad0937a0af8 -r dc279a8fb16c gmp-impl.h
--- a/gmp-impl.h Thu Aug 15 22:45:45 2019 +0200
+++ b/gmp-impl.h Fri Aug 16 08:00:46 2019 +0200
@@ -3968,6 +3968,14 @@
#define PP_FIRST_OMITTED 3
#endif
+typedef struct
+{
+ mp_limb_t d0, d1;
+} mp_double_limb_t;
+
+#define mpn_gcd_22 __MPN (gcd_22)
+__GMP_DECLSPEC mp_double_limb_t mpn_gcd_22 (mp_limb_t, mp_limb_t, mp_limb_t, mp_limb_t);
+
/* BIT1 means a result value in bit 1 (second least significant bit), with a
zero bit representing +1 and a one bit representing -1. Bits other than
bit 1 are garbage. These are meant to be kept in "int"s, and casts are
diff -r fad0937a0af8 -r dc279a8fb16c mpn/generic/gcd.c
--- a/mpn/generic/gcd.c Thu Aug 15 22:45:45 2019 +0200
+++ b/mpn/generic/gcd.c Fri Aug 16 08:00:46 2019 +0200
@@ -76,61 +76,6 @@
ctx->gn = gn;
}
-#if GMP_NAIL_BITS > 0
-/* Nail supports should be easy, replacing the sub_ddmmss with nails
- * logic. */
-#error Nails not supported.
-#endif
-
-/* Use binary algorithm to compute G <-- GCD (U, V) for usize, vsize == 2.
- Both U and V must be odd. */
-static inline mp_size_t
-gcd_2 (mp_ptr gp, mp_srcptr up, mp_srcptr vp)
-{
- mp_limb_t u0, u1, v0, v1;
- mp_size_t gn;
-
- u0 = up[0];
- u1 = up[1];
- v0 = vp[0];
- v1 = vp[1];
-
- ASSERT (u0 & 1);
- ASSERT (v0 & 1);
-
- /* Check for u0 != v0 needed to ensure that argument to
- * count_trailing_zeros is non-zero. */
- while (u1 != v1 && u0 != v0)
- {
- unsigned long int r;
- if (u1 > v1)
- {
- sub_ddmmss (u1, u0, u1, u0, v1, v0);
- count_trailing_zeros (r, u0);
- u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
- u1 >>= r;
- }
- else /* u1 < v1. */
- {
- sub_ddmmss (v1, v0, v1, v0, u1, u0);
- count_trailing_zeros (r, v0);
- v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
- v1 >>= r;
- }
- }
-
- gp[0] = u0, gp[1] = u1, gn = 1 + (u1 != 0);
-
- /* If U == V == GCD, done. Otherwise, compute GCD (V, |U - V|). */
- if (u1 == v1 && u0 == v0)
- return gn;
-
- v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
- gp[0] = mpn_gcd_1 (gp, gn, v0);
-
- return 1;
-}
-
mp_size_t
mpn_gcd (mp_ptr gp, mp_ptr up, mp_size_t usize, mp_ptr vp, mp_size_t n)
{
@@ -301,8 +246,12 @@
vp[1] >>= r;
}
- ctx.gn = gcd_2(gp, up, vp);
-
+ {
+ mp_double_limb_t g = mpn_gcd_22 (up[1], up[0], vp[1], vp[0]);
+ gp[0] = g.d0;
+ gp[1] = g.d1;
+ ctx.gn = 1 + (g.d1 > 0);
+ }
done:
TMP_FREE;
return ctx.gn;
diff -r fad0937a0af8 -r dc279a8fb16c mpn/generic/gcd_22.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/gcd_22.c Fri Aug 16 08:00:46 2019 +0200
@@ -0,0 +1,81 @@
+/* mpn_gcd_22 -- double limb greatest common divisor.
+
+Copyright 1994, 1996, 2000, 2001, 2009, 2012, 2019 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 "gmp-impl.h"
+#include "longlong.h"
+
+#if GMP_NAIL_BITS > 0
+#error Nails not supported.
+#endif
+
+mp_double_limb_t
+mpn_gcd_22 (mp_limb_t u1, mp_limb_t u0, mp_limb_t v1, mp_limb_t v0)
+{
+ mp_double_limb_t g;
+ ASSERT (u0 & 1);
+ ASSERT (v0 & 1);
+
+ /* Check for u0 != v0 needed to ensure that argument to
+ * count_trailing_zeros is non-zero. */
+ while (u1 != v1 && u0 != v0)
+ {
+ unsigned long int r;
+ if (u1 > v1)
+ {
+ sub_ddmmss (u1, u0, u1, u0, v1, v0);
+ count_trailing_zeros (r, u0);
+ u0 = ((u1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (u0 >> r);
+ u1 >>= r;
+ }
+ else /* u1 < v1. */
+ {
+ sub_ddmmss (v1, v0, v1, v0, u1, u0);
+ count_trailing_zeros (r, v0);
+ v0 = ((v1 << (GMP_NUMB_BITS - r)) & GMP_NUMB_MASK) | (v0 >> r);
+ v1 >>= r;
+ }
+ }
+
+ /* If U == V == GCD, done. Otherwise, compute GCD (V, |U - V|). */
+ if (u1 == v1 && u0 == v0)
+ {
+ g.d0 = u0; g.d1 = u1;
+ return g;
+ }
+ else {
+ mp_limb_t t[2];
+ t[0] = u0; t[1] = u1;
+
+ v0 = (u0 == v0) ? ((u1 > v1) ? u1-v1 : v1-u1) : ((u0 > v0) ? u0-v0 : v0-u0);
+ g.d0 = mpn_gcd_1 (t, 1 + (u1 != 0), v0);
+ g.d1 = 0;
+ return g;
+ }
+}
diff -r fad0937a0af8 -r dc279a8fb16c tests/mpn/Makefile.am
--- a/tests/mpn/Makefile.am Thu Aug 15 22:45:45 2019 +0200
+++ b/tests/mpn/Makefile.am Fri Aug 16 08:00:46 2019 +0200
@@ -29,7 +29,7 @@
t-toom2-sqr t-toom3-sqr t-toom4-sqr t-toom6-sqr t-toom8-sqr \
t-div t-mul t-mullo t-sqrlo t-mulmod_bnm1 t-sqrmod_bnm1 t-mulmid \
t-hgcd t-hgcd_appr t-matrix22 t-invert t-bdiv t-fib2m \
- t-broot t-brootinv t-minvert t-sizeinbase t-gcd_11
+ t-broot t-brootinv t-minvert t-sizeinbase t-gcd_11 t-gcd_22
EXTRA_DIST = toom-shared.h toom-sqr-shared.h
diff -r fad0937a0af8 -r dc279a8fb16c tests/mpn/t-gcd_22.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/mpn/t-gcd_22.c Fri Aug 16 08:00:46 2019 +0200
@@ -0,0 +1,83 @@
+/* Test mpn_gcd_22.
+
+Copyright 2019 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library test suite.
+
+The GNU MP Library test suite is free software; you can redistribute it
+and/or modify it under the terms of the GNU General Public License as
+published by the Free Software Foundation; either version 3 of the License,
+or (at your option) any later version.
+
+The GNU MP Library test suite 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 a copy of the GNU General Public License along with
+the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "gmp-impl.h"
+#include "tests.h"
+
+#ifndef COUNT
+#define COUNT 10000
+#endif
+
+static void
+one_test (mpz_srcptr a, mpz_srcptr b, mpz_srcptr ref)
+{
+ mp_double_limb_t r = mpn_gcd_22 (mpz_getlimbn (a, 1), mpz_getlimbn (a, 0),
+ mpz_getlimbn (b, 1), mpz_getlimbn (b, 0));
+ if (r.d0 != mpz_getlimbn (ref, 0) || r.d1 != mpz_getlimbn (ref, 1))
+ {
+ gmp_fprintf (stderr,
+ "gcd_22 (0x%Zx, 0x%Zx) failed, got: g1 = 0x%Mx g0 = %Mx, ref: 0x%Zx\n",
+ a, b, r.d1, r.d0, ref);
+ abort();
+ }
+}
+
+int
+main (int argc, char **argv)
+{
+ mpz_t a, b, ref;
+ int count = COUNT;
+ int test;
+ gmp_randstate_ptr rands;
+
+ TESTS_REPS (count, argv, argc);
+
+ tests_start ();
+ rands = RANDS;
+
+ mpz_init (a);
+ mpz_init (b);
+ mpz_init (ref);
+ for (test = 0; test < count; test++)
+ {
+ mp_bitcnt_t size = 1 + gmp_urandomm_ui(rands, 2*GMP_NUMB_BITS);
+ if (test & 1)
+ {
+ mpz_urandomb (a, rands, size);
+ mpz_urandomb (b, rands, size);
+ }
+ else
+ {
+ mpz_rrandomb (a, rands, size);
More information about the gmp-commit
mailing list