[BACK]Return to Ngcd.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / engine

Diff for /OpenXM_contrib2/asir2000/engine/Ngcd.c between version 1.3 and 1.4

version 1.3, 2001/10/09 01:36:11 version 1.4, 2018/03/29 01:32:51
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/Ngcd.c,v 1.2 2000/12/21 02:56:31 murao Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/Ngcd.c,v 1.3 2001/10/09 01:36:11 noro Exp $
 */  */
 /*  /*
 #include "ca.h"  #include "ca.h"
Line 22 
Line 22 
   
   
 int igcd_algorithm = 0;  int igcd_algorithm = 0;
     /*  == 0 : Euclid,      /*  == 0 : Euclid,
      *  == 1 : binary,       *  == 1 : binary,
      *  == 2 : bmod,       *  == 2 : bmod,
      *  >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,       *  >= 3 : (Weber's accelerated)/(Jebelean's generalized binary) algorithm,
      */       */
 int igcd_thre_inidiv = 50;  int igcd_thre_inidiv = 50;
     /*      /*
Line 49  int igcdacc_thre = 10;
Line 49  int igcdacc_thre = 10;
 #define ShouldCompRemInit(n1,n2) (igcd_thre_inidiv != 0 && PL(n1) >= igcd_thre_inidiv*PL(n2))  #define ShouldCompRemInit(n1,n2) (igcd_thre_inidiv != 0 && PL(n1) >= igcd_thre_inidiv*PL(n2))
   
 #define IniDiv(n1,n2) \  #define IniDiv(n1,n2) \
         if ( ShouldCompRemInit(n1,n2) ) {\    if ( ShouldCompRemInit(n1,n2) ) {\
                 N q, r; int w, b; \      N q, r; int w, b; \
                 divn(n1,n2,&q,&r); \      divn(n1,n2,&q,&r); \
                 if ( !r ) return(n2); \      if ( !r ) return(n2); \
                 b = trailingzerosn( r, &w ); \      b = trailingzerosn( r, &w ); \
                 q = n1;  n1 = n2;  n2 = q; \      q = n1;  n1 = n2;  n2 = q; \
                 rshiftn( r, w, b, n2 ); \      rshiftn( r, w, b, n2 ); \
         }    }
   
 /*  /*
  *      Binary GCD algorithm by J.Stein   *  Binary GCD algorithm by J.Stein
  *      [J. Comp. Phys. Vol. 1 (1967), pp. 397-405)]:   *  [J. Comp. Phys. Vol. 1 (1967), pp. 397-405)]:
  *      The right-shift binary algorithm is used.   *  The right-shift binary algorithm is used.
  */   */
   
   
 /*  /*
  *      subsidiary routines for gcdbinn below.   *  subsidiary routines for gcdbinn below.
  */   */
 static int /* number of bits */ trailingzeros_nbd( /* BD of N */ nbd, pnw )  static int /* number of bits */ trailingzeros_nbd( /* BD of N */ nbd, pnw )
 int *nbd, *pnw /* number of zero words */;  int *nbd, *pnw /* number of zero words */;
 {  {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         unsigned    unsigned
 #endif  #endif
         int nw, nb, w;    int nw, nb, w;
   
         for ( nw = 0; (w = *nbd) == 0; nbd++ ) nw++;    for ( nw = 0; (w = *nbd) == 0; nbd++ ) nw++;
         TRAILINGZEROS(w,nb);    TRAILINGZEROS(w,nb);
         *pnw = nw;    *pnw = nw;
         return nb;    return nb;
 }  }
   
 #define trailingzerosn(n,pnw) trailingzeros_nbd(BD(n),pnw)  #define trailingzerosn(n,pnw) trailingzeros_nbd(BD(n),pnw)
   
 static int /* PL of N */ rshift_nbd( /* BD of N */ nbd, /* PL of N */ nl,  static int /* PL of N */ rshift_nbd( /* BD of N */ nbd, /* PL of N */ nl,
                        /* # words */ shw, /* # bits */ shb, /* BD of N */ p )             /* # words */ shw, /* # bits */ shb, /* BD of N */ p )
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
 unsigned  unsigned
 #endif  #endif
 int *nbd, nl, shw, shb, *p;  int *nbd, nl, shw, shb, *p;
 {  {
         unsigned int i, v, w, lshb;                             /* <---- */    unsigned int i, v, w, lshb;        /* <---- */
   
         nbd += shw,  i = (nl -= shw);    nbd += shw,  i = (nl -= shw);
         if ( shb == 0 ) {    if ( shb == 0 ) {
                 for ( ; nl > 0; nl-- ) *p++ = *nbd++;      for ( ; nl > 0; nl-- ) *p++ = *nbd++;
                 return i;      return i;
         } else if ( nl < 2 ) {    } else if ( nl < 2 ) {
                 *p = (*nbd) >> shb;      *p = (*nbd) >> shb;
                 return 1;      return 1;
         }    }
         for ( lshb = BSH - shb, v = *nbd++; --nl > 0; v = w ) {    for ( lshb = BSH - shb, v = *nbd++; --nl > 0; v = w ) {
                 w = *nbd++;      w = *nbd++;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 *p++ = (v >> shb) | (w << lshb);                /********/      *p++ = (v >> shb) | (w << lshb);    /********/
 #else  #else
                 *p++ = (v >> shb) | ((w << lshb)&BMASK);      *p++ = (v >> shb) | ((w << lshb)&BMASK);
 #endif  #endif
         }    }
         if ( (v >>= shb) == 0 ) return( i-1 );    if ( (v >>= shb) == 0 ) return( i-1 );
         *p = v;    *p = v;
         return i;    return i;
 }  }
   
 #define rshiftn(ns,shw,shb,nd) (PL(nd)=rshift_nbd(BD(ns),PL(ns),shw,shb,BD(nd)))  #define rshiftn(ns,shw,shb,nd) (PL(nd)=rshift_nbd(BD(ns),PL(ns),shw,shb,BD(nd)))
Line 120  int *nbd, nl, shw, shb, *p;
Line 120  int *nbd, nl, shw, shb, *p;
 #ifdef FULLSET  #ifdef FULLSET
 static N N_of_i_lshifted_by_wb( i, gw, gb )  static N N_of_i_lshifted_by_wb( i, gw, gb )
 int i, gw, gb;  int i, gw, gb;
         /*    /*
          *      returns pointer to a new struct (N)(((int)i) >> (gb + gw*BSH))     *  returns pointer to a new struct (N)(((int)i) >> (gb + gw*BSH))
          */     */
 {  {
         unsigned int j, l, *p;                                  /* <---- */    unsigned int j, l, *p;          /* <---- */
         N n;    N n;
   
         j = i >> (BSH - gb);    j = i >> (BSH - gb);
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         i = (i << gb);                                          /********/    i = (i << gb);            /********/
 #else  #else
         i = (i << gb)&BMASK;    i = (i << gb)&BMASK;
 #endif  #endif
         l = j != 0 ? gw + 2 : gw + 1;    l = j != 0 ? gw + 2 : gw + 1;
         n = NALLOC(l);    n = NALLOC(l);
         PL(n) = l;    PL(n) = l;
         for ( p = BD(n); gw-- > 0; ) *p++ = 0;    for ( p = BD(n); gw-- > 0; ) *p++ = 0;
         *p++ = i;    *p++ = i;
         if ( j != 0 ) *p = j;    if ( j != 0 ) *p = j;
         return n;    return n;
 }  }
 #endif /* FULLSET */  #endif /* FULLSET */
   
 /*  /*
  *      routines to make a new struct   *  routines to make a new struct
  *              (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH))   *    (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH))
  */   */
 static N N_of_nbd_lshifted_by_wb( /* BD of N */ b, /* PL of N */ lb, gw, gb )  static N N_of_nbd_lshifted_by_wb( /* BD of N */ b, /* PL of N */ lb, gw, gb )
 int *b, lb, gw, gb;  int *b, lb, gw, gb;
         /*    /*
          *      returns pointer to a new struct     *  returns pointer to a new struct
          *              (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH))     *    (N)(((BD of N)(b[0],...,b[lb-1])) << (gb + gw*BSH))
          */     */
 {  {
         unsigned int rsh, s, t, *p, l;                          /* <---- */    unsigned int rsh, s, t, *p, l;        /* <---- */
         N n;    N n;
   
         l = lb + gw;    l = lb + gw;
         if ( gb == 0 ) {    if ( gb == 0 ) {
                 n = NALLOC(l);      n = NALLOC(l);
                 PL(n) = l;      PL(n) = l;
                 for ( p = BD(n); gw-- > 0; ) *p++ = 0;      for ( p = BD(n); gw-- > 0; ) *p++ = 0;
                 while ( lb-- > 0 ) *p++ = *b++;      while ( lb-- > 0 ) *p++ = *b++;
                 return n;      return n;
         }    }
         rsh = BSH - gb;  s = b[lb-1];    rsh = BSH - gb;  s = b[lb-1];
         if ( (t = s >> rsh) != 0 ) {    if ( (t = s >> rsh) != 0 ) {
                 n = NALLOC(l+1);      n = NALLOC(l+1);
                 PL(n) = l+1;      PL(n) = l+1;
                 (p = BD(n))[l] = t;      (p = BD(n))[l] = t;
         } else {    } else {
                 n = NALLOC(l);      n = NALLOC(l);
                 PL(n) = l;      PL(n) = l;
                 p = BD(n);      p = BD(n);
         }    }
         while ( gw-- > 0 ) *p++ = 0;    while ( gw-- > 0 ) *p++ = 0;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         *p++ = (t = *b++) << gb;                                /********/    *p++ = (t = *b++) << gb;        /********/
         for ( ; --lb > 0; t = s )    for ( ; --lb > 0; t = s )
                 *p++ = (t >> rsh) | ((s = *b++) << gb);         /********/      *p++ = (t >> rsh) | ((s = *b++) << gb);    /********/
 #else  #else
         *p++ = ((t = *b++) << gb)&BMASK;    *p++ = ((t = *b++) << gb)&BMASK;
         for ( ; --lb > 0; t = s )    for ( ; --lb > 0; t = s )
                 *p++ = (t >> rsh) | (((s = *b++) << gb)&BMASK);      *p++ = (t >> rsh) | (((s = *b++) << gb)&BMASK);
 #endif  #endif
         return n;    return n;
 }  }
   
 #define N_of_n_lshifted_by_wb(a,gw,gb) N_of_nbd_lshifted_by_wb(BD(a),PL(a),gw,gb)  #define N_of_n_lshifted_by_wb(a,gw,gb) N_of_nbd_lshifted_by_wb(BD(a),PL(a),gw,gb)
