// file kernel/n/c/gcd.c: greatest common divisor
/*-----------------------------------------------------------------------+
 |  Copyright 2005-2006, Michel Quercia (michel.quercia@prepas.org)      |
 |                                                                       |
 |  This file is part of Numerix. Numerix is free software; you can      |
 |  redistribute it and/or modify it under the terms of the GNU Lesser   |
 |  General Public License as published by the Free Software Foundation; |
 |  either version 2.1 of the License, or (at your option) any later     |
 |  version.                                                             |
 |                                                                       |
 |  The Numerix Library is distributed in the hope that it will be       |
 |  useful, but WITHOUT ANY WARRANTY; without even the implied warranty  |
 |  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU  |
 |  Lesser General Public License for more details.                      |
 |                                                                       |
 |  You should have received a copy of the GNU Lesser General Public     |
 |  License along with the GNU MP Library; see the file COPYING. If not, |
 |  write to the Free Software Foundation, Inc., 59 Temple Place -       |
 |  Suite 330, Boston, MA 02111-1307, USA.                               |
 +-----------------------------------------------------------------------+
 |                                                                       |
 |                                   PGCD                                |
 |                                                                       |
 +-----------------------------------------------------------------------*/

                 /* +--------------------------------------+
                    |  Oprations avec calcul de longueur  |
                    +--------------------------------------+ */

/* addition : c <- a+b, retourne lc */
extern inline long add(chiffre *a, long la, chiffre *b, long lb, chiffre *c) {
    chiffre r;
    long l;
    if (la > lb) {l = la; r = xn(add)(a,la,b,lb,c);}
    else         {l = lb; r = xn(add)(b,lb,a,la,c);}
    if (r) c[l++] = 1;
    return(l);
}

/* addition sur place : a <- a+b, retourne la */
extern inline long inc(chiffre *a, long la, chiffre *b, long lb) {
    chiffre r;
    long l;
    if (la > lb) {l = la; r = xn(inc)(a,la,b,lb);}
    else         {l = lb; r = xn(add)(b,lb,a,la,a);}
    if (r) a[l++] = 1;
    return(l);
}

/* soustraction sur place : a <- a-b, retourne la */
extern inline long dec(chiffre *a, long la, chiffre *b, long lb) {
    xn(dec)(a,la,b,(lb > la) ? la : lb);
    while ((la > 0) && (a[la-1] == 0)) la--;
    return(la);
}

/* addition avec dcalage : a <- a+b*BASE^n, retourne la */
extern inline long join(chiffre *a, long la, chiffre *b, long lb, long n) {
    chiffre r;
    long l;
         if (la > lb+n) {l = la;   r = xn(inc)(a+n,la-n,b,lb);}
    else if (la >= n)   {l = lb+n; r = xn(add)(b,lb,a+n,la-n,a+n);}
    else {l = lb+n; xn(clear)(a+la,n-la); xn(move)(b,lb,a+n); r = 0;}
    if (r) a[l++] = 1;
    return(l);
}

/* multiplication par un chiffre : c <- a*b, retourne lc */
extern inline long mul_1(chiffre *a, long la, unsigned long b, chiffre *c) {
    chiffre r;
    r = xn(mul_1)(a,la,b,c);
    if (r) c[la++] = r;
    return(la);
}

/* multiplication gnrale : c <- a*b, retourne lc */
extern inline long mul(chiffre *a, long la, chiffre *b, long lb, chiffre *c) {
    if ((la == 0) || (lb == 0)) return(0);
    else {
        if (la > lb) xn(fftmul)(a,la,b,lb,c);
        else         xn(fftmul)(b,lb,a,la,c);
        return(la+lb - (c[la+lb-1] == 0));
    }
}


                              /* +------------+
                                 |  Division  |
                                 +------------+ */

/* entre :
   a = naturel de longueur la
   b = naturel de longueur lb
   c = naturel de longueur lc
   r = 0 ou 1

   contrainte : lb >= 2, a >= b+r

   sortie :
   c <- floor(a/(b+r)), lc <- longueur(c)
   a <- a - b*c,        la <- longueur(a)
*/

static void divise
(chiffre *a, long *la, chiffre *b, long lb, chiffre *c, long *lc, chiffre r) {
    chiffre *aa, bb, *cc;
    long n;

    /* si r = 1, b <- b+1, traite  part le cas de dbordement */
    if ((r) && (xn(inc1)(b,lb))) {
        *lc = *la - lb;
        xn(move)(a+lb,*lc,c);
        for (*la = lb; (*la > 0) && (a[*la-1] == 0); (*la)--);
        *la = inc(a,*la,c,*lc);
        xn(dec1)(b,lb);
        return;
    }

    /* mmoire temporaire pour a et c de faon  viter des dbordements
       intermdiaires */
    cc = xn(alloc)(2*(*la)-lb+2);
    aa = cc + *la - lb + 1;

    /* dcalage  appliquer  b+r pour avoir le bit de poids fort  1 */
    for (bb=b[lb-1], n=0; (bb & (BASE_2)) == 0; bb <<=1, n++);
    if (n) {
        aa[(*la)++] = xn(shift_up)(a,*la,aa,n);
        xn(shift_up)(b,lb,b,n);
    } else {
        xn(move)(a,*la,aa);
        if (aa[*la-1] >= bb) aa[(*la)++] = 0;
    }
    *lc = *la - lb;

    xn(karpdiv)(aa,*lc,b,lb,cc,1);

    /* dcale puis recopie les rsultats */
    if (n) {
        xn(shift_down)(aa,lb,aa,n);
        xn(shift_down)(b,lb,b,n);
    }
    for (*la = lb; (*la) && (aa[*la-1] == 0); (*la)--); xn(move)(aa,*la,a);
    for (        ; (*lc) && (cc[*lc-1] == 0); (*lc)--); xn(move)(cc,*lc,c);

    /* b <- b-r, a <- a + r*q */
    if (r) {xn(dec1)(b,lb); *la = inc(a,*la,c,*lc);}

    xn(free)(cc);
}


                     /* +------------------------------+
                        |  Multiplication matricielle  |
                        +------------------------------+ */

