mpz_sqrt_if_perfect_square
Seth Troisi
braintwo at gmail.com
Tue Jun 2 11:54:20 CEST 2026
I'm doing some prime gap searching which would slightly benefit from
knowing the root of a number iff it's a perfect square.
There's a item in https://gmplib.org/tasks related to this
"mpz_sqrt_if_perfect_square: When mpz_perfect_square_p does its tests it
calculates a square root and then discards it. For some applications it
might be useful to return that root. Suggested by Jason Moxham."
In 2012 there was also a discussion of doing this for mpz_perfect_power:
https://gmplib.org/list-archives/gmp-devel/2012-October/002562.html
I coded up a working implementation for myself. I can work on making it
clean enough for upstream if it's a desired feature.
Thanks,
Seth
-------------- next part --------------
# HG changeset patch
# User Seth Troisi <sethtroisi at google.com>
# Date 1780384287 25200
# Node ID e9607e3ec477fab3732cf00850352cfaa8498547
# Parent 1e504d1875715329ac699f1fb4e0f9aca95de4dc
Added mpz_perfect_square and mpn_perfect_square
diff -r 1e504d187571 -r e9607e3ec477 ChangeLog
--- a/ChangeLog Sat Apr 04 14:56:15 2026 +0200
+++ b/ChangeLog Tue Jun 02 00:11:27 2026 -0700
@@ -1,3 +1,8 @@
+2026-06-01 Seth Troisi <sethtroisi at google.com>
+ * mpz/perfpow.c (mpz_perfect_square): New function with root return.
+ * mpn/perfpow.c (mpn_perfect_square): New function with root return.
+ * tests/mpz/t-perfsqr.c: Add tests for mpz_perfect_square
+
2026-04-04 Marc Glisse <marc.glisse at inria.fr>
* gmpxx.h (operator""): Remove deprecated space.
diff -r 1e504d187571 -r e9607e3ec477 doc/gmp.texi
--- a/doc/gmp.texi Sat Apr 04 14:56:15 2026 +0200
+++ b/doc/gmp.texi Tue Jun 02 00:11:27 2026 -0700
@@ -3602,6 +3602,16 @@
be perfect squares.
@end deftypefun
+ at deftypefun int mpz_perfect_square (mpz_t @var{rop}, const mpz_t @var{op})
+ at cindex Perfect square functions
+ at cindex Root testing functions
+Return non-zero and sets @var{rop} to the square root if @var{op} is a perfect
+square, i.e., if the square root of @var{op} is an integer. Under this
+definition both 0 and 1 are considered to be perfect squares and @var{rop} is
+set to @var{op} in both cases. If @var{op} is not a perfect square the value of
+ at var{rop} is undefined.
+ at end deftypefun
+
@need 2000
@node Number Theoretic Functions, Integer Comparisons, Integer Roots, Integer Functions
@@ -5636,7 +5646,7 @@
would have been zero or non-zero.
A return value of zero indicates a perfect square. See also
- at code{mpn_perfect_square_p}.
+ at code{mpn_perfect_square_p} and @code{mpn_perfect_square}.
@end deftypefun
@deftypefun size_t mpn_sizeinbase (const mp_limb_t *@var{xp}, mp_size_t @var{n}, int @var{base})
@@ -5719,8 +5729,15 @@
@deftypefun int mpn_perfect_square_p (const mp_limb_t *@var{s1p}, mp_size_t @var{n})
Return non-zero iff @{@var{s1p}, @var{n}@} is a perfect square.
-The most significant limb of the input @{@var{s1p}, @var{n}@} must be
-non-zero.
+The most significant limb of the input @{@var{s1p}, @var{n}@} must be non-zero.
+ at end deftypefun
+
+ at deftypefun int mpn_perfect_square (mp_limb_t *@var{r1p}, const mp_limb_t *@var{s1p}, mp_size_t @var{n})
+Return non-zero iff @{@var{s1p}, @var{n}@} is a perfect square.
+iff @{@var{s1p}, @var{n}@} is a perfect square @{@var{r1p}, @math{@GMPceil{@var{n} / 2}}@}
+will be set to the square root.
+The most significant limb of @{@var{s1p}, @var{n}@} must be non-zero.
+The areas at @var{r1p} and @var{s1p} must be completely separate.
@end deftypefun
@deftypefun void mpn_and_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
@@ -9250,9 +9267,10 @@
A significant fraction of non-squares can be quickly identified by checking
whether the input is a quadratic residue modulo small integers.
- at code{mpz_perfect_square_p} first tests the input mod 256, which means just
-examining the low byte. Only 44 different values occur for squares mod 256,
-so 82.8% of inputs can be immediately identified as non-squares.
+ at code{mpz_perfect_square}/@code{mpz_perfect_square_p} first tests the input
+mod 256, which means just examining the low byte. Only 44 different values
+occur for squares mod 256, so 82.8% of inputs can be immediately identified
+as non-squares.
On a 32-bit system similar tests are done mod 9, 5, 7, 13 and 17, for a total
99.25% of inputs identified as non-squares. On a 64-bit system 97 is tested
diff -r 1e504d187571 -r e9607e3ec477 doc/tasks.html
--- a/doc/tasks.html Sat Apr 04 14:56:15 2026 +0200
+++ b/doc/tasks.html Tue Jun 02 00:11:27 2026 -0700
@@ -579,10 +579,6 @@
probably simplest to regard this as a compiler compatibility issue, and
leave it to users or sysadmins to ensure application and library code is
built the same.
-<li> <code>mpz_sqrt_if_perfect_square</code>: When
- <code>mpz_perfect_square_p</code> does its tests it calculates a square
- root and then discards it. For some applications it might be useful to
- return that root. Suggested by Jason Moxham.
<li> <code>mpz_get_ull</code>, <code>mpz_set_ull</code>,
<code>mpz_get_sll</code>, <code>mpz_get_sll</code>: Conversions for
<code>long long</code>. These would aid interoperability, though a
diff -r 1e504d187571 -r e9607e3ec477 gmp-h.in
--- a/gmp-h.in Sat Apr 04 14:56:15 2026 +0200
+++ b/gmp-h.in Tue Jun 02 00:11:27 2026 -0700
@@ -972,6 +972,9 @@
__GMP_DECLSPEC int mpz_perfect_square_p (mpz_srcptr) __GMP_ATTRIBUTE_PURE;
#endif
+#define mpz_perfect_square __gmpz_perfect_square
+__GMP_DECLSPEC int mpz_perfect_square (mpz_ptr, mpz_srcptr);
+
#define mpz_popcount __gmpz_popcount
#if __GMP_INLINE_PROTOTYPES || defined (__GMP_FORCE_mpz_popcount)
__GMP_DECLSPEC mp_bitcnt_t mpz_popcount (mpz_srcptr) __GMP_NOTHROW __GMP_ATTRIBUTE_PURE;
@@ -1570,6 +1573,9 @@
#define mpn_perfect_square_p __MPN(perfect_square_p)
__GMP_DECLSPEC int mpn_perfect_square_p (mp_srcptr, mp_size_t) __GMP_ATTRIBUTE_PURE;
+#define mpn_perfect_square __MPN(perfect_square)
+__GMP_DECLSPEC int mpn_perfect_square (mp_ptr, mp_srcptr, mp_size_t);
+
#define mpn_perfect_power_p __MPN(perfect_power_p)
__GMP_DECLSPEC int mpn_perfect_power_p (mp_srcptr, mp_size_t) __GMP_ATTRIBUTE_PURE;
diff -r 1e504d187571 -r e9607e3ec477 mpn/generic/perfsqr.c
--- a/mpn/generic/perfsqr.c Sat Apr 04 14:56:15 2026 +0200
+++ b/mpn/generic/perfsqr.c Tue Jun 02 00:11:27 2026 -0700
@@ -1,8 +1,11 @@
/* mpn_perfect_square_p(u,usize) -- Return non-zero if U is a perfect square,
zero otherwise.
-Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2012 Free Software
-Foundation, Inc.
+mpn_perfect_square(r, u, usize) -- Return non-zero and set r iff U is a
+perfect square, leave r unchanged and return zero otherwise.
+
+Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2012, 2026 Free
+Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -174,13 +177,11 @@
} \
} while (0)
-
-int
-mpn_perfect_square_p (mp_srcptr up, mp_size_t usize)
+int mpn_perfect_square (mp_ptr rp, mp_srcptr up, mp_size_t usize)
{
ASSERT (usize >= 1);
- TRACE (gmp_printf ("mpn_perfect_square_p %Nd\n", up, usize));
+ TRACE (gmp_printf ("mpn_perfect_square %Nd\n", up, usize));
/* The first test excludes 212/256 (82.8%) of the perfect square candidates
in O(1) time. */
@@ -221,18 +222,31 @@
/* For the third and last test, we finally compute the square root,
to make sure we've really got a perfect square. */
- {
- mp_ptr root_ptr;
- int res;
- TMP_DECL;
-
- TMP_MARK;
- root_ptr = TMP_ALLOC_LIMBS ((usize + 1) / 2);
+ if (rp == NULL)
+ {
+ mp_ptr root_ptr;
+ int res;
+ TMP_DECL;
+
+ TMP_MARK;
+ root_ptr = TMP_ALLOC_LIMBS ((usize + 1) / 2);
+
+ /* Iff mpn_sqrtrem returns zero, the square is perfect. */
+ res = ! mpn_sqrtrem (root_ptr, NULL, up, usize);
+ TMP_FREE;
+
+ return res;
+ }
+ else
+ {
+ /* Iff mpn_sqrtrem returns zero, the square is perfect. */
+ int res = ! mpn_sqrtrem (rp, NULL, up, usize);
+ return res;
+ }
+}
- /* Iff mpn_sqrtrem returns zero, the square is perfect. */
- res = ! mpn_sqrtrem (root_ptr, NULL, up, usize);
- TMP_FREE;
-
- return res;
- }
+int
+mpn_perfect_square_p (mp_srcptr up, mp_size_t usize)
+{
+ return mpn_perfect_square(NULL, up, usize);
}
diff -r 1e504d187571 -r e9607e3ec477 mpz/perfsqr.c
--- a/mpz/perfsqr.c Sat Apr 04 14:56:15 2026 +0200
+++ b/mpz/perfsqr.c Tue Jun 02 00:11:27 2026 -0700
@@ -1,7 +1,12 @@
/* mpz_perfect_square_p(arg) -- Return non-zero if ARG is a perfect square,
zero otherwise.
-Copyright 1991, 1993, 1994, 1996, 2000, 2001 Free Software Foundation, Inc.
+mpz_perfect_square(rop, arg) -- Return non-zero if ARG is a perfect square
+ sets rop to root, otherwise returns zero and leaves rop unchanged.
+
+# TODO review does this file need to be split in to perfsqr.c and perfsqrp.c?
+
+Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2026 Free Software Foundation, Inc.
This file is part of the GNU MP Library.
@@ -32,3 +37,36 @@
#define __GMP_FORCE_mpz_perfect_square_p 1
#include "gmp-impl.h"
+
+int
+mpz_perfect_square (mpz_ptr rop, mpz_srcptr u)
+{
+ int res;
+ mp_size_t u_size = SIZ (u);
+
+ /* negatives are not ever squares */
+ if (__GMP_UNLIKELY ( u_size < 0))
+ {
+ return 0;
+ }
+
+ /* zero is a square with rop = 0 */
+ if (__GMP_UNLIKELY ( u_size == 0))
+ {
+ SIZ(rop) = 0;
+ return 1;
+ }
+
+ int rop_size = (SIZ (u) + 1) / 2;
+ if (SIZ(rop) < rop_size)
+ {
+ MPZ_REALLOC(rop, rop_size);
+ }
+
+ res = mpn_perfect_square (PTR(rop), PTR (u), u_size);
+ if (res)
+ {
+ SIZ(rop) = rop_size;
+ }
+ return res;
+}
diff -r 1e504d187571 -r e9607e3ec477 tests/mpz/t-perfsqr.c
--- a/tests/mpz/t-perfsqr.c Sat Apr 04 14:56:15 2026 +0200
+++ b/tests/mpz/t-perfsqr.c Tue Jun 02 00:11:27 2026 -0700
@@ -1,4 +1,4 @@
-/* Test mpz_perfect_square_p.
+/* Test mpz_perfect_square and mpz_perfect_square_p.
Copyright 2000-2002 Free Software Foundation, Inc.
@@ -26,8 +26,8 @@
#include "mpn/perfsqr.h"
-/* check_modulo() exercises mpz_perfect_square_p on squares which cover each
- possible quadratic residue to each divisor used within
+/* check_modulo() exercises mpz_perfect_square/mpz_perfect_square_p on squares
+ which cover each possible quadratic residue to each divisor used within
mpn_perfect_square_p, ensuring those residues aren't incorrectly claimed
to be non-residues.
@@ -45,11 +45,13 @@
static const unsigned long divisor[] = PERFSQR_DIVISORS;
unsigned long i, j;
- mpz_t alldiv, others, n;
+ mpz_t alldiv, others, n, root, rop;
mpz_init (alldiv);
mpz_init (others);
+ mpz_init (root);
mpz_init (n);
+ mpz_init (rop);
/* product of all divisors */
mpz_set_ui (alldiv, 1L);
@@ -67,20 +69,34 @@
for (j = 1; j <= divisor[i]; j++)
{
/* square */
- mpz_mul_ui (n, others, j);
- mpz_mul (n, n, n);
+ mpz_mul_ui (root, others, j);
+ mpz_mul (n, root, root);
if (! mpz_perfect_square_p (n))
{
printf ("mpz_perfect_square_p got 0, want 1\n");
mpz_trace (" n", n);
abort ();
}
+ if (! mpz_perfect_square (rop, n))
+ {
+ // TODO: REVIEW: Why mpz_trace vs gmp_printf?
+ printf ("mpz_perfect_square got 0, want 1\n");
+ mpz_trace (" n", n);
+ abort ();
+ }
+ if (! mpz_cmp (rop, n))
+ {
+ gmp_printf ("mpz_perfect_square rop %Zd, want %Zd\n", rop, n);
+ abort ();
+ }
}
}
mpz_clear (alldiv);
mpz_clear (others);
+ mpz_clear (root);
mpz_clear (n);
+ mpz_clear (rop);
}
@@ -88,50 +104,80 @@
void
check_sqrt (int reps)
{
- mpz_t x2, x2t, x;
+ mpz_t x2, x2t, x, rop;
mp_size_t x2n;
int res;
int i;
- /* int cnt = 0; */
+ int want;
+ int cnt = 0;
gmp_randstate_ptr rands = RANDS;
mpz_t bs;
mpz_init (bs);
mpz_init (x2);
+ mpz_init (x2t);
mpz_init (x);
- mpz_init (x2t);
+ mpz_init (rop);
for (i = 0; i < reps; i++)
{
mpz_urandomb (bs, rands, 9);
x2n = mpz_get_ui (bs);
mpz_rrandomb (x2, rands, x2n);
- /* mpz_out_str (stdout, -16, x2); puts (""); */
+
+ mpz_sqrt (x, x2);
+ mpz_mul (x2t, x, x);
+ want = mpz_cmp(x2, x2t) == 0;
res = mpz_perfect_square_p (x2);
- mpz_sqrt (x, x2);
- mpz_mul (x2t, x, x);
-
- if (res != (mpz_cmp (x2, x2t) == 0))
+ if (res != want)
{
printf ("mpz_perfect_square_p and mpz_sqrt differ\n");
mpz_trace (" x ", x);
mpz_trace (" x2 ", x2);
mpz_trace (" x2t", x2t);
printf (" mpz_perfect_square_p %d\n", res);
- printf (" mpz_sqrt %d\n", mpz_cmp (x2, x2t) == 0);
+ printf (" mpz_sqrt %d\n", want);
abort ();
}
- /* cnt += res != 0; */
+ res = mpz_perfect_square (rop, x2);
+ if (res != want)
+ {
+ printf ("mpz_perfect_square and mpz_sqrt differ\n");
+ mpz_trace (" x ", x);
+ mpz_trace (" x2 ", x2);
+ mpz_trace (" x2t", x2t);
+ printf (" mpz_perfect_square %d\n", res);
+ printf (" mpz_sqrt %d\n", want);
+ abort ();
+ }
+ if (res && mpz_cmp(x, rop) != 0)
+ {
+ printf ("mpz_perfect_square and mpz_sqrt differ\n");
+ mpz_trace (" x ", x);
+ mpz_trace (" x2 ", x2);
+ mpz_trace (" x2t", x2t);
+ mpz_trace (" mpz_perfect_square rop", rop);
+ mpz_trace (" mpz_sqrt", x);
+ abort ();
+ }
+
+ cnt += res != 0;
}
+
+ if (reps > 1000 && cnt == 0) {
+ printf("No perfect squares found in %d reps", reps);
+ abort();
+ }
/* printf ("%d/%d perfect squares\n", cnt, reps); */
mpz_clear (bs);
mpz_clear (x2);
mpz_clear (x);
mpz_clear (x2t);
+ mpz_clear (rop);
}
More information about the gmp-devel
mailing list