Line 196  int *b, lb, gw, gb;
Line 196  int *b, lb, gw, gb;
   
 #ifdef CALL  #ifdef CALL
 static int bw_int32( n )  static int bw_int32( n )
 unsigned int n;                                                 /* <---- */  unsigned int n;              /* <---- */
 {  {
         int w;    int w;
   
         w = 0;    w = 0;
 #if BSH > 32  #if BSH > 32
         if ( n > 0xffffffff ) w += 32, n >>= 32;    if ( n > 0xffffffff ) w += 32, n >>= 32;
 #endif  #endif
         if ( n >= 0x10000 ) w += 16, n >>= 16;    if ( n >= 0x10000 ) w += 16, n >>= 16;
         if ( n >=   0x100 ) w +=  8, n >>=  8;    if ( n >=   0x100 ) w +=  8, n >>=  8;
         if ( n >=    0x10 ) w +=  4, n >>=  4;    if ( n >=    0x10 ) w +=  4, n >>=  4;
         if ( n >=     0x4 ) w +=  2, n >>=  2;    if ( n >=     0x4 ) w +=  2, n >>=  2;
         if ( n >=     0x2 ) w +=  1, n >>=  1;    if ( n >=     0x2 ) w +=  1, n >>=  1;
         if ( n != 0 ) ++w;    if ( n != 0 ) ++w;
         return w;    return w;
 }  }
 #define BitWidth(n,bw) bw = bw_int32( n )  #define BitWidth(n,bw) bw = bw_int32( n )
 #else  #else
   
 #if BSH > 32  #if BSH > 32
 #define BitWidth(n,bw) {\  #define BitWidth(n,bw) {\
         unsigned int k = (n); \    unsigned int k = (n); \
         bw = 0; \    bw = 0; \
         if ( k > 0xffffffff ) bw += 32, k >>= 32; \    if ( k > 0xffffffff ) bw += 32, k >>= 32; \
         if ( k >= 0x10000 ) bw += 16, k >>= 16; \    if ( k >= 0x10000 ) bw += 16, k >>= 16; \
         if ( k >=   0x100 ) bw +=  8, k >>=  8; \    if ( k >=   0x100 ) bw +=  8, k >>=  8; \
         if ( k >=    0x10 ) bw +=  4, k >>=  4; \    if ( k >=    0x10 ) bw +=  4, k >>=  4; \
         if ( k >=     0x4 ) bw +=  2, k >>=  2; \    if ( k >=     0x4 ) bw +=  2, k >>=  2; \
         if ( k >=     0x2 ) bw +=  1, k >>=  1; \    if ( k >=     0x2 ) bw +=  1, k >>=  1; \
         if ( k != 0 ) bw++; \    if ( k != 0 ) bw++; \
 }  }
 #else  #else
 #define BitWidth(n,bw) {\  #define BitWidth(n,bw) {\
         unsigned int k = (n); \    unsigned int k = (n); \
         bw = 0; \    bw = 0; \
         if ( k >= 0x10000 ) bw += 16, k >>= 16; \    if ( k >= 0x10000 ) bw += 16, k >>= 16; \
         if ( k >=   0x100 ) bw +=  8, k >>=  8; \    if ( k >=   0x100 ) bw +=  8, k >>=  8; \
         if ( k >=    0x10 ) bw +=  4, k >>=  4; \    if ( k >=    0x10 ) bw +=  4, k >>=  4; \
         if ( k >=     0x4 ) bw +=  2, k >>=  2; \    if ( k >=     0x4 ) bw +=  2, k >>=  2; \
         if ( k >=     0x2 ) bw +=  1, k >>=  1; \    if ( k >=     0x2 ) bw +=  1, k >>=  1; \
         if ( k != 0 ) bw++; \    if ( k != 0 ) bw++; \
 }  }
 #endif  #endif
 #endif  #endif