/*
  entre :
   _v = tableau de pointeurs vers 6 naturels a,b,p,s,q,r
   _l = tableau de 6 longueurs la,lb,lp,lq,lr,ls
    n = longueur > 0

  contrainte :
    la >= n, lb >= n
    a[n..la-1] >= q, b[n..lb-1] >= r
    p*s - q*r = 1
    
  sortie :
    a <- s*a[0..n-1] - q*b[0..n-1] + a[n..la-1]*BASE^n
    b <- p*b[0..n-1] - r*a[0..n-1] + b[n..la-1]*BASE^n
    met  jour les longueurs
*/

static void mulinv(chiffre **_v, long *_l, long n) {
    chiffre *x,*y,*z, c;
    long l;

#define  a _v[0]
#define  b _v[1]
#define  p _v[2]
#define  s _v[3]
#define  q _v[4]
#define  r _v[5]
#define la _l[0]
#define lb _l[1]
#define lp _l[2]
#define ls _l[3]
#define lq _l[4]
#define lr _l[5]


    if ((lq) && (lr)) { /* q,r non nuls : combinaison gnrale */

        /* mmoire auxilliaire : 3n + max(lp,ls) + lq + lr + 1 */
        l = (lp > ls) ? lp : ls;
        x = xn(alloc)(3*n+l+lq+lr+1);
        y = x + n + l + 1;
        z = y + n + lq;

        /* x <- s*a[0..n-1], y <- q*b[0..n-1], z <- r*a[0..n-1] */
        mul(a,n,s,ls,x);
        mul(b,n,q,lq,y);
        mul(a,n,r,lr,z);

        /* a[0..n-1] <- x[0..n-1] - y[0..n-1], c <- retenue */
        c = xn(sub)(x,n,y,n,a);

        /* a[n..] <- a[n..] + x[n..] - y[n..] - c */
        xn(dec)(a+n,la-n,y+n,lq);
        if (ls <= la-n) {
            if (c) c = xn(dec1)(a+n,la-n);
            c = xn(inc)(a+n,la-n,x+n,ls) - c;
            if (c) a[la++] = 1;
        } else {
            if (c) c = xn(dec1)(x+n,ls);
            c = xn(inc)(x+n,ls,a+n,la-n) - c;
            la = ls+n;
            if (c) x[la++] = 1; else while ((la > n) && (x[la-1] == 0)) la--;
            xn(move)(x+n,la-n,a+n);
        }

        /* x <- p*b[0..n-1] */
        mul(b,n,p,lp,x);

        /* b[0..n-1] <- x[0..n-1] - z[0..n-1], c <- retenue */
        c = xn(sub)(x,n,z,n,b);

        /* b[n..] <- b[n..] + x[n..] - z[n..] - c */
        xn(dec)(b+n,lb-n,z+n,lr);
        if (lp <= lb-n) {
            if (c) c = xn(dec1)(b+n,lb-n);
            c = xn(inc)(b+n,lb-n,x+n,lp) - c;
            if (c) b[lb++] = 1;
        } else {
            if (c) c = xn(dec1)(x+n,lp);
            c = xn(inc)(x+n,lp,b+n,lb-n) - c;
            lb = lp+n;
            if (c) x[lb++] = 1; else while ((lb > n) && (x[lb-1] == 0)) lb--;
            xn(move)(x+n,lb-n,b+n);
        }

        xn(free)(x);

    } else if (lq) { /* r = 0 => p=s=1, a <- a - q*b[0..n-1] */

        x = xn(alloc)(n+lq);
        mul(b,n,q,lq,x);
        xn(dec)(a,la,x,lq+n);
        xn(free)(x);

    } else if (lr) { /* q = 0 => p=s=1, b <- b - r*a[0..n-1] */

        x = xn(alloc)(n+lr);
        mul(a,n,r,lr,x);
        xn(dec)(b,lb,x,lr+n);
        xn(free)(x);
    }

    /* mise  jour de la et lb (ncessaire mme si q = 0 ou r = 0 car il se
       peut qu'une partie basse soit initialement nulle) */
    while ((la > 0) && (a[la-1] == 0)) la--;
    while ((lb > 0) && (b[lb-1] == 0)) lb--;

#undef  a
#undef  b
#undef  p
#undef  s
#undef  q
#undef  r
#undef la
#undef lb
#undef lp
#undef ls
#undef lq
#undef lr
}

/*
  entre :
   _v0,_v1 = tableau de pointeurs vers 6 naturels a,b,p,s,q,r
   _l0,_l1 = tableau de 6 longueurs la,lb,lp,lq,lr,ls

  contraintes : p0*s0 - q0*r0 = p1*s1 - q1*r1 = 1

  sortie :
    p1 <- p0*p1 + q0*r1
    q1 <- p0*q1 + q0*s1
    r1 <- r0*p1 + s0*r1
    s1 <- r0*q1 + s0*s1
    met  jour les longueurs
*/

