Anomaly in mpn_sqrtrem and mpn_rootrem
bodrato at mail.dm.unipi.it
bodrato at mail.dm.unipi.it
Wed Jun 24 08:48:44 UTC 2015
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/
More information about the gmp-devel
mailing list