[Gmp-commit] /var/hg/gmp: 6 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Mon Jun 22 17:22:30 CEST 2026
details: /var/hg/gmp/rev/68fd3295c9ae
changeset: 18523:68fd3295c9ae
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 22 15:18:50 2026 +0200
description:
mpn/generic/bsqrtinv.c: Return the correct result also for bnb==1.
details: /var/hg/gmp/rev/c2b80cc2c5ef
changeset: 18524:c2b80cc2c5ef
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 22 15:26:06 2026 +0200
description:
mpn/generic/bsqrt.c: Update and cleanup.
gmp-impl.h (mpn_bsqrt): Change the prototype, now it returns int.
details: /var/hg/gmp/rev/24fda3c729e6
changeset: 18525:24fda3c729e6
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 22 15:28:06 2026 +0200
description:
tests/mpn/t-bsqrt.c: New File, test the mpn_bsqrt function.
tests/mpn/Makefile.am (check_PROGRAMS): Add t-bsqrt.
details: /var/hg/gmp/rev/4483ffa833aa
changeset: 18526:4483ffa833aa
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 22 15:45:30 2026 +0200
description:
tune/speed.[ch], tune/common.c: Support measuring mpn_bsqrt.
details: /var/hg/gmp/rev/0a1e009b26a7
changeset: 18527:0a1e009b26a7
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 22 17:17:17 2026 +0200
description:
mpz/remove.c: Always use the result memory space if it is large enough.
tests/mpz/t-remove.c: Test that if the result was stored in available space.
details: /var/hg/gmp/rev/2247098521da
changeset: 18528:2247098521da
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jun 22 17:21:44 2026 +0200
description:
ChangeLog
diffstat:
ChangeLog | 10 +++-
gmp-impl.h | 2 +-
mpn/generic/bsqrt.c | 23 ++++++++---
mpn/generic/bsqrtinv.c | 19 ++++++---
mpz/remove.c | 5 +--
tests/mpn/Makefile.am | 2 +-
tests/mpn/t-bsqrt.c | 89 ++++++++++++++++++++++++++++++++++++++++++++++++++
tests/mpz/t-remove.c | 5 +-
tune/common.c | 7 +++-
tune/speed.c | 1 +
tune/speed.h | 3 +-
11 files changed, 139 insertions(+), 27 deletions(-)
diffs (truncated from 349 to 300 lines):
diff -r afed04e71f87 -r 2247098521da ChangeLog
--- a/ChangeLog Sun Jun 21 16:51:47 2026 +0200
+++ b/ChangeLog Mon Jun 22 17:21:44 2026 +0200
@@ -3,11 +3,15 @@
* mpn/generic/perfpow.c (is_kth_power): Check only one possible
square root.
- * tune/speed.c: Recognise mpn_bsqrtinv.
- * tune/speed.h: Declare speed_mpn_bsqrtinv and the macro.
- * tune/common.c (speed_mpn_bsqrtinv): Implement the function.
+ * tune/speed.c: Recognise mpn_bsqrtinv, mpn_bsqrt.
+ * tune/speed.h: Declare speed_mpn_bsqrt{inv,} and the macro.
+ * tune/common.c (speed_mpn_bsqrt{inv,}): Implement the functions.
* mpn/generic/bsqrtinv.c: Rewrite with a differen iteration.
+ * mpn/generic/bsqrt.c: Return 0 on failure.
+
+ * tests/mpn/t-bsqrt.c: New File, test the mpn_bsqrt function.
+ * tests/mpn/Makefile.am (check_PROGRAMS): Add t-bsqrt.
2026-04-04 Marc Glisse <marc.glisse at inria.fr>
diff -r afed04e71f87 -r 2247098521da gmp-impl.h
--- a/gmp-impl.h Sun Jun 21 16:51:47 2026 +0200
+++ b/gmp-impl.h Mon Jun 22 17:21:44 2026 +0200
@@ -1800,7 +1800,7 @@
__GMP_DECLSPEC void mpn_brootinv (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t, mp_ptr);
#define mpn_bsqrt __MPN(bsqrt)
-__GMP_DECLSPEC void mpn_bsqrt (mp_ptr, mp_srcptr, mp_bitcnt_t, mp_ptr);
+__GMP_DECLSPEC int mpn_bsqrt (mp_ptr, mp_srcptr, mp_bitcnt_t, mp_ptr);
#define mpn_bsqrtinv __MPN(bsqrtinv)
__GMP_DECLSPEC int mpn_bsqrtinv (mp_ptr, mp_srcptr, mp_bitcnt_t, mp_ptr);
diff -r afed04e71f87 -r 2247098521da mpn/generic/bsqrt.c
--- a/mpn/generic/bsqrt.c Sun Jun 21 16:51:47 2026 +0200
+++ b/mpn/generic/bsqrt.c Mon Jun 22 17:21:44 2026 +0200
@@ -1,6 +1,6 @@
-/* mpn_bsqrt, a^{1/2} (mod 2^n).
+/* mpn_bsqrt, a^{1/2} (mod 2^n), for odd a.
-Copyright 2009, 2010, 2012, 2015 Free Software Foundation, Inc.
+Copyright 2009, 2010, 2012, 2015, 2026 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -30,8 +30,9 @@
#include "gmp-impl.h"
-
-void
+/* tp needs (1 + bnb / GMP_NUMB_BITS) limbs of space + the scratch
+ used by mpn_bsqrtinv, i.e 3*(1 + bnb / GMP_NUMB_BITS) */
+int
mpn_bsqrt (mp_ptr rp, mp_srcptr ap, mp_bitcnt_t nb, mp_ptr tp)
{
mp_ptr sp;
@@ -39,9 +40,17 @@
ASSERT (nb > 0);
- n = nb / GMP_NUMB_BITS;
+ n = 1 + nb / GMP_NUMB_BITS;
sp = tp + n;
- mpn_bsqrtinv (tp, ap, nb, sp);
- mpn_mullo_n (rp, tp, ap, n);
+ MPN_FILL (tp, n, CNST_LIMB(0));
+ if (! mpn_bsqrtinv (tp, ap, nb, sp))
+ return 0;
+
+ if (n == 1)
+ *rp = *tp * *ap;
+ else
+ mpn_mullo_n (rp, tp, ap, n);
+
+ return 1;
}
diff -r afed04e71f87 -r 2247098521da mpn/generic/bsqrtinv.c
--- a/mpn/generic/bsqrtinv.c Sun Jun 21 16:51:47 2026 +0200
+++ b/mpn/generic/bsqrtinv.c Mon Jun 22 17:21:44 2026 +0200
@@ -75,31 +75,33 @@
32, 74, 51, 249, 216, 162, 68, 30, 239, 122, 60, 182, 200, 45, 75, 206};
#endif
+/* tp needs 2*(1 + bnb / GMP_NUMB_BITS) limbs of space */
int
mpn_bsqrtinv (mp_ptr rp, mp_srcptr yp, mp_bitcnt_t bnb, mp_ptr tp)
{
+ mp_limb_t y0 = *yp;
ASSERT (bnb > 0);
#ifndef BSQRTINV_RP_NOT_ZEROED
ASSERT (mpn_zero_p (rp, 1 + bnb / GMP_NUMB_BITS));
#endif
- if (bnb == 1)
+ if (UNLIKELY (bnb == 1))
{
- if ((yp[0] & 3) != 1)
- return 0;
+ *rp = y0;
+ return (y0 & 3) == 1;
}
else
{
mp_ptr tp2 = tp + 1 + bnb / GMP_NUMB_BITS;
mp_size_t bn, order[GMP_LIMB_BITS + 1];
- mp_limb_t t0, r0, y0 = *yp;
+ mp_limb_t t0, r0;
int i;
if ((y0 & 7) != 1)
return 0;
#ifdef BSQRTINV_DONT_USE_TABLE
- r0 = 33 + ((y0 & 8) * 5 >> 2) - ((y0 & 16) >> 1);
+ r0 = y0 + ((y0 & 8) >> 2) + ((y0 & 16) >> 1);
t0 = r0 * r0 * y0 >> 1;
r0 -= r0 * t0;
@@ -162,6 +164,9 @@
ASSERT ((yrrm1d4 & (GMP_NUMB_MAX >> (GMP_NUMB_BITS - precomputed_bits + 1))) == 0);
mp_limb_t rt = r0 * yrrm1d4 - r0h - 1; /* (r*(r0*r0*y0-1)/2 - r) >>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;
r0 -= r0 * t0;
@@ -197,10 +202,10 @@
t4 = t3 * r0 - 1;
/* [2*t4+1,-r0] <- r0*(r0^2 y-1)/2 - r0 */
if (t4 & GMP_LIMB_HIGHBIT) {
- *rp = r0 & GMP_NUMB_MAX;
+ rp[0] = r0 & GMP_NUMB_MAX;
rp[1] = ~t4 << 1;
} else {
- *rp = -r0 & GMP_NUMB_MAX;
+ rp[0] = -r0 & GMP_NUMB_MAX;
rp[1] = (t4 << 1) ^ 1;
}
#ifdef BSQRTINV_RP_NOT_ZEROED
diff -r afed04e71f87 -r 2247098521da mpz/remove.c
--- a/mpz/remove.c Sun Jun 21 16:51:47 2026 +0200
+++ b/mpz/remove.c Mon Jun 22 17:21:44 2026 +0200
@@ -87,10 +87,9 @@
mpz_t fpow[GMP_LIMB_BITS]; /* Really MP_SIZE_T_BITS */
int p;
-#if WANT_ORIGINAL_DEST
mp_ptr dp;
dp = PTR (dest);
-#endif
+
/* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0). It is an
upper bound of the result we're seeking. We could also shift down the
operands so that they become odd, to make intermediate values
@@ -128,12 +127,10 @@
mpz_clear (fpow[p]);
}
-#if WANT_ORIGINAL_DEST
if (PTR (x) == dp) {
mpz_swap (dest, x);
mpz_set (dest, x);
}
-#endif
}
else
mpz_set (dest, src);
diff -r afed04e71f87 -r 2247098521da tests/mpn/Makefile.am
--- a/tests/mpn/Makefile.am Sun Jun 21 16:51:47 2026 +0200
+++ b/tests/mpn/Makefile.am Mon Jun 22 17:21:44 2026 +0200
@@ -30,7 +30,7 @@
t-div t-mul t-mullo t-sqrlo t-mulmod_bnm1 t-sqrmod_bnm1 t-mulmid \
t-mulmod_bknp1 t-sqrmod_bknp1 \
t-addaddmul 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-gcd_22 t-gcdext_1
+ t-broot t-bsqrt t-brootinv t-minvert t-sizeinbase t-gcd_11 t-gcd_22 t-gcdext_1
EXTRA_DIST = toom-shared.h toom-sqr-shared.h
diff -r afed04e71f87 -r 2247098521da tests/mpn/t-bsqrt.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/mpn/t-bsqrt.c Mon Jun 22 17:21:44 2026 +0200
@@ -0,0 +1,89 @@
+/* Copyright 2012, 2015, 2026 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 <stdlib.h> /* for abort */
+#include <stdio.h> /* for printf */
+
+#include "gmp-impl.h"
+#include "tests/tests.h"
+
+#define MAX_LIMBS 150
+#define COUNT 500
+
+int
+main (int argc, char **argv)
+{
+ gmp_randstate_ptr rands;
+
+ mp_ptr ap, rp, pp, scratch;
+ int count = COUNT;
+ unsigned i;
+ TMP_DECL;
+
+ TMP_MARK;
+
+ TESTS_REPS (count, argv, argc);
+
+ tests_start ();
+ rands = RANDS;
+
+ ap = TMP_ALLOC_LIMBS (MAX_LIMBS);
+ rp = TMP_ALLOC_LIMBS (MAX_LIMBS);
+ pp = TMP_ALLOC_LIMBS (MAX_LIMBS);
+ scratch = TMP_ALLOC_LIMBS (3*MAX_LIMBS);
+
+ for (i = 0; i < count; i++)
+ {
+ mp_size_t n;
+ mp_bitcnt_t bn;
+ int res;
+
+ n = 1 + gmp_urandomm_ui (rands, MAX_LIMBS);
+
+ if (i & 1)
+ mpn_random2 (ap, n);
+ else
+ mpn_random (ap, n);
+
+ if (i & 0xf)
+ ap[0] = (ap[0] | 7) ^ 6;
+
+ bn = 1 + gmp_urandomm_ui (rands, GMP_NUMB_BITS - (n == 1));
+
+ res = mpn_bsqrt (rp, ap, n * GMP_NUMB_BITS - bn, scratch);
+ if (!res && ((*ap & (7 >> ((n == 1) && (bn == GMP_NUMB_BITS - 1)))) != 1))
+ continue;
+
+ mpn_sqrlo (pp, rp, n);
+
+ if ((!res == (!mpn_cmp (pp, ap, n - (bn > 1)) &&
+ !((pp[n - 1] ^ ap[n - 1]) & GMP_NUMB_MAX >> (bn - 1)))))
+ {
+ gmp_fprintf (stderr,
+ "mpn_bsqrt returned bad result: %u limbs, %u bits, res %i [%i]\n",
+ (unsigned) n, (unsigned) bn, res, i);
+ gmp_fprintf (stderr, "a = %Nx\n", ap, n);
+ gmp_fprintf (stderr, "r = %Nx\n", rp, n);
+ gmp_fprintf (stderr, "r^2 = %Nx\n", pp, n);
+ abort ();
+ }
+ }
+ TMP_FREE;
+ tests_end ();
+ return 0;
+}
diff -r afed04e71f87 -r 2247098521da tests/mpz/t-remove.c
--- a/tests/mpz/t-remove.c Sun Jun 21 16:51:47 2026 +0200
+++ b/tests/mpz/t-remove.c Mon Jun 22 17:21:44 2026 +0200
@@ -32,7 +32,7 @@
{
unsigned long int exp;
mpz_t t, dest, refdest, dividend, divisor;
- mp_size_t dividend_size, divisor_size;
+ mp_size_t dividend_size, divisor_size, prev_alloc;
int i;
int reps = 1000;
unsigned long int pwr, refpwr;
@@ -73,9 +73,10 @@
mpz_mul (dividend, dividend, t);
refpwr = mpz_refremove (refdest, dividend, divisor);
+ prev_alloc = ALLOC(dest);
pwr = mpz_remove (dest, dividend, divisor);
- if (refpwr != pwr || mpz_cmp (refdest, dest) != 0)
+ if (refpwr != pwr || mpz_cmp (refdest, dest) != 0 || ALLOC(dest) < prev_alloc)
{
fprintf (stderr, "ERROR after %d tests\n", i);
More information about the gmp-commit
mailing list