static void mulmat(chiffre **_v0, long *_l0, chiffre **_v1, long *_l1) {
#define  p0 _v0[2]
#define  q0 _v0[4]
#define  r0 _v0[5]
#define  s0 _v0[3]
#define  p1 _v1[2]
#define  q1 _v1[4]
#define  r1 _v1[5]
#define  s1 _v1[3]
#define lp0 _l0[2]
#define lq0 _l0[4]
#define lr0 _l0[5]
#define ls0 _l0[3]
#define lp1 _l1[2]
#define lq1 _l1[4]
#define lr1 _l1[5]
#define ls1 _l1[3]

    chiffre *x,*y,*z;
    long    lx,ly,lz;

    /* mmoire ncessaire : lp0 + max(lp1,lq1)          pour x
                            max(lq0,ls0) + max(lr1,ls1) pour y
                            lr0 + max(lp1,lq1)          pour z */
    lx = (lp1 > lq1) ? lp1 : lq1;
    lz = lx + lr0;
    lx += lp0;
    ly = ((lq0 > ls0) ? lq0 : ls0) + ((lr1 > ls1) ? lr1 : ls1);
    x = xn(alloc)(lx+ly+lz); y = x+lx; z = y+ly;

    lx  = mul( p0, lp0,  p1, lp1,  x  );
    ly  = mul( q0, lq0,  r1, lr1,  y  );
    lz  = mul( r0, lr0,  p1, lp1,  z  );
    lp1 = add( x,  lx,   y,  ly,   p1 );
    ly  = mul( s0, ls0,  r1, lr1,  y  );
    lr1 = add( z,  lz,   y,  ly,   r1 );
    lx  = mul( p0, lp0,  q1, lq1,  x  );
    ly  = mul( q0, lq0,  s1, ls1,  y  );
    lz  = mul( r0, lr0,  q1, lq1,  z  );
    lq1 = add( x,  lx,   y,  ly,   q1 );
    ly  = mul( s0, ls0,  s1, ls1,  y  );
    ls1 = add( z,  lz,   y,  ly,   s1 );


    xn(free)(x);

#undef  p0
#undef  q0
#undef  r0
#undef  s0
#undef  p1
#undef  q1
#undef  r1
#undef  s1
#undef lp0
#undef lq0
#undef lr0
#undef ls0
#undef lp1
#undef lq1
#undef lr1
#undef ls1
}

                       /* +------------------------+
                          |  Pgcd  deux chiffres  |
                          +------------------------+ */
        
/*
  entre :
    x = tableau de 8 chiffres [a0,a1,b0,b1,p,s,q,r]
    a = a0 + BASE*a1, b = b0 + BASE*b1

  contrainte : 0 < b < a < BASE*b

  Dveloppe en fraction continue la fraction a/b tant que les coefficients
  tiennent sur un chiffre

  sortie :
   [a0,a1,b0,b1] <- ind.
   [p,s,q,r]     <- coefficients des combinaisons effectues
*/
void xn(gcd_2)(chiffre *x)
#if defined(assembly_sn_gcd_2)
; 
#elif defined(use_clong)
{
    ndouble a,b,p,q,r,s,t;

    a = (ndouble)x[0] + ((ndouble)x[1] << HW);
    b = (ndouble)x[2] + ((ndouble)x[3] << HW);
    x[4] = x[5] = 1; x[6] = x[7] = 0;
    p = s = 1;       q = r = 0;
    do {
        t = a/b; a -= t*b; q += t*p; s += t*r;
        if ((q >= 2*BASE_2) || (s >= 2*BASE_2)) break;
        x[6] = q; x[5] = s;
        if (a == 0) break;
        t = b/a; b -= t*a; p += t*q; r += t*s;
        if ((p >= 2*BASE_2) || (r >= 2*BASE_2)) break;
        x[4] = p; x[7] = r;
        if (b == 0) break;
    } while(1);
}
#else /* mode dlong/slong : division par soustraction/dcalage */
{
#define a0 x[0]
#define a1 x[1]
#define b0 x[2]
#define b1 x[3]
#define p  x[4]
#define q  x[6]
#define r  x[7]
#define s  x[5]

  chiffre t;
  long i;

  /* [p,q,r,s] <- Id */
  p = s = 1; q = r = 0;

  while(1) {

    /* Ici a >= b. Dcale b,p,r tant que a >= 2*b */
    a1 = a1 - b1 - (a0 < b0);
    a0 -= b0;
    for (i=0; ;i++) {
      t = b1 + (a0 < b0); a0 -= b0;
      if (a1 < t) {a1 -= t; break;}
      a1 -= t;
      b1 = 2*b1 + (b0 >= (BASE_2));
      b0 <<= 1;
      if ((p & (BASE_2)) || (r & (BASE_2))) {p >>= i; r >>= i; return;}
      p <<= 1;
      r <<= 1;
    }

    /* Ici a0:a1 contient a - 2^(i+1)*b,
           b0:b1 contient 2^i*b
	   p     contient 2^i*p
	   r     contient 2^i*r
       Calcule a/b par soustractions et dcalages et sort si q ou s dborde */
    a0 += b0; a1 += b1 + (a0 < b0);
    q += p; s += r;
    if ((q < p) || (s < r)) {q -= p; s -= r; p >>= i; r >>= i; return;}
    while (i) {
      b0 >>= 1; if (b1&1) b0 += BASE_2;
      b1 >>= 1;
      p  >>= 1;
      r  >>= 1;
      i--;
      if ((a1 > b1) || ((a1 == b1) && (a0 >= b0))) {
	a1 -= b1 + (a0 < b0); a0 -= b0;
	q += p; s += r;
	if ((q < p) || (s < r)) {q -= p; s -= r; p >>= i; r >>= i; return;}
      }
    }

    /* Fin de la division de a par b.
       Ici a0:a1 contient le reste, b,p,r ont t restaurs et q,s mis  jour.
       si a = 0, c'est termin */
    if ((a0 == 0) && (a1 == 0)) return;

    /* Ici b >= a. Dcale a,s,q tant que b >= 2*a */
    b1 = b1 - a1 - (b0 < a0);
    b0 -= a0;
    for (i=0; ;i++) {
      t = a1 + (b0 < a0); b0 -= a0;
      if (b1 < t) {b1 -= t; break;}
      b1 -= t;
      a1 = 2*a1 + (a0 >= (BASE_2));
      a0 <<= 1;
      if ((s & (BASE_2)) || (q & (BASE_2))) {s >>= i; q >>= i; return;}
      s <<= 1;
      q <<= 1;
    }

    /* Ici b0:b1 contient b - 2^(i+1)*a,
           a0:a1 contient 2^i*a
	   s     contient 2^i*s
	   q     contient 2^i*q
       Calcule b/a par soustractions et dcalages et sort si p ou r dborde */
    b0 += a0; b1 += a1 + (b0 < a0);
    r += s; p += q;
    if ((r < s) || (p < q)) {r -= s; p -= q; s >>= i; q >>= i; return;}
    while(i) {
      a0 >>= 1; if (a1&1) a0 += BASE_2;
      a1 >>= 1;
      s  >>= 1;
      q  >>= 1;
      i--;
      if ((b1 > a1) || ((b1 == a1) && (b0 >= a0))) {
	b1 -= a1 + (b0 < a0); b0 -= a0;
	r += s; p += q;
	if ((r < s) || (p < q)) {r -= s; p -= q; s >>= i; q >>= i; return;}
      }
    }

    /* Fin de la division de b par a.
       Ici b0:b1 contient le reste, a,s,q ont t restaurs et r,p mis  jour.
       si b = 0, c'est termin */
    if ((b0 == 0) && (b1 == 0)) return;

  } /* while(1) */

#undef a0
#undef a1
#undef b0
#undef b1
#undef p
#undef q
#undef r
#undef s
}
#endif /* assembly_sn_gdd_2, use_clong */

                     /* +-----------------------------+
                        |  Demi-pgcd  deux chiffres  |
                        +-----------------------------+ */
 