Line 244  unsigned int n;       /* <---- */
Line 244  unsigned int n;       /* <---- */
 #include "igcdhack.c"  #include "igcdhack.c"
   
 /*  /*
  *      Implementation of the binary GCD algorithm for two oN structs   *  Implementation of the binary GCD algorithm for two oN structs
  *      (big-integers) in risa.   *  (big-integers) in risa.
  *   *
  *      The major operations in the following algorithms are the binary-shifts   *  The major operations in the following algorithms are the binary-shifts
  *      and the updates of (u, v) by (min(u,v), |u-v|), and are to be open-coded   *  and the updates of (u, v) by (min(u,v), |u-v|), and are to be open-coded
  *      without using routines for oN structures just as in addn() or subn().   *  without using routines for oN structures just as in addn() or subn().
  */   */
   
 static int igcd_binary_2w( u, lu, v, lv, pans )  static int igcd_binary_2w( u, lu, v, lv, pans )
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
 unsigned  unsigned
 #endif  #endif
 int *u, lu, *v, lv, *pans;                              /* <---- */  int *u, lu, *v, lv, *pans;        /* <---- */
         /*  both u[0:lu-1] and v[0:lv-1] are assumed to be odd */    /*  both u[0:lu-1] and v[0:lv-1] are assumed to be odd */
 {  {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         unsigned    unsigned
 #endif  #endif
         int i, h1, l1, h2, l2;                                  /* <---- */    int i, h1, l1, h2, l2;          /* <---- */
   
         l1 = u[0],  l2 = v[0];    l1 = u[0],  l2 = v[0];
         h1 = lu <= 1 ? 0 : u[1];    h1 = lu <= 1 ? 0 : u[1];
         h2 = lv <= 1 ? 0 : v[1];    h2 = lv <= 1 ? 0 : v[1];
         /**/    /**/
 loop:   if ( h1 == 0 ) {  loop:  if ( h1 == 0 ) {
         no_hi1: if ( h2 == 0 ) goto one_word;    no_hi1:  if ( h2 == 0 ) goto one_word;
         no_hi1n:if ( l1 == 1 ) return 0;    no_hi1n:if ( l1 == 1 ) return 0;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 if ( l2 == l1 ) {        if ( l2 == l1 ) {
                         for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ;        for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ;
                         goto one_word;        goto one_word;
                 } else if ( l2 < l1 ) {      } else if ( l2 < l1 ) {
                         l2 -= l1,  h2--;        l2 -= l1,  h2--;
                 } else l2 -= l1;      } else l2 -= l1;
                 i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );      i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
                 l2 |= (h2 << (BSH - i));      l2 |= (h2 << (BSH - i));
 #else  #else
                 if ( (l2 -= l1) == 0 ) {      if ( (l2 -= l1) == 0 ) {
                         for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ;        for ( l2 = h2; (l2&1) == 0; l2 >>= 1 ) ;
                         goto one_word;        goto one_word;
                 } else if ( l2 < 0 ) h2--, l2 += BASE;      } else if ( l2 < 0 ) h2--, l2 += BASE;
                 i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );      i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
                 l2 |= ((h2 << (BSH - i)) & BMASK);      l2 |= ((h2 << (BSH - i)) & BMASK);
 #endif  #endif
                 h2 >>= i;      h2 >>= i;
                 goto no_hi1;      goto no_hi1;
         } else if ( h2 == 0 ) {    } else if ( h2 == 0 ) {
         no_hi2: if ( l2 == 1 ) return 0;    no_hi2:  if ( l2 == 1 ) return 0;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 if ( l1 == l2 ) {      if ( l1 == l2 ) {
                         for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ;        for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ;
                         goto one_word;        goto one_word;
                 } else if ( l1 < l2 ) {      } else if ( l1 < l2 ) {
                         l1 -= l2, h1--;        l1 -= l2, h1--;
                 } else l1 -= l2;      } else l1 -= l2;
                 i = 0;  do { l1 >>= 1, i++; } while ( (l1&1) == 0 );      i = 0;  do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
                 l1 |= (h1 << (BSH - i));      l1 |= (h1 << (BSH - i));
 #else  #else
                 if ( (l1 -= l2) == 0 ) {      if ( (l1 -= l2) == 0 ) {
                         for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ;        for ( l1 = h1; (l1&1) == 0; l1 >>= 1 ) ;
                         goto one_word;        goto one_word;
                 } else if ( l1 < 0 ) h1--, l1 += BASE;      } else if ( l1 < 0 ) h1--, l1 += BASE;
                 i = 0;  do { l1 >>= 1, i++; } while ( (l1&1) == 0 );      i = 0;  do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
                 l1 |= ((h1 << (BSH - i)) & BMASK);      l1 |= ((h1 << (BSH - i)) & BMASK);
 #endif  #endif
                 if ( (h1 >>= i) == 0 ) goto one_word;      if ( (h1 >>= i) == 0 ) goto one_word;
                 goto no_hi2;      goto no_hi2;
         } else if ( l1 == l2 ) {    } else if ( l1 == l2 ) {
                 if ( h1 == h2 ) {      if ( h1 == h2 ) {
                         pans[0] = l1, pans[1] = h1;        pans[0] = l1, pans[1] = h1;
                         return 2;        return 2;
                 } else if ( h1 > h2 ) {      } else if ( h1 > h2 ) {
                         for ( l1 = h1 - h2; (l1&1) == 0; l1 >>= 1 ) ;        for ( l1 = h1 - h2; (l1&1) == 0; l1 >>= 1 ) ;
                         goto no_hi1n;        goto no_hi1n;
                 } else {      } else {
                         for ( l2 = h2 - h1; (l2&1) == 0; l2 >>= 1 ) ;        for ( l2 = h2 - h1; (l2&1) == 0; l2 >>= 1 ) ;
                         goto no_hi2;        goto no_hi2;
                 }      }
         } else if ( h1 == h2 ) {    } else if ( h1 == h2 ) {
                 if ( l1 > l2 ) {      if ( l1 > l2 ) {
                         for ( l1 -= l2; (l1&1) == 0; l1 >>= 1 ) ;        for ( l1 -= l2; (l1&1) == 0; l1 >>= 1 ) ;
                         goto no_hi1n;        goto no_hi1n;
                 } else {      } else {
                         for ( l2 -= l1; (l2&1) == 0; l2 >>= 1 ) ;        for ( l2 -= l1; (l2&1) == 0; l2 >>= 1 ) ;
                         goto no_hi2;        goto no_hi2;
                 }      }
         } else if ( h1 > h2 ) {    } else if ( h1 > h2 ) {
                 h1 -= h2;      h1 -= h2;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 if ( l1 < l2 ) l1 -= l2, h1--;      if ( l1 < l2 ) l1 -= l2, h1--;
                 else l1 -= l2;      else l1 -= l2;
                 i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );      i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
                 l1 |= (h1 << (BSH - i));      l1 |= (h1 << (BSH - i));
 #else  #else
                 if ( (l1 -= l2) < 0 ) h1--, l1 += BASE;      if ( (l1 -= l2) < 0 ) h1--, l1 += BASE;
                 i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );      i = 0; do { l1 >>= 1, i++; } while ( (l1&1) == 0 );
                 l1 |= ((h1 << (BSH - i)) & BMASK);      l1 |= ((h1 << (BSH - i)) & BMASK);
 #endif  #endif
                 h1 >>= i;      h1 >>= i;
         } else {    } else {
                 h2 -= h1;      h2 -= h1;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 if ( l2 < l1 ) l2 -= l1, h2--;      if ( l2 < l1 ) l2 -= l1, h2--;
                 else l2 -= l1;      else l2 -= l1;
                 i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );      i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
                 l2 |= (h2 << (BSH - i));      l2 |= (h2 << (BSH - i));
 #else  #else
                 if ( (l2 -= l1) < 0 ) h2--, l2 += BASE;      if ( (l2 -= l1) < 0 ) h2--, l2 += BASE;
                 i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );      i = 0;  do { l2 >>= 1, i++; } while ( (l2&1) == 0 );
                 l2 |= ((h2 << (BSH - i)) & BMASK);      l2 |= ((h2 << (BSH - i)) & BMASK);
 #endif  #endif
                 h2 >>= i;      h2 >>= i;
         }    }
         goto loop;    goto loop;
 one_word:  one_word:
         if ( l1 == 1 || l2 == 1 ) return 0;    if ( l1 == 1 || l2 == 1 ) return 0;
         else if ( l1 == l2 ) {    else if ( l1 == l2 ) {
                 pans[0] = l1;      pans[0] = l1;
                 return 1;      return 1;
         }    }
 one_word_neq:  one_word_neq:
         if ( l1 > l2 ) {    if ( l1 > l2 ) {
                 l1 -= l2;      l1 -= l2;
                 do { l1 >>= 1; } while ( (l1&1) == 0 );      do { l1 >>= 1; } while ( (l1&1) == 0 );
                 goto one_word;      goto one_word;
         } else {    } else {
                 l2 -= l1;      l2 -= l1;
                 do { l2 >>= 1; } while ( (l2&1) == 0 );      do { l2 >>= 1; } while ( (l2&1) == 0 );
                 goto one_word;      goto one_word;
         }    }
 }  }
   
 static N igcd_binary( n1, n2, nt )  static N igcd_binary( n1, n2, nt )
 N n1, n2, nt;  N n1, n2, nt;
         /*  both n1 and n2 are assumed to be odd */    /*  both n1 and n2 are assumed to be odd */
 {  {
         int l1, *b1, l2, *b2, *bt = BD(nt);    int l1, *b1, l2, *b2, *bt = BD(nt);
         int l;    int l;
   
         if ( (l = cmpn( n1, n2 )) == 0 ) return n1;    if ( (l = cmpn( n1, n2 )) == 0 ) return n1;
         else if ( l < 0 ) { SWAP( n1, n2, N ); }    else if ( l < 0 ) { SWAP( n1, n2, N ); }
         IniDiv( n1, n2 );    IniDiv( n1, n2 );
         if ( UNIN(n2) ) return 0;    if ( UNIN(n2) ) return 0;
         l1 = PL(n1), b1 = BD(n1),  l2 = PL(n2), b2 = BD(n2);    l1 = PL(n1), b1 = BD(n1),  l2 = PL(n2), b2 = BD(n2);
 loop:  loop:
 #if 0000  #if 0000
 {  {
Line 400  loop:
Line 400  loop:
   printf( "\n" );    printf( "\n" );
 }  }
 #endif  #endif
         if ( l1 <= 2 && l2 <= 2 ) {    if ( l1 <= 2 && l2 <= 2 ) {
                 l = igcd_binary_2w( b1, l1, b2, l2, bt );      l = igcd_binary_2w( b1, l1, b2, l2, bt );
                 if ( l == 0 ) return 0;      if ( l == 0 ) return 0;
                 PL(nt) = l;      PL(nt) = l;
                 return nt;      return nt;
         }    }
         /**/    /**/
         l = abs_U_V_maxrshift( b1, l1, b2, l2, bt );    l = abs_U_V_maxrshift( b1, l1, b2, l2, bt );
         /**/    /**/
         if ( l == 0 ) {    if ( l == 0 ) {
                 PL(n1) = l1;      PL(n1) = l1;
                 return n1;      return n1;
         } else if ( l > 0 ) {    } else if ( l > 0 ) {
                 l1 = l;      l1 = l;
                 SWAP( b1, bt, int * ); SWAP( n1, nt, N );      SWAP( b1, bt, int * ); SWAP( n1, nt, N );
         } else {    } else {
                 l2 = -l;      l2 = -l;
                 SWAP( b2, bt, int * ); SWAP( n2, nt, N );      SWAP( b2, bt, int * ); SWAP( n2, nt, N );
         }    }
         goto loop;    goto loop;
 }  }
   
 #define RetTrueGCD(p,gw,gb,nr,l0) \  #define RetTrueGCD(p,gw,gb,nr,l0) \
