Anomaly in mpn_sqrtrem and mpn_rottrem

[bodrato at mail.dm.unipi.it]

Ciao,

Il Ven, 12 Giugno 2015 9:04 am, Torbjörn Granlund ha scritto:
> We might want to look into the plain division-free sqrt(A) = A*sqrt(1/A)
> approach before implementing a tricky division sqrt(A).

We can try improving the current implementation, before implementing any
other algorithm :-)

I wrote a simple patch (it touches very few lines) that allows skipping
the final squaring when mpn_sqrtrem is called with a NULL argument and the
approximation computed so far can not change.
It exploits the spurious half-a-limb that current code adds to the result
when the size is odd, i.e. it changes the timings only for odd sizes.

Before the patch:

$ tune/speed -rp100000000 -s21-1100000 -f15 mpn_sqrt mpn_root.2
mpn_rootrem.2 mpn_sqrtrem
overhead 0.000000002 secs, precision 100000000 units of 2.86e-10 secs, CPU
freq 3500.08 MHz
             mpn_sqrt    mpn_root.2 mpn_rootrem.2   mpn_sqrtrem
21        0.000001380        1.6155        1.5569       #0.9996
315      #0.000015374        1.1222        1.5426        1.1521
4725      0.000918521       #0.9000        1.2458        1.0004
70875     0.029574741       #0.9612        1.2192        1.0007
1063125    0.741147157       #0.9718        1.2147        0.9999

Sometimes mpn_root.2 is faster than the other functios in a wide range

With the patch applied:
$ tune/speed.nsqrt -rp100000000 -s21-1100000 -f15 mpn_sqrt mpn_root.2
mpn_rootrem.2 mpn_sqrtrem
overhead 0.000000002 secs, precision 100000000 units of 2.86e-10 secs, CPU
freq 3500.08 MHz
             mpn_sqrt    mpn_root.2 mpn_rootrem.2   mpn_sqrtrem
21       #0.000001339        1.6496        1.6089        1.0244
315      #0.000012924        1.3318        1.8432        1.3701
4725     #0.000799607        1.0349        1.4312        1.1495
70875    #0.026321653        1.0804        1.3698        1.1253
1063125  #0.656389705        1.0982        1.3736        1.1296

mpn_rot.2 is the faster

-- 
http://bodrato.it/