/*
  entre :
    x = tableau de 8 chiffres [a0,a1,b0,b1,p,s,q,r]
    a = a0 + BASE*a1, b = b0 + BASE*b1

  contrainte : 0 < b < a

  Dveloppe en fraction continue les fractions a/(b+1) et (a+1)/b
  tant que les quotients concident et que les coefficients tiennent
  sur un chiffre

  sortie :
    [a0,a1,b0,b1] <- ind.
    [p,s,q,r]     <- coefficients des combinaisons effectues
*/
void xn(hgcd_2)(chiffre *x)
#if defined(assembly_sn_hgcd_2)
;
#elif defined(use_clong)
{
    ndouble a,b,p,q,r,s,t;

    a = (ndouble)x[0] + ((ndouble)x[1] << HW);
    b = (ndouble)x[2] + ((ndouble)x[3] << HW);
    x[4] = x[5] = 1; x[6] = x[7] = 0;
    p = s = 1;       q = r = 0;
    do {
        t = (a-q)/(b+p); if (t == 0) break;
        a -= t*b; q += t*p; s += t*r;
        if ((q >= 2*BASE_2) || (s >= 2*BASE_2)) break;
        x[6] = q; x[5] = s;
        t = (b-r)/(a+s); if (t == 0) break;
        b -= t*a; p += t*q; r += t*s;
        if ((p >= 2*BASE_2) || (r >= 2*BASE_2)) break;
        x[4] = p; x[7] = r;
    } while(1);
}
#else /* mode dlong/slong : division par soustraction/dcalage */
{
#define a0 x[0]
#define a1 x[1]
#define b0 x[2]
#define b1 x[3]
#define p  x[4]
#define q  x[6]
#define r  x[7]
#define s  x[5]

  chiffre t;
  long i;

  /* [p,q,r,s] <- Id */
  p = s = 1; q = r = 0;

  /* a <- a-q, b <- b+p */
  b0++;  b1 += (b0 == 0);

  while(1) {

    /* Ici a-q >= b+p. Dcale b,p,r tant que a-q >= 2*(b+p) */
    a1 = a1 - b1 - (a0 < b0);
    a0 -= b0;
    for (i=0; ;i++) {
      t = b1 + (a0 < b0); a0 -= b0;
      if (a1 < t) {a1 -= t; break;}
      a1 -= t;
      b1 = 2*b1 + (b0 >= (BASE_2));
      b0 <<= 1;
      if ((p & (BASE_2)) || (r & (BASE_2))) {p >>= i; r >>= i; return;}
      p <<= 1;
      r <<= 1;
    }

    /* Ici a0:a1 contient a-q - 2^(i+1)*(b+p),
           b0:b1 contient 2^i*(b+p)
	   p     contient 2^i*p
	   r     contient 2^i*r
       Calcule (a-q)/(b+p) par soustractions et dcalages et sort si q ou s dborde */
    a0 += b0; a1 += b1 + (a0 < b0);
    q += p; s += r;
    if ((q < p) || (s < r)) {q -= p; s -= r; p >>= i; r >>= i; return;}
    while (i) {
      b0 >>= 1; if (b1&1) b0 += BASE_2;
      b1 >>= 1;
      p  >>= 1;
      r  >>= 1;
      i--;
      if ((a1 > b1) || ((a1 == b1) && (a0 >= b0))) {
	a1 -= b1 + (a0 < b0); a0 -= b0;
	q += p; s += r;
	if ((q < p) || (s < r)) {q -= p; s -= r; p >>= i; r >>= i; return;}
      }
    }
    a0 += q; a1 += (a0 < q);
    b1 -= (b0 < p); b0 -= p;

    /* Fin de la division de a-q par b+p.
       Ici a0:a1 contient le reste, b,p,r ont t restaurs et q,s mis  jour.
       si a+s > b-r, c'est termin */
    a0 += s; a1 += (a0 < s);
    b1 -= (b0 < r); b0 -= r;
    if ((a1 > b1) || ((a1 == b1) && (a0 > b0))) return;

    /* Ici b-r >= a+s. Dcale a,s,q tant que b-r >= 2*(a+s) */
    b1 = b1 - a1 - (b0 < a0);
    b0 -= a0;
    for (i=0; ;i++) {
      t = a1 + (b0 < a0); b0 -= a0;
      if (b1 < t) {b1 -= t; break;}
      b1 -= t;
      a1 = 2*a1 + (a0 >= (BASE_2));
      a0 <<= 1;
      if ((s & (BASE_2)) || (q & (BASE_2))) {s >>= i; q >>= i; return;}
      s <<= 1;
      q <<= 1;
    }

    /* Ici b0:b1 contient b-r - 2^(i+1)*(a+s),
           a0:a1 contient 2^i*(a+s)
	   s     contient 2^i*s
	   q     contient 2^i*q
       Calcule (b-r)/(a+s) par soustractions et dcalages et sort si p ou r dborde */
    b0 += a0; b1 += a1 + (b0 < a0);
    r += s; p += q;
    if ((r < s) || (p < q)) {r -= s; p -= q; s >>= i; q >>= i; return;}
    while(i) {
      a0 >>= 1; if (a1&1) a0 += BASE_2;
      a1 >>= 1;
      s  >>= 1;
      q  >>= 1;
      i--;
      if ((b1 > a1) || ((b1 == a1) && (b0 >= a0))) {
	b1 -= a1 + (b0 < a0); b0 -= a0;
	r += s; p += q;
	if ((r < s) || (p < q)) {r -= s; p -= q; s >>= i; q >>= i; return;}
      }
    }
    b0 += r; b1 += (b0 < r);
    a1 -= (a0 < s); a0 -= s;

    /* Fin de la division de b-r par a+s.
       Ici b0:b1 contient le reste, a,s,q ont t restaurs et r,p mis  jour.
       si b+p > a-q, c'est termin */
    b0 += p; b1 += (b0 < p);
    a1 -= (a0 < q); a0 -= q;
    if ((b1 > a1) || ((b1 == a1) && (b0 > a0))) return;

  } /* while(1) */

#undef a0
#undef a1
#undef b0
#undef b1
#undef p
#undef q
#undef r
#undef s
}
#endif /* assembly_sn_hgcd_2, use_clong */

                         /* +--------------------+
                            |  Pgcd quadratique  |
                            +--------------------+ */

