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

Diff for /OpenXM_contrib2/asir2000/engine/upo.c between version 1.5 and 1.6

version 1.5, 2015/08/14 13:51:55 version 1.6, 2018/03/29 01:32:52
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/engine/upo.c,v 1.4 2015/08/06 10:01:52 fujimoto Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/upo.c,v 1.5 2015/08/14 13:51:55 fujimoto Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 99  default:(p3)=(p1)^mulup2_bb((p2),(q));(q3)=(q1)^mulup2
Line 99  default:(p3)=(p1)^mulup2_bb((p2),(q));(q3)=(q1)^mulup2
   
 #define GF2M_MUL_1_HALF(high,low,a,b)\  #define GF2M_MUL_1_HALF(high,low,a,b)\
 {\  {\
         unsigned int _a,_b,_c,_d,_t,_s;\    unsigned int _a,_b,_c,_d,_t,_s;\
         unsigned int _ll,_cc;\    unsigned int _ll,_cc;\
         _a = (a);\    _a = (a);\
         _b = (b);\    _b = (b);\
         _c = _a&0xff00ff00;     /* c = a4 0 a2 0 */\    _c = _a&0xff00ff00;  /* c = a4 0 a2 0 */\
         _a ^= _c;                               /* a = 0 a3 0 a1 */\    _a ^= _c;        /* a = 0 a3 0 a1 */\
         _d = _b&0xff00ff00;     /* d = 0 0 b2 0 */\    _d = _b&0xff00ff00;  /* d = 0 0 b2 0 */\
         _b ^= _d;                               /* b = 0 0 0 b1 */\    _b ^= _d;        /* b = 0 0 0 b1 */\
         _c |= (_d>>8);          /* c = a4 00 a2 b2 */\    _c |= (_d>>8);    /* c = a4 00 a2 b2 */\
         _b |= (_a<<8);          /* b = a3 00 a1 b1 */\    _b |= (_a<<8);    /* b = a3 00 a1 b1 */\
         _a = (_c>>16);          /* a = a4 00 */\    _a = (_c>>16);    /* a = a4 00 */\
         _c &= 0xffff;           /* c = a2 b2 */\    _c &= 0xffff;    /* c = a2 b2 */\
         _d = (_b>>16);          /* d = a3 00 */\    _d = (_b>>16);    /* d = a3 00 */\
         _b &= 0xffff;           /* b = a1 b1 */\    _b &= 0xffff;    /* b = a1 b1 */\
         _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\    _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\
         _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\    _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\
         _a ^= _c; _d^= _b;\    _a ^= _c; _d^= _b;\
         _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\    _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
         _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\    _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
         _cc ^= _ll;\    _cc ^= _ll;\
         (low) = (unsigned int)(_ll^(_cc<<16));\    (low) = (unsigned int)(_ll^(_cc<<16));\
         (high) = (unsigned int)(_cc>>16);\    (high) = (unsigned int)(_cc>>16);\
 }  }
   
 /*  /*
Line 131  default:(p3)=(p1)^mulup2_bb((p2),(q));(q3)=(q1)^mulup2
Line 131  default:(p3)=(p1)^mulup2_bb((p2),(q));(q3)=(q1)^mulup2
   
 #define GF2M_MUL_1(high,low,a,b)\  #define GF2M_MUL_1(high,low,a,b)\
 {\  {\
         unsigned int _a,_b,_c,_d,_t,_s;\    unsigned int _a,_b,_c,_d,_t,_s;\
         unsigned int _uu,_ll,_cc;\    unsigned int _uu,_ll,_cc;\
         _a = (a);\    _a = (a);\
         _b = (b);\    _b = (b);\
         _c = _a&0xff00ff00;     /* _c = a4 0 a2 0 */\    _c = _a&0xff00ff00;  /* _c = a4 0 a2 0 */\
         _a ^= _c;                       /*  a = 0 a3 0 a1 */\    _a ^= _c;      /*  a = 0 a3 0 a1 */\
         _d = _b&0xff00ff00;     /* _d = b4 0 b2 0 */\    _d = _b&0xff00ff00;  /* _d = b4 0 b2 0 */\
         _b ^= _d;                       /*  b = 0 b3 0 b1 */\    _b ^= _d;      /*  b = 0 b3 0 b1 */\
         _c |= (_d>>8);          /* _c = a4 b4 a2 b2 */\    _c |= (_d>>8);    /* _c = a4 b4 a2 b2 */\
         _b |= (_a<<8);          /* b = a3 b3 a1 b1 */\    _b |= (_a<<8);    /* b = a3 b3 a1 b1 */\
         _a = (_c>>16);          /* a = a4 b4 */\    _a = (_c>>16);    /* a = a4 b4 */\
         _c &= 0xffff;           /* _c = a2 b2 */\    _c &= 0xffff;    /* _c = a2 b2 */\
         _d = (_b>>16);          /* _d = a3 b3 */\    _d = (_b>>16);    /* _d = a3 b3 */\
         _b &= 0xffff;           /* b = a1 b1 */\    _b &= 0xffff;    /* b = a1 b1 */\
         _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\    _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
         _uu = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\    _uu = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
         _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\    _t = up2_bbtab[_c]; _s = up2_bbtab[_b];\
         _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\    _ll = ((_t^_s^up2_bbtab[_c^_b])<<8)^(_t<<16)^_s;\
         _a ^= _c; _d^= _b;\    _a ^= _c; _d^= _b;\
         _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\    _t = up2_bbtab[_a]; _s = up2_bbtab[_d];\
         _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\    _cc = ((_t^_s^up2_bbtab[_a^_d])<<8)^(_t<<16)^_s;\
         _cc ^= (_ll^_uu);\    _cc ^= (_ll^_uu);\
         (low) = (unsigned int)(_ll^(_cc<<16));\    (low) = (unsigned int)(_ll^(_cc<<16));\
         (high) = (unsigned int)(_uu^(_cc>>16));\    (high) = (unsigned int)(_uu^(_cc>>16));\
 }  }
 #endif  #endif
   
 #define REMUP2_ONESTEP(a,q,s,n,w)\  #define REMUP2_ONESTEP(a,q,s,n,w)\
 if ( !s ) a[q] ^= w;\  if ( !s ) a[q] ^= w;\
 else {\  else {\
         a[q] ^= (w<<s);\    a[q] ^= (w<<s);\
         if ( q+1 <= n )\    if ( q+1 <= n )\
                 a[q+1] ^= (w>>(32-s));\      a[q+1] ^= (w>>(32-s));\
 }  }
   
 void ptoup2(n,nr)  void ptoup2(n,nr)
 P n;  P n;
 UP2 *nr;  UP2 *nr;
 {  {
         DCP dc;    DCP dc;
         UP2 r,s;    UP2 r,s;
         int d,w,i;    int d,w,i;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else if ( OID(n) == O_N )    else if ( OID(n) == O_N )
                 if ( EVENN(NM((Q)n)) )      if ( EVENN(NM((Q)n)) )
                         *nr = 0;        *nr = 0;
                 else      else
                         *nr = ONEUP2;        *nr = ONEUP2;
         else {    else {
                 d = UDEG(n); w = (d+1)/BSH + ((d+1)%BSH?1:0);      d = UDEG(n); w = (d+1)/BSH + ((d+1)%BSH?1:0);
                 up2_var = VR(n);      up2_var = VR(n);
                 W_NEWUP2(r,w);      W_NEWUP2(r,w);
                 for ( dc = DC(n); dc; dc = NEXT(dc) )      for ( dc = DC(n); dc; dc = NEXT(dc) )
                         if ( !EVENN(NM((Q)COEF(dc))) ) {        if ( !EVENN(NM((Q)COEF(dc))) ) {
                                 i = QTOS(DEG(dc));          i = QTOS(DEG(dc));
                                 r->b[i/BSH] |= 1<<(i%BSH);          r->b[i/BSH] |= 1<<(i%BSH);
                         }        }
                 for ( i = w-1; i >= 0 && !r->b[i]; i-- );      for ( i = w-1; i >= 0 && !r->b[i]; i-- );
                 r->w = i+1;      r->w = i+1;
                 NEWUP2(s,r->w); *nr = s;      NEWUP2(s,r->w); *nr = s;
                 _copyup2(r,s);      _copyup2(r,s);
         }    }
 }  }
   
 void ptoup2_sparse(n,nr)  void ptoup2_sparse(n,nr)
 P n;  P n;
 UP2 *nr;  UP2 *nr;
 {  {
         DCP dc;    DCP dc;
         UP2 s;    UP2 s;
         int d,i;    int d,i;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else if ( OID(n) == O_N )    else if ( OID(n) == O_N )
                 if ( EVENN(NM((Q)n)) )      if ( EVENN(NM((Q)n)) )
                         *nr = 0;        *nr = 0;
                 else {      else {
                         NEWUP2(s,1); s->w = 1; s->b[0] = 0;        NEWUP2(s,1); s->w = 1; s->b[0] = 0;
                         *nr = s;        *nr = s;
                 }      }
         else {    else {
                 d = UDEG(n);      d = UDEG(n);
                 for ( dc = DC(n), i = 0; dc; dc = NEXT(dc) )      for ( dc = DC(n), i = 0; dc; dc = NEXT(dc) )
                         if ( !EVENN(NM((Q)COEF(dc))) )        if ( !EVENN(NM((Q)COEF(dc))) )
                                 i++;          i++;
                 NEWUP2(s,i); s->w = i; *nr = s;      NEWUP2(s,i); s->w = i; *nr = s;
                 for ( dc = DC(n), i = 0; dc; dc = NEXT(dc) )      for ( dc = DC(n), i = 0; dc; dc = NEXT(dc) )
                         if ( !EVENN(NM((Q)COEF(dc))) )        if ( !EVENN(NM((Q)COEF(dc))) )
                                 s->b[i++] = QTOS(DEG(dc));          s->b[i++] = QTOS(DEG(dc));
         }    }
 }  }
   
 void up2top(n,nr)  void up2top(n,nr)
 UP2 n;  UP2 n;
 P *nr;  P *nr;
 {  {
         int i,d;    int i,d;
         DCP dc0,dc;    DCP dc0,dc;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else if ( !(d = degup2(n)) )    else if ( !(d = degup2(n)) )
                 *nr = (P)ONE;      *nr = (P)ONE;
         else {    else {
                 for ( i = d, dc0 = 0; i >= 0; i-- )      for ( i = d, dc0 = 0; i >= 0; i-- )
                         if ( n->b[i/BSH] & (1<<(i%BSH)) ) {        if ( n->b[i/BSH] & (1<<(i%BSH)) ) {
                                 NEXTDC(dc0,dc); STOQ(i,DEG(dc)); COEF(dc) = (P)ONE;          NEXTDC(dc0,dc); STOQ(i,DEG(dc)); COEF(dc) = (P)ONE;
                         }        }
                 if ( !up2_var )      if ( !up2_var )
                         up2_var = CO->v;        up2_var = CO->v;
                 MKP(up2_var,dc0,*nr);      MKP(up2_var,dc0,*nr);
         }    }
 }  }
   
 void up2tovect(n,nr)  void up2tovect(n,nr)
 UP2 n;  UP2 n;
 VECT *nr;  VECT *nr;
 {  {
         int i,d;    int i,d;
         VECT v;    VECT v;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 d = degup2(n);      d = degup2(n);
                 MKVECT(v,d+1); *nr = v;      MKVECT(v,d+1); *nr = v;
                 for ( i = d; i >= 0; i-- )      for ( i = d; i >= 0; i-- )
                         if ( n->b[i/BSH] & (1<<(i%BSH)) )        if ( n->b[i/BSH] & (1<<(i%BSH)) )
                                 v->body[i] = (pointer)ONE;          v->body[i] = (pointer)ONE;
         }    }
 }  }
   
 void up2ton(p,n)  void up2ton(p,n)
 UP2 p;  UP2 p;
 Q *n;  Q *n;
 {  {
         N nm;    N nm;
         int w;    int w;
   
         if ( !p )    if ( !p )
                 *n = 0;      *n = 0;
         else {    else {
                 w = p->w;      w = p->w;
                 nm = NALLOC(w);      nm = NALLOC(w);
                 nm->p = w;      nm->p = w;
                 bcopy(p->b,nm->b,w*sizeof(unsigned int));      bcopy(p->b,nm->b,w*sizeof(unsigned int));
                 NTOQ(nm,1,*n);      NTOQ(nm,1,*n);
         }    }
 }  }
   
 void ntoup2(n,p)  void ntoup2(n,p)
 Q n;  Q n;
 UP2 *p;  UP2 *p;
 {  {
         N nm;    N nm;
         UP2 t;    UP2 t;
         int w;    int w;
   
         if ( !n )    if ( !n )
                 *p = 0;      *p = 0;
         else {    else {
                 nm = NM(n); w = nm->p;      nm = NM(n); w = nm->p;
                 NEWUP2(t,w); *p = t; t->w = w;      NEWUP2(t,w); *p = t; t->w = w;
                 bcopy(nm->b,t->b,w*sizeof(unsigned int));      bcopy(nm->b,t->b,w*sizeof(unsigned int));
         }    }
 }  }
   
 void gen_simpup2(p,m,r)  void gen_simpup2(p,m,r)
Line 306  UP2 p;
Line 306  UP2 p;
 GEN_UP2 m;  GEN_UP2 m;
 UP2 *r;  UP2 *r;
 {  {
         if ( lm_lazy || !m )    if ( lm_lazy || !m )
                 *r = p;      *r = p;
         else if ( m->id == UP2_DENSE )    else if ( m->id == UP2_DENSE )
                 remup2(p,m->dense,r);      remup2(p,m->dense,r);
         else    else
                 remup2_sparse(p,m->sparse,r);      remup2_sparse(p,m->sparse,r);
         if ( *r && !((*r)->w) )    if ( *r && !((*r)->w) )
                 *r = 0;      *r = 0;
 }  }
   
 void gen_simpup2_destructive(p,m)  void gen_simpup2_destructive(p,m)
 UP2 p;  UP2 p;
 GEN_UP2 m;  GEN_UP2 m;
 {  {
         UP2 t;    UP2 t;
   
         if ( lm_lazy || !m )    if ( lm_lazy || !m )
                 return;      return;
         else if ( m->id == UP2_DENSE ) {    else if ( m->id == UP2_DENSE ) {
                 remup2(p,m->dense,&t);      remup2(p,m->dense,&t);
                 _copyup2(t,p);      _copyup2(t,p);
         } else    } else
                 remup2_sparse_destructive(p,m->sparse);      remup2_sparse_destructive(p,m->sparse);
 }  }
   
 void gen_invup2(p,m,r)  void gen_invup2(p,m,r)
Line 336  UP2 p;
Line 336  UP2 p;
 GEN_UP2 m;  GEN_UP2 m;
 UP2 *r;  UP2 *r;
 {  {
         if ( !m )    if ( !m )
                 error("gen_invup2 : invalid modulus");      error("gen_invup2 : invalid modulus");
         else    else
                 invup2(p,m->dense,r);      invup2(p,m->dense,r);
 }  }
   
 void gen_pwrmodup2(a,b,m,c)  void gen_pwrmodup2(a,b,m,c)
Line 348  Q b;
Line 348  Q b;
 GEN_UP2 m;  GEN_UP2 m;
 UP2 *c;  UP2 *c;
 {  {
         if ( !m )    if ( !m )
                 pwrmodup2(a,b,0,c);      pwrmodup2(a,b,0,c);
         else if ( m->id == UP2_DENSE )    else if ( m->id == UP2_DENSE )
                 pwrmodup2(a,b,m->dense,c);      pwrmodup2(a,b,m->dense,c);
         else    else
                 pwrmodup2_sparse(a,b,m->sparse,c);      pwrmodup2_sparse(a,b,m->sparse,c);
 }  }
   
 void simpup2(p,m,r)  void simpup2(p,m,r)
 UP2 p,m;  UP2 p,m;
 UP2 *r;  UP2 *r;
 {  {
         if ( !lm_lazy && m )    if ( !lm_lazy && m )
                 remup2(p,m,r);      remup2(p,m,r);
         else    else
                 *r = p;      *r = p;
 }  }
   
 int degup2(a)  int degup2(a)
 UP2 a;  UP2 a;
 {  {
         unsigned int l,i,t;    unsigned int l,i,t;
   
         if ( !a || !a->w )    if ( !a || !a->w )
                 return -1;      return -1;
         else {    else {
                 l = a->w; t = a->b[l-1];      l = a->w; t = a->b[l-1];
                 for ( i = 0; t; t>>=1, i++);      for ( i = 0; t; t>>=1, i++);
                 return i+(l-1)*BSH-1;      return i+(l-1)*BSH-1;
         }    }
 }  }
   
 int degup2_sparse(a)  int degup2_sparse(a)
 UP2 a;  UP2 a;
 {  {
         if ( !a || !a->w )    if ( !a || !a->w )
                 return -1;      return -1;
         else    else
                 return a->b[0];      return a->b[0];
 }  }
   
 int degup2_1(a)  int degup2_1(a)
 unsigned int a;  unsigned int a;
 {  {
         int i;    int i;
   
         for ( i = 0; a; a>>=1, i++);    for ( i = 0; a; a>>=1, i++);
         return i-1;    return i-1;
 }  }
   
 void addup2(a,b,c)  void addup2(a,b,c)
 UP2 a,b;  UP2 a,b;
 UP2 *c;  UP2 *c;
 {  {
         int i;    int i;
         UP2 t;    UP2 t;
         int w,wa,wb;    int w,wa,wb;
   
         if ( !a )    if ( !a )
                 *c = b;      *c = b;
         else if ( !b )    else if ( !b )
                 *c = a;      *c = a;
         else {    else {
                 wa = a->w; wb = b->w;      wa = a->w; wb = b->w;
                 if ( wa < wb ) {      if ( wa < wb ) {
                         t = a; a = b; b = t;        t = a; a = b; b = t;
                         w = wa; wa = wb; wb = w;        w = wa; wa = wb; wb = w;
                 }      }
                 NEWUP2(t,wa);      NEWUP2(t,wa);
                 for ( i = 0; i < wb; i++ )      for ( i = 0; i < wb; i++ )
                         t->b[i] = a->b[i]^b->b[i];        t->b[i] = a->b[i]^b->b[i];
                 for ( ; i < wa; i++ )      for ( ; i < wa; i++ )
                         t->b[i] = a->b[i];        t->b[i] = a->b[i];
                 for ( i = wa-1; i >= 0 && !t->b[i]; i-- );      for ( i = wa-1; i >= 0 && !t->b[i]; i-- );
                 if ( i < 0 )      if ( i < 0 )
                         *c = 0;        *c = 0;
                 else {      else {
                         t->w = i+1;        t->w = i+1;
                         *c = t;        *c = t;
                 }      }
         }    }
 }  }
   
 void subup2(a,b,c)  void subup2(a,b,c)
 UP2 a,b;  UP2 a,b;
 UP2 *c;  UP2 *c;
 {  {
         addup2(a,b,c);    addup2(a,b,c);
 }  }
   
   
