mpn_sqrtrem{1,2} - patch for pure C implem

[Marco Bodrato]

Ciao,

Il Sab, 1 Aprile 2017 9:02 pm, Adrien Prost-Boucle ha scritto:
> The new diff is at the end of this message.

I got inspired by Adrien's code and by Paul's one, and I wrote another
candidate replacement for mpn_sqrtrem{1,2}, working for both 32 and 64
bits.

My goals:
 - new code should not be slower than the current :-)
 - get rid of all signed shifts,
 - get rid of all mp_limb_signed_t needs,
 - keep-on using the current invsqrttab (optimal IMO).

I decided to use the iteration (for i an approximation of a^{-1/2})
 i' = i*(3-a*i*i)/2
because it does not need handling signed numbers.

This iteration nearly doubles the precision of the approximation.
Some problems arise because of the "nearly"...

invsqrttab allows us to start with an approximation of 8 bits and little
more (I computed an error smaller than 1/355 < 2^-8.47).
A little perturbation on the first iteration
 i' = i*(3+c-a*i*i)/2 , with c = 3/(2*356^2)
(the opposite sign wrt the current MAGIC constant...) give an error
smaller than 1/150000 < 2^-17; the following iterations should be able to
give at least the desired precision. (17 -> 33 -> 65 ...).

For some computations, even for the 64-bits version, a 32x32->32
multiplication is enough. I #define-d the gmp_uint_fast32_t type, but
currently it is not really guaranteed "fast"...

I'll be happy to receive comments or critics :-)

Best regards,
m

-- 
http://bodrato.it/



-------------- next part --------------
A non-text attachment was scrubbed...
Name: sqrtrem.patch
Type: text/x-patch
Size: 5787 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20170606/b5fac1b4/attachment.bin>