/*
  entre :
  _v = tableau de pointeurs vers 6 naturels a,b,p,s,q,r
  _l = tableau de 6 longueurs la,lb,lp,lq,lr,ls
  mode = 0,1,2

  contraintes :
    la > 0, lb > 0, a[la-1] > 0, b[lb-1] > 0
    capacit(p) >= la, capacit(q) >= la,
    capacit(r) >= lb, capacit(s) >= lb.

  dveloppe en fraction continue la fraction a/b (mode = 0 ou 1)
  ou les fractions a/(b+1) et (a+1)/b tant que les quotients sont
  gaux et la prcision des nombres rsiduels est suffisante (mode = 2).

  dans les trois modes, calcule des entiers naturels p',q',r',s',a',b' tels que
    a = p'*a' + q'*b'
    b = r'*a' + s'*b'
    p'*s' - q'*r' = 1

  avec a' = 0 ou b' = 0 (mode 0 ou 1)
       a' >= q' et b' >= r' (mode 2)

  sortie si mode = 0 :
    a,b <- a',b', met  jour les longueurs
    p,q,r,s <- ind.

  sortie si mode = 1 ou 2 :
    a,b,p,s,q,r <- a',b',p',s',q',r', met  jour les longueurs

  remarque :
    en mode 2 on progresse d'au moins une tape (= un chiffre dans la
    premire division) si et seulement si a <> b.

*/