Line 453  unsigned int *s,*r;
Line 453  unsigned int *s,*r;
 int w;  int w;
 unsigned int m;  unsigned int m;
 {  {
         int i;    int i;
         unsigned int _u,_l,t;    unsigned int _u,_l,t;
   
         for ( i = 0; i < w; i++ )    for ( i = 0; i < w; i++ )
                 if ( t = s[i] ) {      if ( t = s[i] ) {
                         GF2M_MUL_1(_u,_l,t,m);        GF2M_MUL_1(_u,_l,t,m);
                         r[i] ^= _l; r[i+1] ^= _u;        r[i] ^= _l; r[i+1] ^= _u;
                 }      }
 }  }
   
 INLINE void mulup2_nh(s,w,m,r)  INLINE void mulup2_nh(s,w,m,r)
Line 468  unsigned int *s,*r;
Line 468  unsigned int *s,*r;
 int w;  int w;
 unsigned int m; /* 0 <= b <= 0xffff */  unsigned int m; /* 0 <= b <= 0xffff */
 {  {
         int i;    int i;
         unsigned int u,l;    unsigned int u,l;
   
         for ( i = 0, r[0] = 0; i < w; i++ ) {    for ( i = 0, r[0] = 0; i < w; i++ ) {
                 GF2M_MUL_1_HALF(u,l,s[i],m);      GF2M_MUL_1_HALF(u,l,s[i],m);
                 r[i] ^= l; r[i+1] = u;      r[i] ^= l; r[i+1] = u;
         }    }
 }  }
   
 void _mulup2_1(a,b,c)  void _mulup2_1(a,b,c)
 UP2 a,c;  UP2 a,c;
 unsigned int b;  unsigned int b;
 {  {
         int w;    int w;
   
         if ( !a || !a->w || !b )    if ( !a || !a->w || !b )
                 c->w = 0;      c->w = 0;
         else if ( b <= 0xffff )    else if ( b <= 0xffff )
                 _mulup2_h(a,b,c);      _mulup2_h(a,b,c);
         else {    else {
                 w = a->w;      w = a->w;
                 bzero((char *)c->b,(w+1)*sizeof(unsigned int));      bzero((char *)c->b,(w+1)*sizeof(unsigned int));
                 mulup2_n1(a->b,w,b,c->b);      mulup2_n1(a->b,w,b,c->b);
                 c->w = c->b[w] ? w+1 : w;      c->w = c->b[w] ? w+1 : w;
         }    }
 }  }
   
 void _mulup2_h(a,b,c)  void _mulup2_h(a,b,c)
 UP2 a,c;  UP2 a,c;
 unsigned int b; /* 0 <= b <= 0xffff */  unsigned int b; /* 0 <= b <= 0xffff */
 {  {
         int w;    int w;
   
         if ( !a || !a->w || !b )    if ( !a || !a->w || !b )
                 c->w = 0;      c->w = 0;
         else {    else {
                 w = a->w;      w = a->w;
                 mulup2_nh(a->b,w,b,c->b);      mulup2_nh(a->b,w,b,c->b);
                 c->w = c->b[w] ? w+1 : w;      c->w = c->b[w] ? w+1 : w;
         }    }
 }  }
   
 void mulup2(a,b,c)  void mulup2(a,b,c)
 UP2 a,b;  UP2 a,b;
 UP2 *c;  UP2 *c;
 {  {
         UP2 t;    UP2 t;
         int wa,wb,w;    int wa,wb,w;
   
         if ( !a || !b )    if ( !a || !b )
                 *c = 0;      *c = 0;
         else if ( a->w==1 && a->b[0]==1 ) {    else if ( a->w==1 && a->b[0]==1 ) {
                 NEWUP2(t,b->w); *c = t; _copyup2(b,t);      NEWUP2(t,b->w); *c = t; _copyup2(b,t);
         } else if ( b->w==1 && b->b[0]==1 ) {    } else if ( b->w==1 && b->b[0]==1 ) {
                 NEWUP2(t,a->w); *c = t; _copyup2(a,t);      NEWUP2(t,a->w); *c = t; _copyup2(a,t);
         } else {    } else {
                 wa = a->w; wb = b->w;      wa = a->w; wb = b->w;
                 if ( wa != wb ) {      if ( wa != wb ) {
                         if ( wb > wa ) {        if ( wb > wa ) {
                                 t = a; a = b; b = t;          t = a; a = b; b = t;
                                 w = wa; wa = wb; wb = w;          w = wa; wa = wb; wb = w;
                         }        }
                         W_NEWUP2(t,wa);        W_NEWUP2(t,wa);
                         _copyup2(b,t);        _copyup2(b,t);
                         bzero((char *)(t->b+wb),(wa-wb)*sizeof(unsigned int));        bzero((char *)(t->b+wb),(wa-wb)*sizeof(unsigned int));
                         t->w = wa;        t->w = wa;
                         b = t;        b = t;
                 }      }
                 NEWUP2(t,2*wa);      NEWUP2(t,2*wa);
                 _kmulup2_(a->b,b->b,wa,t->b);      _kmulup2_(a->b,b->b,wa,t->b);
                 t->w = 2*wa; _adjup2(t);      t->w = 2*wa; _adjup2(t);
                 *c = t;      *c = t;
         }    }
 }  }
   
 void _kmulup2_(a,b,w,c)  void _kmulup2_(a,b,w,c)
 unsigned int *a,*b,*c;  unsigned int *a,*b,*c;
 int w;  int w;
 {  {
         switch ( w ) {    switch ( w ) {
                 case 1: GF2M_MUL_1(c[1],c[0],*a,*b); break;      case 1: GF2M_MUL_1(c[1],c[0],*a,*b); break;
                 case 2: _mulup2_22(a,b,c); break;      case 2: _mulup2_22(a,b,c); break;
                 case 3: _mulup2_33(a,b,c); break;      case 3: _mulup2_33(a,b,c); break;
                 case 4: _mulup2_44(a,b,c); break;      case 4: _mulup2_44(a,b,c); break;
                 case 5: _mulup2_55(a,b,c); break;      case 5: _mulup2_55(a,b,c); break;
                 case 6: _mulup2_66(a,b,c); break;      case 6: _mulup2_66(a,b,c); break;
                 default: _mulup2_nn(a,b,w,c); break;      default: _mulup2_nn(a,b,w,c); break;
         }    }
 }  }
   
 void _mulup2_nn(a,b,w,c)  void _mulup2_nn(a,b,w,c)
 unsigned int *a,*b,*c;  unsigned int *a,*b,*c;
 int w;  int w;
 {  {
         int wlow,whigh;    int wlow,whigh;
         struct _oUP2 ablow,abhigh,alow,ahigh,blow,bhigh,aa,bb,mid,cmid;    struct _oUP2 ablow,abhigh,alow,ahigh,blow,bhigh,aa,bb,mid,cmid;
   
         /* wlow >= whigh */    /* wlow >= whigh */
         wlow = (w+1)/2; whigh = w-wlow;    wlow = (w+1)/2; whigh = w-wlow;
   
         alow.w = wlow; alow.b = a;    alow.w = wlow; alow.b = a;
         ahigh.w = whigh; ahigh.b = a+wlow;    ahigh.w = whigh; ahigh.b = a+wlow;
   
         blow.w = wlow; blow.b = b;    blow.w = wlow; blow.b = b;
         bhigh.w = whigh; bhigh.b = b+wlow;    bhigh.w = whigh; bhigh.b = b+wlow;
   
         ablow.b = c;    ablow.b = c;
         abhigh.b = c+wlow*2;    abhigh.b = c+wlow*2;
   
         _kmulup2_(a,b,wlow,ablow.b); ablow.w = 2*wlow;    _kmulup2_(a,b,wlow,ablow.b); ablow.w = 2*wlow;
         _kmulup2_(a+wlow,b+wlow,whigh,abhigh.b); abhigh.w = 2*whigh;    _kmulup2_(a+wlow,b+wlow,whigh,abhigh.b); abhigh.w = 2*whigh;
   
         W_NEW_UP2(aa,wlow);    W_NEW_UP2(aa,wlow);
         _addup2_(&alow,&ahigh,&aa); aa.w = wlow;    _addup2_(&alow,&ahigh,&aa); aa.w = wlow;
   
         W_NEW_UP2(bb,wlow);    W_NEW_UP2(bb,wlow);
         _addup2_(&blow,&bhigh,&bb); bb.w = wlow;    _addup2_(&blow,&bhigh,&bb); bb.w = wlow;
   
         W_NEW_UP2(mid,2*wlow);    W_NEW_UP2(mid,2*wlow);
         _kmulup2_(aa.b,bb.b,wlow,mid.b); mid.w = 2*wlow;    _kmulup2_(aa.b,bb.b,wlow,mid.b); mid.w = 2*wlow;
   
         _addtoup2_(&ablow,&mid); _addtoup2_(&abhigh,&mid);    _addtoup2_(&ablow,&mid); _addtoup2_(&abhigh,&mid);
   
         cmid.w = 2*wlow; cmid.b = c+wlow;    cmid.w = 2*wlow; cmid.b = c+wlow;
         _addtoup2_(&mid,&cmid);    _addtoup2_(&mid,&cmid);
 }  }
   
 void _mulup2(a,b,c)  void _mulup2(a,b,c)
 UP2 a,b,c;  UP2 a,b,c;
 {  {
         int wa,wb,w;    int wa,wb,w;
         int i;    int i;
         unsigned int *ba,*bb;    unsigned int *ba,*bb;
         unsigned int mul;    unsigned int mul;
   
         if ( !a || !b || !a->w || !b->w ) {    if ( !a || !b || !a->w || !b->w ) {
                 c->w = 0; return;      c->w = 0; return;
         }    }
         wa = a->w; wb = b->w;    wa = a->w; wb = b->w;
         w = wa + wb;    w = wa + wb;
         bzero((char *)c->b,w*sizeof(unsigned int));    bzero((char *)c->b,w*sizeof(unsigned int));
         for ( i = 0, ba = a->b, bb = b->b; i < b->w; i++, bb++ )    for ( i = 0, ba = a->b, bb = b->b; i < b->w; i++, bb++ )
                 if ( mul = *bb )      if ( mul = *bb )
                         mulup2_n1(ba,wa,mul,c->b+i);        mulup2_n1(ba,wa,mul,c->b+i);
         c->w = w;    c->w = w;
         _adjup2(c);    _adjup2(c);
 }  }
   
 void _mulup2_(a,b,c)  void _mulup2_(a,b,c)
 _UP2 a,b,c;  _UP2 a,b,c;
 {  {
         int wa,wb,w;    int wa,wb,w;
         int i;    int i;
         unsigned int *ba,*bb;    unsigned int *ba,*bb;
         unsigned int mul;    unsigned int mul;
   
         if ( !a || !b || !a->w || !b->w ) {    if ( !a || !b || !a->w || !b->w ) {
                 c->w = 0; return;      c->w = 0; return;
         }    }
         wa = a->w; wb = b->w;    wa = a->w; wb = b->w;
         w = wa + wb;    w = wa + wb;
         bzero((char *)c->b,w*sizeof(unsigned int));    bzero((char *)c->b,w*sizeof(unsigned int));
         for ( i = 0, ba = a->b, bb = b->b; i < b->w; i++, bb++ )    for ( i = 0, ba = a->b, bb = b->b; i < b->w; i++, bb++ )
                 if ( mul = *bb )      if ( mul = *bb )
                         mulup2_n1(ba,wa,mul,c->b+i);        mulup2_n1(ba,wa,mul,c->b+i);
         c->w = w;    c->w = w;
         _adjup2_(c);    _adjup2_(c);
 }  }
   
 void squareup2(n,nr)  void squareup2(n,nr)
 UP2 n;  UP2 n;
 UP2 *nr;  UP2 *nr;
 {  {
         int w,w2,i;    int w,w2,i;
         unsigned int s;    unsigned int s;
         unsigned int *tb,*nb;    unsigned int *tb,*nb;
         UP2 t;    UP2 t;
   
         if ( !n )    if ( !n )
                 *nr = 0;      *nr = 0;
         else {    else {
                 w = n->w; w2 = 2*w;      w = n->w; w2 = 2*w;
                 NEWUP2(t,w2); *nr = t;      NEWUP2(t,w2); *nr = t;
                 tb = t->b;      tb = t->b;
                 nb = n->b;      nb = n->b;
                 for ( i = 0; i < w; i++ ) {      for ( i = 0; i < w; i++ ) {
                         s = nb[i];        s = nb[i];
                         tb[2*i] = up2_sqtab[s&0xff]|(up2_sqtab[(s>>8)&0xff]<<16);        tb[2*i] = up2_sqtab[s&0xff]|(up2_sqtab[(s>>8)&0xff]<<16);
                         tb[2*i+1] = up2_sqtab[(s>>16)&0xff]|(up2_sqtab[s>>24]<<16);        tb[2*i+1] = up2_sqtab[(s>>16)&0xff]|(up2_sqtab[s>>24]<<16);
                 }      }
                 t->w = tb[w2-1]?w2:w2-1;      t->w = tb[w2-1]?w2:w2-1;
         }    }
 }  }
   
 void _squareup2(n,nr)  void _squareup2(n,nr)
 UP2 n;  UP2 n;
 UP2 nr;  UP2 nr;
 {  {
         int w,w2,i;    int w,w2,i;
         unsigned int s;    unsigned int s;
         unsigned int *tb,*nb;    unsigned int *tb,*nb;
         UP2 t;    UP2 t;
   
         if ( !n )    if ( !n )
                 nr->w = 0;      nr->w = 0;
         else {    else {
                 w = n->w; w2 = 2*w;      w = n->w; w2 = 2*w;
                 t = nr;      t = nr;
                 tb = t->b;      tb = t->b;
                 nb = n->b;      nb = n->b;
                 for ( i = 0; i < w; i++ ) {      for ( i = 0; i < w; i++ ) {
                         s = nb[i];        s = nb[i];
                         tb[2*i] = up2_sqtab[s&0xff]|(up2_sqtab[(s>>8)&0xff]<<16);        tb[2*i] = up2_sqtab[s&0xff]|(up2_sqtab[(s>>8)&0xff]<<16);
                         tb[2*i+1] = up2_sqtab[(s>>16)&0xff]|(up2_sqtab[s>>24]<<16);        tb[2*i+1] = up2_sqtab[(s>>16)&0xff]|(up2_sqtab[s>>24]<<16);
                 }      }
                 t->w = tb[w2-1]?w2:w2-1;      t->w = tb[w2-1]?w2:w2-1;
         }    }
 }  }
   
 void _adjup2(n)  void _adjup2(n)
 UP2 n;  UP2 n;
 {  {
         int i;    int i;
         unsigned int *nb;    unsigned int *nb;
   
         nb = n->b;    nb = n->b;
         for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );    for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );
         i++;    i++;
         n->w = i;    n->w = i;
 }  }
   
 void _adjup2_(n)  void _adjup2_(n)
 _UP2 n;  _UP2 n;
 {  {
         int i;    int i;
         unsigned int *nb;    unsigned int *nb;
   
         nb = n->b;    nb = n->b;
         for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );    for ( i = n->w - 1; i >= 0 && !nb[i]; i-- );
         i++;    i++;
         n->w = i;    n->w = i;
 }  }
   
 void _addup2(a,b,c)  void _addup2(a,b,c)
 UP2 a,b,c;  UP2 a,b,c;
 {  {
         int i,wa,wb,w;    int i,wa,wb,w;
         UP2 t;    UP2 t;
         unsigned int *ab,*bb,*cb;    unsigned int *ab,*bb,*cb;
   
         if ( !a || !a->w ) {    if ( !a || !a->w ) {
                 _copyup2(b,c); return;      _copyup2(b,c); return;
         } else if ( !b || !b->w ) {    } else if ( !b || !b->w ) {
                 _copyup2(a,c); return;      _copyup2(a,c); return;
         }    }
         wa = a->w; wb = b->w;    wa = a->w; wb = b->w;
         if ( wa < wb ) {    if ( wa < wb ) {
                 t = a; a = b; b = t;      t = a; a = b; b = t;
                 w = wa; wa = wb; wb = w;      w = wa; wa = wb; wb = w;
         }    }
         ab = a->b; bb = b->b; cb = c->b;    ab = a->b; bb = b->b; cb = c->b;
         for ( i = 0; i < wb; i++ )    for ( i = 0; i < wb; i++ )
                 *cb++ = (*ab++)^(*bb++);      *cb++ = (*ab++)^(*bb++);
         for ( ; i < wa; i++ )    for ( ; i < wa; i++ )
                 *cb++ = *ab++;      *cb++ = *ab++;
         c->w = wa;    c->w = wa;
         _adjup2(c);    _adjup2(c);
 }  }
   
 /* a += b */  /* a += b */
