mpz_sqrt_if_perfect_square

[Niels Möller]

Seth Troisi <braintwo at gmail.com> writes:

> New patch using mpz_perfect_square_root name
> + at deftypefun int mpn_perfect_square_trivial_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.
> +The most significant limb of the input @{@var{s1p}, @var{n}@} must be non-zero.
>  @end deftypefun

Another function that needs the right name... I don't think "trivial" is
the right attribute. Maybe mpn_probably_perfect_square_p (even though
it's not really probabilistic) for consistency with probably_prime_p. We
need something to convey that it's fast and cheap (and not totally
accurate). Regardless of name, maybe say clarly that if it returns
false, the input is certainly not a square.

>  @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 +9263,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_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.

I think these docs should be attached to mpn_perfect_square_trivial_p.

> +int
> +mpz_perfect_square_root (mpz_ptr rop, mpz_srcptr u)
> +{
> +  mp_size_t u_size = SIZ (u);
> +
> +  /* Handle <= 0 */
> +  if (__GMP_UNLIKELY ( u_size <= 0))
> +    {
> +      /* Zero is a square with rop = 0 */
> +      if (__GMP_UNLIKELY ( u_size == 0))

Not sure this inner  __GMP_UNLIKELY is helpful, I'd leave to the branch
predictor to find out if zero or negative is more common.

> +        SIZ(rop) = 0;
> +        return 1;
> +      }
> +      /* Negatives are never perfect squares. */
> +      return 0;
> +    }
> +
> +  if (! mpn_perfect_square_trivial_p (PTR (u), u_size))
> +    return 0;
> +
> +  int rop_size = (SIZ (u) + 1) / 2;

Clearer with mp_size_t instead of int. And I think it may be worth a comment
saying that this size is precise (limb rop_size - 1 will always be non-zero).

> +  /* mpn_sqrtrem doesn't allow sp == np */
> +  if (rop == u)
> +    {
> +      TMP_DECL;
> +      TMP_MARK;
> +
> +      mp_ptr rop_ptr = TMP_ALLOC_LIMBS (rop_size);
> +      int res = ! mpn_sqrtrem (rop_ptr, NULL, PTR (u), u_size);
> +
> +      /* It might be better to always copy so the value in rop
> +         is more determistic. */

Unless I'm missing something, the current code leaves rop unchanged when
res == 0? That seems deterministic enough.

> +      if (res)
> +        {
> +          MPN_COPY (PTR(rop), rop_ptr, rop_size);
> +          SIZ(rop) = rop_size;
> +        }
> +
> +      TMP_FREE;
> +      return res;
> +    }
> +
> +  /* Make sure rop has space for the square root. */
> +  MPZ_REALLOC(rop, rop_size);
> +
> +  int res = ! mpn_sqrtrem (PTR (rop), NULL, PTR (u), u_size);

I think MPZ_REALLOC returns the pointer, so you could write this as 

     int res = ! mpn_sqrtrem (MPZ_REALLOC (rop, rop_size), NULL, PTR (u), u_size);

Or do you think it is clearer with MPZ_REALLOC more separated?

> +  if (res)
> +    {
> +      SIZ(rop) = rop_size;
> +    }
> +  return res;

In the non-square case, I think there's a risk that we produce a rop
that isn't normalized. To keep things simple, maybe just set SIZ(rop) =
0 in that case? I.e.,

  SIZ (rop) = res ? rop_size : 0;
  return res;

Regards,
/Niels

-- 
Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677.
Internet email is subject to wholesale government surveillance.