#ifdef debug_gcd_n2
void xn(gcd_n2_buggy)
#else
void xn(gcd_n2)
#endif
(chiffre **_v, long *_l, int mode)
#ifdef assembly_gcd_n2
;
#else
{
    chiffre *a, *b, *p, *q, *r, *s, *u, *v, x[8];
    long    la, lb, lp, lq, lr, ls, lu, lv, n;
    long i,j,l,t;

#define a0 x[0]
#define a1 x[1]
#define b0 x[2]
#define b1 x[3]
#define p0 x[4]
#define q0 x[6]
#define r0 x[7]
#define s0 x[5]

    /* mmoire de travail : 2*max(la,lb)+2 chiffres */
    n = (_l[0] > _l[1]) ? _l[0] : _l[1];
    u = xn(alloc_tmp)(2*n+2); v = u+n+1;

    /* initialisation : [p,q,r,s] <- id */
    if (mode) {_v[2][0] = _v[3][0] = 1; _l[2] = _l[3] = 1; _l[4] = _l[5] = 0;}

    l = 0; /* rang du premier chiffre significatif */

    while (1) {

        /* force a >= b, en ignorant les chiffres de rang < l */
        i = xn(cmp)(_v[0]+l,_l[0]-l,_v[1]+l,_l[1]-l);
             if (i > 0)     i=0;    /* cas a > b, ok      */
        else if (i < 0)     i=1;    /* cas a < b, change */
        else if (mode == 2) break;  /* cas a = b, abandonne en mode 2 */

        a = _v[i]; b = _v[1-i]; la= _l[i]; lb= _l[1-i];
        t = la - lb;

        /* Extrait autant de bits de poids fort de a et b que possible.
	   Soient aa = a0 + BASE*a1 et bb = b0 + BASE*b1 les ttes de a et b.
           La premire division qui sera effectue est aa/bb si on est en
           mode exact, aa/(bb+1) sinon ; il faut garantir que le quotient
           obtenu est non nul et tient sur un chiffre.

           Cas a[la-1] > b[lb-1] :
              on prend autant de chiffres pour a et b

           Cas a[la-1] < b[lb-1] :
              alors t > 0, on prend un chiffre de plus pour a

           Cas a[la-1] = b[lb-1] et t = 0 :
              on prend autant de chiffres pour a et b
              il se peut qu'on obtienne aa = bb en mode inexact, prvoir
              de produire un quotient gal  1 dans ce cas.

           Cas a[la-1] = b[lb-1] et t > 0 :
              en mode exact avec lb=1 on prend un chiffre pour a
              dans les autres cas on prend un chiffre de plus pour a
        */
        if (lb-l == 1) {
	  a0 = a[la-1]; a1 = 0;
          b0 = b[lb-1]; b1 = 0;
	  if ((a0 < b0) || ((a0 == b0) && (t) && (mode == 2))) {
	    a1 = a0; a0 = a[la-2];
	    t--;
	  }
        }
        else {
	  a0 = a[la-2]; a1 = a[la-1];
          b0 = b[lb-2]; b1 = b[lb-1];
	  if ((a1 < b1) || ((a1 == b1) && ((a0 < b0) || ((a0 == b0) && (t))))) {
	    for (j=HW; (a1 & (BASE_2)) == 0; j--, a1 <<= 1);
	    if (j < HW) {
	      a1 += a0 >> j; a0 = (a0 << (HW-j)) + (a[la-3] >> j);
	      b0 = (b1 << (HW-j)) + (b0 >> j); b1 >>= j;
	    } else {
	      b0 = b1; b1 = 0;
	    }
	    t--;
	  }
	  else if (lb-l > 2) {
	    for (j=HW; (a1 & (BASE_2)) == 0; j--, a1 <<= 1);
	    if (j < HW) {
	      a1 += a0 >> j; a0 = (a0 << (HW-j)) + (a[la-3] >> j);
	      b1 = (b1 << (HW-j)) + (b0 >> j);
	      b0 = (b0 << (HW-j)) + (b[lb-3] >> j);
	    }
	  }
        }

        /* si on a charg tous les bits de a et b et mode < 2, effectue autant
           de divisions que possible tant que les coefficients des combinaisons
           tiennent dans un chiffre */
        if ((mode != 2) && (t==0) && (la <= 2)) xn(gcd_2)(x);

	/* si aa = bb, q0 <- 1 */
	else if ((a1 == b1) && (a0 == b0)) {q0 = 1; r0 = 0;}
	
        /* si t > 0, divise aa par bb+1 */
        else if (t) {
#ifdef ndouble
          ndouble aa = ((ndouble)a1 << HW) + (ndouble)a0;
          ndouble bb = ((ndouble)b1 << HW) + (ndouble)b0 + 1;
	  q0 = aa/bb;
#else
	  while (b1) {
	    b0 = (b1 << (HW-1)) + (b0 >> 1); b1 >>= 1;
	    a0 = (a1 << (HW-1)) + (a0 >> 1); a1 >>= 1;
	  }
	  b0++;
	  if (b0) xn(div_0)(a0,a1,b0,&q0,&r0); else q0 = a1;
#endif
	  r0 = 0;
	}

        /* dernier cas : a et b comportent des bits ngligs, t = 0 et aa > bb.
           Dcompose aa/(bb+1) et (aa+1)/bb en fraction continue tant que les
           coefficients tiennent dans un chiffre et les quotients sont gaux */
        else xn(hgcd_2)(x);

        /*
          mise  jour a,b,p,q,r,s.
          Lorsque r0 > 0, on a effectu plusieurs divisions donc t = 0.
          Lorsque r0 = 0, il n'y a eu qu'une division, ce qui simplifie les
          combinaisons  effectuer mais en revanche on doit tenir compte de t.
        */
        if (r0) { /* a <- s0*a - q0*b, b <- p0*b - r0*a */
            lu = mul_1(a,la,r0,u);   xn(mul_1)(a,la,s0,a);
            lv = mul_1(b,lb,q0,v);   xn(mul_1)(b,lb,p0,b);
            _l[i]   = la = dec(a,la,v,lv);
            _l[1-i] = lb = dec(b,lb,u,lu);

        } else { /* a <- a - q0*b*BASE^t */
            lu = mul_1(b,lb,q0,u);
            xn(dec)(a+t,la-t,u,lu);
            while ((la > 0) && (a[la-1] == 0)) la--; _l[i] = la;
        }

        if (mode) {
            p = _v[2+i]; s = _v[3-i]; q = _v[4+i]; r = _v[5-i]; 
            lp= _l[2+i]; ls= _l[3-i]; lq= _l[4+i]; lr= _l[5-i];

            if (r0) {

                if (lq) { /* p <- p0*p + r0*q, q <- q0*p + s0*q */
                    lu = mul_1(p,lp,q0,u);  lv = mul_1(q,lq,r0,v);
                    lp = mul_1(p,lp,p0,p);  lq = mul_1(q,lq,s0,q);
                    lp = inc  (p,lp,v,lv);  lq = inc  (q,lq,u,lu);
                } else if (q0) {p[0] = p0; q[0] = q0; lq = 1;}

                
                if (lr) { /* r <- p0*r + r0*s, s <- q0*r + s0*s */
                    lu = mul_1(r,lr,q0,u);  lv = mul_1(s,ls,r0,v);
                    lr = mul_1(r,lr,p0,r);  ls = mul_1(s,ls,s0,s);
                    lr = inc  (r,lr,v,lv);  ls = inc  (s,ls,u,lu);
                } else if (r0) {s[0] = s0; r[0] = r0; lr = 1;}

                _l[2+i] = lp; _l[4+i] = lq; _l[5-i] = lr; _l[3-i] = ls;

            } else { /* q += q0*p*BASE^t, s += q0*r*BASE^t */

                lu = mul_1(p,lp,q0,u);
                _l[4+i] = lq = join(q,lq,u,lu,t);

                if (lr) {
                    lu = mul_1(r,lr,q0,u);
                    _l[3-i] = ls = join(s,ls,u,lu,t); 
                }
            }
        
            /* cherche le rang du premier chiffre non affect par une retenue */
            if (mode == 2) {

                if (xn(cmp)(a,lq,q,lq) >= 0) l = lq;
                else for (l=lq+1; (l<la) && (a[l-1] == 0); l++);

                if (xn(cmp)(b,lr,r,lr) >= 0) j = lr;
                else for (j=lr+1; (j<lb) && (b[j-1] == 0); j++);
                if (l < j) l = j;

                if ((ls <= la) && (a[ls-1] < (BASE_2)+(BASE_2-1)-s[ls-1])) j = ls;
                else for (j=ls+1; (j<la) && (a[j-1] == (BASE_2)+(BASE_2-1)); j++);
                if (l < j) l = j;

                if ((lp <= lb) && (b[lp-1] < (BASE_2)+(BASE_2-1)-p[lp-1])) j = lp;
                else for (j=lp+1; (j<lp) && (p[j-1] == (BASE_2)+(BASE_2-1)); j++);
                if (l < j) l = j;

            }
        }

        if ((l >= la) || (l >= lb)) break;
    }

    /* termin */
    xn(free_tmp)(u);

#undef a0
#undef a1
#undef b0
#undef b1
#undef p0
#undef q0
#undef r0
#undef s0

}
#endif /* assembly_sn_gcd_n2 */

                    /* +---------------------------------+
                       |  Algorithme de Lehmer rcursif  |
                       +---------------------------------+ */