Line 430  loop:
Line 430  loop:
 void gcdbinn( n1, n2, nr )  void gcdbinn( n1, n2, nr )
 N n1, n2, *nr;  N n1, n2, *nr;
 {  {
         int s1, s2, gw, gb, t1, t2;    int s1, s2, gw, gb, t1, t2;
         int w1, w2;    int w1, w2;
         N tn1, tn2, tnt, p;    N tn1, tn2, tnt, p;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = n2;      *nr = n2;
                 return;      return;
         } else if ( !n2 ) {    } else if ( !n2 ) {
                 *nr = n1;      *nr = n1;
                 return;      return;
         }    }
         s1 = trailingzerosn( n1, &w1 );    s1 = trailingzerosn( n1, &w1 );
         s2 = trailingzerosn( n2, &w2 );    s2 = trailingzerosn( n2, &w2 );
         if ( w1 == w2 ) gw = w1,  gb = s1 <= s2 ? s1 : s2;    if ( w1 == w2 ) gw = w1,  gb = s1 <= s2 ? s1 : s2;
         else if ( w1 < w2 ) gw = w1,  gb = s1;    else if ( w1 < w2 ) gw = w1,  gb = s1;
         else gw = w2,  gb = s2;    else gw = w2,  gb = s2;
         /*    /*
          *      true GCD must be multiplied by 2^{gw*BSH+gb}.     *  true GCD must be multiplied by 2^{gw*BSH+gb}.
          */     */
         t1 = PL(n1) - w1;    t1 = PL(n1) - w1;
         t2 = PL(n2) - w2;    t2 = PL(n2) - w2;
         if ( t1 < t2 ) t1 = t2;    if ( t1 < t2 ) t1 = t2;
         tn1 = W_NALLOC(t1);  tn2 = W_NALLOC(t1);  tnt = W_NALLOC(t1);    tn1 = W_NALLOC(t1);  tn2 = W_NALLOC(t1);  tnt = W_NALLOC(t1);
         rshiftn( n1, w1, s1, tn1 );    rshiftn( n1, w1, s1, tn1 );
         rshiftn( n2, w2, s2, tn2 );    rshiftn( n2, w2, s2, tn2 );
         p = igcd_binary( tn1, tn2, tnt );    p = igcd_binary( tn1, tn2, tnt );
         RetTrueGCD( p, gw, gb, nr, L0 )    RetTrueGCD( p, gw, gb, nr, L0 )
 }  }
   
   
 /*  /*
  *      The bmod gcd algorithm stated briefly in K.Weber's paper   *  The bmod gcd algorithm stated briefly in K.Weber's paper
  *      [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122].   *  [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122].
  *      It replaces the subtraction (n1 - n2) in the binary algorithm   *  It replaces the subtraction (n1 - n2) in the binary algorithm
  *      by (n1 - S*n2) with such an S that (n1 - S*n2) \equiv 0 \bmod 2^BSH,   *  by (n1 - S*n2) with such an S that (n1 - S*n2) \equiv 0 \bmod 2^BSH,
  *      which should improve the efficiency when n1 \gg n2.   *  which should improve the efficiency when n1 \gg n2.
  */   */
   
 /* subsidiary routines */  /* subsidiary routines */
Line 473  N n1, n2, *nr;
Line 473  N n1, n2, *nr;
 #ifdef CALL  #ifdef CALL
 static int u_div_v_mod_2toBSH( u, v )  static int u_div_v_mod_2toBSH( u, v )
 unsigned int u, v;  unsigned int u, v;
         /*    /*
          *    u/v mod 2^BSH.     *    u/v mod 2^BSH.
          */     */
 {  {
         unsigned int i, lsh_i, m;    unsigned int i, lsh_i, m;
   
         lsh_i = (sizeof(int) << 3) - 1;    lsh_i = (sizeof(int) << 3) - 1;
         m = i = 0;    m = i = 0;
         do {    do {
                 if ( u == 0 ) break;      if ( u == 0 ) break;
                 if ( (u << lsh_i) != 0 ) {      if ( (u << lsh_i) != 0 ) {
                         m += (1 << i);        m += (1 << i);
                         u -= (v << i);        u -= (v << i);
                 }      }
                 lsh_i--;      lsh_i--;
         } while ( ++i != BSH );    } while ( ++i != BSH );
         return m;    return m;
 }  }
   
 #define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2toBSH(U,V)  #define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2toBSH(U,V)
 #else /* CALL */  #else /* CALL */
 #define Comp_U_div_V_mod_BASE(U,V,R) {\  #define Comp_U_div_V_mod_BASE(U,V,R) {\
         unsigned int u = (U), v = (V), i, lsh; \    unsigned int u = (U), v = (V), i, lsh; \
         /* U and V are assumed to be odd */ \    /* U and V are assumed to be odd */ \
         i = R = 1, lsh = (sizeof(int) << 3) - 2;  u = (u - v); \    i = R = 1, lsh = (sizeof(int) << 3) - 2;  u = (u - v); \
         do {    if ( u == 0 ) break; \    do {  if ( u == 0 ) break; \
                 if ( (u << lsh) != 0 ) R += (1 << i),  u = (u - (v << i)); \      if ( (u << lsh) != 0 ) R += (1 << i),  u = (u - (v << i)); \
                 i++, lsh--; \      i++, lsh--; \
         } while ( i < BSH ); \    } while ( i < BSH ); \
 }  }
 #endif /* CALL */  #endif /* CALL */
 #else  #else
 #ifdef CALL  #ifdef CALL
 static int u_div_v_mod_2tos( u, v, s )  static int u_div_v_mod_2tos( u, v, s )
 int u, v, s;  int u, v, s;
         /*    /*
          *    u/v mod 2^s.     *    u/v mod 2^s.
          */     */
 {  {
         int i, lsh_i, mask, m;    int i, lsh_i, mask, m;
   
         mask = (1 << s) - 1;    mask = (1 << s) - 1;
         lsh_i = (sizeof(int) << 3) - 1;    lsh_i = (sizeof(int) << 3) - 1;
         m = i = 0;    m = i = 0;
         u &= mask,  v &= mask;    u &= mask,  v &= mask;
         do {    do {
                 if ( u == 0 ) break;      if ( u == 0 ) break;
                 if ( (u << lsh_i) != 0 ) {      if ( (u << lsh_i) != 0 ) {
                         m += (1 << i);        m += (1 << i);
                         u -= (v << i);        u -= (v << i);
                         u &= mask;        u &= mask;
                 }      }
                 lsh_i--;      lsh_i--;
         } while ( ++i != s );    } while ( ++i != s );
         return m;    return m;
 }  }
   
 #define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2tos(U,V,BSH)  #define Comp_U_div_V_mod_BASE(U,V,R) R = u_div_v_mod_2tos(U,V,BSH)
 #else  #else
 #define Comp_U_div_V_mod_BASE(U,V,R) {\  #define Comp_U_div_V_mod_BASE(U,V,R) {\
         int u = (U), v = (V), i, lsh; \    int u = (U), v = (V), i, lsh; \
         /* U and V are assumed to be odd */ \    /* U and V are assumed to be odd */ \
         i = R = 1, lsh = (sizeof(int) << 3) - 2;  u = (u - v) & BMASK; \    i = R = 1, lsh = (sizeof(int) << 3) - 2;  u = (u - v) & BMASK; \
         do {    if ( u == 0 ) break; \    do {  if ( u == 0 ) break; \
                 if ( (u << lsh) != 0 ) R += (1 << i),  u = (u - (v << i)) & BMASK; \      if ( (u << lsh) != 0 ) R += (1 << i),  u = (u - (v << i)) & BMASK; \
                 i++, lsh--; \      i++, lsh--; \
         } while ( i < BSH ); \    } while ( i < BSH ); \
 }  }
 #endif  #endif
 #endif  #endif
