gcdext_1 questions
Paul Zimmermann
Paul.Zimmermann at loria.fr
Fri Dec 4 13:39:01 CET 2009
Niels,
> > Thus one iteration is equivalent to Stein's binary algorithm, except you
> > replace the subtraction by an addmul_si (and doing that you remove more bits).
>
> A slightly different way to look at it is as follows.
> We have u and v odd, and update
>
> u <-- (u + q v) / 2^{j+number of extra trailing zeros}
>
> where q = - u v^{-1} mod 2^j (represented in some symmetric or
> assymmetric range).
yes exactly, except q = - u v^{-1} mod 2^{j+1} in our algorithm.
> In the binary algorithm as described in Stehlé's and your paper, j is
> determined by the number of "extra trailing bits" from the previous
> iteration. But from a practical point of view, the important thing is
> that u and q v are of the same bitsize, so that u + q*v is not much
> larger than u, growing at most by a carry bit.
>
> So I think one could do a little better by using count_leading_zeros, to
> take into account precisely what happens to the bits also in the most
> significant end (instead of just relying on the growth there being
> bounded as a fibonacci sequence), we may by luck have a bit or two
> cancelling there too.
you mean you want to determine the best q modulo 2^{j+1}, i.e., the value
so that |u + q v| is minimal?
> Like,
>
> while (u != v)
> {
> if (u >v)
> {
> j = #u - #v + 1; // I let # denote bitsize
> q = -u v^{-1} mod 2^j;
> u += q v;
> u >>= (trailing zeros, >= j);
> }
> else { ... Analogous update v += q * v ... }
> }
>
> where there then are several options for how to handle signs of u, v and
> q.
I'm not sure to understand your code, since j is not defined that way in our
algorithm. Do you mean a variant of our algorithm?
Paul
More information about the gmp-devel
mailing list