mpz_sqrt_if_perfect_square

[marco.bodrato at tutanota.com]

Ciao Seth,

27 giu 2026, 09:34 da braintwo at gmail.com:

> On Tue, Jun 9, 2026 at 10:29 AM <> marco.bodrato at tutanota.com> > wrote:
>
>> Ciao,
>>
>> 6 giu 2026, 08:45 da >> nisse at lysator.liu.se>> :
>>
>>> Seth Troisi <>>> braintwo at gmail.com>>> > writes:
>>>
>>>> I changed the name.
>>>> Would mpn_probab_perfect_square_p (no ly) better match mpz_probab_prime? 
>>>>
> I've removed it from > gmp-h.in <http://gmp-h.in>>  and the docs.
>

But then it does not compile!

../../mpz/perfsqrt.c: In function '__gmpz_perfect_square_root':
../../mpz/perfsqrt.c:52:9: error: implicit declaration of function 'mpn_probab_perfect_square_p'; did you mean 'mpn_perfect_square_p'? [-Wimplicit-function-declaration]

I suggest:
diff -r 05ed7e126732 gmp-impl.h
--- a/gmp-impl.h        Fri Jun 26 13:34:23 2026 +0200
+++ b/gmp-impl.h        Sat Jun 27 18:06:28 2026 +0200
@@ -1208,6 +1208,9 @@
#define mpn_mulmid __MPN(mulmid)
__GMP_DECLSPEC void mpn_mulmid (mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
 
+#define mpn_probab_perfect_square_p __MPN(probab_perfect_square_p)
+__GMP_DECLSPEC int mpn_probab_perfect_square_p (mp_srcptr, mp_size_t);
+
#define mpn_submul_1c __MPN(submul_1c)
__GMP_DECLSPEC mp_limb_t mpn_submul_1c (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t, mp_limb_t);


On the docs side:

>>> If docs say that value is undefined, then I don't think it's worth an
>>> allocation to make it have some particular value.
>>>
>> I agree.
>>
>>> If I understand your code correctly, then in the non-square case, rop
>>> will either be unchanged or get floor (sqrt(op)), depending on result of
>>>
>> Maybe, but it doesn't really matter.
>> A new implementation using a different algorithm may "undefine" the undefined result value in a different way :-)
>>
The function is expected to return the square root if it exists,
if it does not exist, nothing is expected.

We should not artificially create expectations by documenting the
arbitrary behavior of the current implementation.

After the words
"If @var{op} is not a perfect square the value of @var{rop}"
I'd write
"is undefined"

We may want to add
", but its allocated space will not grow larger than the expected
size of the square root"
If we really want to promise something that any future implementations
should be able to respect.


Always on the doc side, but not related to your proposed patch.

I noticed that the function mpn_perfect_power_p is not documented.
Should we add something like the following?

diff -r 05ed7e126732 doc/gmp.texi
--- a/doc/gmp.texi      Fri Jun 26 13:34:23 2026 +0200
+++ b/doc/gmp.texi      Sat Jun 27 18:30:30 2026 +0200
@@ -5722,7 +5722,8 @@
@end deftypefun
 
@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.
+ at deftypefunx int mpn_perfect_power_p (const mp_limb_t *@var{s1p}, mp_size_t @var{n})
+Return non-zero iff @{@var{s1p}, @var{n}@} is a perfect square or power.
The most significant limb of the input @{@var{s1p}, @var{n}@} must be
non-zero.
@end deftypefun



> Speed.h has this comment for SPEED_ROUTINE_MPN_PERFECT_SQUARE
>
> /* Calculate worst case (larger prime) for perfect_square */
>

Yes, I know, and I'm not sure which is the prime it is talking about...
Unrelated to your proposed patch, but...
...should we substitute "(larger prime)" with "(actually a square)"?


> If I change the loop body in increment r by 1 on each loop I get these results
>
> $ ./tune/speed  -r -s 100-2000 -p 100000000 -f 2 mpn_sqrtrem mpn_perfect_square_p
> overhead 0.000000001 secs, precision 100000000 units of 2.51e-10 secs, CPU freq 3980.51 MHz
>           mpn_sqrtrem mpn_perfect_square_p
> 100       0.000001297       #0.0110
> 200       0.000003405       #0.0075
> 400       0.000010291       #0.0056
> 800       0.000030102       #0.0051
> 1600      0.000090388       #0.0044
>
Yes of course, they are expected. tune/speed is not a static tool,
it is adapted to the needs of the developers...
It is also interesting to ensure that every value will pass the initial trivial test...
diff -r 05ed7e126732 tune/speed.h
--- a/tune/speed.h      Fri Jun 26 13:34:23 2026 +0200
+++ b/tune/speed.h      Sat Jun 27 19:00:21 2026 +0200
@@ -3702,7 +3702,10 @@
     speed_starttime ();                                                        \
     i = s->reps;                                                       \
     do                                                                 \
-      function (PTR(r), SIZ(r));                                       \
+      {                                                                        \
+       function (PTR(r), SIZ(r));                                      \
+       *PTR(r) = ((*PTR(r)|7)^6)+8;                                     \
+      }                                                                        \
     while (--i != 0);                                                  \
     t = speed_endtime ();                                              \
                                                                        \

Anyway, an even more interesting timing,
with the current not patched tune/speed, is the following:

$ build/tune/speed -p10000000 -r -s 100-2000 -f2 mpn_sqrt mpn_perfect_square_p
overhead 0.000000003 secs, precision 10000000 units of 5.28e-10 secs, CPU freq 1895.63 MHz
             mpn_sqrt mpn_perfect_square_p
100      #0.000002323        1.7285
200      #0.000005981        1.7509
400      #0.000017727        1.6575
800      #0.000057803        1.6152
1600     #0.000179633        1.5603

And I realized that mpn_sqrt itself is faster on general numbers (non-squares),
and slower when its timing is measured on squares only...
Ĝis,
m

-- 
DC: https://vado.li/dc-bodrato