Line 547  int u, v, s;
Line 547  int u, v, s;
   
 static int bmod_n( nu, nv, na )  static int bmod_n( nu, nv, na )
 N nu, nv, na;  N nu, nv, na;
         /*    /*
          *      Computes (u[] \bmod v[]) >> (as much as possible) in r[].     *  Computes (u[] \bmod v[]) >> (as much as possible) in r[].
          */     */
 {  {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         unsigned int *u = BD(nu), *v = BD(nv), *r = BD(na);    unsigned int *u = BD(nu), *v = BD(nv), *r = BD(na);
         unsigned int *p, a, t, l, v0, vh, bv, v0r;    unsigned int *p, a, t, l, v0, vh, bv, v0r;
         int lu = PL(nu), lv = PL(nv), z;    int lu = PL(nu), lv = PL(nv), z;
 #else  #else
         int *u = BD(nu), lu = PL(nu), *v = BD(nv), lv = PL(nv),    int *u = BD(nu), lu = PL(nu), *v = BD(nv), lv = PL(nv),
                      *r = BD(na);           *r = BD(na);
         int *p, a, t, l, z, v0, vh, bv, v0r;    int *p, a, t, l, z, v0, vh, bv, v0r;
 #endif  #endif
   
         v0 = v[0];    v0 = v[0];
         if ( lv == 1 ) {    if ( lv == 1 ) {
                 if ( lu == 1 ) a = u[0] % v0;      if ( lu == 1 ) a = u[0] % v0;
                 else {      else {
                         p = &u[--lu];        p = &u[--lu];
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                         a = (*p) % v0, t = (unsigned)(-((int)v0)) % v0;        a = (*p) % v0, t = (unsigned)(-((int)v0)) % v0;
 #else  #else
                         a = (*p) % v0, t = BASE % v0;        a = (*p) % v0, t = BASE % v0;
 #endif  #endif
                         for ( ; --lu >= 0; a = l ) {        for ( ; --lu >= 0; a = l ) {
                                 --p;          --p;
                                 DMAR(a,t,*p,v0,l);          DMAR(a,t,*p,v0,l);
                                 /*  l <= (a*t + p[0])%v0   */          /*  l <= (a*t + p[0])%v0   */
                         }        }
                 }      }
                 if ( a == 0 ) return 0;      if ( a == 0 ) return 0;
                 while ( (a&1) == 0 ) a >>= 1;      while ( (a&1) == 0 ) a >>= 1;
                 *r = a;      *r = a;
                 return( PL(na) = 1 );      return( PL(na) = 1 );
         }    }
         Comp_U_div_V_mod_BASE( 1, v0, v0r );    Comp_U_div_V_mod_BASE( 1, v0, v0r );
         vh = v[lv -1];    vh = v[lv -1];
         BitWidth( vh, bv );    BitWidth( vh, bv );
         bv--;    bv--;
         t = 1 << bv;    t = 1 << bv;
         l = lv + 1;    l = lv + 1;
         for ( z = -1; lu > l || lu == l && u[lu-1] >= t; z = -z ) {    for ( z = -1; lu > l || lu == l && u[lu-1] >= t; z = -z ) {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 a = (v0r*u[0]);      a = (v0r*u[0]);
 #else  #else
                 a = (v0r*u[0])&BMASK;      a = (v0r*u[0])&BMASK;
 #endif  #endif
                 /**/      /**/
 #if 0000  #if 0000
 {  {
   int i;    int i;
Line 604  N nu, nv, na;
Line 604  N nu, nv, na;
   printf( "\n=>\n" );    printf( "\n=>\n" );
 }  }
 #endif  #endif
                 lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r );      lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r );
                 /**/      /**/
 #if 0000  #if 0000
 printf( "***lu=%d\n", lu );  printf( "***lu=%d\n", lu );
 if ( lu != 0 ) {  if ( lu != 0 ) {
Line 614  if ( lu != 0 ) {
Line 614  if ( lu != 0 ) {
   printf( "\n" );    printf( "\n" );
 }  }
 #endif  #endif
                 if ( lu == 0 ) return 0;      if ( lu == 0 ) return 0;
                 p = r;      p = r;
                 r = u;      r = u;
                 u = p;      u = p;
         }    }
         if ( lu < lv ) goto ret;    if ( lu < lv ) goto ret;
         t = u[lu-1];    t = u[lu-1];
         if ( lu > lv ) l = BSH;    if ( lu > lv ) l = BSH;
         else if ( t < vh ) goto ret;    else if ( t < vh ) goto ret;
         else l = 0;    else l = 0;
         BitWidth( t, a );    BitWidth( t, a );
         l += (a - bv);    l += (a - bv);
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         a = (v0r*u[0])&((unsigned)(-1) >> (BSH - l));    a = (v0r*u[0])&((unsigned)(-1) >> (BSH - l));
 #else  #else
         a = (v0r*u[0])&(BMASK >> (BSH - l));    a = (v0r*u[0])&(BMASK >> (BSH - l));
 #endif  #endif
 #if 0000  #if 0000
 {  {
Line 640  if ( lu != 0 ) {
Line 640  if ( lu != 0 ) {
   printf( "\n   =>\n" );    printf( "\n   =>\n" );
 }  }
 #endif  #endif
         /**/    /**/
         lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r );    lu = abs_U_aV_maxrshift( u, lu, a, v, lv, r );
         /**/    /**/
 #if 0000  #if 0000
 printf( "::: lu=%d\n", lu );  printf( "::: lu=%d\n", lu );
 if ( lu != 0 ) {  if ( lu != 0 ) {
Line 651  if ( lu != 0 ) {
Line 651  if ( lu != 0 ) {
   printf( "\n" );    printf( "\n" );
 }  }
 #endif  #endif
         if ( lu == 0 ) return 0;    if ( lu == 0 ) return 0;
         z = -z;    z = -z;
 ret:    if ( z > 0 ) return( PL(na) = lu );  ret:  if ( z > 0 ) return( PL(na) = lu );
         PL(nu) = lu;    PL(nu) = lu;
         return( -lu );    return( -lu );
 }  }
   
   
 static N igcd_bmod( n1, n2, nt )  static N igcd_bmod( n1, n2, nt )
 N n1, n2, nt;  N n1, n2, nt;
         /*  both n1 and n2 are assumed to be odd */    /*  both n1 and n2 are assumed to be odd */
 {  {
         int l1, l2;    int l1, l2;
         int l;    int l;
   
         if ( (l = cmpn( n1, n2 )) == 0 ) return n1;    if ( (l = cmpn( n1, n2 )) == 0 ) return n1;
         else if ( l < 0 ) { SWAP( n1, n2, N ); }    else if ( l < 0 ) { SWAP( n1, n2, N ); }
         IniDiv( n1, n2 );    IniDiv( n1, n2 );
         if ( UNIN(n2) ) return 0;    if ( UNIN(n2) ) return 0;
 loop:   if ( (l1 = PL(n1)) <= 2 && (l2 = PL(n2)) <= 2 ) {  loop:  if ( (l1 = PL(n1)) <= 2 && (l2 = PL(n2)) <= 2 ) {
                 l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) );      l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) );
                 if ( l == 0 ) return 0;      if ( l == 0 ) return 0;
                 PL(nt) = l;      PL(nt) = l;
                 return nt;      return nt;
         }    }
         /**/    /**/
         l = bmod_n( n1, n2, nt );    l = bmod_n( n1, n2, nt );
         /**/    /**/
         if ( l == 0 ) return n2;    if ( l == 0 ) return n2;
         else if ( l > 0 ) {    else if ( l > 0 ) {
                 N tmp = n1;      N tmp = n1;
   
                 n1 = n2;      n1 = n2;
                 n2 = nt;      n2 = nt;
                 nt = tmp;      nt = tmp;
         } else SWAP( n1, n2, N );    } else SWAP( n1, n2, N );
         goto loop;    goto loop;
 }  }
   
 void gcdbmodn( n1, n2, nr )  void gcdbmodn( n1, n2, nr )
 N n1, n2, *nr;  N n1, n2, *nr;
 {  {
         int s1, s2, gw, gb, t1, t2;    int s1, s2, gw, gb, t1, t2;
         int w1, w2;    int w1, w2;
         N tn1, tn2, tnt, p;    N tn1, tn2, tnt, p;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = n2;      *nr = n2;
                 return;      return;
         } else if ( !n2 ) {    } else if ( !n2 ) {
                 *nr = n1;      *nr = n1;
                 return;      return;
         }    }
         s1 = trailingzerosn( n1, &w1 );    s1 = trailingzerosn( n1, &w1 );
         s2 = trailingzerosn( n2, &w2 );    s2 = trailingzerosn( n2, &w2 );
         if ( w1 == w2 ) gw = w1,  gb = s1 <= s2 ? s1 : s2;    if ( w1 == w2 ) gw = w1,  gb = s1 <= s2 ? s1 : s2;
         else if ( w1 < w2 ) gw = w1,  gb = s1;    else if ( w1 < w2 ) gw = w1,  gb = s1;
         else gw = w2,  gb = s2;    else gw = w2,  gb = s2;
         /*    /*
          *      true GCD must be multiplied by 2^{gw*BSH+gs}.     *  true GCD must be multiplied by 2^{gw*BSH+gs}.
          */     */
         t1 = PL(n1) - w1;    t1 = PL(n1) - w1;
         t2 = PL(n2) - w2;    t2 = PL(n2) - w2;
         if ( t1 < t2 ) t1 = t2;    if ( t1 < t2 ) t1 = t2;
         tn1 = W_NALLOC(t1);  tn2 = W_NALLOC(t1);  tnt = W_NALLOC(t1);    tn1 = W_NALLOC(t1);  tn2 = W_NALLOC(t1);  tnt = W_NALLOC(t1);
         rshiftn( n1, w1, s1, tn1 );    rshiftn( n1, w1, s1, tn1 );
         rshiftn( n2, w2, s2, tn2 );    rshiftn( n2, w2, s2, tn2 );
         p = igcd_bmod( tn1, tn2, tnt );    p = igcd_bmod( tn1, tn2, tnt );
         RetTrueGCD( p, gw, gb, nr, L0 )    RetTrueGCD( p, gw, gb, nr, L0 )
 }  }
   
 /*  /*
  *      The accelerated integer GCD algorithm by K.Weber   *  The accelerated integer GCD algorithm by K.Weber
  *      [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]:   *  [ACM TOMS, Vol.21, No. 1 (1995), pp. 111-122]:
  */   */
   
 static int ReducedRatMod( x, y, pn, pd )  static int ReducedRatMod( x, y, pn, pd )