Line 742  UP2 a,b,c;
Line 742  UP2 a,b,c;
 void _addup2_destructive(a,b)  void _addup2_destructive(a,b)
 UP2 a,b;  UP2 a,b;
 {  {
         int i,wa,wb;    int i,wa,wb;
         unsigned int *ab,*bb;    unsigned int *ab,*bb;
   
         if ( !a->w ) {    if ( !a->w ) {
                 _copyup2(b,a); return;      _copyup2(b,a); return;
         } else if ( !b->w )    } else if ( !b->w )
                 return;      return;
         else {    else {
                 wa = a->w; wb = b->w;      wa = a->w; wb = b->w;
                 ab = a->b; bb = b->b;      ab = a->b; bb = b->b;
                 if ( wa >= wb ) {      if ( wa >= wb ) {
                         for ( i = 0; i < wb; i++ )        for ( i = 0; i < wb; i++ )
                                 *ab++ ^= *bb++;          *ab++ ^= *bb++;
                 } else {      } else {
                         for ( i = 0; i < wa; i++ )        for ( i = 0; i < wa; i++ )
                                 *ab++ ^= *bb++;          *ab++ ^= *bb++;
                         for ( ; i < wb; i++ )        for ( ; i < wb; i++ )
                                 *ab++ = *bb++;          *ab++ = *bb++;
                         a->w = wb;        a->w = wb;
                 }      }
                 _adjup2(a);      _adjup2(a);
         }    }
 }  }
   
 void _addup2_(a,b,c)  void _addup2_(a,b,c)
 _UP2 a,b,c;  _UP2 a,b,c;
 {  {
         int i,wa,wb,w;    int i,wa,wb,w;
         _UP2 t;    _UP2 t;
         unsigned int *ab,*bb,*cb;    unsigned int *ab,*bb,*cb;
   
         wa = a->w; wb = b->w;    wa = a->w; wb = b->w;
         if ( wa < wb ) {    if ( wa < wb ) {
                 t = a; a = b; b = t;      t = a; a = b; b = t;
                 w = wa; wa = wb; wb = w;      w = wa; wa = wb; wb = w;
         }    }
         ab = a->b; bb = b->b; cb = c->b;    ab = a->b; bb = b->b; cb = c->b;
         for ( i = 0; i < wb; i++ )    for ( i = 0; i < wb; i++ )
                 *cb++ = (*ab++)^(*bb++);      *cb++ = (*ab++)^(*bb++);
         for ( ; i < wa; i++ )    for ( ; i < wa; i++ )
                 *cb++ = *ab++;      *cb++ = *ab++;
         c->w = wa;    c->w = wa;
 }  }
   
 void _addtoup2_(a,b)  void _addtoup2_(a,b)
 _UP2 a,b;  _UP2 a,b;
 {  {
         int i,wa;    int i,wa;
         unsigned int *ab,*bb;    unsigned int *ab,*bb;
   
         wa = a->w;    wa = a->w;
         ab = a->b; bb = b->b;    ab = a->b; bb = b->b;
         for ( i = 0; i < wa; i++ )    for ( i = 0; i < wa; i++ )
                 *bb++ ^= *ab++;      *bb++ ^= *ab++;
 }  }
   
 /* 8bit x 8bit; also works if deg(a*b) < 32 */  /* 8bit x 8bit; also works if deg(a*b) < 32 */
Line 803  _UP2 a,b;
Line 803  _UP2 a,b;
 unsigned int mulup2_bb(a,b)  unsigned int mulup2_bb(a,b)
 unsigned int a,b;  unsigned int a,b;
 {  {
         unsigned int t;    unsigned int t;
   
         if ( !a || !b )    if ( !a || !b )
                 return 0;      return 0;
         else {    else {
                 t = 0;      t = 0;
                 while ( b ) {      while ( b ) {
                         if ( b & 1 )        if ( b & 1 )
                                 t ^= a;          t ^= a;
                         a <<= 1;        a <<= 1;
                         b >>= 1;        b >>= 1;
                 }      }
                 return t;      return t;
         }    }
 }  }
   
 void init_up2_tab()  void init_up2_tab()
 {  {
         unsigned int i,j;    unsigned int i,j;
   
         for ( i = 0; i < 256; i++ )    for ( i = 0; i < 256; i++ )
                 for ( j = 0; j <= i; j++ ) {      for ( j = 0; j <= i; j++ ) {
                         up2_bbtab[i<<8|j] = mulup2_bb(i,j);        up2_bbtab[i<<8|j] = mulup2_bb(i,j);
                         up2_bbtab[j<<8|i] = up2_bbtab[i<<8|j];        up2_bbtab[j<<8|i] = up2_bbtab[i<<8|j];
                 }      }
         for ( i = 0; i < 256; i++ )    for ( i = 0; i < 256; i++ )
                 up2_sqtab[i] = mulup2_bb(i,i);      up2_sqtab[i] = mulup2_bb(i,i);
 }  }
   
 /*  /*
Line 840  void init_up2_tab()
Line 840  void init_up2_tab()
 INLINE unsigned int quoup2_11(a,b)  INLINE unsigned int quoup2_11(a,b)
 unsigned int a,b;  unsigned int a,b;
 {  {
         unsigned int q,i;    unsigned int q,i;
   
         q = 0;    q = 0;
         for ( i = ((unsigned int)1)<<(BSH-1); i; i >>=1, b >>= 1 )    for ( i = ((unsigned int)1)<<(BSH-1); i; i >>=1, b >>= 1 )
                 if ( i & a ) {      if ( i & a ) {
                         a ^= b;        a ^= b;
                         q |= i;        q |= i;
                 }      }
         return q;    return q;
 }  }
   
 void divup2_1(a1,a2,e1,e2,qp,rp)  void divup2_1(a1,a2,e1,e2,qp,rp)
Line 856  unsigned int a1,a2;
Line 856  unsigned int a1,a2;
 int e1,e2;  int e1,e2;
 unsigned int *qp,*rp;  unsigned int *qp,*rp;
 {  {
         int i;    int i;
         unsigned t,q;    unsigned t,q;
   
         for ( i=e1-e2, t=1<<e1, a2<<=i, q = 0; i>=0; i--, t>>=1, a2>>=1 ) {    for ( i=e1-e2, t=1<<e1, a2<<=i, q = 0; i>=0; i--, t>>=1, a2>>=1 ) {
                 q<<=1;      q<<=1;
                 if ( a1 & t ) {      if ( a1 & t ) {
                         q |= 1;        q |= 1;
                         a1 ^= a2;        a1 ^= a2;
                 }      }
         }    }
         *qp = q; *rp = a1;    *qp = q; *rp = a1;
 }  }
   
 void qrup2(a,b,q,r)  void qrup2(a,b,q,r)
 UP2 a,b;  UP2 a,b;
 UP2 *q,*r;  UP2 *q,*r;
 {  {
         unsigned int msa,msb,t,q0;    unsigned int msa,msb,t,q0;
         int s,i,wq,wb;    int s,i,wq,wb;
         UP2 as,bs,_q;    UP2 as,bs,_q;
   
         if ( !b )    if ( !b )
                 error("qrup2 : division by 0");      error("qrup2 : division by 0");
         else if ( !a ) {    else if ( !a ) {
                 *q = 0; *r = 0; return;      *q = 0; *r = 0; return;
         } else if ( degup2(a) < degup2(b) ) {    } else if ( degup2(a) < degup2(b) ) {
                 *q = 0; *r = a; return;      *q = 0; *r = a; return;
         }    }
         msb = b->b[b->w-1];    msb = b->b[b->w-1];
         for ( t = msb, s = 0; t; t>>=1, s++ );    for ( t = msb, s = 0; t; t>>=1, s++ );
         s = BSH-s;    s = BSH-s;
   
         W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);    W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);
         _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);    _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);
         as->b[as->w] = 0;    as->b[as->w] = 0;
   
         wb = bs->w;    wb = bs->w;
         wq = as->w-wb;    wq = as->w-wb;
         NEWUP2(_q,wq+1); *q = _q;    NEWUP2(_q,wq+1); *q = _q;
   
         msb = bs->b[bs->w-1];    msb = bs->b[bs->w-1];
         for ( i = wq; i >= 0; i-- ) {    for ( i = wq; i >= 0; i-- ) {
                 msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));      msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));
                 if ( msa ) {      if ( msa ) {
                         q0 = quoup2_11(msa,msb);        q0 = quoup2_11(msa,msb);
                         mulup2_n1(bs->b,wb,q0,as->b+i);        mulup2_n1(bs->b,wb,q0,as->b+i);
                 } else      } else
                         q0 = 0;        q0 = 0;
                 _q->b[i] = q0;      _q->b[i] = q0;
         }    }
         for ( i = wq; i >= 0 && !_q->b[i]; i-- );    for ( i = wq; i >= 0 && !_q->b[i]; i-- );
         _q->w = i+1;    _q->w = i+1;
   
         for ( i = wb-1; i >= 0 && !as->b[i]; i-- );    for ( i = wb-1; i >= 0 && !as->b[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 as->w = i+1;      as->w = i+1;
                 NEWUP2(*r,as->w);      NEWUP2(*r,as->w);
                 _bshiftup2(as,s,*r);      _bshiftup2(as,s,*r);
         }    }
 }  }
   
 /* q->w >= a->w-b->w+2, r->w >= b->w */  /* q->w >= a->w-b->w+2, r->w >= b->w */
Line 925  void _qrup2(a,b,q,r)
Line 925  void _qrup2(a,b,q,r)
 UP2 a,b;  UP2 a,b;
 UP2 q,r;  UP2 q,r;
 {  {
         unsigned int msa,msb,t,q0;    unsigned int msa,msb,t,q0;
         int s,i,wq,wb;    int s,i,wq,wb;
   
         if ( !b || !b->w )    if ( !b || !b->w )
                 error("_qrup2 : division by 0");      error("_qrup2 : division by 0");
         else if ( !a || !a->w ) {    else if ( !a || !a->w ) {
                 q->w = 0; r->w = 0; return;      q->w = 0; r->w = 0; return;
         } else if ( degup2(a) < degup2(b) ) {    } else if ( degup2(a) < degup2(b) ) {
                 q->w = 0; _copyup2(a,r); return;      q->w = 0; _copyup2(a,r); return;
         }    }
         msb = b->b[b->w-1];    msb = b->b[b->w-1];
         for ( t = msb, s = BSH; t; t>>=1, s-- );    for ( t = msb, s = BSH; t; t>>=1, s-- );
   
         _copyup2(a,r);    _copyup2(a,r);
         wb = b->w;    wb = b->w;
         wq = r->w-wb;    wq = r->w-wb;
         r->b[r->w] = 0;    r->b[r->w] = 0;
   
         msb = (b->b[b->w-1]<<s)|(b->w==1||!s?0:b->b[b->w-2]>>(BSH-s));    msb = (b->b[b->w-1]<<s)|(b->w==1||!s?0:b->b[b->w-2]>>(BSH-s));
         for ( i = wq; i >= 0; i-- ) {    for ( i = wq; i >= 0; i-- ) {
                 msa = (s==BSH-1?0:r->b[wb+i]<<(s+1))|(wb+i-1<0?0:r->b[wb+i-1]>>(BSH-1-s));      msa = (s==BSH-1?0:r->b[wb+i]<<(s+1))|(wb+i-1<0?0:r->b[wb+i-1]>>(BSH-1-s));
                 if ( msa ) {      if ( msa ) {
                         q0 = quoup2_11(msa,msb);        q0 = quoup2_11(msa,msb);
                         mulup2_n1(b->b,wb,q0,r->b+i);        mulup2_n1(b->b,wb,q0,r->b+i);
                 } else      } else
                         q0 = 0;        q0 = 0;
                 q->b[i] = q0;      q->b[i] = q0;
         }    }
         for ( i = wq; i >= 0 && !q->b[i]; i-- );    for ( i = wq; i >= 0 && !q->b[i]; i-- );
         q->w = i+1;    q->w = i+1;
   
         for ( i = wb-1; i >= 0 && !r->b[i]; i-- );    for ( i = wb-1; i >= 0 && !r->b[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 r->w = 0;      r->w = 0;
         else    else
                 r->w = i+1;      r->w = i+1;
 }  }
   
 void remup2(a,b,c)  void remup2(a,b,c)
 UP2 a,b;  UP2 a,b;
 UP2 *c;  UP2 *c;
 {  {
         unsigned int msa,msb,t,q;    unsigned int msa,msb,t,q;
         int s,i,wq,wb;    int s,i,wq,wb;
         UP2 as,bs,r;    UP2 as,bs,r;
   
         if ( !b || !b->w )    if ( !b || !b->w )
                 error("remup2 : division by 0");      error("remup2 : division by 0");
         else if ( !a || !a->w ) {    else if ( !a || !a->w ) {
                 *c = 0; return;      *c = 0; return;
         } else if ( degup2(a) < degup2(b) ) {    } else if ( degup2(a) < degup2(b) ) {
                 *c = a; return;      *c = a; return;
         }    }
         msb = b->b[b->w-1];    msb = b->b[b->w-1];
         for ( t = msb, s = 0; t; t>>=1, s++ );    for ( t = msb, s = 0; t; t>>=1, s++ );
         s = BSH-s;    s = BSH-s;
   
         W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);    W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);
         _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);    _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);
         as->b[as->w] = 0;    as->b[as->w] = 0;
   
         wb = bs->w;    wb = bs->w;
         wq = as->w-wb;    wq = as->w-wb;
         msb = bs->b[bs->w-1];    msb = bs->b[bs->w-1];
         for ( i = wq; i >= 0; i-- ) {    for ( i = wq; i >= 0; i-- ) {
                 msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));      msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));
                 if ( msa ) {      if ( msa ) {
                         q = quoup2_11(msa,msb);        q = quoup2_11(msa,msb);
                         mulup2_n1(bs->b,wb,q,as->b+i);        mulup2_n1(bs->b,wb,q,as->b+i);
                 }      }
         }    }
         for ( i = wb-1; i >= 0 && !as->b[i]; i-- );    for ( i = wb-1; i >= 0 && !as->b[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 *c = 0;      *c = 0;
         else {    else {
                 as->w = i+1;      as->w = i+1;
                 NEWUP2(r,as->w); *c = r;      NEWUP2(r,as->w); *c = r;
                 _bshiftup2(as,s,r);      _bshiftup2(as,s,r);
         }    }
 }  }
   
 void _remup2(a,b,c)  void _remup2(a,b,c)
 UP2 a,b,c;  UP2 a,b,c;
 {  {
         unsigned int msa,msb,t,q;    unsigned int msa,msb,t,q;
         int s,i,wq,wb;    int s,i,wq,wb;
         UP2 as,bs;    UP2 as,bs;
   
         if ( !b || !b->w )    if ( !b || !b->w )
                 error("_remup2 : division by 0");      error("_remup2 : division by 0");
         else if ( !a || !a->w ) {    else if ( !a || !a->w ) {
                 c->w = 0; return;      c->w = 0; return;
         } else if ( degup2(a) < degup2(b) ) {    } else if ( degup2(a) < degup2(b) ) {
                 _copyup2(a,c); return;      _copyup2(a,c); return;
         }    }
         msb = b->b[b->w-1];    msb = b->b[b->w-1];
         for ( t = msb, s = 0; t; t>>=1, s++ );    for ( t = msb, s = 0; t; t>>=1, s++ );
         s = BSH-s;    s = BSH-s;
   
         W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);    W_NEWUP2(as,a->w+2); W_NEWUP2(bs,b->w+1);
         _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);    _bshiftup2(a,-s,as); _bshiftup2(b,-s,bs);
         as->b[as->w] = 0;    as->b[as->w] = 0;
   
         wb = bs->w;    wb = bs->w;
         wq = as->w-wb;    wq = as->w-wb;
         msb = bs->b[bs->w-1];    msb = bs->b[bs->w-1];
         for ( i = wq; i >= 0; i-- ) {    for ( i = wq; i >= 0; i-- ) {
                 msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));      msa = (as->b[wb+i]<<1) | (as->b[wb+i-1]>>(BSH-1));
                 if ( msa ) {      if ( msa ) {
                         q = quoup2_11(msa,msb);        q = quoup2_11(msa,msb);
                         mulup2_n1(bs->b,wb,q,as->b+i);        mulup2_n1(bs->b,wb,q,as->b+i);
                 }      }
         }    }
         for ( i = wb-1; i >= 0 && !as->b[i]; i-- );    for ( i = wb-1; i >= 0 && !as->b[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 c->w = 0;      c->w = 0;
         else {    else {
                 as->w = i+1;      as->w = i+1;
                 _bshiftup2(as,s,c);      _bshiftup2(as,s,c);
         }    }
 }  }
   
 /* b = b->w|b->b[0]|b->b[1]|...  -> b = x^b->[0]+x^b->[1]+... (b->w terms) */  /* b = b->w|b->b[0]|b->b[1]|...  -> b = x^b->[0]+x^b->[1]+... (b->w terms) */
