Anomaly in mpn_sqrtrem and mpn_rootrem

[Niels Möller]

bodrato at mail.dm.unipi.it writes:

> Not a gain in speed, probably, but we can avoid writing code twice.

Using some variant of sqrtrem is one alternative for the organization. 

> If H=floor(sqrt(A/B^2n)), then the residual we need to compute the lower
> half of the square root is floor(A/B^2n)-H^2 .

I don't follow you here. Consider a simple case: A of size 2n, and we
try to compute an approximation of sqrt(B^{2n-2} A).

Then first compute high half as

  H \approx sqrt(B^{n-2} floor (B^{-n} A))

And the residual is
 
  E = H^2 - A (appr. n limbs, due to cancellation)

So the pieces of A needed for the two computations have very little
overlap, I haven't done the error analysis, but I expect at most one
limb, i.e., E is n+1 limbs, possibly signed.

For other odd/even cases, we also need residuals of the form

  E = B H^2 - A

and

  E = H^2 - B A

> Well, if H=floor(sqrt(floor(A/B^2n))), then H=floor(sqrt(A)/B^n).
> I mean, if the square root of the high half of the number is correctly
> rounded, it will coincide with the high half of the root, the additional
> low libs can not change it.

Hmmm, that's related to sqrt being concave? Computations involving floor
are a bit non-obvious to me. In my notation above, an "exact" H would be

  H = floor (sqrt (floor (A/B^2))) = floor (sqrt(A) / B)

Say we have such an exact H, what kind of rounding and adjustments do we
need after the adjustment step

  X <-- B^{n-1} H - B^{n-1} E / 2H    ; some rounding for the division

If it's possible to get a correct X = floor (sqrt (B^{n-2} A)) without
actually computing X^2, that makes an "exact" algorithm more attractive.

Regards,
/Niels

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