mpz_sqrt_if_perfect_square

[Seth Troisi]

Thanks for the comments I've addressed most of them in this version

> The top limb could be zero, right?
The docs say "The most significant limb of @{@var{s1p}, @var{n}@} must be
non-zero."
I believe this makes the size of sqrt(u) exactly (n+1)/ 2 for any perfect
square.

I'm was unsure how much I should optimize for an extra memory allocation vs
clobbering rop.
1. The fastest code would allow clobbering rop and avoids TMP_ALLOC by
exposing the pre-checks and has mpz/perfsqr.c call MPZ_REALLOC before
calling mpn_sqrtrem.
2. Alternatively I could reimplement MPZ_REALLOC  inside mpn/perfsqr.c when
rop isn't large enough to store the root.
3. To avoiding clobbering rop, I could do a MPZ_REALLOC(rop, (op_size + 1)
/ 2) if rop is smaller, TMP_ALLOC in the non-trivial case, and only do
MPN_COPY if op is a perfect square.

I went with 1 and added mpn_perfect_square_trivial_p along with docs in to
gmp.texi, it might benefit from a better name
(mpn_quick_not_perfect_square_p).
I suspect I'm supposed to mark this as internal only or possibly break it
out of the perfsqr.c file. Let me know if I need to organize the code
differently.

As some extra due-diligence I added a test for -4, 0, 1 to t-perfsqr and
ran speed before and after (both in debug mode) and didn't see any
regression.


On Tue, Jun 2, 2026 at 7:53 AM Niels Möller <nisse at lysator.liu.se> wrote:

> Seth Troisi <braintwo at gmail.com> writes:
>
> > I'm doing some prime gap searching which would slightly benefit from
> > knowing the root of a number iff it's a perfect square.
> ...
> > I coded up a working implementation for myself. I can work on making it
> > clean enough for upstream if it's a desired feature.
>
> Makes sense to me. Below a few initial comments after a quick first look.
>
> > + 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
>
> Would it be difficult to arrange so rop is unchanged in the case that
> input is not a square? "undefined" is a bit scary.
>
> > + 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
>
> Would be good to document if output area is unchanged or clobbered in
> the non-square case. Clobbered would be fine for an mpn function, but
> unchanged would likely be helpful if we want "unchanged" for the mpz
> function.
>
> > +    {
> > +      /* Iff mpn_sqrtrem returns zero, the square is perfect.  */
> > +      int res = ! mpn_sqrtrem (rp, NULL, up, usize);
> > +      return res;
> > +    }
>
> Would be nice to either move declaration of res and the return statement
> out
> of the two if blocks, or write this simply as
>
>   return ! mpn_sqrtrem (rp, NULL, up, usize);.
>
> > +  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 I understand the internal apis correctly, this can be simplified to
>
>   res = mpn_perfect_square (MPZ_REALLOC(rop, rop_size), PTR (u), u_size);
>
> > +  if (res)
> > +    {
> > +      SIZ(rop) = rop_size;
>
> The top limb could be zero, right? Then a normalization step is needed
> here. If it is guaranteed that at most one top limb is zero, something
> like
>
>          rop_size -= PTR(rop)[rop_size-1] == 0;
>
> Regards,
> /Niels
>
> --
> Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677.
> Internet email is subject to wholesale government surveillance.
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mpn_perfect_square.PATCH
Type: application/octet-stream
Size: 15518 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20260602/8a996bcf/attachment.obj>