About sqrt_exact...

[marco.bodrato at tutanota.com]

Ciao,

17 lug 2026, 11:49 da marco.bodrato at tutanota.com:

>> Can this be faster than the current code? I don't know.
>>
> It is, but maybe too slightly faster for such a complex code, is it worth adding?Here are the timings on my laptop:
>
Thanks to an improvement in the function mpn_bsqrt, the last version is even faster:


> $ build/tune/speed -p100000000 -c -s 1-1600 -f1.5 mpn_perfect_square_p
> overhead 5.25 cycles, precision 100000000 units of 5.28e-10 secs, CPU freq 1895.63 MHz
>
                last           after          before
1              68.28           68.43           70.22	    
2             135.65          133.89          128.77	    
3             199.28          189.03          319.57	    
4             191.37          185.58          272.31	    
6             329.87          318.39          414.99	    
9             480.41          575.23          854.72	    
13            647.84          791.94          967.31	    
19            902.47         1203.43         1448.53	    
28           1401.42         2038.94         1763.29	    
42           2176.32         2859.41         2804.96	    
63           3884.83         4504.65         4614.51	    
94           6339.52         6674.98         6774.17	    
141         10145.48        10713.67        12355.94	    
211         16610.33        17956.97        21342.91	    
316         28961.38        31805.80        39641.78	    
474         50867.42        57801.74        74848.89	    
711         98423.81       112779.93       147468.30	    
1066       205105.03       233826.29       280769.99	    
1599       397119.31       452747.42       529568.46


> But a lot of work is still needed, because:
>  - the code passes "make check" and also some deeper testing, but... it is not completely correct, there are some possible corner cases that are not handled yet;
>
Meanwhile, I added a test that my code does not pass ;-)

The problem is: we compute some binary digits from above and some from below.
On both ends there is a small uncertainty,
from above, the number we obtain is "maybe too large, by 1";
from below, the digits we obtain can be negated.
The current code computes (at least) GMP_NUMB_BITS-3 overlapping digits.

That's where the corner cases came from:
assume the above digits  are "0x000...000",
but maybe the estimate is too large and
if we subtract subtract 1, we get "0xFFF..FFF".
on the other side, assume the digits are "0xFFF..FFF",
but we can negate them and they can be "0x000...000".

We don't have enough information to decide,
which sequence should we choose.

Of course other similar combinations are possible.

A possible strategy to solve this is to detect that kind of corner cases (easy) and handle it with some special code (Try one of the two alternatives, square and check, if wrong, try the other alternative, square again, if wrong again, then return "non square"), but I did not implement it yet. The average time should not be affected.

> If someone wants to test this code, and maybe reply with some comment, I'll be happy.Is it worth working on that approach?
> Or it's too risky to add a complex piece of code for a marginal speed-up, maybe?
>

I attach the latest code.

PS: with this patch, if the input is not a perfect square, the value stored in rop by mpz_perfect_square_root is really undefined :-)
Ĝis,
m
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sqrt_exact.diff
Type: text/x-patch
Size: 14892 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20260720/0eec1b45/attachment.bin>