Line 1053  void remup2_sparse(a,b,c)
Line 1053  void remup2_sparse(a,b,c)
 UP2 a,b;  UP2 a,b;
 UP2 *c;  UP2 *c;
 {  {
         int i,j,k,wa,wb,d,ds,db,dr,r;    int i,j,k,wa,wb,d,ds,db,dr,r;
         unsigned int ha,hb;    unsigned int ha,hb;
         unsigned int *tb;    unsigned int *tb;
   
         UP2 t,s;    UP2 t,s;
   
         if ( !b )    if ( !b )
                 error("remup2_sparse : division by 0");      error("remup2_sparse : division by 0");
         else if ( !a ) {    else if ( !a ) {
                 *c = 0; return;      *c = 0; return;
         } else if ( degup2(a) < (db = degup2_sparse(b)) ) {    } else if ( degup2(a) < (db = degup2_sparse(b)) ) {
                 *c = a; return;      *c = a; return;
         } else if ( b->w == (int)(b->b[0])+1 ) {    } else if ( b->w == (int)(b->b[0])+1 ) {
                 NEWUP2(*c,a->w);      NEWUP2(*c,a->w);
                 _copyup2(a,*c);      _copyup2(a,*c);
                 remup2_type1_destructive(*c,b->b[0]);      remup2_type1_destructive(*c,b->b[0]);
                 if ( (*c)->w <= 0 )      if ( (*c)->w <= 0 )
                         *c = 0;        *c = 0;
         } else {    } else {
                 W_NEWUP2(t,a->w); _copyup2(a,t);      W_NEWUP2(t,a->w); _copyup2(a,t);
                 wa = a->w; wb = b->w; tb = t->b;      wa = a->w; wb = b->w; tb = t->b;
                 for ( i = wa-1; (ds = BSH*i-db) >= 0;  ) {      for ( i = wa-1; (ds = BSH*i-db) >= 0;  ) {
                         if ( !(ha = tb[i]) )        if ( !(ha = tb[i]) )
                                 i--;          i--;
                         else {        else {
                                 tb[i] = 0;          tb[i] = 0;
                                 for ( j = 1; j < wb; j++ ) {          for ( j = 1; j < wb; j++ ) {
                                         d = ds+b->b[j];            d = ds+b->b[j];
                                         k = d/BSH; r = d%BSH;            k = d/BSH; r = d%BSH;
                                         if ( !r )            if ( !r )
                                                 tb[k] ^= ha;              tb[k] ^= ha;
                                         else {            else {
                                                 tb[k] ^= (ha<<r);              tb[k] ^= (ha<<r);
                                                 if ( k+1 <= i )              if ( k+1 <= i )
                                                         tb[k+1] ^= (ha>>(BSH-r));                tb[k+1] ^= (ha>>(BSH-r));
                                         }            }
                                 }          }
                                 if ( !tb[i] )          if ( !tb[i] )
                                         i--;            i--;
                         }        }
                 }      }
                 dr = db%BSH;      dr = db%BSH;
                 hb = 1<<dr;      hb = 1<<dr;
                 i = db/BSH;      i = db/BSH;
                 while ( (ha=tb[i]) >= hb ) {      while ( (ha=tb[i]) >= hb ) {
                         ha >>= dr;        ha >>= dr;
                         t->b[i] &= (hb-1);        t->b[i] &= (hb-1);
                         for ( j = 1; j < wb; j++ ) {        for ( j = 1; j < wb; j++ ) {
                                 d = b->b[j];          d = b->b[j];
                                 k = d/BSH; r = d%BSH;          k = d/BSH; r = d%BSH;
                                 if ( !r )          if ( !r )
                                         tb[k] ^= ha;            tb[k] ^= ha;
                                 else {          else {
                                         tb[k] ^= (ha<<r);            tb[k] ^= (ha<<r);
                                         if ( k+1 <= i )            if ( k+1 <= i )
                                                 tb[k+1] ^= (ha>>(BSH-r));              tb[k+1] ^= (ha>>(BSH-r));
                                 }          }
                         }        }
                 }      }
                 t->w = i+1;      t->w = i+1;
                 _adjup2(t);      _adjup2(t);
                 if ( t->w == 0 )      if ( t->w == 0 )
                         *c = 0;        *c = 0;
                 else {      else {
                         NEWUP2(s,t->w); *c = s;        NEWUP2(s,t->w); *c = s;
                         _copyup2(t,s);        _copyup2(t,s);
                 }      }
         }    }
 }  }
   
 void remup2_sparse_destructive(a,b)  void remup2_sparse_destructive(a,b)
 UP2 a,b;  UP2 a,b;
 {  {
         int i,j,k,wb,d,ds,db,dr,r;    int i,j,k,wb,d,ds,db,dr,r;
         unsigned int ha,hb;    unsigned int ha,hb;
         unsigned int *ab,*bb;    unsigned int *ab,*bb;
   
         if ( !b )    if ( !b )
                 error("remup2_sparse_destructive : division by 0");      error("remup2_sparse_destructive : division by 0");
         else if ( !a || !a->w || ( degup2(a) < (db = degup2_sparse(b)) ) )    else if ( !a || !a->w || ( degup2(a) < (db = degup2_sparse(b)) ) )
                 return;      return;
         else if ( b->w == 3 )    else if ( b->w == 3 )
                 remup2_3_destructive(a,b);      remup2_3_destructive(a,b);
         else if ( b->w == 5 )    else if ( b->w == 5 )
                 remup2_5_destructive(a,b);      remup2_5_destructive(a,b);
         else if ( b->w == (int)(b->b[0])+1 )    else if ( b->w == (int)(b->b[0])+1 )
                 remup2_type1_destructive(a,b->b[0]);      remup2_type1_destructive(a,b->b[0]);
         else {    else {
                 wb = b->w; ab = a->b; bb = b->b;      wb = b->w; ab = a->b; bb = b->b;
                 for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {      for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {
                         if ( !(ha = ab[i]) )        if ( !(ha = ab[i]) )
                                 i--;          i--;
                         else {        else {
                                 ab[i] = 0;          ab[i] = 0;
                                 for ( j = 1; j < wb; j++ ) {          for ( j = 1; j < wb; j++ ) {
                                         d = ds+bb[j]; k = d>>5; r = d&31;            d = ds+bb[j]; k = d>>5; r = d&31;
                                         REMUP2_ONESTEP(ab,k,r,i,ha)            REMUP2_ONESTEP(ab,k,r,i,ha)
                                 }          }
                                 if ( !ab[i] )          if ( !ab[i] )
                                         i--;            i--;
                         }        }
                 }      }
                 i = db>>5; dr = db&31;      i = db>>5; dr = db&31;
                 for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {      for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
                         ha >>= dr;        ha >>= dr;
                         ab[i] &= hb;        ab[i] &= hb;
                         for ( j = 1; j < wb; j++ ) {        for ( j = 1; j < wb; j++ ) {
                                 d = bb[j]; k = d>>5; r = d&31;          d = bb[j]; k = d>>5; r = d&31;
                                 REMUP2_ONESTEP(ab,k,r,i,ha)          REMUP2_ONESTEP(ab,k,r,i,ha)
                         }        }
                 }      }
                 a->w = i+1;      a->w = i+1;
                 _adjup2(a);      _adjup2(a);
         }    }
 }  }
   
 /* b = x^d+x^(d-1)+...+1 */  /* b = x^d+x^(d-1)+...+1 */
Line 1175  void remup2_type1_destructive(a,d)
Line 1175  void remup2_type1_destructive(a,d)
 UP2 a;  UP2 a;
 int d;  int d;
 {  {
         int i,k,ds,db,dr;    int i,k,ds,db,dr;
         unsigned int ha,hb,r;    unsigned int ha,hb,r;
         unsigned int *ab;    unsigned int *ab;
   
         ab = a->b; db = d+1;    ab = a->b; db = d+1;
         for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {    for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {
                 if ( !(ha = ab[i]) ) i--;      if ( !(ha = ab[i]) ) i--;
                 else {      else {
                         ab[i] = 0;        ab[i] = 0;
                         k = ds>>5; r = ds&31;        k = ds>>5; r = ds&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         if ( !ab[i] ) i--;        if ( !ab[i] ) i--;
                 }      }
         }    }
         i = db>>5; dr = db&31;    i = db>>5; dr = db&31;
         for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {    for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
                 ha >>= dr; ab[i] &= hb; ab[0] ^= ha;      ha >>= dr; ab[i] &= hb; ab[0] ^= ha;
         }    }
         i = d>>5; hb = 1<<(d&31);    i = d>>5; hb = 1<<(d&31);
         if ( ab[i]&hb ) {    if ( ab[i]&hb ) {
                 ab[i] = (ab[i]^hb)^(hb-1);      ab[i] = (ab[i]^hb)^(hb-1);
                 for ( k = i-1; k >= 0; k-- ) ab[k] ^= 0xffffffff;      for ( k = i-1; k >= 0; k-- ) ab[k] ^= 0xffffffff;
         }    }
         a->w = i+1;    a->w = i+1;
         _adjup2(a);    _adjup2(a);
 }  }
   
 /* b = x^b->b[0]+x^b->b[1]+1 */  /* b = x^b->b[0]+x^b->b[1]+1 */
Line 1207  int d;
Line 1207  int d;
 void remup2_3_destructive(a,b)  void remup2_3_destructive(a,b)
 UP2 a,b;  UP2 a,b;
 {  {
         int i,k,d,ds,db,db1,dr;    int i,k,d,ds,db,db1,dr;
         unsigned int ha,hb,r;    unsigned int ha,hb,r;
         unsigned int *ab,*bb;    unsigned int *ab,*bb;
   
         ab = a->b; bb = b->b;    ab = a->b; bb = b->b;
         db = bb[0]; db1 = bb[1];    db = bb[0]; db1 = bb[1];
         for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {    for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {
                 if ( !(ha = ab[i]) )      if ( !(ha = ab[i]) )
                         i--;        i--;
                 else {      else {
                         ab[i] = 0;        ab[i] = 0;
                         d = ds+db1; k = d>>5; r = d&31;        d = ds+db1; k = d>>5; r = d&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         d = ds; k = d>>5; r = d&31;        d = ds; k = d>>5; r = d&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         if ( !ab[i] )        if ( !ab[i] )
                                 i--;          i--;
                 }      }
         }    }
         i = db>>5; dr = db&31;    i = db>>5; dr = db&31;
         k = db1>>5; r = db1&31;    k = db1>>5; r = db1&31;
         for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {    for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
                 ha >>= dr;      ha >>= dr;
                 ab[i] &= hb;      ab[i] &= hb;
                 REMUP2_ONESTEP(ab,k,r,i,ha)      REMUP2_ONESTEP(ab,k,r,i,ha)
                 ab[0] ^= ha;      ab[0] ^= ha;
         }    }
         a->w = i+1;    a->w = i+1;
         _adjup2(a);    _adjup2(a);
 }  }
   
 /* b = x^b->b[0]+x^b->b[1]+x^b->b[2]+x^b->b[3]+1 */  /* b = x^b->b[0]+x^b->b[1]+x^b->b[2]+x^b->b[3]+1 */
Line 1243  UP2 a,b;
Line 1243  UP2 a,b;
 void remup2_5_destructive(a,b)  void remup2_5_destructive(a,b)
 UP2 a,b;  UP2 a,b;
 {  {
         int i,d,ds,db,db1,db2,db3,dr;    int i,d,ds,db,db1,db2,db3,dr;
         int k,k1,k2,k3;    int k,k1,k2,k3;
         unsigned int ha,hb,r,r1,r2,r3;    unsigned int ha,hb,r,r1,r2,r3;
         unsigned int *ab,*bb;    unsigned int *ab,*bb;
   
         ab = a->b; bb = b->b;    ab = a->b; bb = b->b;
         db = bb[0]; db1 = bb[1]; db2 = bb[2]; db3 = bb[3];    db = bb[0]; db1 = bb[1]; db2 = bb[2]; db3 = bb[3];
         for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {    for ( i = a->w-1; (ds = (i<<5)-db) >= 0;  ) {
                 if ( !(ha = ab[i]) )      if ( !(ha = ab[i]) )
                         i--;        i--;
                 else {      else {
                         ab[i] = 0;        ab[i] = 0;
                         d = ds+db1; k = d>>5; r = d&31;        d = ds+db1; k = d>>5; r = d&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         d = ds+db2; k = d>>5; r = d&31;        d = ds+db2; k = d>>5; r = d&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         d = ds+db3; k = d>>5; r = d&31;        d = ds+db3; k = d>>5; r = d&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         d = ds; k = d>>5; r = d&31;        d = ds; k = d>>5; r = d&31;
                         REMUP2_ONESTEP(ab,k,r,i,ha)        REMUP2_ONESTEP(ab,k,r,i,ha)
                         if ( !ab[i] )        if ( !ab[i] )
                                 i--;          i--;
                 }      }
         }    }
         i = db>>5; dr = db&31;    i = db>>5; dr = db&31;
         k1 = db1>>5; r1 = db1&31;    k1 = db1>>5; r1 = db1&31;
         k2 = db2>>5; r2 = db2&31;    k2 = db2>>5; r2 = db2&31;
         k3 = db3>>5; r3 = db3&31;    k3 = db3>>5; r3 = db3&31;
         for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {    for ( hb = (1<<dr)-1; (ha=ab[i]) > hb; ) {
                 ha >>= dr;      ha >>= dr;
                 ab[i] &= hb;      ab[i] &= hb;
                 REMUP2_ONESTEP(ab,k1,r1,i,ha)      REMUP2_ONESTEP(ab,k1,r1,i,ha)
                 REMUP2_ONESTEP(ab,k2,r2,i,ha)      REMUP2_ONESTEP(ab,k2,r2,i,ha)
                 REMUP2_ONESTEP(ab,k3,r3,i,ha)      REMUP2_ONESTEP(ab,k3,r3,i,ha)
                 ab[0] ^= ha;      ab[0] ^= ha;
         }    }
         a->w = i+1;    a->w = i+1;
         _adjup2(a);    _adjup2(a);
 }  }
   
 void _invup2_1(f1,f2,a1,b1)  void _invup2_1(f1,f2,a1,b1)
 unsigned int f1,f2,*a1,*b1;  unsigned int f1,f2,*a1,*b1;
 {  {
         unsigned int p1,p2,p3,q1,q2,q3,g1,g2,q,r;    unsigned int p1,p2,p3,q1,q2,q3,g1,g2,q,r;
         int d1,d2;    int d1,d2;
   
         p1 = 1; p2 = 0;    p1 = 1; p2 = 0;
         q1 = 0; q2 = 1;    q1 = 0; q2 = 1;
         g1 = f1; g2 = f2;    g1 = f1; g2 = f2;
         d1 = degup2_1(g1); d2 = degup2_1(g2);    d1 = degup2_1(g1); d2 = degup2_1(g2);
         while ( g2 ) {    while ( g2 ) {
                 divup2_1(g1,g2,d1,d2,&q,&r);      divup2_1(g1,g2,d1,d2,&q,&r);
                 g1 = g2; g2 = r;      g1 = g2; g2 = r;
                 GCD_COEF(q,p1,p2,p3,q1,q2,q3)      GCD_COEF(q,p1,p2,p3,q1,q2,q3)
                 p1 = p2; p2 = p3;      p1 = p2; p2 = p3;
                 q1 = q2; q2 = q3;      q1 = q2; q2 = q3;
                 d1 = d2; d2 = degup2_1(g2);      d1 = d2; d2 = degup2_1(g2);
         }    }
         *a1 = p1; *b1 = q1;    *a1 = p1; *b1 = q1;
 }  }
   
 void _gcdup2_1(f1,f2,gcd)  void _gcdup2_1(f1,f2,gcd)
 unsigned int f1,f2,*gcd;  unsigned int f1,f2,*gcd;
 {  {
         unsigned int g1,g2,q,r;    unsigned int g1,g2,q,r;
         int d1,d2;    int d1,d2;
   
         if ( !f1 )    if ( !f1 )
                 *gcd = f2;      *gcd = f2;
         else if ( !f2 )    else if ( !f2 )
                 *gcd = f1;      *gcd = f1;
         else {    else {
                 g1 = f1; g2 = f2;      g1 = f1; g2 = f2;
                 d1 = degup2_1(g1); d2 = degup2_1(g2);      d1 = degup2_1(g1); d2 = degup2_1(g2);
                 while ( 1 ) {      while ( 1 ) {
                         divup2_1(g1,g2,d1,d2,&q,&r);        divup2_1(g1,g2,d1,d2,&q,&r);
                         if ( !r ) {        if ( !r ) {
                                 *gcd = g2;          *gcd = g2;
                                 return;          return;
                         }        }
                         g1 = g2; g2 = r;        g1 = g2; g2 = r;
                         d1 = d2; d2 = degup2_1(g2);        d1 = d2; d2 = degup2_1(g2);
                 }      }
         }    }
 }  }
   
 static struct oEGT eg_a,eg_b;  static struct oEGT eg_a,eg_b;
   
 void up2_init_eg() {  void up2_init_eg() {
         init_eg(&eg_a); init_eg(&eg_b);    init_eg(&eg_a); init_eg(&eg_b);
 }  }
   
 void up2_show_eg() {  void up2_show_eg() {
         print_eg("A",&eg_a);    print_eg("A",&eg_a);
         print_eg("B",&eg_b);    print_eg("B",&eg_b);
         printf("\n");    printf("\n");
 }  }
   
 void invup2(a,m,inv)  void invup2(a,m,inv)
 UP2 a,m;  UP2 a,m;
 UP2 *inv;  UP2 *inv;
 {  {
         int w,e1,e2,d1,d2;    int w,e1,e2,d1,d2;
         UP2 g1,g2,g3,a1,a2,a3,q,r,w1,w2,t;    UP2 g1,g2,g3,a1,a2,a3,q,r,w1,w2,t;
         unsigned int p1,p2,p3,q0,q1,q2,q3,h1,h2,h3;    unsigned int p1,p2,p3,q0,q1,q2,q3,h1,h2,h3;
   
         if ( degup2(a) >= degup2(m) ) {    if ( degup2(a) >= degup2(m) ) {
                 remup2(a,m,&t); a = t;      remup2(a,m,&t); a = t;
         }    }
         w = m->w+2;    w = m->w+2;
         W_NEWUP2(g1,w); W_NEWUP2(g2,w); W_NEWUP2(g3,w);    W_NEWUP2(g1,w); W_NEWUP2(g2,w); W_NEWUP2(g3,w);
         W_NEWUP2(q,w); W_NEWUP2(r,w);    W_NEWUP2(q,w); W_NEWUP2(r,w);
         W_NEWUP2(a1,w); W_NEWUP2(a2,w); W_NEWUP2(a3,w);    W_NEWUP2(a1,w); W_NEWUP2(a2,w); W_NEWUP2(a3,w);
         W_NEWUP2(w1,w); W_NEWUP2(w2,w);    W_NEWUP2(w1,w); W_NEWUP2(w2,w);
   
         a1->w = 1; a1->b[0] = 1;    a1->w = 1; a1->b[0] = 1;
         a2->w = 0;    a2->w = 0;
         _copyup2(a,g1); _copyup2(m,g2);    _copyup2(a,g1); _copyup2(m,g2);
   
         while ( 1 ) {    while ( 1 ) {
                 d1 = degup2(g1); d2 = degup2(g2);      d1 = degup2(g1); d2 = degup2(g2);
                 if ( d1 < d2 ) {      if ( d1 < d2 ) {
                         t = g1; g1 = g2; g2 = t;        t = g1; g1 = g2; g2 = t;
                         t = a1; a1 = a2; a2 = t;        t = a1; a1 = a2; a2 = t;
                 } else if ( d1 < 32 ) {      } else if ( d1 < 32 ) {
                         _invup2_1(g1->b[0],g2->b[0],&p1,&p2);        _invup2_1(g1->b[0],g2->b[0],&p1,&p2);
                         _mulup2_1(a1,p1,w1);        _mulup2_1(a1,p1,w1);
                         _mulup2_1(a2,p2,w2);        _mulup2_1(a2,p2,w2);
                         addup2(w1,w2,inv);        addup2(w1,w2,inv);
                         return;        return;
                 } else if ( d1-d2 >= 16 ) {      } else if ( d1-d2 >= 16 ) {
                         _qrup2(g1,g2,q,g3);        _qrup2(g1,g2,q,g3);
                         t = g1; g1 = g2;        t = g1; g1 = g2;
                         if ( !g3->w ) {        if ( !g3->w ) {
                                 NEWUP2(t,a2->w); *inv = t;          NEWUP2(t,a2->w); *inv = t;
                                 _copyup2(a2,t);          _copyup2(a2,t);
                                 return;          return;
                         } else {        } else {
                                 g2 = g3; g3 = t;          g2 = g3; g3 = t;
                         }        }
                         _mulup2(a2,q,w1); _addup2(a1,w1,a3);        _mulup2(a2,q,w1); _addup2(a1,w1,a3);
                         t = a1; a1 = a2; a2 = a3; a3 = t;        t = a1; a1 = a2; a2 = a3; a3 = t;
                 } else {      } else {
                         _bshiftup2(g1,d1-31,w1); h1 = w1->b[0];        _bshiftup2(g1,d1-31,w1); h1 = w1->b[0];
                         _bshiftup2(g2,d1-31,w2); h2 = w2->b[0];        _bshiftup2(g2,d1-31,w2); h2 = w2->b[0];
                         p1 = 1; p2 = 0;        p1 = 1; p2 = 0;
                         q1 = 0; q2 = 1;        q1 = 0; q2 = 1;
 /*                      get_eg(&eg0); */  /*      get_eg(&eg0); */
                         e1 = degup2_1(h1); /* e1 must be 31 */        e1 = degup2_1(h1); /* e1 must be 31 */
                         while ( 1 ) {        while ( 1 ) {
                                 e2 = degup2_1(h2);          e2 = degup2_1(h2);
                                 if ( e2 <= 15 )          if ( e2 <= 15 )
                                         break;            break;
                                 divup2_1(h1,h2,e1,e2,&q0,&h3);          divup2_1(h1,h2,e1,e2,&q0,&h3);
                                 GCD_COEF(q0,p1,p2,p3,q1,q2,q3)          GCD_COEF(q0,p1,p2,p3,q1,q2,q3)
                                 p1 = p2; p2 = p3;          p1 = p2; p2 = p3;
                                 q1 = q2; q2 = q3;          q1 = q2; q2 = q3;
                                 h1 = h2; h2 = h3;          h1 = h2; h2 = h3;
                                 e1 = e2;          e1 = e2;
                         }        }
 /*                      get_eg(&eg1); */  /*      get_eg(&eg1); */
                         _mulup2_1(g1,p1,w1); _mulup2_1(g2,q1,w2); _addup2(w1,w2,g3);        _mulup2_1(g1,p1,w1); _mulup2_1(g2,q1,w2); _addup2(w1,w2,g3);
                         _mulup2_1(g1,p2,w1); _mulup2_1(g2,q2,w2); _addup2(w1,w2,g2);        _mulup2_1(g1,p2,w1); _mulup2_1(g2,q2,w2); _addup2(w1,w2,g2);
                         t = g1; g1 = g3; g3 = t;        t = g1; g1 = g3; g3 = t;
                         _mulup2_1(a1,p1,w1); _mulup2_1(a2,q1,w2); _addup2(w1,w2,a3);        _mulup2_1(a1,p1,w1); _mulup2_1(a2,q1,w2); _addup2(w1,w2,a3);
                         _mulup2_1(a1,p2,w1); _mulup2_1(a2,q2,w2); _addup2(w1,w2,a2);        _mulup2_1(a1,p2,w1); _mulup2_1(a2,q2,w2); _addup2(w1,w2,a2);
                         t = a1; a1 = a3; a3 = t;        t = a1; a1 = a3; a3 = t;
 /*                      get_eg(&eg2); */  /*      get_eg(&eg2); */
 /*                      add_eg(&eg_a,&eg0,&eg1); */  /*      add_eg(&eg_a,&eg0,&eg1); */
 /*                      add_eg(&eg_b,&eg1,&eg2); */  /*      add_eg(&eg_b,&eg1,&eg2); */
                 }      }
         }    }
 }  }
   
 void gcdup2(a,m,gcd)  void gcdup2(a,m,gcd)
 UP2 a,m;  UP2 a,m;
 UP2 *gcd;  UP2 *gcd;
 {  {
         int w,e1,e2,d1,d2;    int w,e1,e2,d1,d2;
         UP2 g1,g2,g3,q,r,w1,w2,t;    UP2 g1,g2,g3,q,r,w1,w2,t;
         unsigned int p1,p2,p3,q0,q1,q2,q3,h1,h2,h3;    unsigned int p1,p2,p3,q0,q1,q2,q3,h1,h2,h3;
   
         if ( !a ) {    if ( !a ) {
                 *gcd = m; return;      *gcd = m; return;
         } else if ( !m ) {    } else if ( !m ) {
                 *gcd = a; return;      *gcd = a; return;
         }    }
         if ( degup2(a) >= degup2(m) ) {    if ( degup2(a) >= degup2(m) ) {
                 remup2(a,m,&t); a = t;      remup2(a,m,&t); a = t;
                 if ( !a ) {      if ( !a ) {
                         *gcd = m; return;        *gcd = m; return;
                 }      }
         }    }
         w = m->w+2;    w = m->w+2;
         W_NEWUP2(g1,w); W_NEWUP2(g2,w); W_NEWUP2(g3,w);    W_NEWUP2(g1,w); W_NEWUP2(g2,w); W_NEWUP2(g3,w);
         W_NEWUP2(q,w); W_NEWUP2(r,w);    W_NEWUP2(q,w); W_NEWUP2(r,w);
         W_NEWUP2(w1,w); W_NEWUP2(w2,w);    W_NEWUP2(w1,w); W_NEWUP2(w2,w);
   
         _copyup2(a,g1); _copyup2(m,g2);    _copyup2(a,g1); _copyup2(m,g2);
   
         while ( 1 ) {    while ( 1 ) {
                 d1 = degup2(g1); d2 = degup2(g2);      d1 = degup2(g1); d2 = degup2(g2);
                 if ( d1 < d2 ) {      if ( d1 < d2 ) {
                         t = g1; g1 = g2; g2 = t;        t = g1; g1 = g2; g2 = t;
                 } else if ( d2 < 0 ) {      } else if ( d2 < 0 ) {
                         NEWUP2(t,g1->w); *gcd = t;        NEWUP2(t,g1->w); *gcd = t;
                         _copyup2(g1,t);        _copyup2(g1,t);
                         return;        return;
                 } else if ( d1 < 32 ) {      } else if ( d1 < 32 ) {
                         NEWUP2(t,1); t->w = 1; *gcd = t;        NEWUP2(t,1); t->w = 1; *gcd = t;
                         _gcdup2_1(g1->b[0],g2->b[0],&t->b[0]);        _gcdup2_1(g1->b[0],g2->b[0],&t->b[0]);
                         return;        return;
                 } else if ( d1-d2 >= 16 ) {      } else if ( d1-d2 >= 16 ) {
                         _qrup2(g1,g2,q,g3);        _qrup2(g1,g2,q,g3);
                         t = g1; g1 = g2;        t = g1; g1 = g2;
                         if ( !g3->w ) {        if ( !g3->w ) {
                                 NEWUP2(t,g2->w); *gcd = t;          NEWUP2(t,g2->w); *gcd = t;
                                 _copyup2(g2,t);          _copyup2(g2,t);
                                 return;          return;
                         } else {        } else {
                                 g2 = g3; g3 = t;          g2 = g3; g3 = t;
                         }        }
                 } else {      } else {
                         _bshiftup2(g1,d1-31,w1); h1 = w1->b[0];        _bshiftup2(g1,d1-31,w1); h1 = w1->b[0];
                         _bshiftup2(g2,d1-31,w2); h2 = w2->b[0];        _bshiftup2(g2,d1-31,w2); h2 = w2->b[0];
                         p1 = 1; p2 = 0;        p1 = 1; p2 = 0;
                         q1 = 0; q2 = 1;        q1 = 0; q2 = 1;
                         e1 = degup2_1(h1); /* e1 must be 31 */        e1 = degup2_1(h1); /* e1 must be 31 */
                         while ( 1 ) {        while ( 1 ) {
                                 e2 = degup2_1(h2);          e2 = degup2_1(h2);
                                 if ( e2 <= 15 )          if ( e2 <= 15 )
                                         break;            break;
                                 divup2_1(h1,h2,e1,e2,&q0,&h3);          divup2_1(h1,h2,e1,e2,&q0,&h3);
                                 GCD_COEF(q0,p1,p2,p3,q1,q2,q3)          GCD_COEF(q0,p1,p2,p3,q1,q2,q3)
                                 p1 = p2; p2 = p3;          p1 = p2; p2 = p3;
                                 q1 = q2; q2 = q3;          q1 = q2; q2 = q3;
                                 h1 = h2; h2 = h3;          h1 = h2; h2 = h3;
                                 e1 = e2;          e1 = e2;
                         }        }
                         _mulup2_h(g1,p1,w1); _mulup2_h(g2,q1,w2); _addup2(w1,w2,g3);        _mulup2_h(g1,p1,w1); _mulup2_h(g2,q1,w2); _addup2(w1,w2,g3);
                         _mulup2_h(g1,p2,w1); _mulup2_h(g2,q2,w2); _addup2(w1,w2,g2);        _mulup2_h(g1,p2,w1); _mulup2_h(g2,q2,w2); _addup2(w1,w2,g2);
                         t = g1; g1 = g3; g3 = t;        t = g1; g1 = g3; g3 = t;
                 }      }
         }    }
 }  }
   
 void chsgnup2(a,c)  void chsgnup2(a,c)
 UP2 a,*c;  UP2 a,*c;
 {  {
         *c = a;    *c = a;
 }  }
   
 void pwrmodup2(a,b,m,c)  void pwrmodup2(a,b,m,c)
