# HG changeset patch # User Seth Troisi # Date 1780469697 25200 # Tue Jun 02 23:54:57 2026 -0700 # Node ID 1a9d624cfb771f5a9e6e2125386d499e863f8b10 # Parent cf985f9efb5d57a8029bfe4eadf5c652b8d918c8 Added mpz_perfect_square_root Extract initial checks into mpn_perfect_square_trivial_p diff -r cf985f9efb5d -r 1a9d624cfb77 ChangeLog --- a/ChangeLog Wed Jun 03 00:07:28 2026 -0700 +++ b/ChangeLog Tue Jun 02 23:54:57 2026 -0700 @@ -1,3 +1,8 @@ +2026-06-01 Seth Troisi + * mpz/perfsqrt.c (mpz_perfect_square_root): New function with root return. + * mpn/perfsqr.c (mpn_probably_perfect_square_p): Checks used before sqrt. + * tests/mpz/t-perfsqr.c: Add tests for mpz_perfect_square_root + 2026-04-04 Marc Glisse * gmpxx.h (operator""): Remove deprecated space. diff -r cf985f9efb5d -r 1a9d624cfb77 Makefile.am --- a/Makefile.am Wed Jun 03 00:07:28 2026 -0700 +++ b/Makefile.am Tue Jun 02 23:54:57 2026 -0700 @@ -207,7 +207,8 @@ mpz/millerrabin$U.lo mpz/mod$U.lo mpz/mul$U.lo mpz/mul_2exp$U.lo \ mpz/mul_si$U.lo mpz/mul_ui$U.lo \ mpz/n_pow_ui$U.lo mpz/neg$U.lo mpz/nextprime$U.lo \ - mpz/out_raw$U.lo mpz/out_str$U.lo mpz/perfpow$U.lo mpz/perfsqr$U.lo \ + mpz/out_raw$U.lo mpz/out_str$U.lo mpz/perfpow$U.lo \ + mpz/perfsqr$U.lo mpz/perfsqrt$U.lo \ mpz/popcount$U.lo mpz/pow_ui$U.lo mpz/powm$U.lo mpz/powm_sec$U.lo \ mpz/powm_ui$U.lo mpz/primorial_ui$U.lo \ mpz/pprime_p$U.lo mpz/random$U.lo mpz/random2$U.lo \ diff -r cf985f9efb5d -r 1a9d624cfb77 doc/gmp.texi --- a/doc/gmp.texi Wed Jun 03 00:07:28 2026 -0700 +++ b/doc/gmp.texi Tue Jun 02 23:54:57 2026 -0700 @@ -3602,6 +3602,17 @@ be perfect squares. @end deftypefun +@deftypefun int mpz_perfect_square_root (mpz_t @var{rop}, const mpz_t @var{op}) +@cindex Perfect square functions +@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 +@var{rop} will be unchanged or set to @m{\lfloor\sqrt{@var{op}}\rfloor, the +truncated integer part of the square root of @var{op}} +@end deftypefun + @need 2000 @node Number Theoretic Functions, Integer Comparisons, Integer Roots, Integer Functions @@ -5635,8 +5646,7 @@ case the return value is zero or non-zero according to whether the remainder would have been zero or non-zero. -A return value of zero indicates a perfect square. See also -@code{mpn_perfect_square_p}. +A return value of zero indicates a perfect square. See also @code{mpn_perfect_square_p}. @end deftypefun @deftypefun size_t mpn_sizeinbase (const mp_limb_t *@var{xp}, mp_size_t @var{n}, int @var{base}) @@ -5719,8 +5729,13 @@ @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. +@end deftypefun + +@deftypefun int mpn_probably_perfect_square_p (const mp_limb_t *@var{s1p}, mp_size_t @var{n}) +Return non-zero iff @{@var{s1p}, @var{n}@} passes preliminary perfect square checks i.e. +if the return is zero, the input can't be a perfect square. +The most significant limb of the input @{@var{s1p}, @var{n}@} must be non-zero. @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 +9265,10 @@ A significant fraction of non-squares can be quickly identified by checking whether the input is a quadratic residue modulo small integers. -@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. +@code{mpz_perfect_square_p}/@code{mpz_perfect_square_root} 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 cf985f9efb5d -r 1a9d624cfb77 doc/tasks.html --- a/doc/tasks.html Wed Jun 03 00:07:28 2026 -0700 +++ b/doc/tasks.html Tue Jun 02 23:54:57 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. -
  • 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.
  • mpz_get_ull, mpz_set_ull, mpz_get_sll, mpz_get_sll: Conversions for long long. These would aid interoperability, though a diff -r cf985f9efb5d -r 1a9d624cfb77 gmp-h.in --- a/gmp-h.in Wed Jun 03 00:07:28 2026 -0700 +++ b/gmp-h.in Tue Jun 02 23:54:57 2026 -0700 @@ -972,6 +972,9 @@ __GMP_DECLSPEC int mpz_perfect_square_p (mpz_srcptr) __GMP_ATTRIBUTE_PURE; #endif +#define mpz_perfect_square_root __gmpz_perfect_square_root +__GMP_DECLSPEC int mpz_perfect_square_root (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_probably_perfect_square_p __MPN(probably_perfect_square_p) +__GMP_DECLSPEC int mpn_probably_perfect_square_p (mp_srcptr, mp_size_t) __GMP_ATTRIBUTE_PURE; + #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 cf985f9efb5d -r 1a9d624cfb77 mpn/generic/perfsqr.c --- a/mpn/generic/perfsqr.c Wed Jun 03 00:07:28 2026 -0700 +++ b/mpn/generic/perfsqr.c Tue Jun 02 23:54:57 2026 -0700 @@ -1,8 +1,12 @@ /* 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_probably_perfect_square_p(u, usize) -- Returns non-zero if U not trivially +disqualified from being a perfect square checks. Returns zero if impossible +to be a perfect square. + +Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2012, 2026 Free +Software Foundation, Inc. This file is part of the GNU MP Library. @@ -176,12 +180,10 @@ int -mpn_perfect_square_p (mp_srcptr up, mp_size_t usize) +mpn_probably_perfect_square_p (mp_srcptr up, mp_size_t usize) { ASSERT (usize >= 1); - TRACE (gmp_printf ("mpn_perfect_square_p %Nd\n", up, usize)); - /* The first test excludes 212/256 (82.8%) of the perfect square candidates in O(1) time. */ { @@ -217,6 +219,17 @@ according to their residues modulo small primes (or powers of primes). See perfsqr.h. */ PERFSQR_MOD_TEST (up, usize); + return 1; +} + + +int +mpn_perfect_square_p (mp_srcptr up, mp_size_t usize) +{ + TRACE (gmp_printf ("mpn_perfect_square_p %Nd\n", up, usize)); + + if (! mpn_probably_perfect_square_p (up, usize)) + return 0; /* For the third and last test, we finally compute the square root, diff -r cf985f9efb5d -r 1a9d624cfb77 mpz/Makefile.am --- a/mpz/Makefile.am Wed Jun 03 00:07:28 2026 -0700 +++ b/mpz/Makefile.am Tue Jun 02 23:54:57 2026 -0700 @@ -57,7 +57,7 @@ lucnum_ui.c lucnum2_ui.c lucmod.c mfac_uiui.c millerrabin.c \ mod.c mul.c mul_2exp.c mul_si.c mul_ui.c n_pow_ui.c neg.c nextprime.c \ oddfac_1.c \ - out_raw.c out_str.c perfpow.c perfsqr.c popcount.c pow_ui.c powm.c \ + out_raw.c out_str.c perfpow.c perfsqr.c perfsqrt.c popcount.c pow_ui.c powm.c \ powm_sec.c powm_ui.c pprime_p.c prodlimbs.c primorial_ui.c random.c random2.c \ realloc.c realloc2.c remove.c roinit_n.c root.c rootrem.c rrandomb.c \ scan0.c scan1.c set.c set_d.c set_f.c set_q.c set_si.c set_str.c \ diff -r cf985f9efb5d -r 1a9d624cfb77 mpz/perfsqrt.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mpz/perfsqrt.c Tue Jun 02 23:54:57 2026 -0700 @@ -0,0 +1,82 @@ +/* mpz_perfect_square_root(rop, op) -- Return non-zero if op is a perfect + square and sets rop to root, otherwise returns zero. + +Copyright 2026 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" + +int +mpz_perfect_square_root (mpz_ptr rop, mpz_srcptr op) +{ + mp_size_t u_size = SIZ (op); + + /* Handle <= 0 */ + if (__GMP_UNLIKELY ( u_size <= 0)) + { + /* Zero is a square with rop = 0 */ + if ( u_size == 0) + { + SIZ(rop) = 0; + return 1; + } + /* Negatives are never perfect squares. */ + return 0; + } + + if (! mpn_probably_perfect_square_p (PTR (op), u_size)) + return 0; + + /* This is the precise size of sqrt(op) because leading limp is non-zero. */ + mp_size_t rop_size = (SIZ (op) + 1) / 2; + + /* mpn_sqrtrem doesn't allow sp == np */ + if (rop == op) + { + TMP_DECL; + TMP_MARK; + + mp_ptr rop_ptr = TMP_ALLOC_LIMBS (rop_size); + int res = ! mpn_sqrtrem (rop_ptr, NULL, PTR (op), u_size); + + /* It might be better to always copy so the value in rop doesn't + depend on if the arguments are aliased or not. */ + if (res) + { + MPN_COPY (PTR(rop), rop_ptr, rop_size); + SIZ(rop) = rop_size; + } + + TMP_FREE; + return res; + } + + /* Make sure rop has exact space for the square root. */ + SIZ(rop) = rop_size; + return ! mpn_sqrtrem (MPZ_REALLOC (rop, rop_size), NULL, PTR (op), u_size); +} diff -r cf985f9efb5d -r 1a9d624cfb77 tests/mpz/t-perfsqr.c --- a/tests/mpz/t-perfsqr.c Wed Jun 03 00:07:28 2026 -0700 +++ b/tests/mpz/t-perfsqr.c Tue Jun 02 23:54:57 2026 -0700 @@ -1,4 +1,4 @@ -/* Test mpz_perfect_square_p. +/* Test mpz_perfect_square_p and mpz_perfect_square_root. Copyright 2000-2002 Free Software Foundation, Inc. @@ -26,10 +26,10 @@ #include "mpn/perfsqr.h" -/* check_modulo() exercises 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. +/* check_modulo() exercises mpz_perfect_square_p/mpz_perfect_square_root 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. Each divisor is taken separately. It's arranged that n is congruent to 0 modulo the other divisors, 0 of course being a quadratic residue to any @@ -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,71 +69,164 @@ 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_root (rop, n)) + { + printf ("mpz_perfect_square_root got 0, want 1\n"); + mpz_trace (" n", n); + abort (); + } + if (! mpz_cmp (rop, n)) + { + gmp_printf ("mpz_perfect_square_root rop %Zd, want %Zd\n", rop, n); + abort (); + } } } mpz_clear (alldiv); mpz_clear (others); + mpz_clear (root); + mpz_clear (n); + mpz_clear (rop); +} + +/* Check negative, 0, and 1 */ +void +check_edge_cases (void) +{ + mpz_t n, root; + mpz_init (n); + mpz_init (root); + + mpz_set_si(n, -4); + if (mpz_perfect_square_p (n)) + { + printf ("mpz_perfect_square_p -4 not a square\n"); + abort (); + } + if (mpz_perfect_square_root (root, n)) + { + printf ("mpz_perfect_square_root -4 not a square\n"); + abort (); + } + + // 0 and 1 are both perfect squares with respective root 0 and 1. + for (unsigned int m = 0; m < 2; m++) + { + mpz_set_ui(n, m); + if (! mpz_perfect_square_p (n)) + { + printf ("mpz_perfect_square_p %u should be a square\n", m); + abort (); + } + if (! mpz_perfect_square_root (root, n) || mpz_cmp_ui (root, m) != 0) + { + printf ("mpz_perfect_square_root %u should have root = %u\n", m, m); + abort (); + } + } + + mpz_clear (root); mpz_clear (n); } - /* Exercise mpz_perfect_square_p compared to what mpz_sqrt says. */ 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_root (rop, x2); + if (res != want) + { + printf ("mpz_perfect_square_root and mpz_sqrt differ\n"); + mpz_trace (" x ", x); + mpz_trace (" x2 ", x2); + mpz_trace (" x2t", x2t); + printf (" mpz_perfect_square_root %d\n", res); + printf (" mpz_sqrt %d\n", want); + abort (); + } + if (res && mpz_cmp(x, rop) != 0) + { + printf ("mpz_perfect_square_root and mpz_sqrt differ\n"); + mpz_trace (" x ", x); + mpz_trace (" x2 ", x2); + mpz_trace (" x2t", x2t); + mpz_trace (" mpz_perfect_square_root rop", rop); + mpz_trace (" mpz_sqrt", x); + abort (); + } + /* Check that same variable as input and output works */ + mpz_set(rop, x2); + res = mpz_perfect_square_root (rop, rop); + if (res != want || (res && mpz_cmp(x, rop) != 0)) + { + printf ("mpz_perfect_square_root differ when output == input\n"); + mpz_trace (" x2 ", x2); + printf (" mpz_perfect_square_root %d\n", res); + mpz_trace (" mpz_perfect_square_root rop ", rop); + 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); } @@ -146,6 +241,7 @@ if (argc == 2) reps = atoi (argv[1]); + check_edge_cases (); check_modulo (); check_sqrt (reps);