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

[Niels Möller]

"Marco Bodrato" <bodrato at mail.dm.unipi.it> writes:

> The constant (c) I added to the first iteration should be incorporated in
> T1, and of course it is not necessary to really use a constant. A
> different value can be used for each entry in the table...

In Newton iteration, one usually converges from one side, and it's
tempting (to ease analysis and enable use of unsigned arithmetic) to
maintain a known and fixed sign of the error. But one could gain about
one bit of accuracy with error range roughly centered around zero.

For the initial iteration with tabulated values, you'd adjust each table
entry, as you say. In later iterations, one could add a suitably choosen
constant if one is willing to handle signed adjustments cancellation
logic.

When I worked on reciprocal quite a while ago, I aimed for a signed
symmetric error for the initial table lookup, but for the main
iterations, I designed for a non-negative error and I didn't explore
using signed errors.

> I believe that we can define
>
> 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 */

One could save an instruction with using just

  l(a) = (a >> GMP_NUMB_BITS - 20)

and adjusting the tables accordingly. You'd get a larger multiply, but
limiting to 16x16 is probably not a saving on relevant platforms. 

Regards,
/Niels

-- 
Niels Möller. PGP-encrypted email is preferred. Keyid 368C6677.
Internet email is subject to wholesale government surveillance.