Line 733  N x, y;
Line 733  N x, y;
 unsigned  unsigned
 #endif  #endif
 int *pn, *pd;  int *pn, *pd;
         /*    /*
          *    Let m = 2^{2*BSH} = 2*BASE.  We assume x, y > 0 and \gcd(x,m)     *    Let m = 2^{2*BSH} = 2*BASE.  We assume x, y > 0 and \gcd(x,m)
          *    = \gcd(y,m) = 1.  This routine computes n and d (resp. returned     *    = \gcd(y,m) = 1.  This routine computes n and d (resp. returned
          *    in *pn and *pd) such that 0 < n, |d| < \sqrt{m} and     *    in *pn and *pd) such that 0 < n, |d| < \sqrt{m} and
          *    n*y \equiv x*d \bmod m.     *    n*y \equiv x*d \bmod m.
          */     */
 {  {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         unsigned int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l;    unsigned int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l;
         unsigned int th, tl, l1, l2, i, ir;    unsigned int th, tl, l1, l2, i, ir;
         int s1, s2;    int s1, s2;
 #else  #else
         int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l;    int n1h, n1l, d1h, d1l, n2h, n2l, d2h, d2l;
         int th, tl, l1, l2, i, ir;    int th, tl, l1, l2, i, ir;
         int s1, s2;    int s1, s2;
 #endif  #endif
   
         {    {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 unsigned      unsigned
 #endif  #endif
                 int xh, xl, yh, yl, lsh_i;      int xh, xl, yh, yl, lsh_i;
   
                 xl = BD(x)[0];      xl = BD(x)[0];
                 xh = PL(x) > 1 ? BD(x)[1] : 0;      xh = PL(x) > 1 ? BD(x)[1] : 0;
                 yl = BD(y)[0];      yl = BD(y)[0];
                 yh = PL(y) > 1 ? BD(y)[1] : 0;      yh = PL(y) > 1 ? BD(y)[1] : 0;
 #if 0000  #if 0000
 printf( "*** RedRatMod: (0x%08x:0x%08x=%u*2^%d+%u)\n               /(0x%08x:0x%08x=%u*2^%d+%u) mod 2^%d\n",  printf( "*** RedRatMod: (0x%08x:0x%08x=%u*2^%d+%u)\n               /(0x%08x:0x%08x=%u*2^%d+%u) mod 2^%d\n",
         xh,xl, xh,BSH,xl, yh,yl, yh,BSH,yl, BSH );          xh,xl, xh,BSH,xl, yh,yl, yh,BSH,yl, BSH );
 #endif  #endif
                 Comp_U_div_V_mod_BASE( xl, yl, n2l );      Comp_U_div_V_mod_BASE( xl, yl, n2l );
                 DM(n2l,yl,th,tl)  /* n2l*yl = tl+th*BASE, where tl==xl. */;      DM(n2l,yl,th,tl)  /* n2l*yl = tl+th*BASE, where tl==xl. */;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 xh -= th;      xh -= th;
 #else  #else
                 if ( xh > th ) xh -= th;      if ( xh > th ) xh -= th;
                 else xh += (BASE - th);      else xh += (BASE - th);
 #endif  #endif
                 DM(n2l,yh,th,tl)  /* n2l*yh = tl+th*BASE. */;      DM(n2l,yh,th,tl)  /* n2l*yh = tl+th*BASE. */;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                 xh -= tl;      xh -= tl;
 #else  #else
                 if ( xh > tl ) xh -= tl;      if ( xh > tl ) xh -= tl;
                 else xh += (BASE - tl);      else xh += (BASE - tl);
 #endif  #endif
 /*              n2h = i = 0,  lsh_i = 31;*/  /*    n2h = i = 0,  lsh_i = 31;*/
                 n2h = i = 0,  lsh_i = BIT_WIDTH_OF_INT -1;      n2h = i = 0,  lsh_i = BIT_WIDTH_OF_INT -1;
                 do {      do {
                         if ( xh == 0 ) break;        if ( xh == 0 ) break;
                         if ( (xh << lsh_i) != 0 ) {        if ( (xh << lsh_i) != 0 ) {
                                 n2h += (1 << i);          n2h += (1 << i);
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                                 tl = yl << i;          tl = yl << i;
                                 xh -= tl;          xh -= tl;
 #else  #else
                                 tl = (yl << i)&BMASK;          tl = (yl << i)&BMASK;
                                 if ( xh > tl ) xh -= tl;          if ( xh > tl ) xh -= tl;
                                 else xh += (BASE - tl);          else xh += (BASE - tl);
 #endif  #endif
                         }        }
                         lsh_i--;        lsh_i--;
                 } while ( ++i != BSH );      } while ( ++i != BSH );
         }    }
         /*    /*
          *      n2l + n2h*BASE = x/y mod 2^{2*BSH}.     *  n2l + n2h*BASE = x/y mod 2^{2*BSH}.
          */     */
 #if 0000  #if 0000
 printf( "=====> 0x%08x(%u) + 2^%d*0x%08x(%u)\n", n2l, n2l,  BSH, n2h, n2h );  printf( "=====> 0x%08x(%u) + 2^%d*0x%08x(%u)\n", n2l, n2l,  BSH, n2h, n2h );
 #endif  #endif
         d2h = 0, d2l = 1, s2 = 1;    d2h = 0, d2l = 1, s2 = 1;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
         if ( n2h == 0 ) goto done;    if ( n2h == 0 ) goto done;
         BitWidth( n2h, l2 );    BitWidth( n2h, l2 );
         if ( l2 == BSH ) {    if ( l2 == BSH ) {
           d1h = 0, d1l = 1, s1 = -1;      d1h = 0, d1l = 1, s1 = -1;
           th = n2h, tl = n2l;      th = n2h, tl = n2l;
         } else {    } else {
           i = BSH - l2;      i = BSH - l2;
           d1h = 0, d1l = 1 << i, s1 = -1;      d1h = 0, d1l = 1 << i, s1 = -1;
           th = (n2h << i) | (n2l >> l2);      th = (n2h << i) | (n2l >> l2);
           tl = n2l << i;      tl = n2l << i;
         }    }
         if ( tl == 0 ) n1h = - th,  n1l = 0;    if ( tl == 0 ) n1h = - th,  n1l = 0;
         else n1h = BMASK - th,  n1l = - tl;    else n1h = BMASK - th,  n1l = - tl;
         /**/    /**/
         if ( n1h < n2h || (n1h == n2h && n1l < n2l) ) goto swap12;    if ( n1h < n2h || (n1h == n2h && n1l < n2l) ) goto swap12;
         BitWidth( n1h, l1 );    BitWidth( n1h, l1 );
         goto sub12;    goto sub12;
 #else  #else
         n1h = BASE, n1l = 0, l1 = BSH,    n1h = BASE, n1l = 0, l1 = BSH,
         d1h = d1l = 0, s1 = 0;    d1h = d1l = 0, s1 = 0;
 #endif  #endif
         /**/    /**/
         while ( n2h != 0 ) {    while ( n2h != 0 ) {
                 BitWidth( n2h, l2 );      BitWidth( n2h, l2 );
         sub12:  ir = BSH - (i = l1 - l2);    sub12:  ir = BSH - (i = l1 - l2);
                 do {      do {
                         if ( i == 0 ) th = n2h, tl = n2l;        if ( i == 0 ) th = n2h, tl = n2l;
                         else        else
                                 th = (n2h << i) | (n2l >> ir),          th = (n2h << i) | (n2l >> ir),
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                                 tl = n2l << i;          tl = n2l << i;
 #else  #else
                                 tl = (n2l << i) & BMASK;          tl = (n2l << i) & BMASK;
 #endif  #endif
                         if ( th > n1h || (th == n1h && tl > n1l) ) goto next_i;        if ( th > n1h || (th == n1h && tl > n1l) ) goto next_i;
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                         if ( tl > n1l ) n1h--;        if ( tl > n1l ) n1h--;
                         n1l -= tl;        n1l -= tl;
 #else  #else
                         if ( tl <= n1l ) n1l -= tl;        if ( tl <= n1l ) n1l -= tl;
                         else n1l += (BASE - tl), n1h--;        else n1l += (BASE - tl), n1h--;
 #endif  #endif
                         n1h -= th;        n1h -= th;
                         /*  (s1:d1h,d1l) -= ((s2:d2h,d2l) << i);  */        /*  (s1:d1h,d1l) -= ((s2:d2h,d2l) << i);  */
                         if ( s2 != 0 ) {        if ( s2 != 0 ) {
                                 if ( i == 0 ) th = d2h, tl = d2l;          if ( i == 0 ) th = d2h, tl = d2l;
                                 else          else
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                                         th = (d2h << i) | (d2l >> ir),            th = (d2h << i) | (d2l >> ir),
                                         tl = (d2l << i);            tl = (d2l << i);
 #else  #else
                                         th = (d2h << i)&BMASK | (d2l >> ir),            th = (d2h << i)&BMASK | (d2l >> ir),
                                         tl = (d2l << i)&BMASK;            tl = (d2l << i)&BMASK;
 #endif  #endif
                                 if ( s1 == 0 )          if ( s1 == 0 )
                                         s1 = -s2, d1h = th, d1l = tl;            s1 = -s2, d1h = th, d1l = tl;
                                 else if ( s1 != s2 ) {          else if ( s1 != s2 ) {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                                         d1l += tl,            d1l += tl,
                                         d1h += th;            d1h += th;
                                         if ( d1l < tl ) d1h++;            if ( d1l < tl ) d1h++;
 #else  #else
                                         tl += d1l;            tl += d1l;
                                         d1l = tl&BMASK;            d1l = tl&BMASK;
                                         d1h = (d1h + th + (tl >> BSH))&BMASK;            d1h = (d1h + th + (tl >> BSH))&BMASK;
 #endif  #endif
                                         if ( d1h == 0 && d1l == 0 ) s1 = 0;            if ( d1h == 0 && d1l == 0 ) s1 = 0;
                                 } else if ( d1h > th ) {          } else if ( d1h > th ) {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                                         if ( d1l < tl ) d1h--;            if ( d1l < tl ) d1h--;
                                         d1l -= tl;            d1l -= tl;
 #else  #else
                                         if ( d1l >= tl ) d1l -= tl;            if ( d1l >= tl ) d1l -= tl;
                                         else d1l += (BASE - tl), d1h--;            else d1l += (BASE - tl), d1h--;
 #endif  #endif
                                         d1h -= th;            d1h -= th;
                                 } else if ( d1h == th ) {          } else if ( d1h == th ) {
                                         d1h = 0;            d1h = 0;
                                         if ( d1l == tl ) s1 = d2h = 0;            if ( d1l == tl ) s1 = d2h = 0;
                                         else if ( d1l > tl ) d1l -= tl;            else if ( d1l > tl ) d1l -= tl;
                                         else d1l = tl - d1l, s1 = -s1;            else d1l = tl - d1l, s1 = -s1;
                                 } else {          } else {
 #if BSH == BIT_WIDTH_OF_INT  #if BSH == BIT_WIDTH_OF_INT
                                         if ( tl < d1l ) th--;            if ( tl < d1l ) th--;
                                         d1l = tl - d1l;            d1l = tl - d1l;
 #else  #else
                                         if ( tl >= d1l ) d1l = tl - d1l;            if ( tl >= d1l ) d1l = tl - d1l;
                                         else d1l = tl + (BASE - d1l), th--;            else d1l = tl + (BASE - d1l), th--;
 #endif  #endif
                                         d1h = th - d1h;            d1h = th - d1h;
                                         s1 = -s1;            s1 = -s1;
                                 }          }
                         }        }
                 next_i: i--, ir++;      next_i:  i--, ir++;
                 } while ( n1h > n2h || (n1h == n2h && n1l >= n2l) );      } while ( n1h > n2h || (n1h == n2h && n1l >= n2l) );
         swap12: /*  swap 1 and 2  */    swap12:  /*  swap 1 and 2  */
                 th = n1h, tl = n1l;      th = n1h, tl = n1l;
                 n1h = n2h, n1l = n2l;      n1h = n2h, n1l = n2l;
                 n2h = th, n2l = tl;      n2h = th, n2l = tl;
                 l1 = l2;      l1 = l2;
                 th = d1h, tl = d1l, i = s1;      th = d1h, tl = d1l, i = s1;
                 d1h = d2h, d1l = d2l, s1 = s2;      d1h = d2h, d1l = d2l, s1 = s2;
                 d2h = th, d2l = tl, s2 = i;      d2h = th, d2l = tl, s2 = i;
         }    }
         /**/    /**/
 done:   *pn = n2l,  *pd = d2l;  done:  *pn = n2l,  *pd = d2l;
         return s2;    return s2;
 }  }
   
 static int igcd_spurious_factor;  static int igcd_spurious_factor;
   
 #define SaveN(s,d) {\  #define SaveN(s,d) {\
         int i, l; \    int i, l; \
         for ( l = PL(d) = PL(s), i = 0; i < l; i++ ) BD(d)[i] = BD(s)[i]; \    for ( l = PL(d) = PL(s), i = 0; i < l; i++ ) BD(d)[i] = BD(s)[i]; \
 }  }
   
 static N igcd_acc( n1, n2, nt )  static N igcd_acc( n1, n2, nt )
 N n1, n2, nt;  N n1, n2, nt;
         /*  both n1 and n2 are assumed to be odd */    /*  both n1 and n2 are assumed to be odd */
 {  {
         int l1, l2, *b1, *b2, bw1, bw2;    int l1, l2, *b1, *b2, bw1, bw2;
         int l;    int l;
         int n, d;    int n, d;
         N p, s1, s2;    N p, s1, s2;
   
         if ( (l = cmpn( n1, n2 )) == 0 ) return n1;    if ( (l = cmpn( n1, n2 )) == 0 ) return n1;
         else if ( l < 0 ) { SWAP( n1, n2, N ); }    else if ( l < 0 ) { SWAP( n1, n2, N ); }
         if ( ShouldCompRemInit(n1,n2) ) {    if ( ShouldCompRemInit(n1,n2) ) {
                 int w, b;      int w, b;
   
                 divn( n1, n2, &s1, &s2 );      divn( n1, n2, &s1, &s2 );
                 if ( !s2 ) return n2;      if ( !s2 ) return n2;
                 b = trailingzerosn( s2, &w );      b = trailingzerosn( s2, &w );
                 p = n1;  n1 = n2;  n2 = p;      p = n1;  n1 = n2;  n2 = p;
                 rshiftn( s2, w, b, n2 );      rshiftn( s2, w, b, n2 );
                 if ( UNIN(n2) ) return 0;      if ( UNIN(n2) ) return 0;
                 l1 = PL(n1);      l1 = PL(n1);
                 if ( !s1 || PL(s1) < l1 ) s1 = NALLOC(l1);      if ( !s1 || PL(s1) < l1 ) s1 = NALLOC(l1);
         } else if ( UNIN(n2) ) return 0;    } else if ( UNIN(n2) ) return 0;
         else {    else {
                 s1 = NALLOC(PL(n1));      s1 = NALLOC(PL(n1));
                 s2 = NALLOC(PL(n2));      s2 = NALLOC(PL(n2));
         }    }
         SaveN( n1, s1 );    SaveN( n1, s1 );
         SaveN( n2, s2 );    SaveN( n2, s2 );
         igcd_spurious_factor = 0;    igcd_spurious_factor = 0;
 loop:   l1 = PL(n1), l2 = PL(n2);  loop:  l1 = PL(n1), l2 = PL(n2);
         if ( l1 <= 2 && l2 <= 2 ) {    if ( l1 <= 2 && l2 <= 2 ) {
                 l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) );      l = igcd_binary_2w( BD(n1), l1, BD(n2), l2, BD(nt) );
                 if ( l == 0 ) return 0;      if ( l == 0 ) return 0;
                 PL(nt) = l;      PL(nt) = l;
                 SWAP( n2, nt, N );      SWAP( n2, nt, N );
                 goto ret;      goto ret;
         }    }
         /**/    /**/
         b1 = BD(n1), b2 = BD(n2);    b1 = BD(n1), b2 = BD(n2);
         BitWidth( b1[l1 -1], bw1 );    BitWidth( b1[l1 -1], bw1 );
         BitWidth( b2[l2 -1], bw2 );    BitWidth( b2[l2 -1], bw2 );
         if ( (l1*BSH + bw1) - (l2*BSH + bw2) <= igcdacc_thre ) {    if ( (l1*BSH + bw1) - (l2*BSH + bw2) <= igcdacc_thre ) {
                 l = ReducedRatMod( n1, n2, &n, &d );      l = ReducedRatMod( n1, n2, &n, &d );
                 l = l < 0 ? aUplusbV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ) :      l = l < 0 ? aUplusbV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) ) :
                     abs_axU_bxV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) );          abs_axU_bxV_maxrshift( n, b2, l2, d, b1, l1, BD(nt) );
                 igcd_spurious_factor++;      igcd_spurious_factor++;
                 if ( l == 0 ) goto ret;      if ( l == 0 ) goto ret;
                 PL(nt) = l;      PL(nt) = l;
         } else {    } else {
                 l = bmod_n( n1, n2, nt );      l = bmod_n( n1, n2, nt );
                 if ( l == 0 ) goto ret;      if ( l == 0 ) goto ret;
                 else if ( l < 0 ) {      else if ( l < 0 ) {
                         SWAP( n1, n2, N );        SWAP( n1, n2, N );
                         goto loop;        goto loop;
                 }      }
         }    }
         p = n1;    p = n1;
         n1 = n2;    n1 = n2;
         n2 = nt;    n2 = nt;
         nt = p;    nt = p;
         goto loop;    goto loop;
         /**/    /**/
 ret:    if ( igcd_spurious_factor != 0 && !UNIN(n2) ) {  ret:  if ( igcd_spurious_factor != 0 && !UNIN(n2) ) {
                 if ( (p = igcd_bmod( n2, s1, n1 )) == 0 ) return 0;      if ( (p = igcd_bmod( n2, s1, n1 )) == 0 ) return 0;
                 if ( (p = igcd_bmod( p, s2, nt )) == 0 ) return 0;      if ( (p = igcd_bmod( p, s2, nt )) == 0 ) return 0;
                 return p;      return p;
         } else return n2;    } else return n2;
 }  }
   
 void gcdaccn( n1, n2, nr )  void gcdaccn( n1, n2, nr )
 N n1, n2, *nr;  N n1, n2, *nr;
 {  {
         int s1, s2, gw, gb, t1, t2;    int s1, s2, gw, gb, t1, t2;
         int w1, w2;    int w1, w2;
         N tn1, tn2, tnt, p;    N tn1, tn2, tnt, p;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = n2;      *nr = n2;
                 return;      return;
         } else if ( !n2 ) {    } else if ( !n2 ) {
                 *nr = n1;      *nr = n1;
                 return;      return;
         }    }
         s1 = trailingzerosn( n1, &w1 );    s1 = trailingzerosn( n1, &w1 );
         s2 = trailingzerosn( n2, &w2 );    s2 = trailingzerosn( n2, &w2 );
         if ( w1 == w2 ) gw = w1,  gb = s1 <= s2 ? s1 : s2;    if ( w1 == w2 ) gw = w1,  gb = s1 <= s2 ? s1 : s2;
         else if ( w1 < w2 ) gw = w1,  gb = s1;    else if ( w1 < w2 ) gw = w1,  gb = s1;
         else gw = w2,  gb = s2;    else gw = w2,  gb = s2;
         /*    /*
          *      true GCD must be multiplied by 2^{gw*BSH+gs}.     *  true GCD must be multiplied by 2^{gw*BSH+gs}.
          */     */
         t1 = PL(n1) - w1;    t1 = PL(n1) - w1;
         t2 = PL(n2) - w2;    t2 = PL(n2) - w2;
         if ( t1 < t2 ) t1 = t2;    if ( t1 < t2 ) t1 = t2;
         tn1 = W_NALLOC(t1);  tn2 = W_NALLOC(t1);  tnt = W_NALLOC(t1);    tn1 = W_NALLOC(t1);  tn2 = W_NALLOC(t1);  tnt = W_NALLOC(t1);
         rshiftn( n1, w1, s1, tn1 );    rshiftn( n1, w1, s1, tn1 );
         rshiftn( n2, w2, s2, tn2 );    rshiftn( n2, w2, s2, tn2 );
         /**/    /**/
         p = igcd_acc( tn1, tn2, tnt );    p = igcd_acc( tn1, tn2, tnt );
         /**/    /**/
         if ( p == 0 ) goto L0;    if ( p == 0 ) goto L0;
         RetTrueGCD( p, gw, gb, nr, L0 )    RetTrueGCD( p, gw, gb, nr, L0 )
 }  }
   
   