Line 1501  Q b;
Line 1501  Q b;
 UP2 m;  UP2 m;
 UP2 *c;  UP2 *c;
 {  {
         N n;    N n;
         UP2 y,t,t1;    UP2 y,t,t1;
         int k;    int k;
   
         if ( !b )    if ( !b )
                 *c = (UP2)ONEUP2;      *c = (UP2)ONEUP2;
         else if ( !a )    else if ( !a )
                 *c = 0;      *c = 0;
         else {    else {
                 y = (UP2)ONEUP2;      y = (UP2)ONEUP2;
                 n = NM(b);      n = NM(b);
                 if ( !m )      if ( !m )
                         for ( k = n_bits(n)-1; k >= 0; k-- ) {        for ( k = n_bits(n)-1; k >= 0; k-- ) {
                                 squareup2(y,&t);          squareup2(y,&t);
                                 if ( n->b[k/BSH] & (1<<(k%BSH)) )          if ( n->b[k/BSH] & (1<<(k%BSH)) )
                                         mulup2(t,a,&y);            mulup2(t,a,&y);
                                 else          else
                                         y = t;            y = t;
                         }        }
                 else      else
                         for ( k = n_bits(n)-1; k >= 0; k-- ) {        for ( k = n_bits(n)-1; k >= 0; k-- ) {
                                 squareup2(y,&t);          squareup2(y,&t);
                                 remup2(t,m,&t1); t = t1;          remup2(t,m,&t1); t = t1;
                                 if ( n->b[k/BSH] & (1<<(k%BSH)) ) {          if ( n->b[k/BSH] & (1<<(k%BSH)) ) {
                                         mulup2(t,a,&y);            mulup2(t,a,&y);
                                         remup2(y,m,&t1); y = t1;            remup2(y,m,&t1); y = t1;
                                 }          }
                                 else          else
                                         y = t;            y = t;
                         }        }
                 *c = y;      *c = y;
         }    }
 }  }
   
 void pwrmodup2_sparse(a,b,m,c)  void pwrmodup2_sparse(a,b,m,c)
Line 1541  Q b;
Line 1541  Q b;
 UP2 m;  UP2 m;
 UP2 *c;  UP2 *c;
 {  {
         N n;    N n;
         UP2 y,t,t1;    UP2 y,t,t1;
         int k;    int k;
   
         if ( !b )    if ( !b )
                 *c = (UP2)ONEUP2;      *c = (UP2)ONEUP2;
         else if ( !a )    else if ( !a )
                 *c = 0;      *c = 0;
         else {    else {
                 y = (UP2)ONEUP2;      y = (UP2)ONEUP2;
                 n = NM(b);      n = NM(b);
                 if ( !m )      if ( !m )
                         for ( k = n_bits(n)-1; k >= 0; k-- ) {        for ( k = n_bits(n)-1; k >= 0; k-- ) {
                                 squareup2(y,&t);          squareup2(y,&t);
                                 if ( n->b[k/BSH] & (1<<(k%BSH)) )          if ( n->b[k/BSH] & (1<<(k%BSH)) )
                                         mulup2(t,a,&y);            mulup2(t,a,&y);
                                 else          else
                                         y = t;            y = t;
                         }        }
                 else      else
                         for ( k = n_bits(n)-1; k >= 0; k-- ) {        for ( k = n_bits(n)-1; k >= 0; k-- ) {
                                 squareup2(y,&t);          squareup2(y,&t);
                                 remup2_sparse(t,m,&t1); t = t1;          remup2_sparse(t,m,&t1); t = t1;
                                 if ( n->b[k/BSH] & (1<<(k%BSH)) ) {          if ( n->b[k/BSH] & (1<<(k%BSH)) ) {
                                         mulup2(t,a,&y);            mulup2(t,a,&y);
                                         remup2_sparse(y,m,&t1); y = t1;            remup2_sparse(y,m,&t1); y = t1;
                                 }          }
                                 else          else
                                         y = t;            y = t;
                         }        }
                 *c = y;      *c = y;
         }    }
 }  }
   
 int compup2(n1,n2)  int compup2(n1,n2)
 UP2 n1,n2;  UP2 n1,n2;
 {  {
         int i;    int i;
         unsigned int *m1,*m2;    unsigned int *m1,*m2;
   
         if ( !n1 )    if ( !n1 )
                 if ( !n2 )      if ( !n2 )
                         return 0;        return 0;
                 else      else
                         return -1;        return -1;
         else if ( !n2 )    else if ( !n2 )
                 return 1;      return 1;
         else if ( n1->w > n2->w )    else if ( n1->w > n2->w )
                 return 1;      return 1;
         else if ( n1->w < n2->w )    else if ( n1->w < n2->w )
                 return -1;      return -1;
         else {    else {
                 for ( i = n1->w-1, m1 = n1->b+i, m2 = n2->b+i;      for ( i = n1->w-1, m1 = n1->b+i, m2 = n2->b+i;
                         i >= 0; i--, m1--, m2-- )        i >= 0; i--, m1--, m2-- )
                         if ( *m1 > *m2 )        if ( *m1 > *m2 )
                                 return 1;          return 1;
                         else if ( *m1 < *m2 )        else if ( *m1 < *m2 )
                                 return -1;          return -1;
                 return 0;      return 0;
         }    }
 }  }
   
 void _copyup2(n,r)  void _copyup2(n,r)
 UP2 n,r;  UP2 n,r;
 {  {
         r->w = n->w;    r->w = n->w;
         bcopy(n->b,r->b,n->w*sizeof(unsigned int));    bcopy(n->b,r->b,n->w*sizeof(unsigned int));
 }  }
   
 void _bshiftup2(n,b,r)  void _bshiftup2(n,b,r)