/* mmes conventions d'appel que gcd_n2 */

#ifdef debug_lehmer
void xn(lehmer_buggy)
#else
void xn(lehmer)
#endif
(chiffre **_v, long *_l, int mode) {
#define  a _v[0]
#define  b _v[1]
#define  p _v[2]
#define  s _v[3]
#define  q _v[4]
#define  r _v[5]
#define la _l[0]
#define lb _l[1]
#define lp _l[2]
#define ls _l[3]
#define lq _l[4]
#define lr _l[5]

    static long lehmer_tab[] = {lehmer_lim_0,lehmer_lim_1,lehmer_lim_1};
    chiffre *_v0[6];
    long j,n,_l0[6];

    /* petits arguments -> gcd_n2 */
    n = (la > lb) ? lb : la;
    if (n <= lehmer_tab[mode]) {xn(gcd_n2)(_v,_l,mode); return;}

    /* si |la-lb| >= n/4, divise max(a,b) par min(a,b) */
    if (la+lb >= (9*n)/4) {
        j = (la < lb);
        divise(_v[j],_l+j,_v[1-j],_l[1-j],_v[4+j],_l+4+j,(mode == 2));
        p[0] = s[0] = 1;
        lp = ls = 1; _l[5-j] = 0;
    }

    else {
        /* dcoupe le petit argument en deux et calcule le demi-pgcd associ */
        n >>= 1;
        a += n; b += n; la -= n; lb -= n;
        xn(lehmer)(_v,_l,2);
        a -= n; b -= n; la += n; lb += n;

        /*
          si l'appel rcursif retourne l'identit cela signifie que les parties
          hautes sont gales. Dans ce cas on compare les parties basses et
          on produit une division avec quotient = 1 si elles diffrent ou
          si mode != 2, on abandonne sinon.
        */
        if ((lq == 0) && (lr == 0)) switch(xn(cmp)(a,n,b,n)) {
            case  0: if (mode == 2) return;
            case  1: q[0] = 1; lq = 1; la = n; break;
            default: r[0] = 1; lr = 1; lb = n; break;
        }

        /* applique la transformation aux parties basses */
        mulinv(_v,_l,n);

    }

    /* cherche le rang du premier chiffre non affect par une retenue */
    if (mode < 2) n = 0;
    else {
        
        if (xn(cmp)(a,lq,q,lq) >= 0) n = lq;
        else for (n=lq+1; (n<la) && (a[n-1] == 0); n++);

        if (xn(cmp)(b,lr,r,lr) >= 0) j = lr;
        else for (j=lr+1; (j<lb) && (b[j-1] == 0); j++);
        if (n < j) n = j;
        
        if ((ls <= la) && (a[ls-1] < (BASE_2)+(BASE_2-1)-s[ls-1])) j = ls;
        else for (j=ls+1; (j<la) && (a[j-1] == (BASE_2)+(BASE_2-1)); j++);
        if (n < j) n = j;
        
        if ((lp <= lb) && (b[lp-1] < (BASE_2)+(BASE_2-1)-p[lp-1])) j = lp;
        else for (j=lp+1; (j<lp) && (p[j-1] == (BASE_2)+(BASE_2-1)); j++);
        if (n < j) n = j;
        
    }
    if ((la <= n) || (lb <= n)) return;

    /* mode = 0 -> termine rcursivement sans conserver p,q,r,s */
    if (mode == 0) {xn(lehmer)(_v,_l,0); return;}

    /* mode = 1 ou 2 -> recopie p,q,r,s */
    _v0[2] = xn(alloc)(lp+lq+lr+ls); xn(move)(p,lp,_v0[2]); _l0[2] = lp;
    _v0[3] = _v0[2] + lp;            xn(move)(s,ls,_v0[3]); _l0[3] = ls;
    _v0[4] = _v0[3] + ls;            xn(move)(q,lq,_v0[4]); _l0[4] = lq;
    _v0[5] = _v0[4] + lq;            xn(move)(r,lr,_v0[5]); _l0[5] = lr;

    /* 2me appel rcursif avec les parties hautes */
    a += n; b += n; la -= n; lb -= n;
    xn(lehmer)(_v,_l,mode);

    /* mise  jour des parties basses */
    if (mode == 2) {
        a -= n; b -= n; la += n; lb += n;
        mulinv(_v,_l,n);
    }

    /* multiplie les matrices obtenues */
    mulmat(_v0,_l0,_v,_l);
    xn(free)(_v0[2]);

#undef  a
#undef  b
#undef  p
#undef  s
#undef  q
#undef  r
#undef la
#undef lb
#undef lp
#undef ls
#undef lq
#undef lr
}


                              /* +------------+
                                 |  Contrle  |
                                 +------------+ */

