sqrt algorithm (was: Re: Anomaly in mpn_sqrtrem and mpn_rottrem)

[Niels Möller]

nisse at lysator.liu.se (Niels Möller) writes:

> So back to the drawing board.

I had to redo the splitting a bit in sqrt_nm1, and write some special
case code for n = 3. But now I have code that survives some testing.

I'm attaching my current algorithm description and code. The code
currently doesn't really exploit cancellation. Also, it computes
divisions like 

  floor (B^{k-1} E / 2H)

by zero-padding, and calling mpn_tdiv_qr. 

To make this fast, we need some variant of divappr_q which don't require
any of the uninteresting low limbs. Or alternatively, resurrect the
notion of fraction limbs.

Compare to bdiv_q functions, which computes N / D mod B^n, where both
inputs and the output are n limbs. A divappr_q function ought to compute
an approximation of N B^n / D, but possibly with a little mmore
flexibility in the input and output sizes.

The most flexible way is perhaps to mimic the "fraction
limbs"-interface. Say we have inputs N = {np, nn}, D = {dp, dn}, and a
scaling factor k, and compute an approximation of

   Q = floor (N B^k / D)

of qn = nn + dn + 1 + k limbs. Or we could tie things together a bit
more, by requiring that dn = nn = qn + 1 or so.

For the divisions in the sqrt algorithms, I'm not sure of exactly how
the sizes of the numerator E, the denominator H, and the quotient,
relate, but they ought to all be pretty close to k.

Regards,
/Niels

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: sqrt.txt
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20150712/126d84bf/attachment.txt>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sqrt.c
Type: text/x-csrc
Size: 11048 bytes
Desc: not available
URL: <https://gmplib.org/list-archives/gmp-devel/attachments/20150712/126d84bf/attachment.bin>
-------------- next part --------------

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