Anomaly in mpn_sqrtrem and mpn_rootrem

[bodrato at mail.dm.unipi.it]

Ciao,

Il Mer, 24 Giugno 2015 10:25 pm, Niels Möller ha scritto:
> bodrato at mail.dm.unipi.it writes:
>> If we need both the sqrt of the high-half and the residual, we should
>> directly use sqrtrem, don't we?
>
> Maybe, but it's not obvious. For simplicity, consider the case of
> smallish numbers, so we use mullo rather than mulmod_bnm1 for the
> cancellation. Then H (high half of the square root) depends only on the
> high half of A, while the residual is computed from H and the low half
> of A. So I see no obvious gain in computing the residual sooner.

Not a gain in speed, probably, but we can avoid writing code twice.
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 do not know which algorithm is the faster to obtain both H and the
residual, but we should implement it just once, in sqrtrem.

>> Once we compute a residual, correcting the
>> "approximation" to obtain an "exact" result does not cost too much,
>> right?
>
> Exact for the current A input, yes. But one level up in the recursion,
> where additional low limbs are used, it's no longer exact. Nothing is

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.

> very obvious here, but my gut feeling is that "exactness" isn't a very
> useful property for the intermediate results based on truncated versions
> of the top-level input. And it's best to defer the adjustment to
> top-level.

Adjustment is a little bit more complex if the previous intermediate
result was not exact, but you are probably right: a unique top-level
adjustment might be the best strategy.
Actually I opted for a simpler reuse-already-working-code strategy :-)

Best regards,
m

-- 
http://bodrato.it/papers/