#if defined (debug_gcd_n2) || defined(debug_lehmer)

/*
  entre :
  func = pointeur vers la fonction  controler
  name = nom en clair de la fonction
  _v,_l,_mode = paramtres pour f

  sortie :
  mme rsultat que f ou erreur fatale en cas de prcondition ou
  postcondition non remplie.

*/
typedef void gcd_func(chiffre **_v, long *_l, int mode);
static void check(gcd_func *func, char *name, chiffre **_v, long *_l, int mode) {
    long la0,lb0,lx,ly,i,n,ok;
    chiffre *a0,*b0,*x,*y,*z;
    char ss[80];

#define  a  _v[0]
#define  b  _v[1]
#define  p  _v[2]
#define  q  _v[4]
#define  r  _v[5]
#define  s  _v[3]
#define la  _l[0]
#define lb  _l[1]
#define lp  _l[2]
#define lq  _l[4]
#define lr  _l[5]
#define ls  _l[3]

    /* copie les arguments */
    la0 = la; lb0 = lb; n = (la > lb) ? la : lb;
    a0 = xn(alloc_tmp)(la+lb+4*n); b0 = a0 + la; x = b0 + lb; y = x + 2*n;
    xn(move)(a,la,a0);
    xn(move)(b,lb,b0);

    if (mode == 0) {

        /* calcule le pgcd en mode 1, le re calcule en mode 0 puis compare
           les rsultats */
        xn(move)(a,la,x); xn(move)(b,lb,y);
        z = a; a = x; x = z; z = b; b = y; y = z;
        func(_v,_l,1);
        z = a; a = x; x = z; z = b; b = y; y = z;
        lx = la; ly = lb; la = la0; lb = lb0;
        func(_v,_l,0);

        if (((xn(cmp)(a,la,x,lx)) || (xn(cmp)(b,lb,y,ly))) &&
            ((xn(cmp)(a,la,y,ly)) || (xn(cmp)(b,lb,x,lx)))) {
            sprintf(ss,"error in %s",name);
            xn(internal_error)(ss, 4, a0,la0,b0,lb0,a,la,b,lb);
        }
        xn(free_tmp)(a0);
        return;
    }

    /* calcule le pgcd douteux */
    func(_v,_l,mode);

    do {

        /*  priori, c'est bon */
        ok = 1;

        /* vrifie que max(la,lb,lp,lq,lr,ls) <= n */
        for (i=0; (i < 6) && (_l[i] <= n); i++);
        if (i < 6) {ok = 0; break;}

        /* vrifie que a = 0 ou b = 0 (mode 1) */
        if ((mode != 2) && (la) && (lb)) {ok = 0; break;}

        /* vrifie que a >= q et b >= r (mode 2) */
        if ((mode == 2) &&
            ((xn(cmp)(a,la,q,lq) < 0) || (xn(cmp)(b,lb,r,lr) < 0)))
        {ok = 0; break;}

        /* vrifie que ps-qr = 1 */
        lx = mul(p,lp,s,ls,x);
        ly = mul(q,lq,r,lr,y);
        if ((lx < ly) || (xn(dec)(x,lx,y,ly)))    {ok = 0; break;}
        y[0] = 1; ly = 1; if (xn(cmp)(x,lx,y,ly)) {ok = 0; break;}

        /* vrifie que [p,q,r,s]*[a,b] = [a0,b0] */
        lx = mul(p,lp,a,la,x);
        ly = mul(q,lq,b,lb,y);
        lx = inc(x,lx,y,ly);
        if (xn(cmp)(x,lx,a0,la0)) {ok = 0; break;}

        lx = mul(r,lr,a,la,x);
        ly = mul(s,ls,b,lb,y);
        lx = inc(x,lx,y,ly);
        if (xn(cmp)(x,lx,b0,lb0)) {ok = 0; break;}

    } while(0);

    /* erreur ? */
    if (!ok) {
        sprintf(ss,"error in %s",name);
        xn(internal_error)(ss, 8, a0,la0,b0,lb0,
                                  a, la, b, lb,
                                  p, lp, q, lq,
                                  r, lr ,s, ls);
    }

    xn(free_tmp)(a0);

#undef  a
#undef  b
#undef  p
#undef  s
#undef  q
#undef  r
#undef la
#undef lb
#undef lp
#undef ls
#undef lq
#undef lr
}

#ifdef debug_gcd_n2
void xn(gcd_n2)(chiffre **_v, long *_l, int mode) {
    check(xn(gcd_n2_buggy),"gcd_n2",_v,_l,mode);
}
#endif

#ifdef debug_lehmer
void xn(lehmer)(chiffre **_v, long *_l, int mode) {
    check(xn(lehmer_buggy),"lehmer",_v,_l,mode);
}
#endif

#endif /* defined (debug_gcd_n2) || defined(debug_lehmer) */