Line 1615  UP2 n;
Line 1615  UP2 n;
 int b;  int b;
 UP2 r;  UP2 r;
 {  {
         int w,l,nl,i,j;    int w,l,nl,i,j;
         unsigned int msw;    unsigned int msw;
         unsigned int *p,*pr;    unsigned int *p,*pr;
   
         if ( b == 0 ) {    if ( b == 0 ) {
                 _copyup2(n,r); return;      _copyup2(n,r); return;
         }    }
         if ( b > 0 ) { /* >> */    if ( b > 0 ) { /* >> */
                 w = b / BSH; l = n->w-w;      w = b / BSH; l = n->w-w;
                 if ( l <= 0 ) {      if ( l <= 0 ) {
                         r->w = 0; return;        r->w = 0; return;
                 }      }
                 b %= BSH; p = n->b+w;      b %= BSH; p = n->b+w;
                 if ( !b ) {      if ( !b ) {
                         r->w = l;        r->w = l;
                         bcopy(p,r->b,l*sizeof(unsigned int));        bcopy(p,r->b,l*sizeof(unsigned int));
                         return;        return;
                 }      }
                 msw = p[l-1];      msw = p[l-1];
                 for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;      for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
                 if ( b >= i ) {      if ( b >= i ) {
                         l--;        l--;
                         if ( !l ) {        if ( !l ) {
                                 r->w = 0; return;          r->w = 0; return;
                         }        }
                         r->w = l; pr = r->b;        r->w = l; pr = r->b;
                         for ( j = 0; j < l; j++, p++ )        for ( j = 0; j < l; j++, p++ )
                                 *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);          *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
                 } else {      } else {
                         r->w = l; pr = r->b;        r->w = l; pr = r->b;
                         for ( j = 1; j < l; j++, p++ )        for ( j = 1; j < l; j++, p++ )
                                 *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);          *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
                         *pr = *p>>b;        *pr = *p>>b;
                 }      }
         } else { /* << */    } else { /* << */
                 b = -b;      b = -b;
                 w = b / BSH; b %= BSH; l = n->w; p = n->b;      w = b / BSH; b %= BSH; l = n->w; p = n->b;
                 if ( !b ) {      if ( !b ) {
                         nl = l+w; r->w = nl;        nl = l+w; r->w = nl;
                         bzero((char *)r->b,w*sizeof(unsigned int));        bzero((char *)r->b,w*sizeof(unsigned int));
                         bcopy(p,r->b+w,l*sizeof(unsigned int));        bcopy(p,r->b+w,l*sizeof(unsigned int));
                         return;        return;
                 }      }
                 msw = p[l-1];      msw = p[l-1];
                 for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;      for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
                 if ( b + i > BSH ) {      if ( b + i > BSH ) {
                         nl = l+w+1;        nl = l+w+1;
                         r->w = nl; pr = r->b+w;        r->w = nl; pr = r->b+w;
                         bzero((char *)r->b,w*sizeof(unsigned int));        bzero((char *)r->b,w*sizeof(unsigned int));
                         *pr++ = *p++<<b;        *pr++ = *p++<<b;
                         for ( j = 1; j < l; j++, p++ )        for ( j = 1; j < l; j++, p++ )
                                 *pr++ = (*p<<b)|(*(p-1)>>(BSH-b));          *pr++ = (*p<<b)|(*(p-1)>>(BSH-b));
                         *pr = *(p-1)>>(BSH-b);        *pr = *(p-1)>>(BSH-b);
                 } else {      } else {
                         nl = l+w;        nl = l+w;
                         r->w = nl; pr = r->b+w;        r->w = nl; pr = r->b+w;
                         bzero((char *)r->b,w*sizeof(unsigned int));        bzero((char *)r->b,w*sizeof(unsigned int));
                         *pr++ = *p++<<b;        *pr++ = *p++<<b;
                         for ( j = 1; j < l; j++, p++ )        for ( j = 1; j < l; j++, p++ )
                                 *pr++ = (*p<<b)|(*(p-1)>>(BSH-b));          *pr++ = (*p<<b)|(*(p-1)>>(BSH-b));
                 }      }
         }    }
 }  }
   
 void _bshiftup2_destructive(n,b)  void _bshiftup2_destructive(n,b)
 UP2 n;  UP2 n;
 int b;  int b;
 {  {
         int w,l,nl,i,j;    int w,l,nl,i,j;
         unsigned int msw;    unsigned int msw;
         unsigned int *p,*pr;    unsigned int *p,*pr;
   
         if ( b == 0 || !n->w )    if ( b == 0 || !n->w )
                 return;      return;
         if ( b > 0 ) { /* >> */    if ( b > 0 ) { /* >> */
                 w = b / BSH; l = n->w-w;      w = b / BSH; l = n->w-w;
                 if ( l <= 0 ) {      if ( l <= 0 ) {
                         n->w = 0; return;        n->w = 0; return;
                 }      }
                 b %= BSH; p = n->b+w;      b %= BSH; p = n->b+w;
                 if ( !b ) {      if ( !b ) {
                         n->w = l;        n->w = l;
                         bcopy(p,n->b,l*sizeof(unsigned int));        bcopy(p,n->b,l*sizeof(unsigned int));
                         return;        return;
                 }      }
                 msw = p[l-1];      msw = p[l-1];
                 for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;      for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
                 if ( b >= i ) {      if ( b >= i ) {
                         l--;        l--;
                         if ( !l ) {        if ( !l ) {
                                 n->w = 0; return;          n->w = 0; return;
                         }        }
                         n->w = l; pr = n->b;        n->w = l; pr = n->b;
                         for ( j = 0; j < l; j++, p++ )        for ( j = 0; j < l; j++, p++ )
                                 *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);          *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
                 } else {      } else {
                         n->w = l; pr = n->b;        n->w = l; pr = n->b;
                         for ( j = 1; j < l; j++, p++ )        for ( j = 1; j < l; j++, p++ )
                                 *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);          *pr++ = (*(p+1)<<(BSH-b))|(*p>>b);
                         *pr = *p>>b;        *pr = *p>>b;
                 }      }
         } else { /* << */    } else { /* << */
                 b = -b;      b = -b;
                 w = b / BSH; b %= BSH; l = n->w; p = n->b;      w = b / BSH; b %= BSH; l = n->w; p = n->b;
                 if ( !b ) {      if ( !b ) {
                         nl = l+w; n->w = nl;        nl = l+w; n->w = nl;
                         bcopy(p,n->b+w,l*sizeof(unsigned int));        bcopy(p,n->b+w,l*sizeof(unsigned int));
                         bzero((char *)n->b,w*sizeof(unsigned int));        bzero((char *)n->b,w*sizeof(unsigned int));
                         return;        return;
                 }      }
                 msw = p[l-1];      msw = p[l-1];
                 for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;      for ( i = BSH-1; !(msw&(((unsigned int)1)<<i)); i-- ); i++;
                 if ( b + i > BSH ) {      if ( b + i > BSH ) {
                         nl = l+w+1;        nl = l+w+1;
                         n->w = nl; p = n->b+l-1; pr = n->b+w+l;        n->w = nl; p = n->b+l-1; pr = n->b+w+l;
                         *pr-- = *p>>(BSH-b);        *pr-- = *p>>(BSH-b);
                         for ( j = l-1; j >= 1; j--, p-- )        for ( j = l-1; j >= 1; j--, p-- )
                                 *pr-- = (*p<<b)|(*(p-1)>>(BSH-b));          *pr-- = (*p<<b)|(*(p-1)>>(BSH-b));
                         *pr = *p<<b;        *pr = *p<<b;
                         bzero((char *)n->b,w*sizeof(unsigned int));        bzero((char *)n->b,w*sizeof(unsigned int));
                 } else {      } else {
                         nl = l+w;        nl = l+w;
                         n->w = nl; p = n->b+l-1; pr = n->b+w+l-1;        n->w = nl; p = n->b+l-1; pr = n->b+w+l-1;
                         for ( j = l-1; j >= 1; j--, p-- )        for ( j = l-1; j >= 1; j--, p-- )
                                 *pr-- = (*p<<b)|(*(p-1)>>(BSH-b));          *pr-- = (*p<<b)|(*(p-1)>>(BSH-b));
                         *pr = *p<<b;        *pr = *p<<b;
                         bzero((char *)n->b,w*sizeof(unsigned int));        bzero((char *)n->b,w*sizeof(unsigned int));
                 }      }
         }    }
 }  }
   
 void diffup2(f,r)  void diffup2(f,r)
 UP2 f;  UP2 f;
 UP2 *r;  UP2 *r;
 {  {
         int d,i,w;    int d,i,w;
         UP2 t;    UP2 t;
   
         d = degup2(f);    d = degup2(f);
         w = f->w;    w = f->w;
         NEWUP2(t,w);    NEWUP2(t,w);
         _bshiftup2(f,1,t);    _bshiftup2(f,1,t);
         for ( i = 0; i < d; i++ )    for ( i = 0; i < d; i++ )
                 if ( i%2 )      if ( i%2 )
                         t->b[i/BSH] &= ~(1<<(i%BSH));        t->b[i/BSH] &= ~(1<<(i%BSH));
         for ( i = w-1; i >= 0 && !t->b[i]; i-- );    for ( i = w-1; i >= 0 && !t->b[i]; i-- );
         if ( i < 0 )    if ( i < 0 )
                 *r = 0;      *r = 0;
         else {    else {
                 *r = t;      *r = t;
                 t->w = i+1;      t->w = i+1;
         }    }
 }  }
   
 int sqfrcheckup2(f)  int sqfrcheckup2(f)
 UP2 f;  UP2 f;
 {  {
         UP2 df,g;    UP2 df,g;
   
         diffup2(f,&df);    diffup2(f,&df);
         gcdup2(f,df,&g);    gcdup2(f,df,&g);
         if ( degup2(g) >= 1 )    if ( degup2(g) >= 1 )
                 return 0;      return 0;
         else    else
                 return 1;      return 1;
 }  }
   
 int irredcheckup2(f)  int irredcheckup2(f)
 UP2 f;  UP2 f;
 {  {
         int n,w,i,j,k,hcol;    int n,w,i,j,k,hcol;
         unsigned int hbit;    unsigned int hbit;
         unsigned int *u;    unsigned int *u;
         unsigned int **mat;    unsigned int **mat;
         UP2 p,t;    UP2 p,t;
   
         if ( !sqfrcheckup2(f) )    if ( !sqfrcheckup2(f) )
                 return 0;      return 0;
         n = degup2(f);    n = degup2(f);
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         mat = (unsigned int **)almat(n,w);    mat = (unsigned int **)almat(n,w);
         W_NEWUP2(p,w+1);    W_NEWUP2(p,w+1);
         W_NEWUP2(t,w+1);    W_NEWUP2(t,w+1);
         p->w = 1; p->b[0] = 1; /*  p = 1 mod f */    p->w = 1; p->b[0] = 1; /*  p = 1 mod f */
         for ( i = 1; i < n; i++ ) {    for ( i = 1; i < n; i++ ) {
                 _bshiftup2(p,-2,t); _remup2(t,f,p);      _bshiftup2(p,-2,t); _remup2(t,f,p);
                 _copy_up2bits(p,mat,i);      _copy_up2bits(p,mat,i);
   
         }    }
 /*      _print_frobmat(mat,n,w); */  /*  _print_frobmat(mat,n,w); */
         for ( j = 1; j < n; j++ ) {    for ( j = 1; j < n; j++ ) {
                 hcol = j/BSH; hbit = 1<<(j%BSH);      hcol = j/BSH; hbit = 1<<(j%BSH);
                 for ( i = j-1; i < n; i++ )      for ( i = j-1; i < n; i++ )
                         if ( mat[i][hcol] & hbit )        if ( mat[i][hcol] & hbit )
                                 break;          break;
                 if ( i == n )      if ( i == n )
                         return 0;        return 0;
                 if ( i != j-1 ) {      if ( i != j-1 ) {
                         u = mat[i]; mat[i] = mat[j-1]; mat[j-1] = u;        u = mat[i]; mat[i] = mat[j-1]; mat[j-1] = u;
                 }      }
                 for ( i = j; i < n; i++ )      for ( i = j; i < n; i++ )
                         if ( mat[i][hcol] & hbit )        if ( mat[i][hcol] & hbit )
                                 for ( k = hcol; k < w; k++ )          for ( k = hcol; k < w; k++ )
                                         mat[i][k] ^= mat[j-1][k];            mat[i][k] ^= mat[j-1][k];
         }    }
         return 1;    return 1;
 }  }
   
 int irredcheck_dddup2(f)  int irredcheck_dddup2(f)
 UP2 f;  UP2 f;
 {  {
         UP2 x,u,t,s,gcd;    UP2 x,u,t,s,gcd;
         int n,i;    int n,i;
   
         if ( !sqfrcheckup2(f) )    if ( !sqfrcheckup2(f) )
                 return -1;      return -1;
         n = degup2(f);    n = degup2(f);
   
         W_NEWUP2(u,1); u->w = 1; u->b[0] = 2; /* u = x mod f */    W_NEWUP2(u,1); u->w = 1; u->b[0] = 2; /* u = x mod f */
         x = u;    x = u;
         for ( i = 1; 2*i <= n; i++ ) {    for ( i = 1; 2*i <= n; i++ ) {
                 squareup2(u,&t); remup2(t,f,&u); addup2(u,x,&s);      squareup2(u,&t); remup2(t,f,&u); addup2(u,x,&s);
                 gcdup2(f,s,&gcd);      gcdup2(f,s,&gcd);
                 if ( degup2(gcd) > 0 )      if ( degup2(gcd) > 0 )
                         return 0;        return 0;
         }    }
         return 1;    return 1;
 }  }
   
 void _copy_up2bits(p,mat,pos)  void _copy_up2bits(p,mat,pos)
Line 1849  UP2 p;
Line 1849  UP2 p;
 unsigned int **mat;  unsigned int **mat;
 int pos;  int pos;
 {  {
         int d,col,j,jcol,jsh;    int d,col,j,jcol,jsh;
         unsigned int bit;    unsigned int bit;
   
         d = degup2(p);    d = degup2(p);
         col = pos/BSH; bit = 1<<(pos%BSH);    col = pos/BSH; bit = 1<<(pos%BSH);
   
         for ( j = 0; j <= d; j++ ) {    for ( j = 0; j <= d; j++ ) {
                 jcol = j/BSH; jsh = j%BSH;      jcol = j/BSH; jsh = j%BSH;
                 if ( p->b[jcol]&(1<<jsh) )      if ( p->b[jcol]&(1<<jsh) )
                         mat[j][col] |= bit;        mat[j][col] |= bit;
         }    }
         mat[pos][col] ^= bit;    mat[pos][col] ^= bit;
 }  }
   
 int compute_multiplication_matrix(p0,mp)  int compute_multiplication_matrix(p0,mp)
 P p0;  P p0;
 GF2MAT *mp;  GF2MAT *mp;
 {  {
         UP2 p;    UP2 p;
         int n,w,i,j,k,l;    int n,w,i,j,k,l;
         unsigned int **a,**b,**c,**d,**t,**m;    unsigned int **a,**b,**c,**d,**t,**m;
         GF2MAT am,bm;    GF2MAT am,bm;
   
         ptoup2(p0,&p);    ptoup2(p0,&p);
         n = degup2(p);    n = degup2(p);
         w = p->w;    w = p->w;
         if ( !compute_representation_conversion_matrix(p0,&am,&bm) )    if ( !compute_representation_conversion_matrix(p0,&am,&bm) )
                 return 0;      return 0;
         a = am->body; b = bm->body;    a = am->body; b = bm->body;
 /*      _print_frobmat(a,n,w); printf("\n"); */  /*  _print_frobmat(a,n,w); printf("\n"); */
 /*      _print_frobmat(b,n,w); printf("\n"); */  /*  _print_frobmat(b,n,w); printf("\n"); */
         c = (unsigned int **)almat(n,w);    c = (unsigned int **)almat(n,w);
         for ( i = 0; i < n-1; i++ ) {    for ( i = 0; i < n-1; i++ ) {
                 j = i+1;      j = i+1;
                 c[i][j/BSH] |= 1<<(j%BSH);      c[i][j/BSH] |= 1<<(j%BSH);
         }    }
         bcopy(p->b,c[n-1],p->w*sizeof(unsigned int));    bcopy(p->b,c[n-1],p->w*sizeof(unsigned int));
   
 /*      _print_frobmat(b,n,w); printf("\n"); */  /*  _print_frobmat(b,n,w); printf("\n"); */
         t = (unsigned int **)almat(n,w);    t = (unsigned int **)almat(n,w);
         mulgf2mat(n,a,c,t);    mulgf2mat(n,a,c,t);
         d = (unsigned int **)almat(n,w);    d = (unsigned int **)almat(n,w);
         mulgf2mat(n,t,b,d);    mulgf2mat(n,t,b,d);
 /*      _print_frobmat(d,n,w); printf("\n"); */  /*  _print_frobmat(d,n,w); printf("\n"); */
   
         m = (unsigned int **)almat(n,w);    m = (unsigned int **)almat(n,w);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0; j < n; j++ ) {      for ( j = 0; j < n; j++ ) {
                         k = (n-1-(j-i))%n; l = (n-1+i)%n;        k = (n-1-(j-i))%n; l = (n-1+i)%n;
                         if ( d[k][l/BSH] & 1<<(l%BSH) )        if ( d[k][l/BSH] & 1<<(l%BSH) )
                                 m[n-1-i][(n-1-j)/BSH] |= 1<<((n-1-j)%BSH);          m[n-1-i][(n-1-j)/BSH] |= 1<<((n-1-j)%BSH);
                 }      }
 /*      _print_frobmat(m,n,w); printf("\n"); */  /*  _print_frobmat(m,n,w); printf("\n"); */
         TOGF2MAT(n,n,m,*mp);    TOGF2MAT(n,n,m,*mp);
         return 1;    return 1;
 }  }
   
 #define GF2N_PBTOPB 0  #define GF2N_PBTOPB 0
Line 1921  P p0,p1;
Line 1921  P p0,p1;
 int to;  int to;
 GF2MAT *m01,*m10;  GF2MAT *m01,*m10;
 {  {
         UP2 up0;    UP2 up0;
         int n,w;    int n,w;
         unsigned int **p01,**p10;    unsigned int **p01,**p10;
         GF2N root;    GF2N root;
   
         setmod_gf2n(p1);    setmod_gf2n(p1);
   
         ptoup2(p0,&up0);    ptoup2(p0,&up0);
         n = degup2(up0); w = up0->w;    n = degup2(up0); w = up0->w;
         find_root_up2(up0,&root);    find_root_up2(up0,&root);
         compute_change_of_basis_matrix_with_root(p0,p1,to,root,m01,m10);    compute_change_of_basis_matrix_with_root(p0,p1,to,root,m01,m10);
 }  }
   
 void compute_change_of_basis_matrix_with_root(p0,p1,to,root,m01,m10)  void compute_change_of_basis_matrix_with_root(p0,p1,to,root,m01,m10)
