Anomaly in mpn_sqrtrem and mpn_rootrem
Niels Möller
nisse at lysator.liu.se
Wed Jul 8 06:31:24 UTC 2015
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.
More information about the gmp-devel
mailing list