Anomaly in mpn_sqrtrem and mpn_rootrem

[bodrato at mail.dm.unipi.it]

Ciao,

Il Mer, 10 Giugno 2015 9:54 pm, Niels Möller ha scritto:
> sqrt_n (n, A):
>
>   Input: A = {a_{n-1}, ..., a_0}, n >= 1
>   Output: X = {x_{n-1}, ..., x_0} \approx sqrt(B^{n-1} A)
>
>   if (n == 1)
>     x_0 <-- floor (sqrt (a_0))
>     return X = {x_0}
>
>   k <-- floor (n/2)
>   H <-- sqrt_n(k+1, {a_{n-1},...,a_{n-1-k}}) // k+1 limbs
>
>   if (n odd)	// n = 2k+1
>     E <-- H^2 - A
>     X <-- H B^k - floor (B^k E / 2 H)
>
>   else		// n = 2k
[...]

> Note that all the computations of the error/residual E is subject to a
> lot of cancellation, so one can use mullo (for small sizes) or
> wraparound arithmetic using mulmod_bnm1 for larger sizes.

If we need both the sqrt of the high-half and the residual, we should
directly use sqrtrem, don't we? Once we compute a residual, correcting the
"approximation" to obtain an "exact" result does not cost too much, right?
But maybe the use of wraparound compensate the propagation of an updated
residual...

> And the divisions really ought to use divappr, rather than exact
> truncating division.

You are right, I did not try divappr yet. I'll do.

Best regards,
m

-- 
http://bodrato.it/