Line 1032  N n1, n2, *nr;
Line 1032  N n1, n2, *nr;
 void gcdn_HMEXT( n1, n2, nr )  void gcdn_HMEXT( n1, n2, nr )
 N n1, n2, *nr;  N n1, n2, *nr;
 {  {
         int b1, b2, w1, w2, gw, gb;    int b1, b2, w1, w2, gw, gb;
         int l1, l2;    int l1, l2;
         N tn1, tn2, tnt, a;    N tn1, tn2, tnt, a;
   
         if ( !n1 ) {    if ( !n1 ) {
                 *nr = n2; return;      *nr = n2; return;
         } else if ( !n2 ) {    } else if ( !n2 ) {
                 *nr = n1; return;      *nr = n1; return;
         }    }
         b1 = trailingzerosn( n1, &w1 );    b1 = trailingzerosn( n1, &w1 );
         b2 = trailingzerosn( n2, &w2 );    b2 = trailingzerosn( n2, &w2 );
         if ( w1 == w2 ) gw = w1, gb = b1 <= b2 ? b1 : b2;    if ( w1 == w2 ) gw = w1, gb = b1 <= b2 ? b1 : b2;
         else if ( w1 < w2 ) gw = w1, gb = b1;    else if ( w1 < w2 ) gw = w1, gb = b1;
         else gw = w2, gb = b2;    else gw = w2, gb = b2;
         /*    /*
          *      true GCD must be multiplied by 2^{gw*BSH+gb}.     *  true GCD must be multiplied by 2^{gw*BSH+gb}.
          */     */
         l1 = PL(n1) - w1;    l1 = PL(n1) - w1;
         l2 = PL(n2) - w2;    l2 = PL(n2) - w2;
         if ( l1 < l2 ) l1 = l2;    if ( l1 < l2 ) l1 = l2;
         tn1 = W_NALLOC( l1 );  tn2 = W_NALLOC( l1 );  tnt = W_NALLOC( l1 );    tn1 = W_NALLOC( l1 );  tn2 = W_NALLOC( l1 );  tnt = W_NALLOC( l1 );
         rshiftn( n1, w1, b1, tn1 );    rshiftn( n1, w1, b1, tn1 );
         rshiftn( n2, w2, b2, tn2 );    rshiftn( n2, w2, b2, tn2 );
         /**/    /**/
         if ( igcd_algorithm == 1 ) {    if ( igcd_algorithm == 1 ) {
                 a = igcd_binary( tn1, tn2, tnt );      a = igcd_binary( tn1, tn2, tnt );
         } else if ( igcd_algorithm == 2 ) {    } else if ( igcd_algorithm == 2 ) {
                 a = igcd_bmod( tn1, tn2, tnt );      a = igcd_bmod( tn1, tn2, tnt );
         } else {    } else {
                 a = igcd_acc( tn1, tn2, tnt );      a = igcd_acc( tn1, tn2, tnt );
                 if ( igcd_spurious_factor != 0 ) {      if ( igcd_spurious_factor != 0 ) {
                 }      }
         }    }
         RetTrueGCD( a, gw, gb, nr, L0 )    RetTrueGCD( a, gw, gb, nr, L0 )
 }  }
   
   
Line 1075  N maxrshn( n, p )
Line 1075  N maxrshn( n, p )
 N n;  N n;
 int *p;  int *p;
 {  {
         int nw, nb, c, l;    int nw, nb, c, l;
         N new;    N new;
   
         nb = trailingzerosn( n, &nw );    nb = trailingzerosn( n, &nw );
         l = PL(n);    l = PL(n);
         c = BD(n)[l -1];    c = BD(n)[l -1];
         l -= nw;    l -= nw;
         if ( (c >> nb) == 0 ) l--;    if ( (c >> nb) == 0 ) l--;
         new = NALLOC(l);    new = NALLOC(l);
         rshiftn( n, nw, nb, new );    rshiftn( n, nw, nb, new );
         *p = nb + nw*BSH;    *p = nb + nw*BSH;
         return new;    return new;
 }  }
 #endif  #endif
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>