Line 1940  int to;
Line 1940  int to;
 GF2N root;  GF2N root;
 GF2MAT *m01,*m10;  GF2MAT *m01,*m10;
 {  {
         UP2 up0,t,u,s;    UP2 up0,t,u,s;
         int n,w,i;    int n,w,i;
         unsigned int **a,**b,**g,**h;    unsigned int **a,**b,**g,**h;
         unsigned int **p01,**p10;    unsigned int **p01,**p10;
         P tmp;    P tmp;
   
         setmod_gf2n(p1);    setmod_gf2n(p1);
   
         ptoup2(p0,&up0);    ptoup2(p0,&up0);
         n = degup2(up0); w = up0->w;    n = degup2(up0); w = up0->w;
         substp(CO,p0,p0->v,(P)root,&tmp);    substp(CO,p0,p0->v,(P)root,&tmp);
         if ( tmp )    if ( tmp )
                 error("error in find_root_up2");      error("error in find_root_up2");
         u = root->body;    u = root->body;
   
         p01 = (unsigned int **)almat(n,w);    p01 = (unsigned int **)almat(n,w);
         p10 = (unsigned int **)almat(n,w);    p10 = (unsigned int **)almat(n,w);
         if ( to == GF2N_PBTOPB ) {    if ( to == GF2N_PBTOPB ) {
                 t = ONEUP2;      t = ONEUP2;
                 bcopy(t->b,p01[0],t->w*sizeof(unsigned int));      bcopy(t->b,p01[0],t->w*sizeof(unsigned int));
                 for ( i = 1; i < n; i++ ) {      for ( i = 1; i < n; i++ ) {
                         mulup2(t,u,&s); gen_simpup2_destructive(s,current_mod_gf2n); t = s;        mulup2(t,u,&s); gen_simpup2_destructive(s,current_mod_gf2n); t = s;
                         bcopy(t->b,p01[i],t->w*sizeof(unsigned int));        bcopy(t->b,p01[i],t->w*sizeof(unsigned int));
                 }      }
         } else {    } else {
                 t = u;      t = u;
                 bcopy(t->b,p01[n-1],t->w*sizeof(unsigned int));      bcopy(t->b,p01[n-1],t->w*sizeof(unsigned int));
                 for ( i = 1; i < n; i++ ) {      for ( i = 1; i < n; i++ ) {
                         squareup2(t,&s); gen_simpup2_destructive(s,current_mod_gf2n); t = s;        squareup2(t,&s); gen_simpup2_destructive(s,current_mod_gf2n); t = s;
                         bcopy(t->b,p01[n-1-i],t->w*sizeof(unsigned int));        bcopy(t->b,p01[n-1-i],t->w*sizeof(unsigned int));
                 }      }
         }    }
         if ( !invgf2mat(n,p01,p10) )    if ( !invgf2mat(n,p01,p10) )
                 error("compute_change_of_basis_matrix : something is wrong");      error("compute_change_of_basis_matrix : something is wrong");
         TOGF2MAT(n,n,p01,*m01);    TOGF2MAT(n,n,p01,*m01);
         TOGF2MAT(n,n,p10,*m10);    TOGF2MAT(n,n,p10,*m10);
 }  }
   
 /*  /*
Line 1989  int compute_representation_conversion_matrix(p0,np,pn)
Line 1989  int compute_representation_conversion_matrix(p0,np,pn)
 P p0;  P p0;
 GF2MAT *np,*pn;  GF2MAT *np,*pn;
 {  {
         UP2 x,s;    UP2 x,s;
         int n,w,i;    int n,w,i;
         unsigned int **ntop,**pton;    unsigned int **ntop,**pton;
   
         setmod_gf2n(p0);    setmod_gf2n(p0);
   
         n = UDEG(p0); w = n/BSH + (n%BSH?1:0);    n = UDEG(p0); w = n/BSH + (n%BSH?1:0);
   
         /* x = alpha */    /* x = alpha */
         NEWUP2(x,1); x->w = 1; x->b[0]=2;    NEWUP2(x,1); x->w = 1; x->b[0]=2;
         gen_simpup2_destructive(x,current_mod_gf2n);    gen_simpup2_destructive(x,current_mod_gf2n);
   
         ntop = (unsigned int **)almat(n,w);    ntop = (unsigned int **)almat(n,w);
         pton = (unsigned int **)almat(n,w);    pton = (unsigned int **)almat(n,w);
   
         bcopy(x->b,ntop[n-1],x->w*sizeof(unsigned int));    bcopy(x->b,ntop[n-1],x->w*sizeof(unsigned int));
         for ( i = 1; i < n; i++ ) {    for ( i = 1; i < n; i++ ) {
                 squareup2(x,&s); gen_simpup2_destructive(s,current_mod_gf2n); x = s;      squareup2(x,&s); gen_simpup2_destructive(s,current_mod_gf2n); x = s;
                 bcopy(x->b,ntop[n-1-i],x->w*sizeof(unsigned int));      bcopy(x->b,ntop[n-1-i],x->w*sizeof(unsigned int));
         }    }
         if ( !invgf2mat(n,ntop,pton) )    if ( !invgf2mat(n,ntop,pton) )
                 return 0;      return 0;
         TOGF2MAT(n,n,ntop,*np);    TOGF2MAT(n,n,ntop,*np);
         TOGF2MAT(n,n,pton,*pn);    TOGF2MAT(n,n,pton,*pn);
         return 1;    return 1;
 }  }
   
 void mul_nb(mat,a,b,c)  void mul_nb(mat,a,b,c)
 GF2MAT mat;  GF2MAT mat;
 unsigned int *a,*b,*c;  unsigned int *a,*b,*c;
 {  {
         int n,w,i;    int n,w,i;
         unsigned int *wa,*wb,*t;    unsigned int *wa,*wb,*t;
         unsigned int **m;    unsigned int **m;
   
         n = mat->row;    n = mat->row;
         m = mat->body;    m = mat->body;
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         wa = (unsigned int *)ALLOCA(w*sizeof(unsigned int));    wa = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
         wb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));    wb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
         t = (unsigned int *)ALLOCA(w*sizeof(unsigned int));    t = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
         bcopy(a,wa,w*sizeof(unsigned int));    bcopy(a,wa,w*sizeof(unsigned int));
         bcopy(b,wb,w*sizeof(unsigned int));    bcopy(b,wb,w*sizeof(unsigned int));
         bzero((char *)c,w*sizeof(unsigned int));    bzero((char *)c,w*sizeof(unsigned int));
         for ( i = n-1; i >= 0; i-- ) {    for ( i = n-1; i >= 0; i-- ) {
                 mulgf2vectmat(n,wa,m,t);      mulgf2vectmat(n,wa,m,t);
                 if ( mulgf2vectvect(n,t,wb) )      if ( mulgf2vectvect(n,t,wb) )
                         c[i/BSH] |= 1<<(i%BSH);        c[i/BSH] |= 1<<(i%BSH);
                 leftshift(wa,n);      leftshift(wa,n);
                 leftshift(wb,n);      leftshift(wb,n);
         }    }
 }  }
   
   
Line 2051  void leftshift(a,n)
Line 2051  void leftshift(a,n)
 unsigned int *a;  unsigned int *a;
 int n;  int n;
 {  {
         int r,w,i;    int r,w,i;
         unsigned int msb;    unsigned int msb;
   
         r = n%BSH;    r = n%BSH;
         w = n/BSH + (r?1:0);    w = n/BSH + (r?1:0);
         msb = a[(n-1)/BSH]&(1<<((n-1)%BSH));    msb = a[(n-1)/BSH]&(1<<((n-1)%BSH));
         for ( i = w-1; i > 0; i-- )    for ( i = w-1; i > 0; i-- )
                 a[i] = (a[i]<<1)|(a[i-1]>>(BSH-1));      a[i] = (a[i]<<1)|(a[i-1]>>(BSH-1));
         a[0] = (a[0]<<1)|(msb?1:0);    a[0] = (a[0]<<1)|(msb?1:0);
         if ( r )    if ( r )
                 a[w-1] &= (1<<r)-1;      a[w-1] &= (1<<r)-1;
 }  }
   
 void mat_to_gf2mat(a,b)  void mat_to_gf2mat(a,b)
 MAT a;  MAT a;
 unsigned int ***b;  unsigned int ***b;
 {  {
         int n,w,i,j;    int n,w,i,j;
         unsigned int **m;    unsigned int **m;
   
         n = a->row;    n = a->row;
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         *b = m = (unsigned int **)almat(n,w);    *b = m = (unsigned int **)almat(n,w);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0; j < n; j++ )      for ( j = 0; j < n; j++ )
                         if ( a->body[i][j] )        if ( a->body[i][j] )
                                 m[i][j/BSH] |= 1<<(j%BSH);          m[i][j/BSH] |= 1<<(j%BSH);
 }  }
   
 void gf2mat_to_mat(a,n,b)  void gf2mat_to_mat(a,n,b)
 unsigned int **a;  unsigned int **a;
 MAT *b;  MAT *b;
 {  {
         int w,i,j;    int w,i,j;
         MAT m;    MAT m;
   
         MKMAT(m,n,n); *b = m;    MKMAT(m,n,n); *b = m;
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0; j < n; j++ )      for ( j = 0; j < n; j++ )
                         if ( a[i][j/BSH] & 1<<(j%BSH) )        if ( a[i][j/BSH] & 1<<(j%BSH) )
                                 m->body[i][j] = (pointer)ONE;          m->body[i][j] = (pointer)ONE;
 }  }
   
 void mulgf2mat(n,a,b,c)  void mulgf2mat(n,a,b,c)
 int n;  int n;
 unsigned int **a,**b,**c;  unsigned int **a,**b,**c;
 {  {
         int i,j,k,w;    int i,j,k,w;
   
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 bzero((char *)c[i],w*sizeof(unsigned int));      bzero((char *)c[i],w*sizeof(unsigned int));
                 for ( j = 0; j < n; j++ )      for ( j = 0; j < n; j++ )
                         if ( a[i][j/BSH] & 1<<(j%BSH) )        if ( a[i][j/BSH] & 1<<(j%BSH) )
                                 for ( k = 0; k < w; k++ )          for ( k = 0; k < w; k++ )
                                         c[i][k] ^= b[j][k];            c[i][k] ^= b[j][k];
         }    }
 }  }
   
 /* c = a*b; where a, c are row vectors */  /* c = a*b; where a, c are row vectors */
Line 2119  unsigned int *a;
Line 2119  unsigned int *a;
 unsigned int **b;  unsigned int **b;
 unsigned int *c;  unsigned int *c;
 {  {
         int j,k,w;    int j,k,w;
   
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         bzero((char *)c,w*sizeof(unsigned int));    bzero((char *)c,w*sizeof(unsigned int));
         for ( j = 0; j < n; j++ )    for ( j = 0; j < n; j++ )
                 if ( a[j/BSH] & 1<<(j%BSH) )      if ( a[j/BSH] & 1<<(j%BSH) )
                         for ( k = 0; k < w; k++ )        for ( k = 0; k < w; k++ )
                                 c[k] ^= b[j][k];          c[k] ^= b[j][k];
 }  }
   
 int mulgf2vectvect(n,a,b)  int mulgf2vectvect(n,a,b)
 int n;  int n;
 unsigned int *a,*b;  unsigned int *a,*b;
 {  {
         unsigned int t,r;    unsigned int t,r;
         int i,w;    int i,w;
   
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         for ( i = 0, r = 0; i < w; i++ )    for ( i = 0, r = 0; i < w; i++ )
                 for ( t = a[i]&b[i]; t; t >>=1 )      for ( t = a[i]&b[i]; t; t >>=1 )
                         r ^= (t&1);        r ^= (t&1);
         return r;    return r;
 }  }
   
 int invgf2mat(n,a,b)  int invgf2mat(n,a,b)
 int n;  int n;
 unsigned int **a,**b;  unsigned int **a,**b;
 {  {
         int i,j,k,hcol,hbit,w;    int i,j,k,hcol,hbit,w;
         unsigned int *u;    unsigned int *u;
         unsigned int **s;    unsigned int **s;
   
         w = n/BSH + (n%BSH?1:0);    w = n/BSH + (n%BSH?1:0);
         s = (unsigned int **)almat(n,w);    s = (unsigned int **)almat(n,w);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 bcopy(a[i],s[i],w*sizeof(unsigned int));      bcopy(a[i],s[i],w*sizeof(unsigned int));
   
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 b[i][i/BSH] |= 1<<(i%BSH);      b[i][i/BSH] |= 1<<(i%BSH);
   
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 hcol = j/BSH; hbit = 1<<(j%BSH);      hcol = j/BSH; hbit = 1<<(j%BSH);
                 for ( i = j; i < n; i++ )      for ( i = j; i < n; i++ )
                         if ( s[i][hcol] & hbit )        if ( s[i][hcol] & hbit )
                                 break;          break;
                 if ( i == n )      if ( i == n )
                         return 0;        return 0;
                 if ( i != j ) {      if ( i != j ) {
                         u = s[i]; s[i] = s[j]; s[j] = u;        u = s[i]; s[i] = s[j]; s[j] = u;
                         u = b[i]; b[i] = b[j]; b[j] = u;        u = b[i]; b[i] = b[j]; b[j] = u;
                 }      }
                 for ( i = 0; i < n; i++ ) {      for ( i = 0; i < n; i++ ) {
                         if ( i == j )        if ( i == j )
                                 continue;          continue;
                         if ( s[i][hcol] & hbit ) {        if ( s[i][hcol] & hbit ) {
                                 for ( k = hcol; k < w; k++ )          for ( k = hcol; k < w; k++ )
                                         s[i][k] ^= s[j][k];            s[i][k] ^= s[j][k];
                                 for ( k = 0; k < w; k++ )          for ( k = 0; k < w; k++ )
                                         b[i][k] ^= b[j][k];            b[i][k] ^= b[j][k];
                         }        }
                 }      }
         }    }
         return 1;    return 1;
 }  }
   
 INLINE void _mulup2_11(a1,a2,ar)  INLINE void _mulup2_11(a1,a2,ar)
 unsigned int a1,a2;  unsigned int a1,a2;
 unsigned int *ar;  unsigned int *ar;
 {  {
         GF2M_MUL_1(ar[1],ar[0],a1,a2);    GF2M_MUL_1(ar[1],ar[0],a1,a2);
 }  }
   
 void _mulup2_22(a1,a2,ar)  void _mulup2_22(a1,a2,ar)
 unsigned int *a1,*a2,*ar;  unsigned int *a1,*a2,*ar;
 {  {
         unsigned int m[2];    unsigned int m[2];
   
         _mulup2_11(*a1,*a2,ar);    _mulup2_11(*a1,*a2,ar);
         _mulup2_11(*(a1+1),*(a2+1),ar+2);    _mulup2_11(*(a1+1),*(a2+1),ar+2);
         _mulup2_11(a1[0]^a1[1],a2[0]^a2[1],m);    _mulup2_11(a1[0]^a1[1],a2[0]^a2[1],m);
         m[0] ^= ar[0]^ar[2]; m[1] ^= ar[1]^ar[3];    m[0] ^= ar[0]^ar[2]; m[1] ^= ar[1]^ar[3];
   
         ar[1] ^= m[0]; ar[2] ^= m[1];    ar[1] ^= m[0]; ar[2] ^= m[1];
 }  }
   
 #if 0  #if 0
 void _mulup2_33(a1,a2,ar)  void _mulup2_33(a1,a2,ar)
 unsigned int *a1,*a2,*ar;  unsigned int *a1,*a2,*ar;
 {  {
         unsigned int m[4];    unsigned int m[4];
         unsigned int c1[2],c2[2];    unsigned int c1[2],c2[2];
   
         c1[0] = a1[2]^a1[0]; c1[1] = a1[1];    c1[0] = a1[2]^a1[0]; c1[1] = a1[1];
         c2[0] = a2[2]^a2[0]; c2[1] = a2[1];    c2[0] = a2[2]^a2[0]; c2[1] = a2[1];
   
         _mulup2_22(a1,a2,ar);    _mulup2_22(a1,a2,ar);
         _mulup2_11(*(a1+2),*(a2+2),ar+4);    _mulup2_11(*(a1+2),*(a2+2),ar+4);
         _mulup2_22(c1,c2,m);    _mulup2_22(c1,c2,m);
   
         m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];    m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
         m[2] ^= ar[2]; m[3] ^= ar[3];    m[2] ^= ar[2]; m[3] ^= ar[3];
   
         ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];    ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
 }  }
 #else  #else
 /* (ar[5]...ar[0]) = (a1[2]a1[1]a1[0])*(a2[2]a2[1]a2[0])  */  /* (ar[5]...ar[0]) = (a1[2]a1[1]a1[0])*(a2[2]a2[1]a2[0])  */
