CACUL (RAPIDE) DE PGCD

//==========================================================================
//                                ENTIERS 64 BITS
//==========================================================================
typedef unsigned __int64 U64;
//==========================================================================
//                                      PGCD
//==========================================================================
U64 pgcdStupid(U64 a,U64 b)
{
while(1) // boucle infinie
  {
  if(b > a) // on veut toujours avoir <a> plus grand que <b>
    {
    U64 tmp;
    tmp = a;
    a   = b;
    b   = tmp;
    }

  if(0==b) break; // bingo ! PGCD(a,0) = a

  a -= b;
  }

return a;
} // pgcdStupid()
//--------------------------------------------------------------------------
U64 pgcdSchool(U64 a,U64 b)
{
if(b==0) return a; // cas trivial ... PGCD(a,0)=a

while(1) // boucle infinie
  {
  U64 r;

  r = a%b;

  if(0==r) break; // bingo ! PGCD(a=k.b,b)=b

  // PGCD(a,b)=PGCD(b,r=a%b)
  a = b;
  b = r;
  }

return b;
} // pgcdSchool()
//--------------------------------------------------------------------------
U64 pgcdFast(U64 a,U64 b)
{
U64 pow2;

// cas trivial ... PGCD(a,0)=a ou PGCD(0,b)=b
// dans les deux cas : PGCD(a,b)=a+b si l'un des deux nombres est nul
if(0==a || b==0) return a+b;

// on s'occupe de toutes les puissances de deux contenues dans le PGCD(a,b)
// PGCD( a=2^n.a' , b=2^n.b' )=2^n.PGCD(a',b')
pow2 = 1u;
while(!(a&1)) // tant que a est pair
  {
  if(!(b&1)) // si b est de plus pair
    {
    // PGCD( a=2.a' , b=2.b' )=2.PGCD(a',b')
    a     >>= 1;
    b     >>= 1;
    pow2  <<= 1;
    }
  else
    {
    // on rend <a> pair
    do a>>=1;while(!(a&1));
    break; // toutes les puissances de 2 du PGCD(a,b) ont ete traitees
    }
  }

// on rend <b> pair
while(!(b&1)) b>>=1;

// ici les deux nombres sont pairs, le temps de calcul precedent est negligeable
// devant ce qui suit...
// information : la difference de deux nombres impairs est pair
do
  {
  if(a==b)
    {
    break; // bingo ! PGCD(a,b)=a
    }
  else if(a > b)
    {
    a -= b; // PGCD(a,b)=PGCD(a-b,b)
    do a>>=1;while(!(a&1)); // on peut rendre <a> pair car on a plus de puissance de deux dans le PGCD
    }
  else // b > a
    {
    b -= a; // PGCD(a,b)=PGCD(a,b-a)
    do b>>=1;while(!(b&1)); // on peut rendre <b> pair car on a plus de puissance de deux dans le PGCD
    }
  }while(1);

return pow2*a;
} // pgcdFast()