mpz_sqrt_if_perfect_square

[Niels Möller]

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.