[Gmp-commit] /var/hg/gmp: 5 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Mon Jul 20 14:30:46 CEST 2026
details: /var/hg/gmp/rev/2b662fbd0183
changeset: 18562:2b662fbd0183
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Tue Jul 14 20:48:57 2026 +0200
description:
mpn/generic/sqrtrem.c (mpn_sqrtrem1): Swap roles of sqrt and rem.
details: /var/hg/gmp/rev/6bc54872f4b6
changeset: 18563:6bc54872f4b6
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jul 20 14:12:49 2026 +0200
description:
mpn/generic/bsqrt.c: Use a last iteration to increase precision
tests/mpn/t-bsqrt.c: Check exactly the required bits.
details: /var/hg/gmp/rev/aee4e58bc6f0
changeset: 18564:aee4e58bc6f0
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jul 20 14:15:53 2026 +0200
description:
tests/mpz/t-perfsqr.c (tests/mpz/t-perfsqr.c): New tests on large squares and almost-squares
details: /var/hg/gmp/rev/3a97c649ad6f
changeset: 18565:3a97c649ad6f
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jul 20 14:16:11 2026 +0200
description:
ChangeLog
details: /var/hg/gmp/rev/d2888762e87d
changeset: 18566:d2888762e87d
user: Marco Bodrato <bodrato at mail.dm.unipi.it>
date: Mon Jul 20 14:30:32 2026 +0200
description:
tests/mpz/t-perfsqr.c: Balance the number of repetitions, to avoid a too long test
diffstat:
ChangeLog | 5 ++
mpn/generic/bsqrt.c | 107 +++++++++++++++++++++++++++++++++++++++++++------
mpn/generic/sqrtrem.c | 22 +++++-----
tests/mpn/t-bsqrt.c | 36 +++++++++++++---
tests/mpz/t-perfsqr.c | 76 +++++++++++++++++++++++++++++++++++-
5 files changed, 212 insertions(+), 34 deletions(-)
diffs (truncated from 375 to 300 lines):
diff -r 29913bdd5d7c -r d2888762e87d ChangeLog
--- a/ChangeLog Sun Jul 12 08:19:13 2026 +0200
+++ b/ChangeLog Mon Jul 20 14:30:32 2026 +0200
@@ -1,3 +1,8 @@
+2026-07-20 Marco Bodrato <marco.bodrato at protonmail.com>
+
+ * mpn/generic/bsqrt.c: Use a last iteration to increase precision
+ * tests/mpn/t-bsqrt.c: Check exactly the required bits.
+
2026-07-12 Marco Bodrato <marco.bodrato at protonmail.com>
* tests/mpz/t-perfsqr.c: Check the format of return values.
diff -r 29913bdd5d7c -r d2888762e87d mpn/generic/bsqrt.c
--- a/mpn/generic/bsqrt.c Sun Jul 12 08:19:13 2026 +0200
+++ b/mpn/generic/bsqrt.c Mon Jul 20 14:30:32 2026 +0200
@@ -30,27 +30,106 @@
#include "gmp-impl.h"
-/* tp needs (1 + bnb / GMP_NUMB_BITS) limbs of space + the scratch
- used by mpn_bsqrtinv, i.e 3*(1 + bnb / GMP_NUMB_BITS) */
+/* tp needs (1 + nb / GMP_NUMB_BITS) limbs of space + the scratch
+ used by mpn_bsqrtinv, i.e 3*(1 + nb / GMP_NUMB_BITS)
+
+ For large enough sizes, it calls mpn_bsqrtinv for just half the
+ needed precision, the last half is computed here.
+
+ If T is the result of mpn_bsqrtinv, and A the input,
+ R = A * T;
+ result = R + T * (A-R^2)/2.
+ In the computation of A-R^2, the lowest half is zero.
+
+ WARNING:
+ it reads 1 + (nb - 1) / GMP_NUMB_BITS limbs, and
+ it writes 1 + (nb - 2) / GMP_NUMB_BITS limbs.
+*/
int
mpn_bsqrt (mp_ptr rp, mp_srcptr ap, mp_bitcnt_t nb, mp_ptr tp)
{
- mp_ptr sp;
- mp_size_t n;
+ ASSERT (nb > 1);
+
+ --nb;
+ if (nb <= GMP_NUMB_BITS)
+ {
+ if (! mpn_bsqrtinv (tp, ap, nb, NULL))
+ return 0;
+ *rp = *tp * *ap;
+ }
+ else
+ {
+ mp_size_t nn= 1 + nb / GMP_NUMB_BITS;
+ if (nn <= 4) /* THINK: should we tune this? */
+ {
+ mp_ptr sp = tp + nn;
- ASSERT (nb > 0);
+ MPN_FILL (tp, nn, CNST_LIMB(0));
+ if (! mpn_bsqrtinv (tp, ap, nb, sp))
+ return 0;
- n = 1 + nb / GMP_NUMB_BITS;
- sp = tp + n;
+ nn -= nb % GMP_NUMB_BITS == 0;
+ mpn_mullo_n (rp, tp, ap, nn);
+ }
+ else
+ {
+ mp_ptr sp;
+ mp_size_t n;
+ mp_bitcnt_t bnb = nb + 2 >> 1;
+
+ n = 1 + bnb / GMP_NUMB_BITS;
+ sp = tp + n;
+
+ MPN_FILL (tp, n, CNST_LIMB(0));
+ if (! mpn_bsqrtinv (tp, ap, bnb, sp))
+ return 0;
+
+ /* S = A * T; This is the low part of the result. */
+ 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 < 64) /* THINK: should we tune this? */
+ {
+ /* R^2; */
+ mpn_sqr (sp, rp, n);
+ --n;
+
+ /* (A-R^2)/2 */
+ ASSERT (mpn_cmp (sp, ap, n) == 0);
+ sp += n;
+#if HAVE_NATIVE_mpn_rsh1sub_n
+ mpn_rsh1sub_n (sp, ap + n, sp, nn - n);
+#else
+ mpn_sub_n (sp, ap + n, sp, nn - n);
+ ASSERT_NOCARRY (mpn_rshift (sp, sp, nn - n, 1));
+#endif
+ }
+ else
+ {
+ mp_size_t rn = mpn_sqrmod_bnm1_next_size (n + 1);
+ TMP_DECL;
- if (n == 1)
- *rp = *tp * *ap;
- else
- mpn_mullo_n (rp, tp, ap, n);
+ TMP_MARK;
+ mp_ptr scratch = TMP_ALLOC_LIMBS (mpn_sqrmod_bnm1_itch (rn, n));
+ /* S^2; */
+ mpn_sqrmod_bnm1 (sp + rn, rn, rp, n, scratch);
+ TMP_FREE;
+
+ --n;
+ mp_limb_t bw = mpn_sub (sp + rn, sp + rn, rn, ap, n);
+ MPN_DECR_U (sp + rn, rn, bw);
+ /* (A-R^2)/2 */
+ bw = mpn_sub_n (sp, ap + n, sp + n + rn, rn - n);
+ mpn_sub_nc (sp + rn - n, ap + rn, sp + rn, nn - rn, bw);
+ ASSERT_NOCARRY (mpn_rshift (sp, sp, nn - n, 1));
+ }
+
+ nn -= nb % GMP_NUMB_BITS == 0;
+ /* result = R + T * (A-R^2)/2. */
+ mp_limb_t saved_limb = rp [n];
+ mpn_mullo_n (rp + n, sp, tp, nn - n);
+ mpn_add_1 (rp + n, rp + n, nn - n, saved_limb);
+ }
+ }
return 1;
}
diff -r 29913bdd5d7c -r d2888762e87d mpn/generic/sqrtrem.c
--- a/mpn/generic/sqrtrem.c Sun Jul 12 08:19:13 2026 +0200
+++ b/mpn/generic/sqrtrem.c Mon Jul 20 14:30:32 2026 +0200
@@ -8,8 +8,8 @@
IN FACT, IT IS ALMOST GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A
FUTURE GMP RELEASE.
-Copyright 1999-2002, 2004, 2005, 2008, 2010, 2012, 2015, 2017 Free Software
-Foundation, Inc.
+Copyright 1999-2002, 2004, 2005, 2008, 2010, 2012, 2015, 2017, 2026
+Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -101,7 +101,7 @@
0x02,0x02,0x01,0x01,0x01,0x01,0x00,0x00 /* sqrt(1/1f8)..sqrt(1/1ff) */
};
-/* Compute s = floor(sqrt(a0)), and *rp = a0 - s^2. */
+/* Compute *sp = floor(sqrt(a0)), and returns a0 - (*sp)^2. */
#if GMP_NUMB_BITS > 32
#define MAGIC CNST_LIMB(0x10000000000) /* 0xffe7debbfc < MAGIC < 0x232b1850f410 */
@@ -110,7 +110,7 @@
#endif
static mp_limb_t
-mpn_sqrtrem1 (mp_ptr rp, mp_limb_t a0)
+mpn_sqrtrem1 (mp_ptr sp, mp_limb_t a0)
{
#if GMP_NUMB_BITS > 32
mp_limb_t a1;
@@ -160,8 +160,8 @@
x0++;
}
- *rp = a0 - x2;
- return x0;
+ *sp = x0;
+ return a0 - x2;
}
@@ -190,8 +190,8 @@
ASSERT (np[1] >= GMP_NUMB_HIGHBIT / 2);
np0 = np[0];
- sp0 = mpn_sqrtrem1 (rp, np[1]);
- rp0 = rp[0];
+ rp0 = mpn_sqrtrem1 (sp, np[1]);
+ sp0 = sp[0];
/* rp0 <= 2*sp0 < 2^(Prec + 1) */
rp0 = (rp0 << (Prec - 1)) + (np0 >> (Prec + 1));
q = rp0 / sp0;
@@ -455,14 +455,14 @@
if (nn == 1) {
if (c == 0)
{
- sp[0] = mpn_sqrtrem1 (&rl, high);
+ rl = mpn_sqrtrem1 (sp, high);
if (rp != NULL)
rp[0] = rl;
}
else
{
- cc = mpn_sqrtrem1 (&rl, high << (2*c)) >> c;
- sp[0] = cc;
+ rl = mpn_sqrtrem1 (sp, high << (2*c));
+ cc = sp[0] >>= c;
if (rp != NULL)
rp[0] = rl = high - cc*cc;
}
diff -r 29913bdd5d7c -r d2888762e87d tests/mpn/t-bsqrt.c
--- a/tests/mpn/t-bsqrt.c Sun Jul 12 08:19:13 2026 +0200
+++ b/tests/mpn/t-bsqrt.c Mon Jul 20 14:30:32 2026 +0200
@@ -30,7 +30,8 @@
{
gmp_randstate_ptr rands;
- mp_ptr ap, rp, pp, scratch;
+ mp_ptr ap, rp, pp, sp;
+ mp_limb_t before_rp, before_sp;
int count = COUNT;
unsigned i;
TMP_DECL;
@@ -43,14 +44,17 @@
rands = RANDS;
ap = TMP_ALLOC_LIMBS (MAX_LIMBS);
- rp = TMP_ALLOC_LIMBS (MAX_LIMBS);
+ rp = TMP_ALLOC_LIMBS (MAX_LIMBS + 2) + 1;
pp = TMP_ALLOC_LIMBS (MAX_LIMBS);
- scratch = TMP_ALLOC_LIMBS (3*MAX_LIMBS);
+ sp = TMP_ALLOC_LIMBS (3*MAX_LIMBS + 2) + 1;
+ before_rp = rp [-1] = gmp_urandomm_ui (rands, GMP_NUMB_MAX);
+ before_sp = sp [-1] = gmp_urandomm_ui (rands, GMP_NUMB_MAX);
for (i = 0; i < count; i++)
{
mp_size_t n;
mp_bitcnt_t bn;
+ mp_limb_t after_rp, after_sp;
int res;
n = 1 + gmp_urandomm_ui (rands, MAX_LIMBS);
@@ -63,19 +67,35 @@
if (i & 0xf)
ap[0] = (ap[0] | 7) ^ 6;
- bn = 1 + gmp_urandomm_ui (rands, GMP_NUMB_BITS - (n == 1));
+ bn = 1 + gmp_urandomm_ui (rands, GMP_NUMB_BITS - 2*(n == 1));
- res = mpn_bsqrt (rp, ap, n * GMP_NUMB_BITS - bn, scratch);
+ after_rp = rp [n - (bn >= GMP_NUMB_BITS - 1)] = gmp_urandomm_ui (rands, GMP_NUMB_MAX);
+ after_sp = sp [3 * n] = gmp_urandomm_ui (rands, GMP_NUMB_MAX);
+ res = mpn_bsqrt (rp, ap, n * GMP_NUMB_BITS - bn, sp);
+ if (rp [n - (bn >= GMP_NUMB_BITS - 1)] != after_rp || rp [-1] != before_rp ||
+ sp [3 * n] != after_sp || sp [-1] != before_sp)
+ {
+ gmp_fprintf (stderr,
+ "mpn_bsqrt memoty bounds violated: %u limbs, - %u bits, res %i [%i]\n",
+ (unsigned) n, (unsigned) bn, res, i);
+ gmp_fprintf (stderr, "before_rp: %Mx <> %Mx, \t", before_rp, rp[-1]);
+ gmp_fprintf (stderr, " after_rp: %Mx <> %Mx\n", after_rp, rp[n - (bn >= GMP_NUMB_BITS - 1)]);
+ gmp_fprintf (stderr, "before_sp: %Mx <> %Mx, \t", before_sp, sp[-1]);
+ gmp_fprintf (stderr, " after_sp: %Mx <> %Mx\n", after_sp, sp[3 * n]);
+ gmp_fprintf (stderr, "a = %Nx\n", ap, n);
+ abort ();
+ }
+
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)))))
+ if (!res || ((n!=1) && (mpn_cmp (pp, ap, n - 1) != 0)) ||
+ ((bn != GMP_NUMB_BITS) && ((pp[n - 1] ^ ap[n - 1]) & GMP_NUMB_MAX >> bn != 0)))
{
gmp_fprintf (stderr,
- "mpn_bsqrt returned bad result: %u limbs, %u bits, res %i [%i]\n",
+ "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);
diff -r 29913bdd5d7c -r d2888762e87d tests/mpz/t-perfsqr.c
--- a/tests/mpz/t-perfsqr.c Sun Jul 12 08:19:13 2026 +0200
+++ b/tests/mpz/t-perfsqr.c Mon Jul 20 14:30:32 2026 +0200
@@ -231,6 +231,79 @@
mpz_clear (rop);
}
+/* Exercise mpz_perfect_square_p on large squares. */
+void
+check_sqares (int reps)
+{
+ mpz_t x2, x, rop;
+ mp_size_t xn;
+ int res;
+ int i;
+ int want;
+ mp_bitcnt_t l, h;
+ gmp_randstate_ptr rands = RANDS;
More information about the gmp-commit
mailing list