Line 2229  unsigned int *a1,*a2,*ar;
Line 2229  unsigned int *a1,*a2,*ar;
 void _mulup2_33(a1,a2,ar)  void _mulup2_33(a1,a2,ar)
 unsigned int *a1,*a2,*ar;  unsigned int *a1,*a2,*ar;
 {  {
         unsigned int m[4];    unsigned int m[4];
         unsigned int c[2];    unsigned int c[2];
         unsigned int t,s;    unsigned int t,s;
   
         /* (ar[1],ar[0]) = a1[0]*a2[0] */    /* (ar[1],ar[0]) = a1[0]*a2[0] */
         _mulup2_11(a1[0],a2[0],ar);    _mulup2_11(a1[0],a2[0],ar);
   
         /* (ar[3],ar[2]) = hi(m) = a1[1]*a1[1] */    /* (ar[3],ar[2]) = hi(m) = a1[1]*a1[1] */
         _mulup2_11(a1[1],a2[1],ar+2); m[2] = ar[2]; m[3] = ar[3];    _mulup2_11(a1[1],a2[1],ar+2); m[2] = ar[2]; m[3] = ar[3];
   
         /* (ar[5],ar[4]) = a1[2]*a2[2] */    /* (ar[5],ar[4]) = a1[2]*a2[2] */
         _mulup2_11(a1[2],a2[2],ar+4);    _mulup2_11(a1[2],a2[2],ar+4);
   
         /*    /*
          * (ar[3],ar[2],ar[1],ar[0]) = (a1[1],a1[0])*(a2[1],a2[0])     * (ar[3],ar[2],ar[1],ar[0]) = (a1[1],a1[0])*(a2[1],a2[0])
          * a1[1]*a2[1] is already in (ar[3],ar[2])     * a1[1]*a2[1] is already in (ar[3],ar[2])
          */     */
         /* c = (a1[1]+a1[0])*(a2[1]+a2[0]) */    /* c = (a1[1]+a1[0])*(a2[1]+a2[0]) */
         t = a1[1]^a1[0]; s = a2[1]^a2[0]; _mulup2_11(t,s,c);    t = a1[1]^a1[0]; s = a2[1]^a2[0]; _mulup2_11(t,s,c);
         /* c += (ar[1],ar[0])+(ar[3],ar[2]) */    /* c += (ar[1],ar[0])+(ar[3],ar[2]) */
         c[0] ^= ar[2]^ar[0]; c[1] ^= ar[3]^ar[1];    c[0] ^= ar[2]^ar[0]; c[1] ^= ar[3]^ar[1];
         /* (ar[2],ar[1]) += c */    /* (ar[2],ar[1]) += c */
         ar[1] ^= c[0]; ar[2] ^= c[1];    ar[1] ^= c[0]; ar[2] ^= c[1];
   
         /*    /*
          * m = (a1[1],a1[2]+a1[0])*(a2[1],a2[2]+a2[0])     * m = (a1[1],a1[2]+a1[0])*(a2[1],a2[2]+a2[0])
          * a1[1]*a2[1] is already in hi(m)     * a1[1]*a2[1] is already in hi(m)
          *     *
          */     */
         /* low(m) = (a1[0]+a1[2])*(a2[0]+a2[2]) */    /* low(m) = (a1[0]+a1[2])*(a2[0]+a2[2]) */
         t = a1[2]^a1[0]; s = a2[2]^a2[0]; _mulup2_11(t,s,m);    t = a1[2]^a1[0]; s = a2[2]^a2[0]; _mulup2_11(t,s,m);
         /* c = (a1[1]+(a1[0]+a1[2]))*(a2[1]+(a2[0]+a2[2])) */    /* c = (a1[1]+(a1[0]+a1[2]))*(a2[1]+(a2[0]+a2[2])) */
         t ^= a1[1]; s ^= a2[1]; _mulup2_11(t,s,c);    t ^= a1[1]; s ^= a2[1]; _mulup2_11(t,s,c);
         /* c += low(m)+hi(m) */    /* c += low(m)+hi(m) */
         c[0] ^= m[2]^m[0]; c[1] ^= m[3]^m[1];    c[0] ^= m[2]^m[0]; c[1] ^= m[3]^m[1];
         /* mid(m) += c */    /* mid(m) += c */
         m[1] ^= c[0]; m[2] ^= c[1];    m[1] ^= c[0]; m[2] ^= c[1];
   
         /* m += (ar[3],ar[2],ar[1],ar[0])+(ar[5],ar[4]) */    /* m += (ar[3],ar[2],ar[1],ar[0])+(ar[5],ar[4]) */
         m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];    m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
         m[2] ^= ar[2]; m[3] ^= ar[3];    m[2] ^= ar[2]; m[3] ^= ar[3];
   
         /* (ar[5],ar[4],ar[3],ar[2]) += m */    /* (ar[5],ar[4],ar[3],ar[2]) += m */
         ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];    ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
 }  }
 #endif  #endif
   
 void _mulup2_44(a1,a2,ar)  void _mulup2_44(a1,a2,ar)
 unsigned int *a1,*a2,*ar;  unsigned int *a1,*a2,*ar;
 {  {
         unsigned int m[4];    unsigned int m[4];
         unsigned int c1[2],c2[2];    unsigned int c1[2],c2[2];
   
         c1[0] = a1[2]^a1[0]; c1[1] = a1[3]^a1[1];    c1[0] = a1[2]^a1[0]; c1[1] = a1[3]^a1[1];
         c2[0] = a2[2]^a2[0]; c2[1] = a2[3]^a2[1];    c2[0] = a2[2]^a2[0]; c2[1] = a2[3]^a2[1];
   
         _mulup2_22(a1,a2,ar);    _mulup2_22(a1,a2,ar);
         _mulup2_22(a1+2,a2+2,ar+4);    _mulup2_22(a1+2,a2+2,ar+4);
         _mulup2_22(c1,c2,m);    _mulup2_22(c1,c2,m);
   
         m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];    m[0] ^= ar[0]^ar[4]; m[1] ^= ar[1]^ar[5];
         m[2] ^= ar[2]^ar[6]; m[3] ^= ar[3]^ar[7];    m[2] ^= ar[2]^ar[6]; m[3] ^= ar[3]^ar[7];
   
         ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];    ar[2] ^= m[0]; ar[3] ^= m[1]; ar[4] ^= m[2]; ar[5] ^= m[3];
 }  }
   
 void _mulup2_55(a1,a2,ar)  void _mulup2_55(a1,a2,ar)
 unsigned int *a1,*a2,*ar;  unsigned int *a1,*a2,*ar;
 {  {
         unsigned int m[6];    unsigned int m[6];
         unsigned int c1[3],c2[3];    unsigned int c1[3],c2[3];
   
         c1[0]=a1[3]^a1[0]; c1[1]=a1[4]^a1[1]; c1[2]=a1[2];    c1[0]=a1[3]^a1[0]; c1[1]=a1[4]^a1[1]; c1[2]=a1[2];
         c2[0]=a2[3]^a2[0]; c2[1]=a2[4]^a2[1]; c2[2]=a2[2];    c2[0]=a2[3]^a2[0]; c2[1]=a2[4]^a2[1]; c2[2]=a2[2];
   
         _mulup2_33(a1,a2,ar);    _mulup2_33(a1,a2,ar);
         _mulup2_22(a1+3,a2+3,ar+6);    _mulup2_22(a1+3,a2+3,ar+6);
         _mulup2_33(c1,c2,m);    _mulup2_33(c1,c2,m);
   
         m[0] ^= ar[0]^ar[6]; m[1] ^= ar[1]^ar[7]; m[2] ^= ar[2]^ar[8];    m[0] ^= ar[0]^ar[6]; m[1] ^= ar[1]^ar[7]; m[2] ^= ar[2]^ar[8];
         m[3] ^= ar[3]^ar[9]; m[4] ^= ar[4]; m[5] ^= ar[5];    m[3] ^= ar[3]^ar[9]; m[4] ^= ar[4]; m[5] ^= ar[5];
   
         ar[3] ^= m[0]; ar[4] ^= m[1]; ar[5] ^= m[2];    ar[3] ^= m[0]; ar[4] ^= m[1]; ar[5] ^= m[2];
         ar[6] ^= m[3]; ar[7] ^= m[4]; ar[8] ^= m[5];    ar[6] ^= m[3]; ar[7] ^= m[4]; ar[8] ^= m[5];
 }  }
   
 void _mulup2_66(a1,a2,ar)  void _mulup2_66(a1,a2,ar)
 unsigned int *a1,*a2,*ar;  unsigned int *a1,*a2,*ar;
 {  {
         unsigned int m[6];    unsigned int m[6];
         unsigned int c1[3],c2[3];    unsigned int c1[3],c2[3];
   
         c1[0]=a1[3]^a1[0]; c1[1]=a1[4]^a1[1]; c1[2]=a1[5]^a1[2];    c1[0]=a1[3]^a1[0]; c1[1]=a1[4]^a1[1]; c1[2]=a1[5]^a1[2];
         c2[0]=a2[3]^a2[0]; c2[1]=a2[4]^a2[1]; c2[2]=a2[5]^a2[2];    c2[0]=a2[3]^a2[0]; c2[1]=a2[4]^a2[1]; c2[2]=a2[5]^a2[2];
   
         _mulup2_33(a1,a2,ar);    _mulup2_33(a1,a2,ar);
         _mulup2_33(a1+3,a2+3,ar+6);    _mulup2_33(a1+3,a2+3,ar+6);
         _mulup2_33(c1,c2,m);    _mulup2_33(c1,c2,m);
   
         m[0] ^= ar[0]^ar[6]; m[1] ^= ar[1]^ar[7]; m[2] ^= ar[2]^ar[8];    m[0] ^= ar[0]^ar[6]; m[1] ^= ar[1]^ar[7]; m[2] ^= ar[2]^ar[8];
         m[3] ^= ar[3]^ar[9]; m[4] ^= ar[4]^ar[10]; m[5] ^= ar[5]^ar[11];    m[3] ^= ar[3]^ar[9]; m[4] ^= ar[4]^ar[10]; m[5] ^= ar[5]^ar[11];
   
         ar[3] ^= m[0]; ar[4] ^= m[1]; ar[5] ^= m[2];    ar[3] ^= m[0]; ar[4] ^= m[1]; ar[5] ^= m[2];
         ar[6] ^= m[3]; ar[7] ^= m[4]; ar[8] ^= m[5];    ar[6] ^= m[3]; ar[7] ^= m[4]; ar[8] ^= m[5];
 }  }
   
 #if 0  #if 0
Line 2340  void printup2_(f,l)
Line 2340  void printup2_(f,l)
 unsigned int *f;  unsigned int *f;
 int l;  int l;
 {  {
         int i;    int i;
   
         for ( i = 32*l-1; i >= 0; i-- ) {    for ( i = 32*l-1; i >= 0; i-- ) {
                 if ( f[i>>5]&(1<<(i&31)) )      if ( f[i>>5]&(1<<(i&31)) )
                         fprintf(stderr,"+x^%d",i);        fprintf(stderr,"+x^%d",i);
         }    }
         fprintf(stderr,"\n");    fprintf(stderr,"\n");
 }  }
 #endif  #endif
   
Line 2355  UP2 a;
Line 2355  UP2 a;
 int n;  int n;
 UP2 *inv;  UP2 *inv;
 {  {
         int lf,lg,i,j,k,lg2,df,dg,l,w;    int lf,lg,i,j,k,lg2,df,dg,l,w;
         unsigned int r;    unsigned int r;
         UP2 wf,wg,wb,wc,ret,t;    UP2 wf,wg,wb,wc,ret,t;
         unsigned int *p;    unsigned int *p;
   
         lf = a->w;    lf = a->w;
         lg = (n>>5)+1;    lg = (n>>5)+1;
   
         W_NEWUP2(wf,lf+1); _copyup2(a,wf);    W_NEWUP2(wf,lf+1); _copyup2(a,wf);
         W_NEWUP2(wg,lg+1); wg->w = lg;    W_NEWUP2(wg,lg+1); wg->w = lg;
         for ( i = lg-1, p = wg->b; i>=0; i-- )    for ( i = lg-1, p = wg->b; i>=0; i-- )
                 p[i] = 0xffffffff;      p[i] = 0xffffffff;
         if ( r = (n+1)&31 )    if ( r = (n+1)&31 )
                 p[lg-1] &= (1<<r)-1;      p[lg-1] &= (1<<r)-1;
         lg2 = lg*2;    lg2 = lg*2;
         W_NEWUP2(wb,lg2+2); wb->w = 1; wb->b[0] = 1;    W_NEWUP2(wb,lg2+2); wb->w = 1; wb->b[0] = 1;
         W_NEWUP2(wc,lg2+2); wc->w = 0;    W_NEWUP2(wc,lg2+2); wc->w = 0;
         k = 0;    k = 0;
         df = degup2(wf); dg = degup2(wg);    df = degup2(wf); dg = degup2(wg);
         while ( 1 ) {    while ( 1 ) {
                 lf = wf->w;      lf = wf->w;
                 p = wf->b;      p = wf->b;
                 for ( j = 0; j < lf && !p[j]; j++ );      for ( j = 0; j < lf && !p[j]; j++ );
                 for ( l = j<<5, w = p[j]; !(w&1); w >>=1, l++ );      for ( l = j<<5, w = p[j]; !(w&1); w >>=1, l++ );
                 _bshiftup2_destructive(wf,l);      _bshiftup2_destructive(wf,l);
                 _bshiftup2_destructive(wc,-l);      _bshiftup2_destructive(wc,-l);
                 k += l;      k += l;
                 df -= l;      df -= l;
                 if ( wf->w == 1 && wf->b[0] == 1 ) {      if ( wf->w == 1 && wf->b[0] == 1 ) {
                         k %= (n+1);        k %= (n+1);
                         if ( k != 1 ) {        if ( k != 1 ) {
                                 _bshiftup2_destructive(wb,k-(n+1));          _bshiftup2_destructive(wb,k-(n+1));
                                 remup2_type1_destructive(wb,n);          remup2_type1_destructive(wb,n);
                         }        }
                         NEWUP2(ret,wb->w);        NEWUP2(ret,wb->w);
                         _copyup2(wb,ret);        _copyup2(wb,ret);
                         *inv = ret;        *inv = ret;
                         return;        return;
                 }      }
                 if ( df < dg ) {      if ( df < dg ) {
                         i = df; df = dg; dg = i;        i = df; df = dg; dg = i;
                         t = wf; wf = wg; wg = t;        t = wf; wf = wg; wg = t;
                         t = wb; wb = wc; wc = t;        t = wb; wb = wc; wc = t;
                 }      }
                 /* df >= dg */      /* df >= dg */
                 _addup2_destructive(wf,wg);      _addup2_destructive(wf,wg);
                 df = degup2(wf);      df = degup2(wf);
                 _addup2_destructive(wb,wc);      _addup2_destructive(wb,wc);
         }    }
 }  }
   
 UP2 *compute_tab_gf2n(f)  UP2 *compute_tab_gf2n(f)
 UP2 f;  UP2 f;
 {  {
         GEN_UP2 mod;    GEN_UP2 mod;
         int m,n,w,i;    int m,n,w,i;
         UP2 t;    UP2 t;
         UP2 *tab;    UP2 *tab;
   
         mod = current_mod_gf2n;    mod = current_mod_gf2n;
         m = degup2(mod->dense);    m = degup2(mod->dense);
         n = degup2(f);    n = degup2(f);
         w = f->w;    w = f->w;
         tab = (UP2 *)CALLOC(m,sizeof(UP2));    tab = (UP2 *)CALLOC(m,sizeof(UP2));
         NEWUP2(tab[0],w); tab[0]->w = 1; tab[0]->b[0] = 2;    NEWUP2(tab[0],w); tab[0]->w = 1; tab[0]->b[0] = 2;
         for ( i = 1; i < m; i++ ) {    for ( i = 1; i < m; i++ ) {
                 squareup2(tab[i-1],&t); remup2(t,f,&tab[i]);      squareup2(tab[i-1],&t); remup2(t,f,&tab[i]);
         }    }
         return tab;    return tab;
 }  }
   
 UP compute_trace_gf2n(tab,c,n)  UP compute_trace_gf2n(tab,c,n)
Line 2431  UP2 *tab;
Line 2431  UP2 *tab;
 GF2N c;  GF2N c;
 int n;  int n;
 {  {
         GEN_UP2 mod;    GEN_UP2 mod;
         int w,m,i,j;    int w,m,i,j;
         UP2 *trace;    UP2 *trace;
         UP2 c1,c2,s,t;    UP2 c1,c2,s,t;
         UP r;    UP r;
         GF2N a;    GF2N a;
   
         mod = current_mod_gf2n;    mod = current_mod_gf2n;
         w = mod->dense->w;    w = mod->dense->w;
         m = degup2(mod->dense);    m = degup2(mod->dense);
         trace = (UP2 *)ALLOCA(n*sizeof(UP2));    trace = (UP2 *)ALLOCA(n*sizeof(UP2));
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 W_NEWUP2(trace[i],w); trace[i]->w = 0;      W_NEWUP2(trace[i],w); trace[i]->w = 0;
         }    }
         W_NEWUP2(c1,2*w); _copyup2(c->body,c1);    W_NEWUP2(c1,2*w); _copyup2(c->body,c1);
         W_NEWUP2(c2,2*w);    W_NEWUP2(c2,2*w);
   
         for ( i = 0; i < m; i++ ) {    for ( i = 0; i < m; i++ ) {
                 t = tab[i];      t = tab[i];
                 /* trace += c^(2^i)*tab[i] */      /* trace += c^(2^i)*tab[i] */
                 for ( j = degup2(t); j >= 0; j-- )      for ( j = degup2(t); j >= 0; j-- )
                         if ( t->b[j/BSH] & (1<<(j%BSH)) )        if ( t->b[j/BSH] & (1<<(j%BSH)) )
                                 _addup2_destructive(trace[j],c1);          _addup2_destructive(trace[j],c1);
                 _squareup2(c1,c2); gen_simpup2_destructive(c2,mod);      _squareup2(c1,c2); gen_simpup2_destructive(c2,mod);
                 t = c1; c1 = c2; c2 = t;      t = c1; c1 = c2; c2 = t;
         }    }
         for ( i = n-1; i >= 0 && !trace[i]->w; i-- );    for ( i = n-1; i >= 0 && !trace[i]->w; i-- );
         if ( i < 0 )    if ( i < 0 )
                 r = 0;      r = 0;
         else {    else {
                 r = UPALLOC(i); r->d = i;      r = UPALLOC(i); r->d = i;
                 for ( j = 0; j <= i; j++ )      for ( j = 0; j <= i; j++ )
                         if ( t = trace[j] ) {        if ( t = trace[j] ) {
                                 NEWUP2(s,t->w); _copyup2(t,s);          NEWUP2(s,t->w); _copyup2(t,s);
                                 MKGF2N(s,a); r->c[j] = (Num)a;          MKGF2N(s,a); r->c[j] = (Num)a;
                         }        }
         }    }
         return r;    return r;
 }  }
   
 void up2toup(f,r)  void up2toup(f,r)
 UP2 f;  UP2 f;
 UP *r;  UP *r;
 {  {
         int d,i;    int d,i;
         UP2 c;    UP2 c;
         UP t;    UP t;
         GF2N a;    GF2N a;
   
         if ( !f || !f->w )    if ( !f || !f->w )
                 *r = 0;      *r = 0;
         else {    else {
                 d = degup2(f);      d = degup2(f);
                 t = UPALLOC(d); *r = t;      t = UPALLOC(d); *r = t;
                 t->d = d;      t->d = d;
                 for ( i = 0; i <= d; i++ ) {      for ( i = 0; i <= d; i++ ) {
                         if ( f->b[i/BSH] & (1<<(i%BSH)) ) {        if ( f->b[i/BSH] & (1<<(i%BSH)) ) {
                                 NEWUP2(c,1);          NEWUP2(c,1);
                                 c->w = 1; c->b[0] = 1;          c->w = 1; c->b[0] = 1;
                                 MKGF2N(c,a);          MKGF2N(c,a);
                                 t->c[i] = (Num)a;          t->c[i] = (Num)a;
                         }        }
                 }      }
         }    }
 }  }
   
 void find_root_up2(f,r)  void find_root_up2(f,r)
 UP2 f;  UP2 f;
 GF2N *r;  GF2N *r;
 {  {
         int n;    int n;
         UP2 *tab;    UP2 *tab;
         UP uf,trace,gcd,quo,rem;    UP uf,trace,gcd,quo,rem;
         GF2N c;    GF2N c;
         int i;    int i;
   
         /* computeation of tab : tab[i] = t^(2^i) mod f (i=0,...,n-1) */    /* computeation of tab : tab[i] = t^(2^i) mod f (i=0,...,n-1) */
         n = degup2(f);    n = degup2(f);
         tab = compute_tab_gf2n(f);    tab = compute_tab_gf2n(f);
         up2toup(f,&uf);    up2toup(f,&uf);
         while ( uf->d > 1 ) {    while ( uf->d > 1 ) {
                 randomgf2n(&c);      randomgf2n(&c);
                 if ( !c )      if ( !c )
                         continue;        continue;
                 /* trace = (ct)+(ct)^2+...+(ct)^(2^(n-1)) */      /* trace = (ct)+(ct)^2+...+(ct)^(2^(n-1)) */
                 trace = compute_trace_gf2n(tab,c,n);      trace = compute_trace_gf2n(tab,c,n);
                 gcdup(trace,uf,&gcd);      gcdup(trace,uf,&gcd);
                 if ( gcd->d > 0 && gcd->d < uf->d ) {      if ( gcd->d > 0 && gcd->d < uf->d ) {
 /*                      fprintf(stderr,"deg(GCD)=%d\n",gcd->d); */  /*      fprintf(stderr,"deg(GCD)=%d\n",gcd->d); */
                         if ( 2*gcd->d > uf->d ) {        if ( 2*gcd->d > uf->d ) {
                                 qrup(uf,gcd,&quo,&rem);          qrup(uf,gcd,&quo,&rem);
                                 uf = quo;          uf = quo;
                         } else {        } else {
                                 uf = gcd;          uf = gcd;
                         }        }
                 }      }
         }    }
         /* now deg(uf) = 1. The root of uf is ug->c[0]/uf->c[1] */    /* now deg(uf) = 1. The root of uf is ug->c[0]/uf->c[1] */
         divgf2n((GF2N)uf->c[0],(GF2N)uf->c[1],r);    divgf2n((GF2N)uf->c[0],(GF2N)uf->c[1],r);
 }  }
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>