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