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

[Marco Bodrato]

Ciao,

Il Ven, 9 Giugno 2017 1:41 pm, paul zimmermann ha scritto:

>> h(a) = a >> GMP_NUMB_BITS - 9 /* the highest 9 bits */
>> l(a) = (a >> GMP_NUMB_BITS - 20) & ((1<<12)-1) /* the next 11 bits */
>> Finally we can totally skip the first approximation and directly compute
>>
>> i' = T1[h(a)] + l(a) * T2[h(a)]
>>
>> using a piece-wise linear function.

> Not sure l(a) is better than using the upper (say) 16 bits of a, in case
> you want a 16-bit approximation i'.

I'd prefer a 17 or 18-bit approximation ;-)

Niels too suggested to use all the bits and not only the lower ones.
He is right, this approach saves an operation.

But there are some side effects... I'd like to use a single table for
 - sqrt(32-bits) -> 16-bits
 - sqrt(2x32-bits) -> 32-bits
 - sqrt(64-bits) -> 32-bits
 - sqrt(2x64-bits) -> 64-bits
And I'd like to store both T1[x] and T2[x] into a single 32-bit entry...

I'll try which strategy gives better results.

> To compute T1 and T2 you have two choices:
>
> 1) either try to minimize the maximal error after the first Newton
> iteration
>    (i.e., on i')
>
> 2) or try to minimize the maximal error after k iterations. See the
> article
>    by Jean-Michel Muller and Peter Kornerup:
>    https://hal.archives-ouvertes.fr/inria-00071899/document

The article is interesting, but not applicable here, I believe...

We do not know in advance the number of iteration. I don't think it is
worth storing a different table for each function.
 - sqrt(32-bits) -> 16-bits
will directly use the value

 - sqrt(2x32-bits) -> 32-bits
 - sqrt(64-bits) -> 32-bits
use an iteration, not on the inverse square root, but on the root;

 - sqrt(2x64-bits) -> 64-bits
need one iteration on the inverse, and one on the root.

Moreover the piece-wise linear approximation approach means choosing two
parameters to obtain a reasonable approximation for a full bunch of
intervals, this is not a careful choice for each one of the intervals.

Last but not least, there is a third choice, the one Niels describe with
"maintain a known and fixed sign of the error". Generating an estimate
that is NOT IN the interval of possible final results, but outside... on
the wanted side :-)

The tables in the current code of MPFR probably lead to this third goal.

Regards,
m

-- 
http://bodrato.it/papers/