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

Niels Möller nisse at
Fri Jun 9 10:04:01 UTC 2017

"Marco Bodrato" <bodrato at> 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

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. 


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

More information about the gmp-devel mailing list