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

Diff for /OpenXM_contrib2/asir2000/builtin/array.c between version 1.75 and 1.76

version 1.75, 2017/09/17 02:34:02 version 1.76, 2018/03/29 01:32:50
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/builtin/array.c,v 1.74 2017/09/15 01:52:51 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/array.c,v 1.75 2017/09/17 02:34:02 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
Line 104  void Plusolve_prep();
Line 104  void Plusolve_prep();
 void Plusolve_main();  void Plusolve_main();
   
 struct ftab array_tab[] = {  struct ftab array_tab[] = {
         {"lu_mat",Plu_mat,1},    {"lu_mat",Plu_mat,1},
         {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},    {"solve_by_lu_gfmmat",Psolve_by_lu_gfmmat,4},
         {"lu_gfmmat",Plu_gfmmat,2},    {"lu_gfmmat",Plu_gfmmat,2},
         {"mat_to_gfmmat",Pmat_to_gfmmat,2},    {"mat_to_gfmmat",Pmat_to_gfmmat,2},
         {"generic_gauss_elim",Pgeneric_gauss_elim,1},    {"generic_gauss_elim",Pgeneric_gauss_elim,1},
         {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},    {"generic_gauss_elim_mod",Pgeneric_gauss_elim_mod,2},
         {"indep_rows_mod",Pindep_rows_mod,2},    {"indep_rows_mod",Pindep_rows_mod,2},
         {"newvect",Pnewvect,-2},    {"newvect",Pnewvect,-2},
         {"vect",Pvect,-99999999},    {"vect",Pvect,-99999999},
         {"vector",Pnewvect,-2},    {"vector",Pnewvect,-2},
         {"exponent_vector",Pexponent_vector,-99999999},    {"exponent_vector",Pexponent_vector,-99999999},
         {"newmat",Pnewmat,-3},    {"newmat",Pnewmat,-3},
         {"matrix",Pnewmat,-3},    {"matrix",Pnewmat,-3},
         {"mat",Pmat,-99999999},    {"mat",Pmat,-99999999},
         {"matr",Pmat,-99999999},    {"matr",Pmat,-99999999},
         {"matc",Pmatc,-99999999},    {"matc",Pmatc,-99999999},
         {"newbytearray",Pnewbytearray,-2},    {"newbytearray",Pnewbytearray,-2},
         {"memoryplot_to_coord",Pmemoryplot_to_coord,1},    {"memoryplot_to_coord",Pmemoryplot_to_coord,1},
         {"sepmat_destructive",Psepmat_destructive,2},    {"sepmat_destructive",Psepmat_destructive,2},
         {"sepvect",Psepvect,2},    {"sepvect",Psepvect,2},
         {"qsort",Pqsort,-2},    {"qsort",Pqsort,-2},
         {"vtol",Pvtol,1},    {"vtol",Pvtol,1},
         {"ltov",Pltov,1},    {"ltov",Pltov,1},
         {"size",Psize,1},    {"size",Psize,1},
         {"det",Pdet,-2},    {"det",Pdet,-2},
         {"nd_det",Pnd_det,-2},    {"nd_det",Pnd_det,-2},
         {"invmat",Pinvmat,-2},    {"invmat",Pinvmat,-2},
         {"leqm",Pleqm,2},    {"leqm",Pleqm,2},
         {"leqm1",Pleqm1,2},    {"leqm1",Pleqm1,2},
         {"geninvm",Pgeninvm,2},    {"geninvm",Pgeninvm,2},
         {"geninvm_swap",Pgeninvm_swap,2},    {"geninvm_swap",Pgeninvm_swap,2},
         {"geninv_sf_swap",Pgeninv_sf_swap,1},    {"geninv_sf_swap",Pgeninv_sf_swap,1},
         {"remainder",Premainder,2},    {"remainder",Premainder,2},
         {"sremainder",Psremainder,2},    {"sremainder",Psremainder,2},
         {"mulmat_gf2n",Pmulmat_gf2n,1},    {"mulmat_gf2n",Pmulmat_gf2n,1},
         {"bconvmat_gf2n",Pbconvmat_gf2n,-4},    {"bconvmat_gf2n",Pbconvmat_gf2n,-4},
         {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},    {"mul_vect_mat_gf2n",Pmul_vect_mat_gf2n,2},
         {"mul_mat_vect_int",Pmul_mat_vect_int,2},    {"mul_mat_vect_int",Pmul_mat_vect_int,2},
         {"nbmul_gf2n",PNBmul_gf2n,3},    {"nbmul_gf2n",PNBmul_gf2n,3},
         {"x962_irredpoly_up2",Px962_irredpoly_up2,2},    {"x962_irredpoly_up2",Px962_irredpoly_up2,2},
         {"irredpoly_up2",Pirredpoly_up2,2},    {"irredpoly_up2",Pirredpoly_up2,2},
         {"nbpoly_up2",Pnbpoly_up2,2},    {"nbpoly_up2",Pnbpoly_up2,2},
         {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},    {"mat_swap_row_destructive",Pmat_swap_row_destructive,3},
         {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},    {"mat_swap_col_destructive",Pmat_swap_col_destructive,3},
         {"mat_col",Pmat_col,2},    {"mat_col",Pmat_col,2},
         {"lusolve_prep",Plusolve_prep,1},    {"lusolve_prep",Plusolve_prep,1},
         {"lusolve_main",Plusolve_main,1},    {"lusolve_main",Plusolve_main,1},
         {"triangleq",Ptriangleq,1},    {"triangleq",Ptriangleq,1},
         {0,0,0},    {0,0,0},
 };  };
   
 typedef struct _ent { int j; unsigned int e; } ent;  typedef struct _ent { int j; unsigned int e; } ent;
Line 171  static int modulus;
Line 171  static int modulus;
   
 void Plusolve_prep(NODE arg,Q *rp)  void Plusolve_prep(NODE arg,Q *rp)
 {  {
         char *fname;    char *fname;
         FILE *in;    FILE *in;
         int len,i,rank;    int len,i,rank;
         int *rhs;    int *rhs;
   
         fname = BDY((STRING)ARG0(arg));    fname = BDY((STRING)ARG0(arg));
         in = fopen(fname,"r");    in = fopen(fname,"r");
         modulus = getw(in);    modulus = getw(in);
         len = getw(in);    len = getw(in);
         ul = (int *)MALLOC_ATOMIC(len*sizeof(int));    ul = (int *)MALLOC_ATOMIC(len*sizeof(int));
         u = (ent **)MALLOC(len*sizeof(ent *));    u = (ent **)MALLOC(len*sizeof(ent *));
         ll = (int *)MALLOC_ATOMIC(len*sizeof(int));    ll = (int *)MALLOC_ATOMIC(len*sizeof(int));
         l = (ent **)MALLOC(len*sizeof(ent *));    l = (ent **)MALLOC(len*sizeof(ent *));
         for ( i = 0; i < len; i++ ) {    for ( i = 0; i < len; i++ ) {
                 u[i] = get_row(in,&ul[i]);      u[i] = get_row(in,&ul[i]);
         }    }
         for ( i = 0; i < len; i++ ) {    for ( i = 0; i < len; i++ ) {
                 l[i] = get_row(in,&ll[i]);      l[i] = get_row(in,&ll[i]);
         }    }
         fclose(in);    fclose(in);
         *rp = ONE;    *rp = ONE;
 }  }
   
 void Plusolve_main(NODE arg,VECT *rp)  void Plusolve_main(NODE arg,VECT *rp)
 {  {
         Q *d,*p;    Q *d,*p;
         VECT v,r;    VECT v,r;
         int len,i;    int len,i;
         int *rhs;    int *rhs;
   
         v = (VECT)ARG0(arg); len = v->len;    v = (VECT)ARG0(arg); len = v->len;
         d = (Q *)BDY(v);    d = (Q *)BDY(v);
         rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));    rhs = (int *)MALLOC_ATOMIC(len*sizeof(int));
         for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]);    for ( i = 0; i < len; i++ ) rhs[i] = QTOS(d[i]);
         solve_l(ll,l,len,rhs,modulus);    solve_l(ll,l,len,rhs,modulus);
         solve_u(ul,u,len,rhs,modulus);    solve_u(ul,u,len,rhs,modulus);
         NEWVECT(r); r->len = len;    NEWVECT(r); r->len = len;
         r->body = (pointer *)MALLOC(len*sizeof(pointer));    r->body = (pointer *)MALLOC(len*sizeof(pointer));
         p = (Q *)r->body;    p = (Q *)r->body;
         for ( i = 0; i < len; i++ )    for ( i = 0; i < len; i++ )
                 STOQ(rhs[i],p[i]);      STOQ(rhs[i],p[i]);
         *rp = r;    *rp = r;
 }  }
   
 ent *get_row(FILE *in,int *l)  ent *get_row(FILE *in,int *l)
 {  {
         int len,i;    int len,i;
         ent *a;    ent *a;
   
         *l = len = getw(in);    *l = len = getw(in);
         a = (ent *)MALLOC_ATOMIC(len*sizeof(ent));    a = (ent *)MALLOC_ATOMIC(len*sizeof(ent));
         for ( i = 0; i < len; i++ ) {    for ( i = 0; i < len; i++ ) {
                 a[i].j = getw(in);      a[i].j = getw(in);
                 a[i].e = getw(in);      a[i].e = getw(in);
         }    }
         return a;    return a;
 }  }
   
 void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)  void lu_gauss(int *ul,ent **u,int *ll,ent **l,int n,int mod)
 {  {
         int i,j,k,s,mul;    int i,j,k,s,mul;
         unsigned int inv;    unsigned int inv;
         int *ll2;    int *ll2;
   
         ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int));    ll2 = (int *)MALLOC_ATOMIC(n*sizeof(int));
         for ( i = 0; i < n; i++ ) ll2[i] = 0;    for ( i = 0; i < n; i++ ) ll2[i] = 0;
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 fprintf(stderr,"i=%d\n",i);      fprintf(stderr,"i=%d\n",i);
                 inv = invm(u[i][0].e,mod);      inv = invm(u[i][0].e,mod);
                 for ( k = i+1; k < n; k++ )      for ( k = i+1; k < n; k++ )
                         if ( u[k][0].j == n-i ) {        if ( u[k][0].j == n-i ) {
                                 s = u[k][0].e;          s = u[k][0].e;
                                 DMAR(s,inv,0,mod,mul);          DMAR(s,inv,0,mod,mul);
                                 lu_elim(ul,u,k,i,mul,mod);          lu_elim(ul,u,k,i,mul,mod);
                                 lu_append(ll,l,ll2,k,i,mul);          lu_append(ll,l,ll2,k,i,mul);
                         }        }
         }    }
 }  }
   
 #define INITLEN 10  #define INITLEN 10
   
 void lu_append(int *l,ent **a,int *l2,int k,int i,int mul)  void lu_append(int *l,ent **a,int *l2,int k,int i,int mul)
 {  {
         int len;    int len;
         ent *p;    ent *p;
   
         len = l[k];    len = l[k];
         if ( !len ) {    if ( !len ) {
                 a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent));      a[k] = p = (ent *)MALLOC_ATOMIC(INITLEN*sizeof(ent));
                 p[0].j = i; p[0].e = mul;      p[0].j = i; p[0].e = mul;
                 l[k] = 1; l2[k] = INITLEN;      l[k] = 1; l2[k] = INITLEN;
         } else {    } else {
                 if ( l2[k] == l[k] ) {      if ( l2[k] == l[k] ) {
                         l2[k] *= 2;        l2[k] *= 2;
                         a[k] = REALLOC(a[k],l2[k]*sizeof(ent));        a[k] = REALLOC(a[k],l2[k]*sizeof(ent));
                 }      }
                 p =a[k];      p =a[k];
                 p[l[k]].j = i; p[l[k]].e = mul;      p[l[k]].j = i; p[l[k]].e = mul;
                 l[k]++;      l[k]++;
         }    }
 }  }
   
 /* a[k] = a[k]-mul*a[i] */  /* a[k] = a[k]-mul*a[i] */
   
 void lu_elim(int *l,ent **a,int k,int i,int mul,int mod)  void lu_elim(int *l,ent **a,int k,int i,int mul,int mod)
 {  {
         ent *ak,*ai,*w;    ent *ak,*ai,*w;
         int lk,li,j,m,p,q,r,s,t,j0;    int lk,li,j,m,p,q,r,s,t,j0;
   
         ak = a[k]; ai = a[i]; lk = l[k]; li = l[i];    ak = a[k]; ai = a[i]; lk = l[k]; li = l[i];
         w = (ent *)alloca((lk+li)*sizeof(ent));    w = (ent *)alloca((lk+li)*sizeof(ent));
         p = 0; q = 0; j = 0;    p = 0; q = 0; j = 0;
         mul = mod-mul;    mul = mod-mul;
         while ( p < lk && q < li ) {    while ( p < lk && q < li ) {
                 if ( ak[p].j > ai[q].j ) {      if ( ak[p].j > ai[q].j ) {
                         w[j] = ak[p]; j++; p++;        w[j] = ak[p]; j++; p++;
                 } else if ( ak[p].j < ai[q].j ) {      } else if ( ak[p].j < ai[q].j ) {
                         w[j].j = ai[q].j;        w[j].j = ai[q].j;
                         t = ai[q].e;        t = ai[q].e;
                         DMAR(t,mul,0,mod,r);        DMAR(t,mul,0,mod,r);
                         w[j].e = r;        w[j].e = r;
                         j++; q++;        j++; q++;
                 } else {      } else {
                         t = ai[q].e; s = ak[p].e;        t = ai[q].e; s = ak[p].e;
                         DMAR(t,mul,s,mod,r);        DMAR(t,mul,s,mod,r);
                         if ( r ) {        if ( r ) {
                                 w[j].j = ai[q].j; w[j].e = r; j++;          w[j].j = ai[q].j; w[j].e = r; j++;
                         }        }
                         p++; q++;        p++; q++;
                 }      }
         }    }
         if ( q == li )    if ( q == li )
                 while ( p < lk ) {      while ( p < lk ) {
                         w[j] = ak[p]; j++; p++;        w[j] = ak[p]; j++; p++;
                 }      }
         else if ( p == lk )    else if ( p == lk )
                 while ( q < li ) {      while ( q < li ) {
                         w[j].j = ai[q].j;        w[j].j = ai[q].j;
                         t = ai[q].e;        t = ai[q].e;
                         DMAR(t,mul,0,mod,r);        DMAR(t,mul,0,mod,r);
                         w[j].e = r;        w[j].e = r;
                         j++; q++;        j++; q++;
                 }      }
         if ( j <= lk ) {    if ( j <= lk ) {
                 for ( m = 0; m < j; m++ ) ak[m] = w[m];      for ( m = 0; m < j; m++ ) ak[m] = w[m];
         } else {    } else {
                 a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent));      a[k] = ak = (ent *)MALLOC_ATOMIC(j*sizeof(ent));
                 for ( m = 0; m < j; m++ ) ak[m] = w[m];      for ( m = 0; m < j; m++ ) ak[m] = w[m];
         }    }
         l[k] = j;    l[k] = j;
 }  }
   
 void solve_l(int *ll,ent **l,int n,int *rhs,int mod)  void solve_l(int *ll,ent **l,int n,int *rhs,int mod)
 {  {
         int j,k,s,len;    int j,k,s,len;
         ent *p;    ent *p;
   
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 len = ll[j]; p = l[j];      len = ll[j]; p = l[j];
                 for ( k = 0, s = 0; k < len; k++ )      for ( k = 0, s = 0; k < len; k++ )
                         s = dmar(p[k].e,rhs[p[k].j],s,mod);        s = dmar(p[k].e,rhs[p[k].j],s,mod);
                 rhs[j] -=  s;      rhs[j] -=  s;
                 if ( rhs[j] < 0 ) rhs[j] += mod;      if ( rhs[j] < 0 ) rhs[j] += mod;
         }    }
 }  }
   
 void solve_u(int *ul,ent **u,int n,int *rhs,int mod)  void solve_u(int *ul,ent **u,int n,int *rhs,int mod)
 {  {
         int j,k,s,len,inv;    int j,k,s,len,inv;
         ent *p;    ent *p;
   
         for ( j = n-1; j >= 0; j-- ) {    for ( j = n-1; j >= 0; j-- ) {
                 len = ul[j]; p = u[j];      len = ul[j]; p = u[j];
                 for ( k = 1, s = 0; k < len; k++ )      for ( k = 1, s = 0; k < len; k++ )
                         s = dmar(p[k].e,rhs[p[k].j],s,mod);        s = dmar(p[k].e,rhs[p[k].j],s,mod);
                 rhs[j] -=  s;      rhs[j] -=  s;
                 if ( rhs[j] < 0 ) rhs[j] += mod;      if ( rhs[j] < 0 ) rhs[j] += mod;
                 inv = invm((unsigned int)p[0].e,mod);      inv = invm((unsigned int)p[0].e,mod);
                 rhs[j] = dmar(rhs[j],inv,0,mod);      rhs[j] = dmar(rhs[j],inv,0,mod);
         }    }
 }  }
   
 int comp_obj(Obj *a,Obj *b)  int comp_obj(Obj *a,Obj *b)
 {  {
         return arf_comp(CO,*a,*b);    return arf_comp(CO,*a,*b);
 }  }
   
 static FUNC generic_comp_obj_func;  static FUNC generic_comp_obj_func;
Line 364  static NODE generic_comp_obj_option;
Line 364  static NODE generic_comp_obj_option;
   
 int generic_comp_obj(Obj *a,Obj *b)  int generic_comp_obj(Obj *a,Obj *b)
 {  {
         Q r;    Q r;
   
         BDY(generic_comp_obj_arg)=(pointer)(*a);    BDY(generic_comp_obj_arg)=(pointer)(*a);
         BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);    BDY(NEXT(generic_comp_obj_arg))=(pointer)(*b);
         r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option);    r = (Q)bevalf_with_opts(generic_comp_obj_func,generic_comp_obj_arg,generic_comp_obj_option);
         if ( !r )    if ( !r )
                 return 0;      return 0;
         else    else
                 return SGN(r)>0?1:-1;      return SGN(r)>0?1:-1;
 }  }
   
   
 void Pqsort(NODE arg,LIST *rp)  void Pqsort(NODE arg,LIST *rp)
 {  {
         VECT vect;    VECT vect;
         NODE n,n1;    NODE n,n1;
         P p;    P p;
         V v;    V v;
         FUNC func;    FUNC func;
         int len,i;    int len,i;
         pointer *a;    pointer *a;
         Obj t;    Obj t;
   
         t = ARG0(arg);    t = ARG0(arg);
     if (OID(t) == O_LIST) {      if (OID(t) == O_LIST) {
         n = (NODE)BDY((LIST)t);          n = (NODE)BDY((LIST)t);
         len = length(n);          len = length(n);
Line 401  void Pqsort(NODE arg,LIST *rp)
Line 401  void Pqsort(NODE arg,LIST *rp)
     }else {      }else {
         vect = (VECT)t;          vect = (VECT)t;
     }      }
         if ( argc(arg) == 1 )    if ( argc(arg) == 1 )
                 qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);      qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))comp_obj);
         else {    else {
                 p = (P)ARG1(arg);      p = (P)ARG1(arg);
                 if ( !p || OID(p)!=2 )      if ( !p || OID(p)!=2 )
                         error("qsort : invalid argument");        error("qsort : invalid argument");
                 v = VR(p);      v = VR(p);
                 gen_searchf(NAME(v),&func);      gen_searchf(NAME(v),&func);
                 if ( !func ) {      if ( !func ) {
                         if ( (int)v->attr != V_SR )        if ( (int)v->attr != V_SR )
                                 error("qsort : no such function");          error("qsort : no such function");
                         func = (FUNC)v->priv;        func = (FUNC)v->priv;
                 }      }
                 generic_comp_obj_func = func;      generic_comp_obj_func = func;
                 MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);      MKNODE(n,0,0); MKNODE(generic_comp_obj_arg,0,n);
                 generic_comp_obj_option = current_option;      generic_comp_obj_option = current_option;
                 qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);      qsort(BDY(vect),vect->len,sizeof(Obj),(int (*)(const void *,const void *))generic_comp_obj);
         }    }
     if (OID(t) == O_LIST) {      if (OID(t) == O_LIST) {
         a = BDY(vect);          a = BDY(vect);
         for ( i = len - 1, n = 0; i >= 0; i-- ) {          for ( i = len - 1, n = 0; i >= 0; i-- ) {
Line 432  void Pqsort(NODE arg,LIST *rp)
Line 432  void Pqsort(NODE arg,LIST *rp)
   
 void PNBmul_gf2n(NODE arg,GF2N *rp)  void PNBmul_gf2n(NODE arg,GF2N *rp)
 {  {
         GF2N a,b;    GF2N a,b;
         GF2MAT mat;    GF2MAT mat;
         int n,w;    int n,w;
         unsigned int *ab,*bb;    unsigned int *ab,*bb;
         UP2 r;    UP2 r;
   
         a = (GF2N)ARG0(arg);    a = (GF2N)ARG0(arg);
         b = (GF2N)ARG1(arg);    b = (GF2N)ARG1(arg);
         mat = (GF2MAT)ARG2(arg);    mat = (GF2MAT)ARG2(arg);
         if ( !a || !b )    if ( !a || !b )
                 *rp = 0;      *rp = 0;
         else {    else {
                 n = mat->row;      n = mat->row;
                 w = (n+BSH-1)/BSH;      w = (n+BSH-1)/BSH;
   
                 ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));      ab = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
                 bzero((char *)ab,w*sizeof(unsigned int));      bzero((char *)ab,w*sizeof(unsigned int));
                 bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));      bcopy(a->body->b,ab,(a->body->w)*sizeof(unsigned int));
   
                 bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));      bb = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
                 bzero((char *)bb,w*sizeof(unsigned int));      bzero((char *)bb,w*sizeof(unsigned int));
                 bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));      bcopy(b->body->b,bb,(b->body->w)*sizeof(unsigned int));
   
                 NEWUP2(r,w);      NEWUP2(r,w);
                 bzero((char *)r->b,w*sizeof(unsigned int));      bzero((char *)r->b,w*sizeof(unsigned int));
                 mul_nb(mat,ab,bb,r->b);      mul_nb(mat,ab,bb,r->b);
                 r->w = w;      r->w = w;
                 _adjup2(r);      _adjup2(r);
                 if ( !r->w )      if ( !r->w )
                         *rp = 0;        *rp = 0;
                 else      else
                         MKGF2N(r,*rp);        MKGF2N(r,*rp);
         }    }
 }  }
   
 void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)  void Pmul_vect_mat_gf2n(NODE arg,GF2N *rp)
 {  {
         GF2N a;    GF2N a;
         GF2MAT mat;    GF2MAT mat;
         int n,w;    int n,w;
         unsigned int *b;    unsigned int *b;
         UP2 r;    UP2 r;
   
         a = (GF2N)ARG0(arg);    a = (GF2N)ARG0(arg);
         mat = (GF2MAT)ARG1(arg);    mat = (GF2MAT)ARG1(arg);
         if ( !a )    if ( !a )
                 *rp = 0;      *rp = 0;
         else {    else {
                 n = mat->row;      n = mat->row;
                 w = (n+BSH-1)/BSH;      w = (n+BSH-1)/BSH;
                 b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));      b = (unsigned int *)ALLOCA(w*sizeof(unsigned int));
                 bzero((char *)b,w*sizeof(unsigned int));      bzero((char *)b,w*sizeof(unsigned int));
                 bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));      bcopy(a->body->b,b,(a->body->w)*sizeof(unsigned int));
                 NEWUP2(r,w);      NEWUP2(r,w);
                 bzero((char *)r->b,w*sizeof(unsigned int));      bzero((char *)r->b,w*sizeof(unsigned int));
                 mulgf2vectmat(mat->row,b,mat->body,r->b);      mulgf2vectmat(mat->row,b,mat->body,r->b);
                 r->w = w;      r->w = w;
                 _adjup2(r);      _adjup2(r);
                 if ( !r->w )      if ( !r->w )
                         *rp = 0;        *rp = 0;
                 else {      else {
                         MKGF2N(r,*rp);        MKGF2N(r,*rp);
                 }      }
         }    }
 }  }
   
 void Pbconvmat_gf2n(NODE arg,LIST *rp)  void Pbconvmat_gf2n(NODE arg,LIST *rp)
 {  {
         P p0,p1;    P p0,p1;
         int to;    int to;
         GF2MAT p01,p10;    GF2MAT p01,p10;
         GF2N root;    GF2N root;
         NODE n0,n1;    NODE n0,n1;
   
         p0 = (P)ARG0(arg);    p0 = (P)ARG0(arg);
         p1 = (P)ARG1(arg);    p1 = (P)ARG1(arg);
         to = ARG2(arg)?1:0;    to = ARG2(arg)?1:0;
         if ( argc(arg) == 4 ) {    if ( argc(arg) == 4 ) {
                 root = (GF2N)ARG3(arg);      root = (GF2N)ARG3(arg);
                 compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);      compute_change_of_basis_matrix_with_root(p0,p1,to,root,&p01,&p10);
         } else    } else
                 compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);      compute_change_of_basis_matrix(p0,p1,to,&p01,&p10);
         MKNODE(n1,p10,0); MKNODE(n0,p01,n1);    MKNODE(n1,p10,0); MKNODE(n0,p01,n1);
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void Pmulmat_gf2n(NODE arg,GF2MAT *rp)  void Pmulmat_gf2n(NODE arg,GF2MAT *rp)
 {  {
         GF2MAT m;    GF2MAT m;
   
         if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )    if ( !compute_multiplication_matrix((P)ARG0(arg),&m) )
                 error("mulmat_gf2n : input is not a normal polynomial");      error("mulmat_gf2n : input is not a normal polynomial");
         *rp = m;    *rp = m;
 }  }
   
 void Psepmat_destructive(NODE arg,LIST *rp)  void Psepmat_destructive(NODE arg,LIST *rp)
 {  {
         MAT mat,mat1;    MAT mat,mat1;
         int i,j,row,col;    int i,j,row,col;
         Q **a,**a1;    Q **a,**a1;
         Q ent;    Q ent;
         N nm,mod,rem,quo;    N nm,mod,rem,quo;
         int sgn;    int sgn;
         NODE n0,n1;    NODE n0,n1;
   
         mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));    mat = (MAT)ARG0(arg); mod = NM((Q)ARG1(arg));
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         MKMAT(mat1,row,col);    MKMAT(mat1,row,col);
         a = (Q **)mat->body; a1 = (Q **)mat1->body;    a = (Q **)mat->body; a1 = (Q **)mat1->body;
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0; j < col; j++ ) {      for ( j = 0; j < col; j++ ) {
                         ent = a[i][j];        ent = a[i][j];
                         if ( !ent )        if ( !ent )
                                 continue;          continue;
                         nm = NM(ent);        nm = NM(ent);
                         sgn = SGN(ent);        sgn = SGN(ent);
                         divn(nm,mod,&quo,&rem);        divn(nm,mod,&quo,&rem);
 /*                      if ( quo != nm && rem != nm ) */  /*      if ( quo != nm && rem != nm ) */
 /*                              GCFREE(nm); */  /*        GCFREE(nm); */
 /*                      GCFREE(ent); */  /*      GCFREE(ent); */
                         NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);        NTOQ(rem,sgn,a[i][j]); NTOQ(quo,sgn,a1[i][j]);
                 }      }
         MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);    MKNODE(n1,mat1,0); MKNODE(n0,mat,n1);
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void Psepvect(NODE arg,VECT *rp)  void Psepvect(NODE arg,VECT *rp)
 {  {
         sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);    sepvect((VECT)ARG0(arg),QTOS((Q)ARG1(arg)),rp);
 }  }
   
 void sepvect(VECT v,int d,VECT *rp)  void sepvect(VECT v,int d,VECT *rp)
 {  {
         int i,j,k,n,q,q1,r;    int i,j,k,n,q,q1,r;
         pointer *pv,*pw,*pu;    pointer *pv,*pw,*pu;
         VECT w,u;    VECT w,u;
   
         n = v->len;    n = v->len;
         if ( d > n )    if ( d > n )
                 d = n;      d = n;
         q = n/d; r = n%d; q1 = q+1;    q = n/d; r = n%d; q1 = q+1;
         MKVECT(w,d); *rp = w;    MKVECT(w,d); *rp = w;
         pv = BDY(v); pw = BDY(w); k = 0;    pv = BDY(v); pw = BDY(w); k = 0;
         for ( i = 0; i < r; i++ ) {    for ( i = 0; i < r; i++ ) {
                 MKVECT(u,q1); pw[i] = (pointer)u;      MKVECT(u,q1); pw[i] = (pointer)u;
                 for ( pu = BDY(u), j = 0; j < q1; j++, k++ )      for ( pu = BDY(u), j = 0; j < q1; j++, k++ )
                         pu[j] = pv[k];        pu[j] = pv[k];
         }    }
         for ( ; i < d; i++ ) {    for ( ; i < d; i++ ) {
                 MKVECT(u,q); pw[i] = (pointer)u;      MKVECT(u,q); pw[i] = (pointer)u;
                 for ( pu = BDY(u), j = 0; j < q; j++, k++ )      for ( pu = BDY(u), j = 0; j < q; j++, k++ )
                         pu[j] = pv[k];        pu[j] = pv[k];
         }    }
 }  }
   
 void Pnewvect(NODE arg,VECT *rp)  void Pnewvect(NODE arg,VECT *rp)
 {  {
         int len,i,r;    int len,i,r;
         VECT vect;    VECT vect;
         pointer *vb;    pointer *vb;
         LIST list;    LIST list;
         NODE tn;    NODE tn;
   
         asir_assert(ARG0(arg),O_N,"newvect");    asir_assert(ARG0(arg),O_N,"newvect");
         len = QTOS((Q)ARG0(arg));    len = QTOS((Q)ARG0(arg));
         if ( len < 0 )    if ( len < 0 )
                 error("newvect : invalid size");      error("newvect : invalid size");
         MKVECT(vect,len);    MKVECT(vect,len);
         if ( argc(arg) == 2 ) {    if ( argc(arg) == 2 ) {
                 list = (LIST)ARG1(arg);      list = (LIST)ARG1(arg);
                 asir_assert(list,O_LIST,"newvect");      asir_assert(list,O_LIST,"newvect");
 #if 0  #if 0
                 for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );      for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
                 if ( r > len ) {      if ( r > len ) {
                         *rp = vect;        *rp = vect;
                         return;        return;
                 }      }
 #endif  #endif
                 for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )      for ( i = 0, tn = BDY(list), vb = BDY(vect); tn; i++, tn = NEXT(tn) )
                         vb[i] = (pointer)BDY(tn);        vb[i] = (pointer)BDY(tn);
         }    }
         *rp = vect;    *rp = vect;
 }  }
   
 void Pvect(NODE arg,VECT *rp) {  void Pvect(NODE arg,VECT *rp) {
         int len,i;    int len,i;
         VECT vect;    VECT vect;
         pointer *vb;    pointer *vb;
         NODE tn;    NODE tn;
   
         if ( !arg ) {    if ( !arg ) {
                 *rp =0;      *rp =0;
                 return;      return;
         }    }
   
         for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);    for (len = 0, tn = arg; tn; tn = NEXT(tn), len++);
         if ( len == 1 ) {    if ( len == 1 ) {
                 if ( ARG0(arg) != 0 ) {      if ( ARG0(arg) != 0 ) {
                         switch ( OID(ARG0(arg)) ) {        switch ( OID(ARG0(arg)) ) {
                                 case O_VECT:          case O_VECT:
                                         *rp = ARG0(arg);            *rp = ARG0(arg);
                                         return;            return;
                                 case O_LIST:          case O_LIST:
                                         for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );            for ( len = 0, tn = ARG0(arg); tn; tn = NEXT(tn), len++ );
                                         MKVECT(vect,len-1);            MKVECT(vect,len-1);
                                         for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);            for ( i = 0, tn = BDY((LIST)ARG0(arg)), vb =BDY(vect);
                                                         tn; i++, tn = NEXT(tn) )                tn; i++, tn = NEXT(tn) )
                                                 vb[i] = (pointer)BDY(tn);              vb[i] = (pointer)BDY(tn);
                                         *rp=vect;            *rp=vect;
                                         return;            return;
                         }        }
                 }      }
         }    }
         MKVECT(vect,len);    MKVECT(vect,len);
         for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )    for ( i = 0, tn = arg, vb = BDY(vect); tn; i++, tn = NEXT(tn) )
                 vb[i] = (pointer)BDY(tn);      vb[i] = (pointer)BDY(tn);
         *rp = vect;    *rp = vect;
 }  }
   
 void Pexponent_vector(NODE arg,DP *rp)  void Pexponent_vector(NODE arg,DP *rp)
 {  {
         nodetod(arg,rp);    nodetod(arg,rp);
 }  }
   
 void Pnewbytearray(NODE arg,BYTEARRAY *rp)  void Pnewbytearray(NODE arg,BYTEARRAY *rp)
 {  {
         int len,i,r;    int len,i,r;
         BYTEARRAY array;    BYTEARRAY array;
         unsigned char *vb;    unsigned char *vb;
         char *str;    char *str;
         LIST list;    LIST list;
         NODE tn;    NODE tn;
         int ac;    int ac;
         struct stat sbuf;    struct stat sbuf;
         char *fname;    char *fname;
         FILE *fp;    FILE *fp;
   
         ac = argc(arg);    ac = argc(arg);
         if ( ac == 1 ) {    if ( ac == 1 ) {
                 if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument");      if ( !OID((Obj)ARG0(arg)) ) error("newbytearray : invalid argument");
                 switch ( OID((Obj)ARG0(arg)) ) {      switch ( OID((Obj)ARG0(arg)) ) {
                         case O_STR:        case O_STR:
                                 fname = BDY((STRING)ARG0(arg));          fname = BDY((STRING)ARG0(arg));
                                 fp = fopen(fname,"rb");          fp = fopen(fname,"rb");
                                 if ( !fp ) error("newbytearray : fopen failed");          if ( !fp ) error("newbytearray : fopen failed");
                                 if ( stat(fname,&sbuf) < 0 )          if ( stat(fname,&sbuf) < 0 )
                                         error("newbytearray : stat failed");            error("newbytearray : stat failed");
                                 len = sbuf.st_size;          len = sbuf.st_size;
                                 MKBYTEARRAY(array,len);          MKBYTEARRAY(array,len);
                                 fread(BDY(array),len,sizeof(char),fp);          fread(BDY(array),len,sizeof(char),fp);
                                 break;          break;
                         case O_N:        case O_N:
                                 if ( !RATN(ARG0(arg)) )          if ( !RATN(ARG0(arg)) )
                                         error("newbytearray : invalid argument");            error("newbytearray : invalid argument");
                                 len = QTOS((Q)ARG0(arg));          len = QTOS((Q)ARG0(arg));
                                 if ( len < 0 )          if ( len < 0 )
                                         error("newbytearray : invalid size");            error("newbytearray : invalid size");
                                 MKBYTEARRAY(array,len);          MKBYTEARRAY(array,len);
                                 break;          break;
                         default:        default:
                                 error("newbytearray : invalid argument");          error("newbytearray : invalid argument");
                 }      }
         } else if ( ac == 2 ) {    } else if ( ac == 2 ) {
                 asir_assert(ARG0(arg),O_N,"newbytearray");      asir_assert(ARG0(arg),O_N,"newbytearray");
                 len = QTOS((Q)ARG0(arg));      len = QTOS((Q)ARG0(arg));
                 if ( len < 0 )      if ( len < 0 )
                         error("newbytearray : invalid size");        error("newbytearray : invalid size");
                 MKBYTEARRAY(array,len);      MKBYTEARRAY(array,len);
                 if ( !ARG1(arg) )      if ( !ARG1(arg) )
                         error("newbytearray : invalid initialization");        error("newbytearray : invalid initialization");
                 switch ( OID((Obj)ARG1(arg)) ) {      switch ( OID((Obj)ARG1(arg)) ) {
                         case O_LIST:        case O_LIST:
                                 list = (LIST)ARG1(arg);          list = (LIST)ARG1(arg);
                                 asir_assert(list,O_LIST,"newbytearray");          asir_assert(list,O_LIST,"newbytearray");
                                 for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );          for ( r = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) );
                                 if ( r <= len ) {          if ( r <= len ) {
                                         for ( i = 0, tn = BDY(list), vb = BDY(array); tn;            for ( i = 0, tn = BDY(list), vb = BDY(array); tn;
                                                 i++, tn = NEXT(tn) )              i++, tn = NEXT(tn) )
                                                 vb[i] = (unsigned char)QTOS((Q)BDY(tn));              vb[i] = (unsigned char)QTOS((Q)BDY(tn));
                                 }          }
                                 break;          break;
                         case O_STR:        case O_STR:
                                 str = BDY((STRING)ARG1(arg));          str = BDY((STRING)ARG1(arg));
                                 r = strlen(str);          r = strlen(str);
                                 if ( r <= len )          if ( r <= len )
                                         bcopy(str,BDY(array),r);            bcopy(str,BDY(array),r);
                                 break;          break;
                         default:        default:
                                 if ( !ARG1(arg) )          if ( !ARG1(arg) )
                                         error("newbytearray : invalid initialization");            error("newbytearray : invalid initialization");
                 }      }
         } else    } else
                 error("newbytearray : invalid argument");      error("newbytearray : invalid argument");
         *rp = array;    *rp = array;
 }  }
   
 #define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7)))  #define MEMORY_GETPOINT(a,len,x,y) (((a)[(len)*(y)+((x)>>3)])&(1<<((x)&7)))
   
 void Pmemoryplot_to_coord(NODE arg,LIST *rp)  void Pmemoryplot_to_coord(NODE arg,LIST *rp)
 {  {
         int len,blen,y,i,j;    int len,blen,y,i,j;
         unsigned char *a;    unsigned char *a;
         NODE r0,r,n;    NODE r0,r,n;
         LIST l;    LIST l;
         BYTEARRAY ba;    BYTEARRAY ba;
         Q iq,jq;    Q iq,jq;
   
         asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");    asir_assert(ARG0(arg),O_LIST,"memoryplot_to_coord");
         arg = BDY((LIST)ARG0(arg));    arg = BDY((LIST)ARG0(arg));
         len = QTOS((Q)ARG0(arg));    len = QTOS((Q)ARG0(arg));
         blen = (len+7)/8;    blen = (len+7)/8;
         y = QTOS((Q)ARG1(arg));    y = QTOS((Q)ARG1(arg));
         ba = (BYTEARRAY)ARG2(arg); a = ba->body;    ba = (BYTEARRAY)ARG2(arg); a = ba->body;
         r0 = 0;    r0 = 0;
         for ( j = 0; j < y; j++ )    for ( j = 0; j < y; j++ )
                 for ( i = 0; i < len; i++ )      for ( i = 0; i < len; i++ )
                         if ( MEMORY_GETPOINT(a,blen,i,j) ) {        if ( MEMORY_GETPOINT(a,blen,i,j) ) {
                                 NEXTNODE(r0,r);          NEXTNODE(r0,r);
                                 STOQ(i,iq); STOQ(j,jq);          STOQ(i,iq); STOQ(j,jq);
                                 n = mknode(2,iq,jq);          n = mknode(2,iq,jq);
                                 MKLIST(l,n);          MKLIST(l,n);
                                 BDY(r) = l;          BDY(r) = l;
                         }        }
         if ( r0 ) NEXT(r) = 0;    if ( r0 ) NEXT(r) = 0;
         MKLIST(*rp,r0);    MKLIST(*rp,r0);
 }  }
   
 void Pnewmat(NODE arg,MAT *rp)  void Pnewmat(NODE arg,MAT *rp)
 {  {
         int row,col;    int row,col;
         int i,j,r,c;    int i,j,r,c;
         NODE tn,sn;    NODE tn,sn;
         MAT m;    MAT m;
         pointer **mb;    pointer **mb;
         LIST list;    LIST list;
   
         asir_assert(ARG0(arg),O_N,"newmat");    asir_assert(ARG0(arg),O_N,"newmat");
         asir_assert(ARG1(arg),O_N,"newmat");    asir_assert(ARG1(arg),O_N,"newmat");
         row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));    row = QTOS((Q)ARG0(arg)); col = QTOS((Q)ARG1(arg));
         if ( row < 0 || col < 0 )    if ( row < 0 || col < 0 )
                 error("newmat : invalid size");      error("newmat : invalid size");
         MKMAT(m,row,col);    MKMAT(m,row,col);
         if ( argc(arg) == 3 ) {    if ( argc(arg) == 3 ) {
                 list = (LIST)ARG2(arg);      list = (LIST)ARG2(arg);
                 asir_assert(list,O_LIST,"newmat");      asir_assert(list,O_LIST,"newmat");
                 for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {      for ( r = 0, c = 0, tn = BDY(list); tn; r++, tn = NEXT(tn) ) {
                         for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );        for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) );
                         c = MAX(c,j);        c = MAX(c,j);
                 }      }
                 if ( (r > row) || (c > col) ) {      if ( (r > row) || (c > col) ) {
                         *rp = m;        *rp = m;
                         return;        return;
                 }      }
                 for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {      for ( i = 0, tn = BDY(list), mb = BDY(m); tn; i++, tn = NEXT(tn) ) {
                         asir_assert(BDY(tn),O_LIST,"newmat");        asir_assert(BDY(tn),O_LIST,"newmat");
                         for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )        for ( j = 0, sn = BDY((LIST)BDY(tn)); sn; j++, sn = NEXT(sn) )
                                 mb[i][j] = (pointer)BDY(sn);          mb[i][j] = (pointer)BDY(sn);
                 }      }
         }    }
         *rp = m;    *rp = m;
 }  }
   
 void Pmat(NODE arg, MAT *rp)  void Pmat(NODE arg, MAT *rp)
 {  {
         int row,col;    int row,col;
         int i;    int i;
         MAT m;    MAT m;
         pointer **mb;    pointer **mb;
         pointer *ent;    pointer *ent;
         NODE tn, sn;    NODE tn, sn;
         VECT v;    VECT v;
   
         if ( !arg ) {    if ( !arg ) {
                 *rp =0;      *rp =0;
                 return;      return;
         }    }
   
         for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);    for (row = 0, tn = arg; tn; tn = NEXT(tn), row++);
         if ( row == 1 ) {    if ( row == 1 ) {
                 if ( OID(ARG0(arg)) == O_MAT ) {      if ( OID(ARG0(arg)) == O_MAT ) {
                         *rp=ARG0(arg);        *rp=ARG0(arg);
                         return;        return;
                 } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {      } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
                         error("mat : invalid argument");        error("mat : invalid argument");
                 }      }
         }    }
         if ( OID(ARG0(arg)) == O_VECT ) {    if ( OID(ARG0(arg)) == O_VECT ) {
                 v = ARG0(arg);      v = ARG0(arg);
                 col = v->len;      col = v->len;
         } else if ( OID(ARG0(arg)) == O_LIST ) {    } else if ( OID(ARG0(arg)) == O_LIST ) {
                 for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);      for (col = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), col++);
         } else {    } else {
                 error("mat : invalid argument");      error("mat : invalid argument");
         }    }
   
         MKMAT(m,row,col);    MKMAT(m,row,col);
         for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {    for (row = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), row++) {
                 if ( BDY(tn) == 0 ) {      if ( BDY(tn) == 0 ) {
                         error("mat : invalid argument");        error("mat : invalid argument");
                 } else if ( OID(BDY(tn)) == O_VECT ) {      } else if ( OID(BDY(tn)) == O_VECT ) {
                         v = tn->body;        v = tn->body;
                         ent = BDY(v);        ent = BDY(v);
                         for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];        for (i = 0; i < v->len; i++ ) mb[row][i] = (Obj)ent[i];
                 } else if ( OID(BDY(tn)) == O_LIST ) {      } else if ( OID(BDY(tn)) == O_LIST ) {
                         for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )        for (col = 0, sn = BDY((LIST)BDY(tn)); sn; col++, sn = NEXT(sn) )
                                 mb[row][col] = (pointer)BDY(sn);          mb[row][col] = (pointer)BDY(sn);
                 } else {      } else {
                         error("mat : invalid argument");        error("mat : invalid argument");
                 }      }
         }    }
         *rp = m;    *rp = m;
 }  }
   
 void Pmatc(NODE arg, MAT *rp)  void Pmatc(NODE arg, MAT *rp)
 {  {
         int row,col;    int row,col;
         int i;    int i;
         MAT m;    MAT m;
         pointer **mb;    pointer **mb;
         pointer *ent;    pointer *ent;
         NODE tn, sn;    NODE tn, sn;
         VECT v;    VECT v;
   
         if ( !arg ) {    if ( !arg ) {
                 *rp =0;      *rp =0;
                 return;      return;
         }    }
   
         for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);    for (col = 0, tn = arg; tn; tn = NEXT(tn), col++);
         if ( col == 1 ) {    if ( col == 1 ) {
                 if ( OID(ARG0(arg)) == O_MAT ) {      if ( OID(ARG0(arg)) == O_MAT ) {
                         *rp=ARG0(arg);        *rp=ARG0(arg);
                         return;        return;
                 } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {      } else if ( !(OID(ARG0(arg)) == O_LIST || OID(ARG0(arg)) == O_VECT)) {
                         error("matc : invalid argument");        error("matc : invalid argument");
                 }      }
         }    }
         if ( OID(ARG0(arg)) == O_VECT ) {    if ( OID(ARG0(arg)) == O_VECT ) {
                 v = ARG0(arg);      v = ARG0(arg);
                 row = v->len;      row = v->len;
         } else if ( OID(ARG0(arg)) == O_LIST ) {    } else if ( OID(ARG0(arg)) == O_LIST ) {
                 for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);      for (row = 0, tn = BDY((LIST)ARG0(arg)); tn ; tn = NEXT(tn), row++);
         } else {    } else {
                 error("matc : invalid argument");      error("matc : invalid argument");
         }    }
   
         MKMAT(m,row,col);    MKMAT(m,row,col);
         for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {    for (col = 0, tn = arg, mb = BDY(m); tn; tn = NEXT(tn), col++) {
                 if ( BDY(tn) == 0 ) {      if ( BDY(tn) == 0 ) {
                         error("matc : invalid argument");        error("matc : invalid argument");
                 } else if ( OID(BDY(tn)) == O_VECT ) {      } else if ( OID(BDY(tn)) == O_VECT ) {
                         v = tn->body;        v = tn->body;
                         ent = BDY(v);        ent = BDY(v);
                         for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];        for (i = 0; i < v->len; i++ ) mb[i][col] = (Obj)ent[i];
                 } else if ( OID(BDY(tn)) == O_LIST ) {      } else if ( OID(BDY(tn)) == O_LIST ) {
                         for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )        for (row = 0, sn = BDY((LIST)BDY(tn)); sn; row++, sn = NEXT(sn) )
                                 mb[row][col] = (pointer)BDY(sn);          mb[row][col] = (pointer)BDY(sn);
                 } else {      } else {
                         error("matc : invalid argument");        error("matc : invalid argument");
                 }      }
         }    }
         *rp = m;    *rp = m;
 }  }
   
 void Pvtol(NODE arg,LIST *rp)  void Pvtol(NODE arg,LIST *rp)
 {  {
         NODE n,n1;    NODE n,n1;
         VECT v;    VECT v;
         pointer *a;    pointer *a;
         int len,i;    int len,i;
   
         if ( OID(ARG0(arg)) == O_LIST ) {    if ( OID(ARG0(arg)) == O_LIST ) {
                 *rp = ARG0(arg);      *rp = ARG0(arg);
                 return;      return;
         }    }
         asir_assert(ARG0(arg),O_VECT,"vtol");    asir_assert(ARG0(arg),O_VECT,"vtol");
         v = (VECT)ARG0(arg); len = v->len; a = BDY(v);    v = (VECT)ARG0(arg); len = v->len; a = BDY(v);
         for ( i = len - 1, n = 0; i >= 0; i-- ) {    for ( i = len - 1, n = 0; i >= 0; i-- ) {
                 MKNODE(n1,a[i],n); n = n1;      MKNODE(n1,a[i],n); n = n1;
         }    }
         MKLIST(*rp,n);    MKLIST(*rp,n);
 }  }
   
 void Pltov(NODE arg,VECT *rp)  void Pltov(NODE arg,VECT *rp)
 {  {
         NODE n;    NODE n;
         VECT v,v0;    VECT v,v0;
         int len,i;    int len,i;
   
         if ( OID(ARG0(arg)) == O_VECT ) {    if ( OID(ARG0(arg)) == O_VECT ) {
                 v0 = (VECT)ARG0(arg); len = v0->len;      v0 = (VECT)ARG0(arg); len = v0->len;
                 MKVECT(v,len);      MKVECT(v,len);
                 for ( i = 0; i < len; i++ ) {      for ( i = 0; i < len; i++ ) {
                         BDY(v)[i] = BDY(v0)[i];        BDY(v)[i] = BDY(v0)[i];
                 }      }
                 *rp = v;      *rp = v;
                 return;      return;
         }    }
         asir_assert(ARG0(arg),O_LIST,"ltov");    asir_assert(ARG0(arg),O_LIST,"ltov");
         n = (NODE)BDY((LIST)ARG0(arg));    n = (NODE)BDY((LIST)ARG0(arg));
         len = length(n);    len = length(n);
         MKVECT(v,len);    MKVECT(v,len);
         for ( i = 0; i < len; i++, n = NEXT(n) )    for ( i = 0; i < len; i++, n = NEXT(n) )
                 BDY(v)[i] = BDY(n);      BDY(v)[i] = BDY(n);
         *rp = v;    *rp = v;
 }  }
   
 void Premainder(NODE arg,Obj *rp)  void Premainder(NODE arg,Obj *rp)
 {  {
         Obj a;    Obj a;
         VECT v,w;    VECT v,w;
         MAT m,l;    MAT m,l;
         pointer *vb,*wb;    pointer *vb,*wb;
         pointer **mb,**lb;    pointer **mb,**lb;
         int id,i,j,n,row,col,t,smd,sgn;    int id,i,j,n,row,col,t,smd,sgn;
         Q md,q;    Q md,q;
   
         a = (Obj)ARG0(arg); md = (Q)ARG1(arg);    a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
         if ( !a )    if ( !a )
                 *rp = 0;      *rp = 0;
         else {    else {
                 id = OID(a);      id = OID(a);
                 switch ( id ) {      switch ( id ) {
                         case O_N:        case O_N:
                         case O_P:        case O_P:
                                 cmp(md,(P)a,(P *)rp); break;          cmp(md,(P)a,(P *)rp); break;
                         case O_VECT:        case O_VECT:
                                 smd = QTOS(md);          smd = QTOS(md);
                                 v = (VECT)a; n = v->len; vb = v->body;          v = (VECT)a; n = v->len; vb = v->body;
                                 MKVECT(w,n); wb = w->body;          MKVECT(w,n); wb = w->body;
                                 for ( i = 0; i < n; i++ ) {          for ( i = 0; i < n; i++ ) {
                                         if ( q = (Q)vb[i] ) {            if ( q = (Q)vb[i] ) {
                                                 sgn = SGN(q); t = rem(NM(q),smd);              sgn = SGN(q); t = rem(NM(q),smd);
                                                 STOQ(t,q);              STOQ(t,q);
                                                 if ( q )              if ( q )
                                                         SGN(q) = sgn;                SGN(q) = sgn;
                                         }            }
                                         wb[i] = (pointer)q;            wb[i] = (pointer)q;
                                 }          }
                                 *rp = (Obj)w;          *rp = (Obj)w;
                                 break;          break;
                         case O_MAT:        case O_MAT:
                                 m = (MAT)a; row = m->row; col = m->col; mb = m->body;          m = (MAT)a; row = m->row; col = m->col; mb = m->body;
                                 MKMAT(l,row,col); lb = l->body;          MKMAT(l,row,col); lb = l->body;
                                 for ( i = 0; i < row; i++ )          for ( i = 0; i < row; i++ )
                                         for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )            for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
                                                 cmp(md,(P)vb[j],(P *)&wb[j]);              cmp(md,(P)vb[j],(P *)&wb[j]);
                                 *rp = (Obj)l;          *rp = (Obj)l;
                                 break;          break;
                         default:        default:
                                 error("remainder : invalid argument");          error("remainder : invalid argument");
                 }      }
         }    }
 }  }
   
 void Psremainder(NODE arg,Obj *rp)  void Psremainder(NODE arg,Obj *rp)
 {  {
         Obj a;    Obj a;
         VECT v,w;    VECT v,w;
         MAT m,l;    MAT m,l;
         pointer *vb,*wb;    pointer *vb,*wb;
         pointer **mb,**lb;    pointer **mb,**lb;
         unsigned int t,smd;    unsigned int t,smd;
         int id,i,j,n,row,col;    int id,i,j,n,row,col;
         Q md,q;    Q md,q;
   
         a = (Obj)ARG0(arg); md = (Q)ARG1(arg);    a = (Obj)ARG0(arg); md = (Q)ARG1(arg);
         if ( !a )    if ( !a )
                 *rp = 0;      *rp = 0;
         else {    else {
                 id = OID(a);      id = OID(a);
                 switch ( id ) {      switch ( id ) {
                         case O_N:        case O_N:
                         case O_P:        case O_P:
                                 cmp(md,(P)a,(P *)rp); break;          cmp(md,(P)a,(P *)rp); break;
                         case O_VECT:        case O_VECT:
                                 smd = QTOS(md);          smd = QTOS(md);
                                 v = (VECT)a; n = v->len; vb = v->body;          v = (VECT)a; n = v->len; vb = v->body;
                                 MKVECT(w,n); wb = w->body;          MKVECT(w,n); wb = w->body;
                                 for ( i = 0; i < n; i++ ) {          for ( i = 0; i < n; i++ ) {
                                         if ( q = (Q)vb[i] ) {            if ( q = (Q)vb[i] ) {
                                                 t = (unsigned int)rem(NM(q),smd);              t = (unsigned int)rem(NM(q),smd);
                                                 if ( SGN(q) < 0 )              if ( SGN(q) < 0 )
                                                         t = (smd - t) % smd;                t = (smd - t) % smd;
                                                 UTOQ(t,q);              UTOQ(t,q);
                                         }            }
                                         wb[i] = (pointer)q;            wb[i] = (pointer)q;
                                 }          }
                                 *rp = (Obj)w;          *rp = (Obj)w;
                                 break;          break;
                         case O_MAT:        case O_MAT:
                                 m = (MAT)a; row = m->row; col = m->col; mb = m->body;          m = (MAT)a; row = m->row; col = m->col; mb = m->body;
                                 MKMAT(l,row,col); lb = l->body;          MKMAT(l,row,col); lb = l->body;
                                 for ( i = 0; i < row; i++ )          for ( i = 0; i < row; i++ )
                                         for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )            for ( j = 0, vb = mb[i], wb = lb[i]; j < col; j++ )
                                                 cmp(md,(P)vb[j],(P *)&wb[j]);              cmp(md,(P)vb[j],(P *)&wb[j]);
                                 *rp = (Obj)l;          *rp = (Obj)l;
                                 break;          break;
                         default:        default:
                                 error("remainder : invalid argument");          error("remainder : invalid argument");
                 }      }
         }    }
 }  }
   
 void Psize(NODE arg,LIST *rp)  void Psize(NODE arg,LIST *rp)
 {  {
   
         int n,m;    int n,m;
         Q q;    Q q;
         NODE t,s;    NODE t,s;
   
         if ( !ARG0(arg) )    if ( !ARG0(arg) )
                  t = 0;       t = 0;
         else {    else {
                 switch (OID(ARG0(arg))) {      switch (OID(ARG0(arg))) {
                         case O_VECT:        case O_VECT:
                                 n = ((VECT)ARG0(arg))->len;          n = ((VECT)ARG0(arg))->len;
                                 STOQ(n,q); MKNODE(t,q,0);          STOQ(n,q); MKNODE(t,q,0);
                                 break;          break;
                         case O_MAT:        case O_MAT:
                                 n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;          n = ((MAT)ARG0(arg))->row; m = ((MAT)ARG0(arg))->col;
                                 STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);          STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
                                 break;          break;
                         case O_IMAT:        case O_IMAT:
                                 n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;          n = ((IMAT)ARG0(arg))->row; m = ((IMAT)ARG0(arg))->col;
                                 STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);          STOQ(m,q); MKNODE(s,q,0); STOQ(n,q); MKNODE(t,q,s);
                                 break;          break;
                         default:        default:
                                 error("size : invalid argument"); break;          error("size : invalid argument"); break;
                 }      }
         }    }
         MKLIST(*rp,t);    MKLIST(*rp,t);
 }  }
   
 void Pdet(NODE arg,P *rp)  void Pdet(NODE arg,P *rp)
 {  {
         MAT m;    MAT m;
         int n,i,j,mod;    int n,i,j,mod;
         P d;    P d;
         P **mat,**w;    P **mat,**w;
   
         m = (MAT)ARG0(arg);    m = (MAT)ARG0(arg);
         asir_assert(m,O_MAT,"det");    asir_assert(m,O_MAT,"det");
         if ( m->row != m->col )    if ( m->row != m->col )
                 error("det : non-square matrix");      error("det : non-square matrix");
         else if ( argc(arg) == 1 )    else if ( argc(arg) == 1 )
                 detp(CO,(P **)BDY(m),m->row,rp);      detp(CO,(P **)BDY(m),m->row,rp);
         else {    else {
                 n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);      n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
                 w = (P **)almat_pointer(n,n);      w = (P **)almat_pointer(n,n);
                 for ( i = 0; i < n; i++ )      for ( i = 0; i < n; i++ )
                         for ( j = 0; j < n; j++ )        for ( j = 0; j < n; j++ )
                                 ptomp(mod,mat[i][j],&w[i][j]);          ptomp(mod,mat[i][j],&w[i][j]);
                 detmp(CO,mod,w,n,&d);      detmp(CO,mod,w,n,&d);
                 mptop(d,rp);      mptop(d,rp);
         }    }
 }  }
   
 void Pinvmat(NODE arg,LIST *rp)  void Pinvmat(NODE arg,LIST *rp)
 {  {
         MAT m,r;    MAT m,r;
         int n,i,j,mod;    int n,i,j,mod;
         P dn;    P dn;
         P **mat,**imat,**w;    P **mat,**imat,**w;
         NODE nd;    NODE nd;
   
         m = (MAT)ARG0(arg);    m = (MAT)ARG0(arg);
         asir_assert(m,O_MAT,"invmat");    asir_assert(m,O_MAT,"invmat");
         if ( m->row != m->col )    if ( m->row != m->col )
                 error("invmat : non-square matrix");      error("invmat : non-square matrix");
         else if ( argc(arg) == 1 ) {    else if ( argc(arg) == 1 ) {
                 n = m->row;      n = m->row;
                 invmatp(CO,(P **)BDY(m),n,&imat,&dn);      invmatp(CO,(P **)BDY(m),n,&imat,&dn);
                 NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;      NEWMAT(r); r->row = n; r->col = n; r->body = (pointer **)imat;
                 nd = mknode(2,r,dn);      nd = mknode(2,r,dn);
                 MKLIST(*rp,nd);      MKLIST(*rp,nd);
         } else {    } else {
                 n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);      n = m->row; mod = QTOS((Q)ARG1(arg)); mat = (P **)BDY(m);
                 w = (P **)almat_pointer(n,n);      w = (P **)almat_pointer(n,n);
                 for ( i = 0; i < n; i++ )      for ( i = 0; i < n; i++ )
                         for ( j = 0; j < n; j++ )        for ( j = 0; j < n; j++ )
                                 ptomp(mod,mat[i][j],&w[i][j]);          ptomp(mod,mat[i][j],&w[i][j]);
 #if 0  #if 0
                 detmp(CO,mod,w,n,&d);      detmp(CO,mod,w,n,&d);
                 mptop(d,rp);      mptop(d,rp);
 #else  #else
                 error("not implemented yet");      error("not implemented yet");
 #endif  #endif
         }    }
 }  }
   
 /*  /*
         input : a row x col matrix A    input : a row x col matrix A
                 A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...      A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
   
         output : [B,D,R,C]    output : [B,D,R,C]
                 B : a rank(A) x col-rank(A) matrix      B : a rank(A) x col-rank(A) matrix
                 D : the denominator      D : the denominator
                 R : a vector of length rank(A)      R : a vector of length rank(A)
                 C : a vector of length col-rank(A)      C : a vector of length col-rank(A)
                 B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...      B[I] <-> D*x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */  */
   
 void Pgeneric_gauss_elim(NODE arg,LIST *rp)  void Pgeneric_gauss_elim(NODE arg,LIST *rp)
 {  {
         NODE n0,opt,p;    NODE n0,opt,p;
         MAT m,nm;    MAT m,nm;
         int *ri,*ci;    int *ri,*ci;
         VECT rind,cind;    VECT rind,cind;
         Q dn,q;    Q dn,q;
         int i,row,col,t,rank;    int i,row,col,t,rank;
         int is_hensel = 0;    int is_hensel = 0;
         char *key;    char *key;
         Obj value;    Obj value;
   
         if ( current_option ) {    if ( current_option ) {
                 for ( opt = current_option; opt; opt = NEXT(opt) ) {      for ( opt = current_option; opt; opt = NEXT(opt) ) {
                         p = BDY((LIST)BDY(opt));        p = BDY((LIST)BDY(opt));
                         key = BDY((STRING)BDY(p));        key = BDY((STRING)BDY(p));
                         value = (Obj)BDY(NEXT(p));        value = (Obj)BDY(NEXT(p));
                         if ( !strcmp(key,"hensel") && value ) {        if ( !strcmp(key,"hensel") && value ) {
                                 is_hensel = value ? 1 : 0;          is_hensel = value ? 1 : 0;
                                 break;          break;
                         }        }
                 }      }
         }    }
         asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");    asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim");
         m = (MAT)ARG0(arg);    m = (MAT)ARG0(arg);
         row = m->row; col = m->col;    row = m->row; col = m->col;
         if ( is_hensel )    if ( is_hensel )
                 rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);      rank = generic_gauss_elim_hensel(m,&nm,&dn,&ri,&ci);
         else    else
                 rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);      rank = generic_gauss_elim(m,&nm,&dn,&ri,&ci);
         t = col-rank;    t = col-rank;
         MKVECT(rind,rank);    MKVECT(rind,rank);
         MKVECT(cind,t);    MKVECT(cind,t);
         for ( i = 0; i < rank; i++ ) {    for ( i = 0; i < rank; i++ ) {
                 STOQ(ri[i],q);      STOQ(ri[i],q);
                 BDY(rind)[i] = (pointer)q;      BDY(rind)[i] = (pointer)q;
         }    }
         for ( i = 0; i < t; i++ ) {    for ( i = 0; i < t; i++ ) {
                 STOQ(ci[i],q);      STOQ(ci[i],q);
                 BDY(cind)[i] = (pointer)q;      BDY(cind)[i] = (pointer)q;
         }    }
         n0 = mknode(4,nm,dn,rind,cind);    n0 = mknode(4,nm,dn,rind,cind);
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void Pindep_rows_mod(NODE arg,VECT *rp)  void Pindep_rows_mod(NODE arg,VECT *rp)
 {  {
         MAT m,mat;    MAT m,mat;
         VECT rind;    VECT rind;
         Q **tmat;    Q **tmat;
         int **wmat,**row0;    int **wmat,**row0;
         Q *rib;    Q *rib;
         int *rowstat,*p;    int *rowstat,*p;
         Q q;    Q q;
         int md,i,j,k,l,row,col,t,rank;    int md,i,j,k,l,row,col,t,rank;
   
         asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");    asir_assert(ARG0(arg),O_MAT,"indep_rows_mod");
         asir_assert(ARG1(arg),O_N,"indep_rows_mod");    asir_assert(ARG1(arg),O_N,"indep_rows_mod");
         m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
         row = m->row; col = m->col; tmat = (Q **)m->body;    row = m->row; col = m->col; tmat = (Q **)m->body;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
   
         row0 = (int **)ALLOCA(row*sizeof(int *));    row0 = (int **)ALLOCA(row*sizeof(int *));
         for ( i = 0; i < row; i++ ) row0[i] = wmat[i];    for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
   
         rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));    rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (Q)tmat[i][j] ) {        if ( q = (Q)tmat[i][j] ) {
                                 t = rem(NM(q),md);          t = rem(NM(q),md);
                                 if ( t && SGN(q) < 0 )          if ( t && SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         } else        } else
                                 wmat[i][j] = 0;          wmat[i][j] = 0;
         rank = indep_rows_mod(wmat,row,col,md,rowstat);    rank = indep_rows_mod(wmat,row,col,md,rowstat);
   
         MKVECT(rind,rank);    MKVECT(rind,rank);
         rib = (Q *)rind->body;    rib = (Q *)rind->body;
         for ( j = 0; j < rank; j++ ) {    for ( j = 0; j < rank; j++ ) {
                 STOQ(rowstat[j],rib[j]);      STOQ(rowstat[j],rib[j]);
         }    }
     *rp = rind;      *rp = rind;
 }  }
   
 /*  /*
         input : a row x col matrix A    input : a row x col matrix A
                 A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...      A[I] <-> A[I][0]*x_0+A[I][1]*x_1+...
   
         output : [B,R,C]    output : [B,R,C]
                 B : a rank(A) x col-rank(A) matrix      B : a rank(A) x col-rank(A) matrix
                 R : a vector of length rank(A)      R : a vector of length rank(A)
                 C : a vector of length col-rank(A)      C : a vector of length col-rank(A)
                 RN : a vector of length rank(A) indicating useful rows      RN : a vector of length rank(A) indicating useful rows
   
                 B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...      B[I] <-> x_{R[I]}+B[I][0]x_{C[0]}+B[I][1]x_{C[1]}+...
 */  */
   
 void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)  void Pgeneric_gauss_elim_mod(NODE arg,LIST *rp)
 {  {
         NODE n0;    NODE n0;
         MAT m,mat;    MAT m,mat;
         VECT rind,cind,rnum;    VECT rind,cind,rnum;
         Q **tmat;    Q **tmat;
         int **wmat,**row0;    int **wmat,**row0;
         Q *rib,*cib,*rnb;    Q *rib,*cib,*rnb;
         int *colstat,*p;    int *colstat,*p;
         Q q;    Q q;
         int md,i,j,k,l,row,col,t,rank;    int md,i,j,k,l,row,col,t,rank;
   
         asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");    asir_assert(ARG0(arg),O_MAT,"generic_gauss_elim_mod");
         asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");    asir_assert(ARG1(arg),O_N,"generic_gauss_elim_mod");
         m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
         row = m->row; col = m->col; tmat = (Q **)m->body;    row = m->row; col = m->col; tmat = (Q **)m->body;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
   
         row0 = (int **)ALLOCA(row*sizeof(int *));    row0 = (int **)ALLOCA(row*sizeof(int *));
         for ( i = 0; i < row; i++ ) row0[i] = wmat[i];    for ( i = 0; i < row; i++ ) row0[i] = wmat[i];
   
         colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (Q)tmat[i][j] ) {        if ( q = (Q)tmat[i][j] ) {
                                 t = rem(NM(q),md);          t = rem(NM(q),md);
                                 if ( t && SGN(q) < 0 )          if ( t && SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         } else        } else
                                 wmat[i][j] = 0;          wmat[i][j] = 0;
         rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);    rank = generic_gauss_elim_mod(wmat,row,col,md,colstat);
   
         MKVECT(rnum,rank);    MKVECT(rnum,rank);
         rnb = (Q *)rnum->body;    rnb = (Q *)rnum->body;
         for ( i = 0; i < rank; i++ )    for ( i = 0; i < rank; i++ )
                 for ( j = 0, p = wmat[i]; j < row; j++ )      for ( j = 0, p = wmat[i]; j < row; j++ )
                         if ( p == row0[j] )        if ( p == row0[j] )
                                 STOQ(j,rnb[i]);          STOQ(j,rnb[i]);
   
         MKMAT(mat,rank,col-rank);    MKMAT(mat,rank,col-rank);
         tmat = (Q **)mat->body;    tmat = (Q **)mat->body;
         for ( i = 0; i < rank; i++ )    for ( i = 0; i < rank; i++ )
                 for ( j = k = 0; j < col; j++ )      for ( j = k = 0; j < col; j++ )
                         if ( !colstat[j] ) {        if ( !colstat[j] ) {
                                 UTOQ(wmat[i][j],tmat[i][k]); k++;          UTOQ(wmat[i][j],tmat[i][k]); k++;
                         }        }
   
         MKVECT(rind,rank);    MKVECT(rind,rank);
         MKVECT(cind,col-rank);    MKVECT(cind,col-rank);
         rib = (Q *)rind->body; cib = (Q *)cind->body;    rib = (Q *)rind->body; cib = (Q *)cind->body;
         for ( j = k = l = 0; j < col; j++ )    for ( j = k = l = 0; j < col; j++ )
                 if ( colstat[j] ) {      if ( colstat[j] ) {
                         STOQ(j,rib[k]); k++;        STOQ(j,rib[k]); k++;
                 } else {      } else {
                         STOQ(j,cib[l]); l++;        STOQ(j,cib[l]); l++;
                 }      }
         n0 = mknode(4,mat,rind,cind,rnum);    n0 = mknode(4,mat,rind,cind,rnum);
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void Pleqm(NODE arg,VECT *rp)  void Pleqm(NODE arg,VECT *rp)
 {  {
         MAT m;    MAT m;
         VECT vect;    VECT vect;
         pointer **mat;    pointer **mat;
         Q *v;    Q *v;
         Q q;    Q q;
         int **wmat;    int **wmat;
         int md,i,j,row,col,t,n,status;    int md,i,j,row,col,t,n,status;
   
         asir_assert(ARG0(arg),O_MAT,"leqm");    asir_assert(ARG0(arg),O_MAT,"leqm");
         asir_assert(ARG1(arg),O_N,"leqm");    asir_assert(ARG1(arg),O_N,"leqm");
         m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
         row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (Q)mat[i][j] ) {        if ( q = (Q)mat[i][j] ) {
                                 t = rem(NM(q),md);          t = rem(NM(q),md);
                                 if ( SGN(q) < 0 )          if ( SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         } else        } else
                                 wmat[i][j] = 0;          wmat[i][j] = 0;
         status = gauss_elim_mod(wmat,row,col,md);    status = gauss_elim_mod(wmat,row,col,md);
         if ( status < 0 )    if ( status < 0 )
                 *rp = 0;      *rp = 0;
         else if ( status > 0 )    else if ( status > 0 )
                 *rp = (VECT)ONE;      *rp = (VECT)ONE;
         else {    else {
                 n = col - 1;      n = col - 1;
                 MKVECT(vect,n);      MKVECT(vect,n);
                 for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {      for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
                         t = (md-wmat[i][n])%md; STOQ(t,v[i]);        t = (md-wmat[i][n])%md; STOQ(t,v[i]);
                 }      }
                 *rp = vect;      *rp = vect;
         }    }
 }  }
   
 int gauss_elim_mod(int **mat,int row,int col,int md)  int gauss_elim_mod(int **mat,int row,int col,int md)
 {  {
         int i,j,k,inv,a,n;    int i,j,k,inv,a,n;
         int *t,*pivot;    int *t,*pivot;
   
         n = col - 1;    n = col - 1;
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 for ( i = j; i < row && !mat[i][j]; i++ );      for ( i = j; i < row && !mat[i][j]; i++ );
                 if ( i == row )      if ( i == row )
                         return 1;        return 1;
                 if ( i != j ) {      if ( i != j ) {
                         t = mat[i]; mat[i] = mat[j]; mat[j] = t;        t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                 }      }
                 pivot = mat[j];      pivot = mat[j];
                 inv = invm(pivot[j],md);      inv = invm(pivot[j],md);
                 for ( k = j; k <= n; k++ ) {      for ( k = j; k <= n; k++ ) {
 /*                      pivot[k] = dmar(pivot[k],inv,0,md); */  /*      pivot[k] = dmar(pivot[k],inv,0,md); */
                         DMAR(pivot[k],inv,0,md,pivot[k])        DMAR(pivot[k],inv,0,md,pivot[k])
                 }      }
                 for ( i = 0; i < row; i++ ) {      for ( i = 0; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( i != j && (a = t[j]) )        if ( i != j && (a = t[j]) )
                                 for ( k = j, a = md - a; k <= n; k++ ) {          for ( k = j, a = md - a; k <= n; k++ ) {
                                         unsigned int tk;            unsigned int tk;
 /*                                      t[k] = dmar(pivot[k],a,t[k],md); */  /*          t[k] = dmar(pivot[k],a,t[k],md); */
                                         DMAR(pivot[k],a,t[k],md,tk)            DMAR(pivot[k],a,t[k],md,tk)
                                         t[k] = tk;            t[k] = tk;
                                 }          }
                 }      }
         }    }
         for ( i = n; i < row && !mat[i][n]; i++ );    for ( i = n; i < row && !mat[i][n]; i++ );
         if ( i == row )    if ( i == row )
                 return 0;      return 0;
         else    else
                 return -1;      return -1;
 }  }
   
 struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;  struct oEGT eg_mod,eg_elim,eg_elim1,eg_elim2,eg_chrem,eg_gschk,eg_intrat,eg_symb;
Line 1375  struct oEGT eg_conv;
Line 1375  struct oEGT eg_conv;
   
 int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)  int generic_gauss_elim(MAT mat,MAT *nm,Q *dn,int **rindp,int **cindp)
 {  {
         int **wmat;    int **wmat;
         Q **bmat;    Q **bmat;
         N **tmat;    N **tmat;
         Q *bmi;    Q *bmi;
         N *tmi;    N *tmi;
         Q q;    Q q;
         int *wmi;    int *wmi;
         int *colstat,*wcolstat,*rind,*cind;    int *colstat,*wcolstat,*rind,*cind;
         int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;    int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
         N m1,m2,m3,s,u;    N m1,m2,m3,s,u;
         MAT r,crmat;    MAT r,crmat;
         struct oEGT tmp0,tmp1;    struct oEGT tmp0,tmp1;
         struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;    struct oEGT eg_mod_split,eg_elim_split,eg_chrem_split;
         struct oEGT eg_intrat_split,eg_gschk_split;    struct oEGT eg_intrat_split,eg_gschk_split;
         int ret;    int ret;
   
         init_eg(&eg_mod_split); init_eg(&eg_chrem_split);    init_eg(&eg_mod_split); init_eg(&eg_chrem_split);
         init_eg(&eg_elim_split); init_eg(&eg_intrat_split);    init_eg(&eg_elim_split); init_eg(&eg_intrat_split);
         init_eg(&eg_gschk_split);    init_eg(&eg_gschk_split);
         bmat = (Q **)mat->body;    bmat = (Q **)mat->body;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
         colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));    wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
         for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
                 if ( DP_Print ) {      if ( DP_Print ) {
                         fprintf(asir_out,"."); fflush(asir_out);        fprintf(asir_out,"."); fflush(asir_out);
                 }      }
                 md = get_lprime(ind);      md = get_lprime(ind);
                 get_eg(&tmp0);      get_eg(&tmp0);
                 for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
                         for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )        for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
                                 if ( q = (Q)bmi[j] ) {          if ( q = (Q)bmi[j] ) {
                                         t = rem(NM(q),md);            t = rem(NM(q),md);
                                         if ( t && SGN(q) < 0 )            if ( t && SGN(q) < 0 )
                                                 t = (md - t) % md;              t = (md - t) % md;
                                         wmi[j] = t;            wmi[j] = t;
                                 } else          } else
                                         wmi[j] = 0;            wmi[j] = 0;
                 get_eg(&tmp1);      get_eg(&tmp1);
                 add_eg(&eg_mod,&tmp0,&tmp1);      add_eg(&eg_mod,&tmp0,&tmp1);
                 add_eg(&eg_mod_split,&tmp0,&tmp1);      add_eg(&eg_mod_split,&tmp0,&tmp1);
                 get_eg(&tmp0);      get_eg(&tmp0);
                 rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);      rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
                 get_eg(&tmp1);      get_eg(&tmp1);
                 add_eg(&eg_elim,&tmp0,&tmp1);      add_eg(&eg_elim,&tmp0,&tmp1);
                 add_eg(&eg_elim_split,&tmp0,&tmp1);      add_eg(&eg_elim_split,&tmp0,&tmp1);
                 if ( !ind ) {      if ( !ind ) {
 RESET:  RESET:
                         UTON(md,m1);        UTON(md,m1);
                         rank0 = rank;        rank0 = rank;
                         bcopy(wcolstat,colstat,col*sizeof(int));        bcopy(wcolstat,colstat,col*sizeof(int));
                         MKMAT(crmat,rank,col-rank);        MKMAT(crmat,rank,col-rank);
                         MKMAT(r,rank,col-rank); *nm = r;        MKMAT(r,rank,col-rank); *nm = r;
                         tmat = (N **)crmat->body;        tmat = (N **)crmat->body;
                         for ( i = 0; i < rank; i++ )        for ( i = 0; i < rank; i++ )
                                 for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )          for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                                         if ( !colstat[j] ) {            if ( !colstat[j] ) {
                                                 UTON(wmi[j],tmi[k]); k++;              UTON(wmi[j],tmi[k]); k++;
                                         }            }
                 } else {      } else {
                         if ( rank < rank0 ) {        if ( rank < rank0 ) {
                                 if ( DP_Print ) {          if ( DP_Print ) {
                                         fprintf(asir_out,"lower rank matrix; continuing...\n");            fprintf(asir_out,"lower rank matrix; continuing...\n");
                                         fflush(asir_out);            fflush(asir_out);
                                 }          }
                                 continue;          continue;
                         } else if ( rank > rank0 ) {        } else if ( rank > rank0 ) {
                                 if ( DP_Print ) {          if ( DP_Print ) {
                                         fprintf(asir_out,"higher rank matrix; resetting...\n");            fprintf(asir_out,"higher rank matrix; resetting...\n");
                                         fflush(asir_out);            fflush(asir_out);
                                 }          }
                                 goto RESET;          goto RESET;
                         } else {        } else {
                                 for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );          for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
                                 if ( j < col ) {          if ( j < col ) {
                                         if ( DP_Print ) {            if ( DP_Print ) {
                                                 fprintf(asir_out,"inconsitent colstat; resetting...\n");              fprintf(asir_out,"inconsitent colstat; resetting...\n");
                                                 fflush(asir_out);              fflush(asir_out);
                                         }            }
                                         goto RESET;            goto RESET;
                                 }          }
                         }        }
   
                         get_eg(&tmp0);        get_eg(&tmp0);
                         inv = invm(rem(m1,md),md);        inv = invm(rem(m1,md),md);
                         UTON(md,m2); muln(m1,m2,&m3);        UTON(md,m2); muln(m1,m2,&m3);
                         for ( i = 0; i < rank; i++ )        for ( i = 0; i < rank; i++ )
                                 for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )          for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
                                         if ( !colstat[j] ) {            if ( !colstat[j] ) {
                                                 if ( tmi[k] ) {              if ( tmi[k] ) {
                                                 /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */              /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                                                         t = rem(tmi[k],md);                t = rem(tmi[k],md);
                                                         if ( wmi[j] >= t )                if ( wmi[j] >= t )
                                                                 t = wmi[j]-t;                  t = wmi[j]-t;
                                                         else                else
                                                                 t = md-(t-wmi[j]);                  t = md-(t-wmi[j]);
                                                         DMAR(t,inv,0,md,t1)                DMAR(t,inv,0,md,t1)
                                                         UTON(t1,u);                UTON(t1,u);
                                                         muln(m1,u,&s);                muln(m1,u,&s);
                                                         addn(tmi[k],s,&u); tmi[k] = u;                addn(tmi[k],s,&u); tmi[k] = u;
                                                 } else if ( wmi[j] ) {              } else if ( wmi[j] ) {
                                                 /* f3 = m1*(m1 mod m2)^(-1)*f2 */              /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                                                         DMAR(wmi[j],inv,0,md,t)                DMAR(wmi[j],inv,0,md,t)
                                                         UTON(t,u);                UTON(t,u);
                                                         muln(m1,u,&s); tmi[k] = s;                muln(m1,u,&s); tmi[k] = s;
                                                 }              }
                                                 k++;              k++;
                                         }            }
                         m1 = m3;        m1 = m3;
                         get_eg(&tmp1);        get_eg(&tmp1);
                         add_eg(&eg_chrem,&tmp0,&tmp1);        add_eg(&eg_chrem,&tmp0,&tmp1);
                         add_eg(&eg_chrem_split,&tmp0,&tmp1);        add_eg(&eg_chrem_split,&tmp0,&tmp1);
   
                         get_eg(&tmp0);        get_eg(&tmp0);
                         if ( ind % F4_INTRAT_PERIOD )        if ( ind % F4_INTRAT_PERIOD )
                                 ret = 0;          ret = 0;
                         else        else
                                 ret = intmtoratm(crmat,m1,*nm,dn);          ret = intmtoratm(crmat,m1,*nm,dn);
                         get_eg(&tmp1);        get_eg(&tmp1);
                         add_eg(&eg_intrat,&tmp0,&tmp1);        add_eg(&eg_intrat,&tmp0,&tmp1);
                         add_eg(&eg_intrat_split,&tmp0,&tmp1);        add_eg(&eg_intrat_split,&tmp0,&tmp1);
                         if ( ret ) {        if ( ret ) {
                                 *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));          *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                                 *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));          *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                                 for ( j = k = l = 0; j < col; j++ )          for ( j = k = l = 0; j < col; j++ )
                                         if ( colstat[j] )            if ( colstat[j] )
                                                 rind[k++] = j;              rind[k++] = j;
                                         else            else
                                                 cind[l++] = j;              cind[l++] = j;
                                 get_eg(&tmp0);          get_eg(&tmp0);
                                 if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {          if ( gensolve_check(mat,*nm,*dn,rind,cind) ) {
                                         get_eg(&tmp1);            get_eg(&tmp1);
                                         add_eg(&eg_gschk,&tmp0,&tmp1);            add_eg(&eg_gschk,&tmp0,&tmp1);
                                         add_eg(&eg_gschk_split,&tmp0,&tmp1);            add_eg(&eg_gschk_split,&tmp0,&tmp1);
                                         if ( DP_Print ) {            if ( DP_Print ) {
                                                 print_eg("Mod",&eg_mod_split);              print_eg("Mod",&eg_mod_split);
                                                 print_eg("Elim",&eg_elim_split);              print_eg("Elim",&eg_elim_split);
                                                 print_eg("ChRem",&eg_chrem_split);              print_eg("ChRem",&eg_chrem_split);
                                                 print_eg("IntRat",&eg_intrat_split);              print_eg("IntRat",&eg_intrat_split);
                                                 print_eg("Check",&eg_gschk_split);              print_eg("Check",&eg_gschk_split);
                                                 fflush(asir_out);              fflush(asir_out);
                                         }            }
                                         return rank;            return rank;
                                 }          }
                         }        }
                 }      }
         }    }
 }  }
   
 void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
Line 1530  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
Line 1530  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm);
 /* XXX broken */  /* XXX broken */
 void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)  void lu_dec_cr(MAT mat,MAT lu,Q *dn,int **perm)
 {  {
         Q **a0,**b;    Q **a0,**b;
         Q *aiq;    Q *aiq;
         N **a;    N **a;
         N *ai;    N *ai;
         Q q,q1,dn2,a1,q0,bik;    Q q,q1,dn2,a1,q0,bik;
         MAT m;    MAT m;
         unsigned int md;    unsigned int md;
         int n,ind,i,j,rank,t,inv,t1,ret,min,k;    int n,ind,i,j,rank,t,inv,t1,ret,min,k;
         int **w;    int **w;
         int *wi,*rinfo0,*rinfo;    int *wi,*rinfo0,*rinfo;
         N m1,m2,m3,u,s;    N m1,m2,m3,u,s;
   
         a0 = (Q **)mat->body;    a0 = (Q **)mat->body;
         n = mat->row;    n = mat->row;
         if ( n != mat->col )    if ( n != mat->col )
                 error("lu_dec_cr : non-square matrix");      error("lu_dec_cr : non-square matrix");
         w = (int **)almat(n,n);    w = (int **)almat(n,n);
         MKMAT(m,n,n);    MKMAT(m,n,n);
         a = (N **)m->body;    a = (N **)m->body;
         UTON(1,m1);    UTON(1,m1);
         rinfo0 = 0;    rinfo0 = 0;
         ind = 0;    ind = 0;
         while ( 1 ) {    while ( 1 ) {
                 md = get_lprime(ind);      md = get_lprime(ind);
                 /* mat mod md */      /* mat mod md */
                 for ( i = 0; i < n; i++ )      for ( i = 0; i < n; i++ )
                         for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )        for ( j = 0, aiq = a0[i], wi = w[i]; j < n; j++ )
                                 if ( q = aiq[j] ) {          if ( q = aiq[j] ) {
                                         t = rem(NM(q),md);            t = rem(NM(q),md);
                                         if ( t && SGN(q) < 0 )            if ( t && SGN(q) < 0 )
                                                 t = (md - t) % md;              t = (md - t) % md;
                                         wi[j] = t;            wi[j] = t;
                                 } else          } else
                                         wi[j] = 0;            wi[j] = 0;
   
                 if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;      if ( !lu_mod((unsigned int **)w,n,md,&rinfo) ) continue;
                 printf("."); fflush(stdout);      printf("."); fflush(stdout);
                 if ( !rinfo0 )      if ( !rinfo0 )
                         *perm = rinfo0 = rinfo;        *perm = rinfo0 = rinfo;
                 else {      else {
                         for ( i = 0; i < n; i++ )        for ( i = 0; i < n; i++ )
                                 if ( rinfo[i] != rinfo0[i] ) break;          if ( rinfo[i] != rinfo0[i] ) break;
                         if ( i < n ) continue;        if ( i < n ) continue;
                 }      }
                 if ( UNIN(m1) ) {      if ( UNIN(m1) ) {
                         for ( i = 0; i < n; i++ )        for ( i = 0; i < n; i++ )
                                 for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {          for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ ) {
                                         UTON(wi[j],u); ai[j] = u;            UTON(wi[j],u); ai[j] = u;
                                 }          }
                         UTON(md,m1);        UTON(md,m1);
                 } else {      } else {
                         inv = invm(rem(m1,md),md);        inv = invm(rem(m1,md),md);
                         UTON(md,m2); muln(m1,m2,&m3);        UTON(md,m2); muln(m1,m2,&m3);
                         for ( i = 0; i < n; i++ )        for ( i = 0; i < n; i++ )
                                 for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )          for ( j = 0, ai = a[i], wi = w[i]; j < n; j++ )
                                         if ( ai[i] ) {            if ( ai[i] ) {
                                         /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */            /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
                                                 t = rem(ai[j],md);              t = rem(ai[j],md);
                                                 if ( wi[j] >= t )              if ( wi[j] >= t )
                                                         t = wi[j]-t;                t = wi[j]-t;
                                                 else              else
                                                         t = md-(t-wi[j]);                t = md-(t-wi[j]);
                                                 DMAR(t,inv,0,md,t1)              DMAR(t,inv,0,md,t1)
                                                 UTON(t1,u);              UTON(t1,u);
                                                 muln(m1,u,&s);              muln(m1,u,&s);
                                                 addn(ai[j],s,&u); ai[j] = u;              addn(ai[j],s,&u); ai[j] = u;
                                         } else if ( wi[j] ) {            } else if ( wi[j] ) {
                                                 /* f3 = m1*(m1 mod m2)^(-1)*f2 */              /* f3 = m1*(m1 mod m2)^(-1)*f2 */
                                                 DMAR(wi[j],inv,0,md,t)              DMAR(wi[j],inv,0,md,t)
                                                 UTON(t,u);              UTON(t,u);
                                                 muln(m1,u,&s); ai[j] = s;              muln(m1,u,&s); ai[j] = s;
                                         }            }
                         m1 = m3;        m1 = m3;
                 }      }
                 if ( (++ind%8) == 0 ) {      if ( (++ind%8) == 0 ) {
                         ret = intmtoratm(m,m1,lu,dn);        ret = intmtoratm(m,m1,lu,dn);
                         if ( ret ) {        if ( ret ) {
                                 b = (Q **)lu->body;          b = (Q **)lu->body;
                                 mulq(*dn,*dn,&dn2);          mulq(*dn,*dn,&dn2);
                                 for ( i = 0; i < n; i++ ) {          for ( i = 0; i < n; i++ ) {
                                         for ( j = 0; j < n; j++ ) {            for ( j = 0; j < n; j++ ) {
                                                 q = 0;              q = 0;
                                                 min = MIN(i,j);              min = MIN(i,j);
                                                 for ( k = 0; k <= min; k++ ) {              for ( k = 0; k <= min; k++ ) {
                                                         bik = k==i ? *dn : b[i][k];                bik = k==i ? *dn : b[i][k];
                                                         mulq(bik,b[k][j],&q0);                mulq(bik,b[k][j],&q0);
                                                         addq(q,q0,&q1); q = q1;                addq(q,q0,&q1); q = q1;
                                                 }              }
                                                 mulq(a0[rinfo0[i]][j],dn2,&q1);              mulq(a0[rinfo0[i]][j],dn2,&q1);
                                                 if ( cmpq(q,q1) ) break;              if ( cmpq(q,q1) ) break;
                                         }            }
                                         if ( j < n ) break;            if ( j < n ) break;
                                 }          }
                                 if ( i == n )          if ( i == n )
                                         return;            return;
                         }        }
                 }      }
         }    }
 }  }
   
 void nmat(N **m,int n)  void nmat(N **m,int n)
 {  {
         int i,j;    int i,j;
   
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 for ( j = 0; j < n; j++ ) {      for ( j = 0; j < n; j++ ) {
                         printn(m[i][j]); printf(" ");        printn(m[i][j]); printf(" ");
                 }      }
                 printf("\n");      printf("\n");
         }    }
 }  }
   
 int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)  int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Q *dn,int **rindp,int **cindp)
 {  {
         MAT bmat,xmat;    MAT bmat,xmat;
         Q **a0,**a,**b,**x,**nm;    Q **a0,**a,**b,**x,**nm;
         Q *ai,*bi,*xi;    Q *ai,*bi,*xi;
         int row,col;    int row,col;
         int **w;    int **w;
         int *wi;    int *wi;
         int **wc;    int **wc;
         Q mdq,q,s,u;    Q mdq,q,s,u;
         N tn;    N tn;
         int ind,md,i,j,k,l,li,ri,rank;    int ind,md,i,j,k,l,li,ri,rank;
         unsigned int t;    unsigned int t;
         int *cinfo,*rinfo;    int *cinfo,*rinfo;
         int *rind,*cind;    int *rind,*cind;
         int count;    int count;
         int ret;    int ret;
         struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;    struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
         int period;    int period;
         int *wx,*ptr;    int *wx,*ptr;
         int wxsize,nsize;    int wxsize,nsize;
         N wn;    N wn;
         Q wq;    Q wq;
   
         a0 = (Q **)mat->body;    a0 = (Q **)mat->body;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         w = (int **)almat(row,col);    w = (int **)almat(row,col);
         for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
                 md = get_lprime(ind);      md = get_lprime(ind);
                 STOQ(md,mdq);      STOQ(md,mdq);
                 for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
                         for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )        for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
                                 if ( q = (Q)ai[j] ) {          if ( q = (Q)ai[j] ) {
                                         t = rem(NM(q),md);            t = rem(NM(q),md);
                                         if ( t && SGN(q) < 0 )            if ( t && SGN(q) < 0 )
                                                 t = (md - t) % md;              t = (md - t) % md;
                                         wi[j] = t;            wi[j] = t;
                                 } else          } else
                                         wi[j] = 0;            wi[j] = 0;
   
                 if ( DP_Print > 3 ) {      if ( DP_Print > 3 ) {
                         fprintf(asir_out,"LU decomposition.."); fflush(asir_out);        fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
                 }      }
                 rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);      rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
                 if ( DP_Print > 3 ) {      if ( DP_Print > 3 ) {
                         fprintf(asir_out,"done.\n"); fflush(asir_out);        fprintf(asir_out,"done.\n"); fflush(asir_out);
                 }      }
                 a = (Q **)almat_pointer(rank,rank); /* lhs mat */      a = (Q **)almat_pointer(rank,rank); /* lhs mat */
                 MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */      MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
                 for ( j = li = ri = 0; j < col; j++ )      for ( j = li = ri = 0; j < col; j++ )
                         if ( cinfo[j] ) {        if ( cinfo[j] ) {
                                 /* the column is in lhs */          /* the column is in lhs */
                                 for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
                                         w[i][li] = w[i][j];            w[i][li] = w[i][j];
                                         a[i][li] = a0[rinfo[i]][j];            a[i][li] = a0[rinfo[i]][j];
                                 }          }
                                 li++;          li++;
                         } else {        } else {
                                 /* the column is in rhs */          /* the column is in rhs */
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         b[i][ri] = a0[rinfo[i]][j];            b[i][ri] = a0[rinfo[i]][j];
                                 ri++;          ri++;
                         }        }
   
                         /* solve Ax+B=0; A: rank x rank, B: rank x ri */        /* solve Ax+B=0; A: rank x rank, B: rank x ri */
                         MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;        MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
                         MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;        MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
                         /* use the right part of w as work area */        /* use the right part of w as work area */
                         /* ri = col - rank */        /* ri = col - rank */
                         wc = (int **)almat(rank,ri);        wc = (int **)almat(rank,ri);
                         for ( i = 0; i < rank; i++ )        for ( i = 0; i < rank; i++ )
                                 wc[i] = w[i]+rank;          wc[i] = w[i]+rank;
                         *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));        *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                         *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));        *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
   
                         init_eg(&eg_mul); init_eg(&eg_inv);        init_eg(&eg_mul); init_eg(&eg_inv);
                         init_eg(&eg_check); init_eg(&eg_intrat);        init_eg(&eg_check); init_eg(&eg_intrat);
                         period = F4_INTRAT_PERIOD;        period = F4_INTRAT_PERIOD;
                         nsize = period;        nsize = period;
                         wxsize = rank*ri*nsize;        wxsize = rank*ri*nsize;
                         wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));        wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
                         for ( i = 0; i < wxsize; i++ ) wx[i] = 0;        for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                         for ( q = ONE, count = 0; ; ) {        for ( q = ONE, count = 0; ; ) {
                                 if ( DP_Print > 3 )          if ( DP_Print > 3 )
                                         fprintf(stderr,"o");            fprintf(stderr,"o");
                                 /* wc = -b mod md */          /* wc = -b mod md */
                                 get_eg(&tmp0);          get_eg(&tmp0);
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )            for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
                                                 if ( u = (Q)bi[j] ) {              if ( u = (Q)bi[j] ) {
                                                         t = rem(NM(u),md);                t = rem(NM(u),md);
                                                         if ( t && SGN(u) > 0 )                if ( t && SGN(u) > 0 )
                                                                 t = (md - t) % md;                  t = (md - t) % md;
                                                         wi[j] = t;                wi[j] = t;
                                                 } else              } else
                                                         wi[j] = 0;                wi[j] = 0;
                                 /* wc = A^(-1)wc; wc is not normalized */          /* wc = A^(-1)wc; wc is not normalized */
                                 solve_by_lu_mod(w,rank,md,wc,ri,0);          solve_by_lu_mod(w,rank,md,wc,ri,0);
                                 /* wx += q*wc */          /* wx += q*wc */
                                 ptr = wx;          ptr = wx;
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         for ( j = 0, wi = wc[i]; j < ri; j++ ) {            for ( j = 0, wi = wc[i]; j < ri; j++ ) {
                                                 if ( wi[j] )              if ( wi[j] )
                                                         muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);                muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
                                                 ptr += nsize;              ptr += nsize;
                                         }            }
                                 count++;          count++;
                                 get_eg(&tmp1);          get_eg(&tmp1);
                                 add_eg(&eg_inv,&tmp0,&tmp1);          add_eg(&eg_inv,&tmp0,&tmp1);
                                 get_eg(&tmp0);          get_eg(&tmp0);
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         for ( j = 0; j < ri; j++ ) {            for ( j = 0; j < ri; j++ ) {
                                                 inner_product_mat_int_mod(a,wc,rank,i,j,&u);              inner_product_mat_int_mod(a,wc,rank,i,j,&u);
                                                 addq(b[i][j],u,&s);              addq(b[i][j],u,&s);
                                                 if ( s ) {              if ( s ) {
                                                         t = divin(NM(s),md,&tn);                t = divin(NM(s),md,&tn);
                                                         if ( t )                if ( t )
                                                                 error("generic_gauss_elim_hensel:incosistent");                  error("generic_gauss_elim_hensel:incosistent");
                                                         NTOQ(tn,SGN(s),b[i][j]);                NTOQ(tn,SGN(s),b[i][j]);
                                                 } else              } else
                                                         b[i][j] = 0;                b[i][j] = 0;
                                         }            }
                                 get_eg(&tmp1);          get_eg(&tmp1);
                                 add_eg(&eg_mul,&tmp0,&tmp1);          add_eg(&eg_mul,&tmp0,&tmp1);
                                 /* q = q*md */          /* q = q*md */
                                 mulq(q,mdq,&u); q = u;          mulq(q,mdq,&u); q = u;
                                 if ( count == period ) {          if ( count == period ) {
                                         get_eg(&tmp0);            get_eg(&tmp0);
                                         ptr = wx;            ptr = wx;
                                         for ( i = 0; i < rank; i++ )            for ( i = 0; i < rank; i++ )
                                                 for ( j = 0, xi = x[i]; j < ri;              for ( j = 0, xi = x[i]; j < ri;
                                                         j++, ptr += nsize ) {                j++, ptr += nsize ) {
                                                         for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );                for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
                                                         if ( k >= 0 ) {                if ( k >= 0 ) {
                                                                 wn = NALLOC(k+1);                  wn = NALLOC(k+1);
                                                                 PL(wn) = k+1;                  PL(wn) = k+1;
                                                                 for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];                  for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
                                                                 NTOQ(wn,1,wq);                  NTOQ(wn,1,wq);
                                                                 subq(xi[j],wq,&u); xi[j] = u;                  subq(xi[j],wq,&u); xi[j] = u;
                                                         }                }
                                                 }              }
                                         ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);            ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
                                         get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);            get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
                                         if ( ret ) {            if ( ret ) {
                                                 rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));              rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                                                 cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));              cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                                                 for ( j = k = l = 0; j < col; j++ )              for ( j = k = l = 0; j < col; j++ )
                                                         if ( cinfo[j] )                if ( cinfo[j] )
                                                                 rind[k++] = j;                  rind[k++] = j;
                                                         else                else
                                                                 cind[l++] = j;                  cind[l++] = j;
                                                 get_eg(&tmp0);              get_eg(&tmp0);
                                                 ret = gensolve_check(mat,*nmmat,*dn,rind,cind);              ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
                                                 get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);              get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
                                                 if ( ret ) {              if ( ret ) {
                                                         if ( DP_Print > 3 ) {                if ( DP_Print > 3 ) {
                                                                 fprintf(stderr,"\n");                  fprintf(stderr,"\n");
                                                                 print_eg("INV",&eg_inv);                  print_eg("INV",&eg_inv);
                                                                 print_eg("MUL",&eg_mul);                  print_eg("MUL",&eg_mul);
                                                                 print_eg("INTRAT",&eg_intrat);                  print_eg("INTRAT",&eg_intrat);
                                                                 print_eg("CHECK",&eg_check);                  print_eg("CHECK",&eg_check);
                                                                 fflush(asir_out);                  fflush(asir_out);
                                                         }                }
                                                         *rindp = rind;                *rindp = rind;
                                                         *cindp = cind;                *cindp = cind;
                                                         for ( j = k = 0; j < col; j++ )                for ( j = k = 0; j < col; j++ )
                                                                 if ( !cinfo[j] )                  if ( !cinfo[j] )
                                                                         cind[k++] = j;                    cind[k++] = j;
                                                         return rank;                return rank;
                                                 }              }
                                         } else {            } else {
                                                 period = period*3/2;              period = period*3/2;
                                                 count = 0;              count = 0;
                                                 nsize += period;              nsize += period;
                                                 wxsize += rank*ri*nsize;              wxsize += rank*ri*nsize;
                                                 wx = (int *)REALLOC(wx,wxsize*sizeof(int));              wx = (int *)REALLOC(wx,wxsize*sizeof(int));
                                                 for ( i = 0; i < wxsize; i++ ) wx[i] = 0;              for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                                         }            }
                                 }          }
                         }        }
         }    }
 }  }
   
 int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp)  int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Q *dn,int **rindp,int **cindp)
 {  {
         MAT bmat,xmat;    MAT bmat,xmat;
         Q **a0,**a,**b,**x,**nm;    Q **a0,**a,**b,**x,**nm;
         Q *ai,*bi,*xi;    Q *ai,*bi,*xi;
         int row,col;    int row,col;
         int **w;    int **w;
         int *wi;    int *wi;
         int **wc;    int **wc;
         Q mdq,q,s,u;    Q mdq,q,s,u;
         N tn;    N tn;
         int ind,md,i,j,k,l,li,ri,rank;    int ind,md,i,j,k,l,li,ri,rank;
         unsigned int t;    unsigned int t;
         int *cinfo,*rinfo;    int *cinfo,*rinfo;
         int *rind,*cind;    int *rind,*cind;
         int count;    int count;
         int ret;    int ret;
         struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;    struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
         int period;    int period;
         int *wx,*ptr;    int *wx,*ptr;
         int wxsize,nsize;    int wxsize,nsize;
         N wn;    N wn;
         Q wq;    Q wq;
         NumberField nf;    NumberField nf;
         DP m;    DP m;
         int col1;    int col1;
   
         a0 = (Q **)mat->body;    a0 = (Q **)mat->body;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         w = (int **)almat(row,col);    w = (int **)almat(row,col);
         for ( ind = 0; ; ind++ ) {    for ( ind = 0; ; ind++ ) {
                 md = get_lprime(ind);      md = get_lprime(ind);
                 STOQ(md,mdq);      STOQ(md,mdq);
                 for ( i = 0; i < row; i++ )      for ( i = 0; i < row; i++ )
                         for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )        for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
                                 if ( q = (Q)ai[j] ) {          if ( q = (Q)ai[j] ) {
                                         t = rem(NM(q),md);            t = rem(NM(q),md);
                                         if ( t && SGN(q) < 0 )            if ( t && SGN(q) < 0 )
                                                 t = (md - t) % md;              t = (md - t) % md;
                                         wi[j] = t;            wi[j] = t;
                                 } else          } else
                                         wi[j] = 0;            wi[j] = 0;
   
                 if ( DP_Print ) {      if ( DP_Print ) {
                         fprintf(asir_out,"LU decomposition.."); fflush(asir_out);        fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
                 }      }
                 rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);      rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
                 if ( DP_Print ) {      if ( DP_Print ) {
                         fprintf(asir_out,"done.\n"); fflush(asir_out);        fprintf(asir_out,"done.\n"); fflush(asir_out);
                 }      }
                 for ( i = 0; i < col-1; i++ ) {      for ( i = 0; i < col-1; i++ ) {
                         if ( !cinfo[i] ) {        if ( !cinfo[i] ) {
                                 m = mb[i];          m = mb[i];
                                 for ( j = i+1; j < col-1; j++ )          for ( j = i+1; j < col-1; j++ )
                                         if ( dp_redble(mb[j],m) )            if ( dp_redble(mb[j],m) )
                                                 cinfo[j] = -1;              cinfo[j] = -1;
                         }        }
                 }      }
                 a = (Q **)almat_pointer(rank,rank); /* lhs mat */      a = (Q **)almat_pointer(rank,rank); /* lhs mat */
                 MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */      MKMAT(bmat,rank,col-rank); b = (Q **)bmat->body; /* lhs mat */
                 for ( j = li = ri = 0; j < col; j++ )      for ( j = li = ri = 0; j < col; j++ )
                         if ( cinfo[j] > 0 ) {        if ( cinfo[j] > 0 ) {
                                 /* the column is in lhs */          /* the column is in lhs */
                                 for ( i = 0; i < rank; i++ ) {          for ( i = 0; i < rank; i++ ) {
                                         w[i][li] = w[i][j];            w[i][li] = w[i][j];
                                         a[i][li] = a0[rinfo[i]][j];            a[i][li] = a0[rinfo[i]][j];
                                 }          }
                                 li++;          li++;
                         } else if ( !cinfo[j] ) {        } else if ( !cinfo[j] ) {
                                 /* the column is in rhs */          /* the column is in rhs */
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         b[i][ri] = a0[rinfo[i]][j];            b[i][ri] = a0[rinfo[i]][j];
                                 ri++;          ri++;
                         }        }
   
                         /* solve Ax+B=0; A: rank x rank, B: rank x ri */        /* solve Ax+B=0; A: rank x rank, B: rank x ri */
                         MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;        MKMAT(xmat,rank,ri); x = (Q **)(xmat)->body;
                         MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;        MKMAT(*nmmat,rank,ri); nm = (Q **)(*nmmat)->body;
                         /* use the right part of w as work area */        /* use the right part of w as work area */
                         wc = (int **)almat(rank,ri);        wc = (int **)almat(rank,ri);
                         for ( i = 0; i < rank; i++ )        for ( i = 0; i < rank; i++ )
                                 wc[i] = w[i]+rank;          wc[i] = w[i]+rank;
                         *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));        *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
                         *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));        *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
                         init_eg(&eg_mul); init_eg(&eg_inv);        init_eg(&eg_mul); init_eg(&eg_inv);
                         init_eg(&eg_check); init_eg(&eg_intrat);        init_eg(&eg_check); init_eg(&eg_intrat);
                         period = F4_INTRAT_PERIOD;        period = F4_INTRAT_PERIOD;
                         nsize = period;        nsize = period;
                         wxsize = rank*ri*nsize;        wxsize = rank*ri*nsize;
                         wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));        wx = (int *)MALLOC_ATOMIC(wxsize*sizeof(int));
                         for ( i = 0; i < wxsize; i++ ) wx[i] = 0;        for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                         for ( q = ONE, count = 0; ; ) {        for ( q = ONE, count = 0; ; ) {
                                 if ( DP_Print )          if ( DP_Print )
                                         fprintf(stderr,"o");            fprintf(stderr,"o");
                                 /* wc = -b mod md */          /* wc = -b mod md */
                                 get_eg(&tmp0);          get_eg(&tmp0);
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )            for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
                                                 if ( u = (Q)bi[j] ) {              if ( u = (Q)bi[j] ) {
                                                         t = rem(NM(u),md);                t = rem(NM(u),md);
                                                         if ( t && SGN(u) > 0 )                if ( t && SGN(u) > 0 )
                                                                 t = (md - t) % md;                  t = (md - t) % md;
                                                         wi[j] = t;                wi[j] = t;
                                                 } else              } else
                                                         wi[j] = 0;                wi[j] = 0;
                                 /* wc = A^(-1)wc; wc is not normalized */          /* wc = A^(-1)wc; wc is not normalized */
                                 solve_by_lu_mod(w,rank,md,wc,ri,0);          solve_by_lu_mod(w,rank,md,wc,ri,0);
                                 /* wx += q*wc */          /* wx += q*wc */
                                 ptr = wx;          ptr = wx;
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         for ( j = 0, wi = wc[i]; j < ri; j++ ) {            for ( j = 0, wi = wc[i]; j < ri; j++ ) {
                                                 if ( wi[j] )              if ( wi[j] )
                                                         muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);                muln_1(BD(NM(q)),PL(NM(q)),wi[j],ptr);
                                                 ptr += nsize;              ptr += nsize;
                                         }            }
                                 count++;          count++;
                                 get_eg(&tmp1);          get_eg(&tmp1);
                                 add_eg(&eg_inv,&tmp0,&tmp1);          add_eg(&eg_inv,&tmp0,&tmp1);
                                 get_eg(&tmp0);          get_eg(&tmp0);
                                 for ( i = 0; i < rank; i++ )          for ( i = 0; i < rank; i++ )
                                         for ( j = 0; j < ri; j++ ) {            for ( j = 0; j < ri; j++ ) {
                                                 inner_product_mat_int_mod(a,wc,rank,i,j,&u);              inner_product_mat_int_mod(a,wc,rank,i,j,&u);
                                                 addq(b[i][j],u,&s);              addq(b[i][j],u,&s);
                                                 if ( s ) {              if ( s ) {
                                                         t = divin(NM(s),md,&tn);                t = divin(NM(s),md,&tn);
                                                         if ( t )                if ( t )
                                                                 error("generic_gauss_elim_hensel:incosistent");                  error("generic_gauss_elim_hensel:incosistent");
                                                         NTOQ(tn,SGN(s),b[i][j]);                NTOQ(tn,SGN(s),b[i][j]);
                                                 } else              } else
                                                         b[i][j] = 0;                b[i][j] = 0;
                                         }            }
                                 get_eg(&tmp1);          get_eg(&tmp1);
                                 add_eg(&eg_mul,&tmp0,&tmp1);          add_eg(&eg_mul,&tmp0,&tmp1);
                                 /* q = q*md */          /* q = q*md */
                                 mulq(q,mdq,&u); q = u;          mulq(q,mdq,&u); q = u;
                                 if ( count == period ) {          if ( count == period ) {
                                         get_eg(&tmp0);            get_eg(&tmp0);
                                         ptr = wx;            ptr = wx;
                                         for ( i = 0; i < rank; i++ )            for ( i = 0; i < rank; i++ )
                                                 for ( j = 0, xi = x[i]; j < ri;              for ( j = 0, xi = x[i]; j < ri;
                                                         j++, ptr += nsize ) {                j++, ptr += nsize ) {
                                                         for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );                for ( k = nsize-1; k >= 0 && !ptr[k]; k-- );
                                                         if ( k >= 0 ) {                if ( k >= 0 ) {
                                                                 wn = NALLOC(k+1);                  wn = NALLOC(k+1);
                                                                 PL(wn) = k+1;                  PL(wn) = k+1;
                                                                 for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];                  for ( l = 0; l <= k; l++ ) BD(wn)[l] = (unsigned int)ptr[l];
                                                                 NTOQ(wn,1,wq);                  NTOQ(wn,1,wq);
                                                                 subq(xi[j],wq,&u); xi[j] = u;                  subq(xi[j],wq,&u); xi[j] = u;
                                                         }                }
                                                 }              }
                                         ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);            ret = intmtoratm_q(xmat,NM(q),*nmmat,dn);
                                         get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);            get_eg(&tmp1); add_eg(&eg_intrat,&tmp0,&tmp1);
                                         if ( ret ) {            if ( ret ) {
                                                 for ( j = k = l = 0; j < col; j++ )              for ( j = k = l = 0; j < col; j++ )
                                                         if ( cinfo[j] > 0 )                if ( cinfo[j] > 0 )
                                                                 rind[k++] = j;                  rind[k++] = j;
                                                         else if ( !cinfo[j] )                else if ( !cinfo[j] )
                                                                 cind[l++] = j;                  cind[l++] = j;
                                                 get_eg(&tmp0);              get_eg(&tmp0);
                                                 ret = gensolve_check(mat,*nmmat,*dn,rind,cind);              ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
                                                 get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);              get_eg(&tmp1); add_eg(&eg_check,&tmp0,&tmp1);
                                                 if ( ret ) {              if ( ret ) {
                                                         if ( DP_Print > 3 ) {                if ( DP_Print > 3 ) {
                                                                 fprintf(stderr,"\n");                  fprintf(stderr,"\n");
                                                                 print_eg("INV",&eg_inv);                  print_eg("INV",&eg_inv);
                                                                 print_eg("MUL",&eg_mul);                  print_eg("MUL",&eg_mul);
                                                                 print_eg("INTRAT",&eg_intrat);                  print_eg("INTRAT",&eg_intrat);
                                                                 print_eg("CHECK",&eg_check);                  print_eg("CHECK",&eg_check);
                                                                 fflush(asir_out);                  fflush(asir_out);
                                                         }                }
                                                         return rank;                return rank;
                                                 }              }
                                         } else {            } else {
                                                 period = period*3/2;              period = period*3/2;
                                                 count = 0;              count = 0;
                                                 nsize += period;              nsize += period;
                                                 wxsize += rank*ri*nsize;              wxsize += rank*ri*nsize;
                                                 wx = (int *)REALLOC(wx,wxsize*sizeof(int));              wx = (int *)REALLOC(wx,wxsize*sizeof(int));
                                                 for ( i = 0; i < wxsize; i++ ) wx[i] = 0;              for ( i = 0; i < wxsize; i++ ) wx[i] = 0;
                                         }            }
                                 }          }
                         }        }
         }    }
 }  }
   
 int f4_nocheck;  int f4_nocheck;
   
 int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)  int gensolve_check(MAT mat,MAT nm,Q dn,int *rind,int *cind)
 {  {
         int row,col,rank,clen,i,j,k,l;    int row,col,rank,clen,i,j,k,l;
         Q s,t;    Q s,t;
         Q *w;    Q *w;
         Q *mati,*nmk;    Q *mati,*nmk;
   
         if ( f4_nocheck )    if ( f4_nocheck )
                 return 1;      return 1;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         rank = nm->row; clen = nm->col;    rank = nm->row; clen = nm->col;
         w = (Q *)MALLOC(clen*sizeof(Q));    w = (Q *)MALLOC(clen*sizeof(Q));
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 mati = (Q *)mat->body[i];      mati = (Q *)mat->body[i];
 #if 1  #if 1
                 bzero(w,clen*sizeof(Q));      bzero(w,clen*sizeof(Q));
                 for ( k = 0; k < rank; k++ )      for ( k = 0; k < rank; k++ )
                         for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {        for ( l = 0, nmk = (Q *)nm->body[k]; l < clen; l++ ) {
                                 mulq(mati[rind[k]],nmk[l],&t);          mulq(mati[rind[k]],nmk[l],&t);
                                 addq(w[l],t,&s); w[l] = s;          addq(w[l],t,&s); w[l] = s;
                         }        }
                 for ( j = 0; j < clen; j++ ) {      for ( j = 0; j < clen; j++ ) {
                         mulq(dn,mati[cind[j]],&t);        mulq(dn,mati[cind[j]],&t);
                         if ( cmpq(w[j],t) )        if ( cmpq(w[j],t) )
                                 break;          break;
                 }      }
 #else  #else
                 for ( j = 0; j < clen; j++ ) {      for ( j = 0; j < clen; j++ ) {
                         for ( k = 0, s = 0; k < rank; k++ ) {        for ( k = 0, s = 0; k < rank; k++ ) {
                                 mulq(mati[rind[k]],nm->body[k][j],&t);          mulq(mati[rind[k]],nm->body[k][j],&t);
                                 addq(s,t,&u); s = u;          addq(s,t,&u); s = u;
                         }        }
                         mulq(dn,mati[cind[j]],&t);        mulq(dn,mati[cind[j]],&t);
                         if ( cmpq(s,t) )        if ( cmpq(s,t) )
                                 break;          break;
                 }      }
 #endif  #endif
                 if ( j != clen )      if ( j != clen )
                         break;        break;
         }    }
         if ( i != row )    if ( i != row )
                 return 0;      return 0;
         else    else
                 return 1;      return 1;
 }  }
   
 /* assuming 0 < c < m */  /* assuming 0 < c < m */
   
 int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)  int inttorat(N c,N m,N b,int *sgnp,N *nmp,N *dnp)
 {  {
         Q qq,t,u1,v1,r1;    Q qq,t,u1,v1,r1;
         N q,u2,v2,r2;    N q,u2,v2,r2;
   
         u1 = 0; v1 = ONE; u2 = m; v2 = c;    u1 = 0; v1 = ONE; u2 = m; v2 = c;
         while ( cmpn(v2,b) >= 0 ) {    while ( cmpn(v2,b) >= 0 ) {
                 divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;      divn(u2,v2,&q,&r2); u2 = v2; v2 = r2;
                 NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;      NTOQ(q,1,qq); mulq(qq,v1,&t); subq(u1,t,&r1); u1 = v1; v1 = r1;
         }    }
         if ( cmpn(NM(v1),b) >= 0 )    if ( cmpn(NM(v1),b) >= 0 )
                 return 0;      return 0;
         else {    else {
                 *nmp = v2;      *nmp = v2;
                 *dnp = NM(v1);      *dnp = NM(v1);
                 *sgnp = SGN(v1);      *sgnp = SGN(v1);
                 return 1;      return 1;
         }    }
 }  }
   
 /* mat->body = N ** */  /* mat->body = N ** */
   
 int intmtoratm(MAT mat,N md,MAT nm,Q *dn)  int intmtoratm(MAT mat,N md,MAT nm,Q *dn)
 {  {
         N t,s,b;    N t,s,b;
         Q dn0,dn1,nm1,q;    Q dn0,dn1,nm1,q;
         int i,j,k,l,row,col;    int i,j,k,l,row,col;
         Q **rmat;    Q **rmat;
         N **tmat;    N **tmat;
         N *tmi;    N *tmi;
         Q *nmk;    Q *nmk;
         N u,unm,udn;    N u,unm,udn;
         int sgn,ret;    int sgn,ret;
   
         if ( UNIN(md) )    if ( UNIN(md) )
                 return 0;      return 0;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         bshiftn(md,1,&t);    bshiftn(md,1,&t);
         isqrt(t,&s);    isqrt(t,&s);
         bshiftn(s,64,&b);    bshiftn(s,64,&b);
         if ( !b )    if ( !b )
                 b = ONEN;      b = ONEN;
         dn0 = ONE;    dn0 = ONE;
         tmat = (N **)mat->body;    tmat = (N **)mat->body;
         rmat = (Q **)nm->body;    rmat = (Q **)nm->body;
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0, tmi = tmat[i]; j < col; j++ )      for ( j = 0, tmi = tmat[i]; j < col; j++ )
                         if ( tmi[j] ) {        if ( tmi[j] ) {
                                 muln(tmi[j],NM(dn0),&s);          muln(tmi[j],NM(dn0),&s);
                                 remn(s,md,&u);          remn(s,md,&u);
                                 ret = inttorat(u,md,b,&sgn,&unm,&udn);          ret = inttorat(u,md,b,&sgn,&unm,&udn);
                                 if ( !ret )          if ( !ret )
                                         return 0;            return 0;
                                 else {          else {
                                         NTOQ(unm,sgn,nm1);            NTOQ(unm,sgn,nm1);
                                         NTOQ(udn,1,dn1);            NTOQ(udn,1,dn1);
                                         if ( !UNIQ(dn1) ) {            if ( !UNIQ(dn1) ) {
                                                 for ( k = 0; k < i; k++ )              for ( k = 0; k < i; k++ )
                                                         for ( l = 0, nmk = rmat[k]; l < col; l++ ) {                for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
                                                                 mulq(nmk[l],dn1,&q); nmk[l] = q;                  mulq(nmk[l],dn1,&q); nmk[l] = q;
                                                         }                }
                                                 for ( l = 0, nmk = rmat[i]; l < j; l++ ) {              for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
                                                         mulq(nmk[l],dn1,&q); nmk[l] = q;                mulq(nmk[l],dn1,&q); nmk[l] = q;
                                                 }              }
                                         }            }
                                         rmat[i][j] = nm1;            rmat[i][j] = nm1;
                                         mulq(dn0,dn1,&q); dn0 = q;            mulq(dn0,dn1,&q); dn0 = q;
                                 }          }
                         }        }
         *dn = dn0;    *dn = dn0;
         return 1;    return 1;
 }  }
   
 /* mat->body = Q ** */  /* mat->body = Q ** */
   
 int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)  int intmtoratm_q(MAT mat,N md,MAT nm,Q *dn)
 {  {
         N t,s,b;    N t,s,b;
         Q dn0,dn1,nm1,q;    Q dn0,dn1,nm1,q;
         int i,j,k,l,row,col;    int i,j,k,l,row,col;
         Q **rmat;    Q **rmat;
         Q **tmat;    Q **tmat;
         Q *tmi;    Q *tmi;
         Q *nmk;    Q *nmk;
         N u,unm,udn;    N u,unm,udn;
         int sgn,ret;    int sgn,ret;
   
         if ( UNIN(md) )    if ( UNIN(md) )
                 return 0;      return 0;
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         bshiftn(md,1,&t);    bshiftn(md,1,&t);
         isqrt(t,&s);    isqrt(t,&s);
         bshiftn(s,64,&b);    bshiftn(s,64,&b);
         if ( !b )    if ( !b )
                 b = ONEN;      b = ONEN;
         dn0 = ONE;    dn0 = ONE;
         tmat = (Q **)mat->body;    tmat = (Q **)mat->body;
         rmat = (Q **)nm->body;    rmat = (Q **)nm->body;
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0, tmi = tmat[i]; j < col; j++ )      for ( j = 0, tmi = tmat[i]; j < col; j++ )
                         if ( tmi[j] ) {        if ( tmi[j] ) {
                                 muln(NM(tmi[j]),NM(dn0),&s);          muln(NM(tmi[j]),NM(dn0),&s);
                                 remn(s,md,&u);          remn(s,md,&u);
                                 ret = inttorat(u,md,b,&sgn,&unm,&udn);          ret = inttorat(u,md,b,&sgn,&unm,&udn);
                                 if ( !ret )          if ( !ret )
                                         return 0;            return 0;
                                 else {          else {
                                         if ( SGN(tmi[j])<0 )            if ( SGN(tmi[j])<0 )
                                                 sgn = -sgn;              sgn = -sgn;
                                         NTOQ(unm,sgn,nm1);            NTOQ(unm,sgn,nm1);
                                         NTOQ(udn,1,dn1);            NTOQ(udn,1,dn1);
                                         if ( !UNIQ(dn1) ) {            if ( !UNIQ(dn1) ) {
                                                 for ( k = 0; k < i; k++ )              for ( k = 0; k < i; k++ )
                                                         for ( l = 0, nmk = rmat[k]; l < col; l++ ) {                for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
                                                                 mulq(nmk[l],dn1,&q); nmk[l] = q;                  mulq(nmk[l],dn1,&q); nmk[l] = q;
                                                         }                }
                                                 for ( l = 0, nmk = rmat[i]; l < j; l++ ) {              for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
                                                         mulq(nmk[l],dn1,&q); nmk[l] = q;                mulq(nmk[l],dn1,&q); nmk[l] = q;
                                                 }              }
                                         }            }
                                         rmat[i][j] = nm1;            rmat[i][j] = nm1;
                                         mulq(dn0,dn1,&q); dn0 = q;            mulq(dn0,dn1,&q); dn0 = q;
                                 }          }
                         }        }
         *dn = dn0;    *dn = dn0;
         return 1;    return 1;
 }  }
   
 #define ONE_STEP1  if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;  #define ONE_STEP1  if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
   
 void reduce_reducers_mod(int **mat,int row,int col,int md)  void reduce_reducers_mod(int **mat,int row,int col,int md)
 {  {
         int i,j,k,l,hc,zzz;    int i,j,k,l,hc,zzz;
         int *t,*s,*tj,*ind;    int *t,*s,*tj,*ind;
   
         /* reduce the reducers */    /* reduce the reducers */
         ind = (int *)ALLOCA(row*sizeof(int));    ind = (int *)ALLOCA(row*sizeof(int));
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 t = mat[i];      t = mat[i];
                 for ( j = 0; j < col && !t[j]; j++ );      for ( j = 0; j < col && !t[j]; j++ );
                 /* register the position of the head term */      /* register the position of the head term */
                 ind[i] = j;      ind[i] = j;
                 for ( l = i-1; l >= 0; l-- ) {      for ( l = i-1; l >= 0; l-- ) {
                         /* reduce mat[i] by mat[l] */        /* reduce mat[i] by mat[l] */
                         if ( hc = t[ind[l]] ) {        if ( hc = t[ind[l]] ) {
                                 /* mat[i] = mat[i]-hc*mat[l] */          /* mat[i] = mat[i]-hc*mat[l] */
                                 j = ind[l];          j = ind[l];
                                 s = mat[l]+j;          s = mat[l]+j;
                                 tj = t+j;          tj = t+j;
                                 hc = md-hc;          hc = md-hc;
                                 k = col-j;          k = col-j;
                                 for ( ; k >= 64; k -= 64 ) {          for ( ; k >= 64; k -= 64 ) {
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                         ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1            ONE_STEP1 ONE_STEP1 ONE_STEP1 ONE_STEP1
                                 }          }
                                 for ( ; k > 0; k-- ) {          for ( ; k > 0; k-- ) {
                                         if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;            if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
                                 }          }
                         }        }
                 }      }
         }    }
 }  }
   
 /*  /*
         mat[i] : reducers (i=0,...,nred-1)    mat[i] : reducers (i=0,...,nred-1)
                  spolys (i=nred,...,row-1)             spolys (i=nred,...,row-1)
         mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order    mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
         1. reduce the reducers    1. reduce the reducers
         2. reduce spolys by the reduced reducers    2. reduce spolys by the reduced reducers
 */  */
   
 void pre_reduce_mod(int **mat,int row,int col,int nred,int md)  void pre_reduce_mod(int **mat,int row,int col,int nred,int md)
 {  {
         int i,j,k,l,hc,inv;    int i,j,k,l,hc,inv;
         int *t,*s,*tk,*ind;    int *t,*s,*tk,*ind;
   
 #if 1  #if 1
         /* reduce the reducers */    /* reduce the reducers */
         ind = (int *)ALLOCA(row*sizeof(int));    ind = (int *)ALLOCA(row*sizeof(int));
         for ( i = 0; i < nred; i++ ) {    for ( i = 0; i < nred; i++ ) {
                 /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */      /* make mat[i] monic and mat[i] by mat[0],...,mat[i-1] */
                 t = mat[i];      t = mat[i];
                 for ( j = 0; j < col && !t[j]; j++ );      for ( j = 0; j < col && !t[j]; j++ );
                 /* register the position of the head term */      /* register the position of the head term */
                 ind[i] = j;      ind[i] = j;
                 inv = invm(t[j],md);      inv = invm(t[j],md);
                 for ( k = j; k < col; k++ )      for ( k = j; k < col; k++ )
                         if ( t[k] )        if ( t[k] )
                                 DMAR(t[k],inv,0,md,t[k])          DMAR(t[k],inv,0,md,t[k])
                 for ( l = i-1; l >= 0; l-- ) {      for ( l = i-1; l >= 0; l-- ) {
                         /* reduce mat[i] by mat[l] */        /* reduce mat[i] by mat[l] */
                         if ( hc = t[ind[l]] ) {        if ( hc = t[ind[l]] ) {
                                 /* mat[i] = mat[i]-hc*mat[l] */          /* mat[i] = mat[i]-hc*mat[l] */
                                 for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;          for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
                                         k < col; k++, tk++, s++ )            k < col; k++, tk++, s++ )
                                         if ( *s )            if ( *s )
                                                 DMAR(*s,hc,*tk,md,*tk)              DMAR(*s,hc,*tk,md,*tk)
                         }        }
                 }      }
         }    }
         /* reduce the spolys */    /* reduce the spolys */
         for ( i = nred; i < row; i++ ) {    for ( i = nred; i < row; i++ ) {
                 t = mat[i];      t = mat[i];
                 for ( l = nred-1; l >= 0; l-- ) {      for ( l = nred-1; l >= 0; l-- ) {
                         /* reduce mat[i] by mat[l] */        /* reduce mat[i] by mat[l] */
                         if ( hc = t[ind[l]] ) {        if ( hc = t[ind[l]] ) {
                                 /* mat[i] = mat[i]-hc*mat[l] */          /* mat[i] = mat[i]-hc*mat[l] */
                                 for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;          for ( k = ind[l], hc = md-hc, s = mat[l]+k, tk = t+k;
                                         k < col; k++, tk++, s++ )            k < col; k++, tk++, s++ )
                                         if ( *s )            if ( *s )
                                                 DMAR(*s,hc,*tk,md,*tk)              DMAR(*s,hc,*tk,md,*tk)
                         }        }
                 }      }
         }    }
 #endif  #endif
 }  }
 /*  /*
         mat[i] : reducers (i=0,...,nred-1)    mat[i] : reducers (i=0,...,nred-1)
         mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order    mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
 */  */
   
 void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)  void reduce_sp_by_red_mod(int *sp,int **redmat,int *ind,int nred,int col,int md)
 {  {
         int i,j,k,hc,zzz;    int i,j,k,hc,zzz;
         int *s,*tj;    int *s,*tj;
   
         /* reduce the spolys by redmat */    /* reduce the spolys by redmat */
         for ( i = nred-1; i >= 0; i-- ) {    for ( i = nred-1; i >= 0; i-- ) {
                 /* reduce sp by redmat[i] */      /* reduce sp by redmat[i] */
                 if ( hc = sp[ind[i]] ) {      if ( hc = sp[ind[i]] ) {
                         /* sp = sp-hc*redmat[i] */        /* sp = sp-hc*redmat[i] */
                         j = ind[i];        j = ind[i];
                         hc = md-hc;        hc = md-hc;
                         s = redmat[i]+j;        s = redmat[i]+j;
                         tj = sp+j;        tj = sp+j;
                         for ( k = col-j; k > 0; k-- ) {        for ( k = col-j; k > 0; k-- ) {
                                 if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;          if ( zzz = *s ) { DMAR(zzz,hc,*tj,md,*tj) } tj++; s++;
                         }        }
                 }      }
         }    }
 }  }
   
 /*  /*
         mat[i] : compressed reducers (i=0,...,nred-1)    mat[i] : compressed reducers (i=0,...,nred-1)
         mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order    mat[0] < mat[1] < ... < mat[nred-1] w.r.t the term order
 */  */
   
 void red_by_compress(int m,unsigned int *p,unsigned int *r,  void red_by_compress(int m,unsigned int *p,unsigned int *r,
         unsigned int *ri,unsigned int hc,int len)    unsigned int *ri,unsigned int hc,int len)
 {  {
         unsigned int up,lo;    unsigned int up,lo;
         unsigned int dmy;    unsigned int dmy;
         unsigned int *pj;    unsigned int *pj;
   
         p[*ri] = 0; r++; ri++;    p[*ri] = 0; r++; ri++;
         for ( len--; len; len--, r++, ri++ ) {    for ( len--; len; len--, r++, ri++ ) {
                 pj = p+ *ri;      pj = p+ *ri;
                 DMA(*r,hc,*pj,up,lo);      DMA(*r,hc,*pj,up,lo);
                 if ( up ) {      if ( up ) {
                         DSAB(m,up,lo,dmy,*pj);        DSAB(m,up,lo,dmy,*pj);
                 } else      } else
                         *pj = lo;        *pj = lo;
         }    }
 }  }
   
 /* p -= hc*r */  /* p -= hc*r */
   
 void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)  void red_by_vect(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
 {  {
         unsigned int up,lo,dmy;    unsigned int up,lo,dmy;
   
         *p++ = 0; r++; len--;    *p++ = 0; r++; len--;
         for ( ; len; len--, r++, p++ )    for ( ; len; len--, r++, p++ )
                 if ( *r ) {      if ( *r ) {
                         DMA(*r,hc,*p,up,lo);        DMA(*r,hc,*p,up,lo);
                         if ( up ) {        if ( up ) {
                                 DSAB(m,up,lo,dmy,*p);          DSAB(m,up,lo,dmy,*p);
                         } else        } else
                                 *p = lo;          *p = lo;
                 }      }
 }  }
   
 #if defined(__GNUC__) && SIZEOF_LONG==8  #if defined(__GNUC__) && SIZEOF_LONG==8
Line 2376  void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *
Line 2376  void red_by_vect64(int m, U64 *p,unsigned int *c,U64 *
   
 void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)  void red_by_vect_sf(int m,unsigned int *p,unsigned int *r,unsigned int hc,int len)
 {  {
         *p++ = 0; r++; len--;    *p++ = 0; r++; len--;
         for ( ; len; len--, r++, p++ )    for ( ; len; len--, r++, p++ )
                 if ( *r )      if ( *r )
                         *p = _addsf(_mulsf(*r,hc),*p);        *p = _addsf(_mulsf(*r,hc),*p);
 }  }
   
 extern GZ current_mod_lf;  extern GZ current_mod_lf;
Line 2387  extern int current_mod_lf_size;
Line 2387  extern int current_mod_lf_size;
   
 void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)  void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len)
 {  {
         mpz_set_ui(*p++,0); r++; len--;    mpz_set_ui(*p++,0); r++; len--;
         for ( ; len; len--, r++, p++ ) {    for ( ; len; len--, r++, p++ ) {
        mpz_addmul(*p,*r,hc);         mpz_addmul(*p,*r,hc);
 #if 0  #if 0
        if ( mpz_size(*p) > current_mod_lf_size )         if ( mpz_size(*p) > current_mod_lf_size )
Line 2401  void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len
Line 2401  void red_by_vect_lf(mpz_t *p,mpz_t *r,mpz_t hc,int len
 extern unsigned int **psca;  extern unsigned int **psca;
   
 void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,  void reduce_sp_by_red_mod_compress (int *sp,CDP *redmat,int *ind,
         int nred,int col,int md)    int nred,int col,int md)
 {  {
         int i,len;    int i,len;
         CDP ri;    CDP ri;
         unsigned int hc;    unsigned int hc;
         unsigned int *usp;    unsigned int *usp;
   
         usp = (unsigned int *)sp;    usp = (unsigned int *)sp;
         /* reduce the spolys by redmat */    /* reduce the spolys by redmat */
         for ( i = nred-1; i >= 0; i-- ) {    for ( i = nred-1; i >= 0; i-- ) {
                 /* reduce sp by redmat[i] */      /* reduce sp by redmat[i] */
                 usp[ind[i]] %= md;      usp[ind[i]] %= md;
                 if ( hc = usp[ind[i]] ) {      if ( hc = usp[ind[i]] ) {
                         /* sp = sp-hc*redmat[i] */        /* sp = sp-hc*redmat[i] */
                         hc = md-hc;        hc = md-hc;
                         ri = redmat[i];        ri = redmat[i];
                         len = ri->len;        len = ri->len;
                         red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);        red_by_compress(md,usp,psca[ri->psindex],ri->body,hc,len);
                 }      }
         }    }
         for ( i = 0; i < col; i++ )    for ( i = 0; i < col; i++ )
                 if ( usp[i] >= (unsigned int)md )      if ( usp[i] >= (unsigned int)md )
                         usp[i] %= md;        usp[i] %= md;
 }  }
   
 #define ONE_STEP2  if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;  #define ONE_STEP2  if ( zzz = *pk ) { DMAR(zzz,a,*tk,md,*tk) } pk++; tk++;
   
 int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)  int generic_gauss_elim_mod(int **mat0,int row,int col,int md,int *colstat)
 {  {
         int i,j,k,l,inv,a,rank;    int i,j,k,l,inv,a,rank;
         unsigned int *t,*pivot,*pk;    unsigned int *t,*pivot,*pk;
         unsigned int **mat;    unsigned int **mat;
   
         mat = (unsigned int **)mat0;    mat = (unsigned int **)mat0;
         for ( rank = 0, j = 0; j < col; j++ ) {    for ( rank = 0, j = 0; j < col; j++ ) {
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         mat[i][j] %= md;        mat[i][j] %= md;
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         if ( mat[i][j] )        if ( mat[i][j] )
                                 break;          break;
                 if ( i == row ) {      if ( i == row ) {
                         colstat[j] = 0;        colstat[j] = 0;
                         continue;        continue;
                 } else      } else
                         colstat[j] = 1;        colstat[j] = 1;
                 if ( i != rank ) {      if ( i != rank ) {
                         t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;        t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                 }      }
                 pivot = mat[rank];      pivot = mat[rank];
                 inv = invm(pivot[j],md);      inv = invm(pivot[j],md);
                 for ( k = j, pk = pivot+k; k < col; k++, pk++ )      for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                         if ( *pk ) {        if ( *pk ) {
                                 if ( *pk >= (unsigned int)md )          if ( *pk >= (unsigned int)md )
                                         *pk %= md;            *pk %= md;
                                 DMAR(*pk,inv,0,md,*pk)          DMAR(*pk,inv,0,md,*pk)
                         }        }
                 for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 red_by_vect(md,t+j,pivot+j,md-a,col-j);          red_by_vect(md,t+j,pivot+j,md-a,col-j);
                 }      }
                 rank++;      rank++;
         }    }
         for ( j = col-1, l = rank-1; j >= 0; j-- )    for ( j = col-1, l = rank-1; j >= 0; j-- )
                 if ( colstat[j] ) {      if ( colstat[j] ) {
                         pivot = mat[l];        pivot = mat[l];
                         for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
                                 t = mat[i];          t = mat[i];
                                 t[j] %= md;          t[j] %= md;
                                 if ( a = t[j] )          if ( a = t[j] )
                                         red_by_vect(md,t+j,pivot+j,md-a,col-j);            red_by_vect(md,t+j,pivot+j,md-a,col-j);
                         }        }
                         l--;        l--;
                 }      }
         for ( j = 0, l = 0; l < rank; j++ )    for ( j = 0, l = 0; l < rank; j++ )
                 if ( colstat[j] ) {      if ( colstat[j] ) {
                         t = mat[l];        t = mat[l];
                         for ( k = j; k < col; k++ )        for ( k = j; k < col; k++ )
                                 if ( t[k] >= (unsigned int)md )          if ( t[k] >= (unsigned int)md )
                                         t[k] %= md;            t[k] %= md;
                         l++;        l++;
                 }      }
         return rank;    return rank;
 }  }
   
 int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)  int generic_gauss_elim_mod2(int **mat0,int row,int col,int md,int *colstat,int *rowstat)
 {  {
         int i,j,k,l,inv,a,rank;    int i,j,k,l,inv,a,rank;
         unsigned int *t,*pivot,*pk;    unsigned int *t,*pivot,*pk;
         unsigned int **mat;    unsigned int **mat;
   
         for ( i = 0; i < row; i++ ) rowstat[i] = i;    for ( i = 0; i < row; i++ ) rowstat[i] = i;
         mat = (unsigned int **)mat0;    mat = (unsigned int **)mat0;
         for ( rank = 0, j = 0; j < col; j++ ) {    for ( rank = 0, j = 0; j < col; j++ ) {
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         mat[i][j] %= md;        mat[i][j] %= md;
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         if ( mat[i][j] )        if ( mat[i][j] )
                                 break;          break;
                 if ( i == row ) {      if ( i == row ) {
                         colstat[j] = 0;        colstat[j] = 0;
                         continue;        continue;
                 } else      } else
                         colstat[j] = 1;        colstat[j] = 1;
                 if ( i != rank ) {      if ( i != rank ) {
                         t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;        t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                         k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;        k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
                 }      }
                 pivot = mat[rank];      pivot = mat[rank];
                 inv = invm(pivot[j],md);      inv = invm(pivot[j],md);
                 for ( k = j, pk = pivot+k; k < col; k++, pk++ )      for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                         if ( *pk ) {        if ( *pk ) {
                                 if ( *pk >= (unsigned int)md )          if ( *pk >= (unsigned int)md )
                                         *pk %= md;            *pk %= md;
                                 DMAR(*pk,inv,0,md,*pk)          DMAR(*pk,inv,0,md,*pk)
                         }        }
                 for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 red_by_vect(md,t+j,pivot+j,md-a,col-j);          red_by_vect(md,t+j,pivot+j,md-a,col-j);
                 }      }
                 rank++;      rank++;
         }    }
         for ( j = col-1, l = rank-1; j >= 0; j-- )    for ( j = col-1, l = rank-1; j >= 0; j-- )
                 if ( colstat[j] ) {      if ( colstat[j] ) {
                         pivot = mat[l];        pivot = mat[l];
                         for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
                                 t = mat[i];          t = mat[i];
                                 t[j] %= md;          t[j] %= md;
                                 if ( a = t[j] )          if ( a = t[j] )
                                         red_by_vect(md,t+j,pivot+j,md-a,col-j);            red_by_vect(md,t+j,pivot+j,md-a,col-j);
                         }        }
                         l--;        l--;
                 }      }
         for ( j = 0, l = 0; l < rank; j++ )    for ( j = 0, l = 0; l < rank; j++ )
                 if ( colstat[j] ) {      if ( colstat[j] ) {
                         t = mat[l];        t = mat[l];
                         for ( k = j; k < col; k++ )        for ( k = j; k < col; k++ )
                                 if ( t[k] >= (unsigned int)md )          if ( t[k] >= (unsigned int)md )
                                         t[k] %= md;            t[k] %= md;
                         l++;        l++;
                 }      }
         return rank;    return rank;
 }  }
   
 int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)  int indep_rows_mod(int **mat0,int row,int col,int md,int *rowstat)
 {  {
         int i,j,k,l,inv,a,rank;    int i,j,k,l,inv,a,rank;
         unsigned int *t,*pivot,*pk;    unsigned int *t,*pivot,*pk;
         unsigned int **mat;    unsigned int **mat;
   
         for ( i = 0; i < row; i++ ) rowstat[i] = i;    for ( i = 0; i < row; i++ ) rowstat[i] = i;
         mat = (unsigned int **)mat0;    mat = (unsigned int **)mat0;
         for ( rank = 0, j = 0; j < col; j++ ) {    for ( rank = 0, j = 0; j < col; j++ ) {
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         mat[i][j] %= md;        mat[i][j] %= md;
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         if ( mat[i][j] )        if ( mat[i][j] )
                                 break;          break;
                 if ( i == row ) continue;      if ( i == row ) continue;
                 if ( i != rank ) {      if ( i != rank ) {
                         t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;        t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                         k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;        k = rowstat[i]; rowstat[i] = rowstat[rank]; rowstat[rank] = k;
                 }      }
                 pivot = mat[rank];      pivot = mat[rank];
                 inv = invm(pivot[j],md);      inv = invm(pivot[j],md);
                 for ( k = j, pk = pivot+k; k < col; k++, pk++ )      for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                         if ( *pk ) {        if ( *pk ) {
                                 if ( *pk >= (unsigned int)md )          if ( *pk >= (unsigned int)md )
                                         *pk %= md;            *pk %= md;
                                 DMAR(*pk,inv,0,md,*pk)          DMAR(*pk,inv,0,md,*pk)
                         }        }
                 for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 red_by_vect(md,t+j,pivot+j,md-a,col-j);          red_by_vect(md,t+j,pivot+j,md-a,col-j);
                 }      }
                 rank++;      rank++;
         }    }
         return rank;    return rank;
 }  }
   
 int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)  int generic_gauss_elim_sf(int **mat0,int row,int col,int md,int *colstat)
 {  {
         int i,j,k,l,inv,a,rank;    int i,j,k,l,inv,a,rank;
         unsigned int *t,*pivot,*pk;    unsigned int *t,*pivot,*pk;
         unsigned int **mat;    unsigned int **mat;
   
         mat = (unsigned int **)mat0;    mat = (unsigned int **)mat0;
         for ( rank = 0, j = 0; j < col; j++ ) {    for ( rank = 0, j = 0; j < col; j++ ) {
                 for ( i = rank; i < row; i++ )      for ( i = rank; i < row; i++ )
                         if ( mat[i][j] )        if ( mat[i][j] )
                                 break;          break;
                 if ( i == row ) {      if ( i == row ) {
                         colstat[j] = 0;        colstat[j] = 0;
                         continue;        continue;
                 } else      } else
                         colstat[j] = 1;        colstat[j] = 1;
                 if ( i != rank ) {      if ( i != rank ) {
                         t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;        t = mat[i]; mat[i] = mat[rank]; mat[rank] = t;
                 }      }
                 pivot = mat[rank];      pivot = mat[rank];
                 inv = _invsf(pivot[j]);      inv = _invsf(pivot[j]);
                 for ( k = j, pk = pivot+k; k < col; k++, pk++ )      for ( k = j, pk = pivot+k; k < col; k++, pk++ )
                         if ( *pk )        if ( *pk )
                                 *pk = _mulsf(*pk,inv);          *pk = _mulsf(*pk,inv);
                 for ( i = rank+1; i < row; i++ ) {      for ( i = rank+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);          red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
                 }      }
                 rank++;      rank++;
         }    }
         for ( j = col-1, l = rank-1; j >= 0; j-- )    for ( j = col-1, l = rank-1; j >= 0; j-- )
                 if ( colstat[j] ) {      if ( colstat[j] ) {
                         pivot = mat[l];        pivot = mat[l];
                         for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
                                 t = mat[i];          t = mat[i];
                                 if ( a = t[j] )          if ( a = t[j] )
                                         red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);            red_by_vect_sf(md,t+j,pivot+j,_chsgnsf(a),col-j);
                         }        }
                         l--;        l--;
                 }      }
         return rank;    return rank;
 }  }
   
 /* LU decomposition; a[i][i] = 1/U[i][i] */  /* LU decomposition; a[i][i] = 1/U[i][i] */
   
 int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)  int lu_gfmmat(GFMMAT mat,unsigned int md,int *perm)
 {  {
         int row,col;    int row,col;
         int i,j,k;    int i,j,k;
         unsigned int *t,*pivot;    unsigned int *t,*pivot;
         unsigned int **a;    unsigned int **a;
         unsigned int inv,m;    unsigned int inv,m;
   
         row = mat->row; col = mat->col;    row = mat->row; col = mat->col;
         a = mat->body;    a = mat->body;
         bzero(perm,row*sizeof(int));    bzero(perm,row*sizeof(int));
   
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 perm[i] = i;      perm[i] = i;
         for ( k = 0; k < col; k++ ) {    for ( k = 0; k < col; k++ ) {
                 for ( i = k; i < row && !a[i][k]; i++ );      for ( i = k; i < row && !a[i][k]; i++ );
                 if ( i == row )      if ( i == row )
                         return 0;        return 0;
                 if ( i != k ) {      if ( i != k ) {
                         j = perm[i]; perm[i] = perm[k]; perm[k] = j;        j = perm[i]; perm[i] = perm[k]; perm[k] = j;
                         t = a[i]; a[i] = a[k]; a[k] = t;        t = a[i]; a[i] = a[k]; a[k] = t;
                 }      }
                 pivot = a[k];      pivot = a[k];
                 pivot[k] = inv = invm(pivot[k],md);      pivot[k] = inv = invm(pivot[k],md);
                 for ( i = k+1; i < row; i++ ) {      for ( i = k+1; i < row; i++ ) {
                         t = a[i];        t = a[i];
                         if ( m = t[k] ) {        if ( m = t[k] ) {
                                 DMAR(inv,m,0,md,t[k])          DMAR(inv,m,0,md,t[k])
                                 for ( j = k+1, m = md - t[k]; j < col; j++ )          for ( j = k+1, m = md - t[k]; j < col; j++ )
                                         if ( pivot[j] ) {            if ( pivot[j] ) {
                                                 unsigned int tj;              unsigned int tj;
   
                                                 DMAR(m,pivot[j],t[j],md,tj)              DMAR(m,pivot[j],t[j],md,tj)
                                                 t[j] = tj;              t[j] = tj;
                                         }            }
                         }        }
                 }      }
         }    }
         return 1;    return 1;
 }  }
   
 /*  /*
  Input   Input
         a: a row x col matrix    a: a row x col matrix
         md : a modulus    md : a modulus
   
  Output:   Output:
         return : d = the rank of mat    return : d = the rank of mat
         a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])    a[0..(d-1)][0..(d-1)] : LU decomposition (a[i][i] = 1/U[i][i])
         rinfo: array of length row    rinfo: array of length row
         cinfo: array of length col    cinfo: array of length col
     i-th row in new a <-> rinfo[i]-th row in old a      i-th row in new a <-> rinfo[i]-th row in old a
         cinfo[j]=1 <=> j-th column is contained in the LU decomp.    cinfo[j]=1 <=> j-th column is contained in the LU decomp.
 */  */
   
 int find_lhs_and_lu_mod(unsigned int **a,int row,int col,  int find_lhs_and_lu_mod(unsigned int **a,int row,int col,
         unsigned int md,int **rinfo,int **cinfo)    unsigned int md,int **rinfo,int **cinfo)
 {  {
         int i,j,k,d;    int i,j,k,d;
         int *rp,*cp;    int *rp,*cp;
         unsigned int *t,*pivot;    unsigned int *t,*pivot;
         unsigned int inv,m;    unsigned int inv,m;
   
         *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));    *rinfo = rp = (int *)MALLOC_ATOMIC(row*sizeof(int));
         *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));    *cinfo = cp = (int *)MALLOC_ATOMIC(col*sizeof(int));
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 rp[i] = i;      rp[i] = i;
         for ( k = 0, d = 0; k < col; k++ ) {    for ( k = 0, d = 0; k < col; k++ ) {
                 for ( i = d; i < row && !a[i][k]; i++ );      for ( i = d; i < row && !a[i][k]; i++ );
                 if ( i == row ) {      if ( i == row ) {
                         cp[k] = 0;        cp[k] = 0;
                         continue;        continue;
                 } else      } else
                         cp[k] = 1;        cp[k] = 1;
                 if ( i != d ) {      if ( i != d ) {
                         j = rp[i]; rp[i] = rp[d]; rp[d] = j;        j = rp[i]; rp[i] = rp[d]; rp[d] = j;
                         t = a[i]; a[i] = a[d]; a[d] = t;        t = a[i]; a[i] = a[d]; a[d] = t;
                 }      }
                 pivot = a[d];      pivot = a[d];
                 pivot[k] = inv = invm(pivot[k],md);      pivot[k] = inv = invm(pivot[k],md);
                 for ( i = d+1; i < row; i++ ) {      for ( i = d+1; i < row; i++ ) {
                         t = a[i];        t = a[i];
                         if ( m = t[k] ) {        if ( m = t[k] ) {
                                 DMAR(inv,m,0,md,t[k])          DMAR(inv,m,0,md,t[k])
                                 for ( j = k+1, m = md - t[k]; j < col; j++ )          for ( j = k+1, m = md - t[k]; j < col; j++ )
                                         if ( pivot[j] ) {            if ( pivot[j] ) {
                                                 unsigned int tj;              unsigned int tj;
                                                 DMAR(m,pivot[j],t[j],md,tj)              DMAR(m,pivot[j],t[j],md,tj)
                                                 t[j] = tj;              t[j] = tj;
                                         }            }
                         }        }
                 }      }
                 d++;      d++;
         }    }
         return d;    return d;
 }  }
   
 int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)  int lu_mod(unsigned int **a,int n,unsigned int md,int **rinfo)
 {  {
         int i,j,k;    int i,j,k;
         int *rp;    int *rp;
         unsigned int *t,*pivot;    unsigned int *t,*pivot;
         unsigned int inv,m;    unsigned int inv,m;
   
         *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));    *rinfo = rp = (int *)MALLOC_ATOMIC(n*sizeof(int));
         for ( i = 0; i < n; i++ ) rp[i] = i;    for ( i = 0; i < n; i++ ) rp[i] = i;
         for ( k = 0; k < n; k++ ) {    for ( k = 0; k < n; k++ ) {
                 for ( i = k; i < n && !a[i][k]; i++ );      for ( i = k; i < n && !a[i][k]; i++ );
                 if ( i == n ) return 0;      if ( i == n ) return 0;
                 if ( i != k ) {      if ( i != k ) {
                         j = rp[i]; rp[i] = rp[k]; rp[k] = j;        j = rp[i]; rp[i] = rp[k]; rp[k] = j;
                         t = a[i]; a[i] = a[k]; a[k] = t;        t = a[i]; a[i] = a[k]; a[k] = t;
                 }      }
                 pivot = a[k];      pivot = a[k];
                 inv = invm(pivot[k],md);      inv = invm(pivot[k],md);
                 for ( i = k+1; i < n; i++ ) {      for ( i = k+1; i < n; i++ ) {
                         t = a[i];        t = a[i];
                         if ( m = t[k] ) {        if ( m = t[k] ) {
                                 DMAR(inv,m,0,md,t[k])          DMAR(inv,m,0,md,t[k])
                                 for ( j = k+1, m = md - t[k]; j < n; j++ )          for ( j = k+1, m = md - t[k]; j < n; j++ )
                                         if ( pivot[j] ) {            if ( pivot[j] ) {
                                                 unsigned int tj;              unsigned int tj;
                                                 DMAR(m,pivot[j],t[j],md,tj)              DMAR(m,pivot[j],t[j],md,tj)
                                                 t[j] = tj;              t[j] = tj;
                                         }            }
                         }        }
                 }      }
         }    }
         return 1;    return 1;
 }  }
   
 /*  /*
   Input    Input
         a : n x n matrix; a result of LU-decomposition    a : n x n matrix; a result of LU-decomposition
         md : modulus    md : modulus
         b : n x l matrix    b : n x l matrix
  Output   Output
         b = a^(-1)b    b = a^(-1)b
  */   */
   
 void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)  void solve_by_lu_mod(int **a,int n,int md,int **b,int l,int normalize)
 {  {
         unsigned int *y,*c;    unsigned int *y,*c;
         int i,j,k;    int i,j,k;
         unsigned int t,m,m2;    unsigned int t,m,m2;
   
         y = (int *)MALLOC_ATOMIC(n*sizeof(int));    y = (int *)MALLOC_ATOMIC(n*sizeof(int));
         c = (int *)MALLOC_ATOMIC(n*sizeof(int));    c = (int *)MALLOC_ATOMIC(n*sizeof(int));
         m2 = md>>1;    m2 = md>>1;
         for ( k = 0; k < l; k++ ) {    for ( k = 0; k < l; k++ ) {
                 /* copy b[.][k] to c */      /* copy b[.][k] to c */
                 for ( i = 0; i < n; i++ )      for ( i = 0; i < n; i++ )
                         c[i] = (unsigned int)b[i][k];        c[i] = (unsigned int)b[i][k];
                 /* solve Ly=c */      /* solve Ly=c */
                 for ( i = 0; i < n; i++ ) {      for ( i = 0; i < n; i++ ) {
                         for ( t = c[i], j = 0; j < i; j++ )        for ( t = c[i], j = 0; j < i; j++ )
                                 if ( a[i][j] ) {          if ( a[i][j] ) {
                                         m = md - a[i][j];            m = md - a[i][j];
                                         DMAR(m,y[j],t,md,t)            DMAR(m,y[j],t,md,t)
                                 }          }
                         y[i] = t;        y[i] = t;
                 }      }
                 /* solve Uc=y */      /* solve Uc=y */
                 for ( i = n-1; i >= 0; i-- ) {      for ( i = n-1; i >= 0; i-- ) {
                         for ( t = y[i], j =i+1; j < n; j++ )        for ( t = y[i], j =i+1; j < n; j++ )
                                 if ( a[i][j] ) {          if ( a[i][j] ) {
                                         m = md - a[i][j];            m = md - a[i][j];
                                         DMAR(m,c[j],t,md,t)            DMAR(m,c[j],t,md,t)
                                 }          }
                         /* a[i][i] = 1/U[i][i] */        /* a[i][i] = 1/U[i][i] */
                         DMAR(t,a[i][i],0,md,c[i])        DMAR(t,a[i][i],0,md,c[i])
                 }      }
                 /* copy c to b[.][k] with normalization */      /* copy c to b[.][k] with normalization */
                 if ( normalize )      if ( normalize )
                         for ( i = 0; i < n; i++ )        for ( i = 0; i < n; i++ )
                                 b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);          b[i][k] = (int)(c[i]>m2 ? c[i]-md : c[i]);
                 else      else
                         for ( i = 0; i < n; i++ )        for ( i = 0; i < n; i++ )
                                 b[i][k] = c[i];          b[i][k] = c[i];
         }    }
 }  }
   
 void Pleqm1(NODE arg,VECT *rp)  void Pleqm1(NODE arg,VECT *rp)
 {  {
         MAT m;    MAT m;
         VECT vect;    VECT vect;
         pointer **mat;    pointer **mat;
         Q *v;    Q *v;
         Q q;    Q q;
         int **wmat;    int **wmat;
         int md,i,j,row,col,t,n,status;    int md,i,j,row,col,t,n,status;
   
         asir_assert(ARG0(arg),O_MAT,"leqm1");    asir_assert(ARG0(arg),O_MAT,"leqm1");
         asir_assert(ARG1(arg),O_N,"leqm1");    asir_assert(ARG1(arg),O_N,"leqm1");
         m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
         row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
         wmat = (int **)almat(row,col);    wmat = (int **)almat(row,col);
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (Q)mat[i][j] ) {        if ( q = (Q)mat[i][j] ) {
                                 t = rem(NM(q),md);          t = rem(NM(q),md);
                                 if ( SGN(q) < 0 )          if ( SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         } else        } else
                                 wmat[i][j] = 0;          wmat[i][j] = 0;
         status = gauss_elim_mod1(wmat,row,col,md);    status = gauss_elim_mod1(wmat,row,col,md);
         if ( status < 0 )    if ( status < 0 )
                 *rp = 0;      *rp = 0;
         else if ( status > 0 )    else if ( status > 0 )
                 *rp = (VECT)ONE;      *rp = (VECT)ONE;
         else {    else {
                 n = col - 1;      n = col - 1;
                 MKVECT(vect,n);      MKVECT(vect,n);
                 for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {      for ( i = 0, v = (Q *)vect->body; i < n; i++ ) {
                         t = (md-wmat[i][n])%md; STOQ(t,v[i]);        t = (md-wmat[i][n])%md; STOQ(t,v[i]);
                 }      }
                 *rp = vect;      *rp = vect;
         }    }
 }  }
   
 int gauss_elim_mod1(int **mat,int row,int col,int md)  int gauss_elim_mod1(int **mat,int row,int col,int md)
 {  {
         int i,j,k,inv,a,n;    int i,j,k,inv,a,n;
         int *t,*pivot;    int *t,*pivot;
   
         n = col - 1;    n = col - 1;
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 for ( i = j; i < row && !mat[i][j]; i++ );      for ( i = j; i < row && !mat[i][j]; i++ );
                 if ( i == row )      if ( i == row )
                         return 1;        return 1;
                 if ( i != j ) {      if ( i != j ) {
                         t = mat[i]; mat[i] = mat[j]; mat[j] = t;        t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                 }      }
                 pivot = mat[j];      pivot = mat[j];
                 inv = invm(pivot[j],md);      inv = invm(pivot[j],md);
                 for ( k = j; k <= n; k++ )      for ( k = j; k <= n; k++ )
                         pivot[k] = dmar(pivot[k],inv,0,md);        pivot[k] = dmar(pivot[k],inv,0,md);
                 for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( i != j && (a = t[j]) )        if ( i != j && (a = t[j]) )
                                 for ( k = j, a = md - a; k <= n; k++ )          for ( k = j, a = md - a; k <= n; k++ )
                                         t[k] = dmar(pivot[k],a,t[k],md);            t[k] = dmar(pivot[k],a,t[k],md);
                 }      }
         }    }
         for ( i = n; i < row && !mat[i][n]; i++ );    for ( i = n; i < row && !mat[i][n]; i++ );
         if ( i == row ) {    if ( i == row ) {
                 for ( j = n-1; j >= 0; j-- ) {      for ( j = n-1; j >= 0; j-- ) {
                         for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {        for ( i = j-1, a = (md-mat[j][n])%md; i >= 0; i-- ) {
                                 mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);          mat[i][n] = dmar(mat[i][j],a,mat[i][n],md);
                                 mat[i][j] = 0;          mat[i][j] = 0;
                         }        }
                 }      }
                 return 0;      return 0;
         } else    } else
                 return -1;      return -1;
 }  }
   
 void Pgeninvm(NODE arg,LIST *rp)  void Pgeninvm(NODE arg,LIST *rp)
 {  {
         MAT m;    MAT m;
         pointer **mat;    pointer **mat;
         Q **tmat;    Q **tmat;
         Q q;    Q q;
         unsigned int **wmat;    unsigned int **wmat;
         int md,i,j,row,col,t,status;    int md,i,j,row,col,t,status;
         MAT mat1,mat2;    MAT mat1,mat2;
         NODE node1,node2;    NODE node1,node2;
   
         asir_assert(ARG0(arg),O_MAT,"leqm1");    asir_assert(ARG0(arg),O_MAT,"leqm1");
         asir_assert(ARG1(arg),O_N,"leqm1");    asir_assert(ARG1(arg),O_N,"leqm1");
         m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
         row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
         wmat = (unsigned int **)almat(row,col+row);    wmat = (unsigned int **)almat(row,col+row);
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 bzero((char *)wmat[i],(col+row)*sizeof(int));      bzero((char *)wmat[i],(col+row)*sizeof(int));
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (Q)mat[i][j] ) {        if ( q = (Q)mat[i][j] ) {
                                 t = rem(NM(q),md);          t = rem(NM(q),md);
                                 if ( SGN(q) < 0 )          if ( SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         }        }
                 wmat[i][col+i] = 1;      wmat[i][col+i] = 1;
         }    }
         status = gauss_elim_geninv_mod(wmat,row,col,md);    status = gauss_elim_geninv_mod(wmat,row,col,md);
         if ( status > 0 )    if ( status > 0 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);      MKMAT(mat1,col,row); MKMAT(mat2,row-col,row);
                 for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )      for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
                         for ( j = 0; j < row; j++ )        for ( j = 0; j < row; j++ )
                                 UTOQ(wmat[i][j+col],tmat[i][j]);          UTOQ(wmat[i][j+col],tmat[i][j]);
                 for ( tmat = (Q **)mat2->body; i < row; i++ )      for ( tmat = (Q **)mat2->body; i < row; i++ )
                         for ( j = 0; j < row; j++ )        for ( j = 0; j < row; j++ )
                                 UTOQ(wmat[i][j+col],tmat[i-col][j]);          UTOQ(wmat[i][j+col],tmat[i-col][j]);
                 MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);       MKNODE(node2,mat2,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
         }    }
 }  }
   
 int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)  int gauss_elim_geninv_mod(unsigned int **mat,int row,int col,int md)
 {  {
         int i,j,k,inv,a,n,m;    int i,j,k,inv,a,n,m;
         unsigned int *t,*pivot;    unsigned int *t,*pivot;
   
         n = col; m = row+col;    n = col; m = row+col;
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 for ( i = j; i < row && !mat[i][j]; i++ );      for ( i = j; i < row && !mat[i][j]; i++ );
                 if ( i == row )      if ( i == row )
                         return 1;        return 1;
                 if ( i != j ) {      if ( i != j ) {
                         t = mat[i]; mat[i] = mat[j]; mat[j] = t;        t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                 }      }
                 pivot = mat[j];      pivot = mat[j];
                 inv = invm(pivot[j],md);      inv = invm(pivot[j],md);
                 for ( k = j; k < m; k++ )      for ( k = j; k < m; k++ )
                         pivot[k] = dmar(pivot[k],inv,0,md);        pivot[k] = dmar(pivot[k],inv,0,md);
                 for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
                                         t[k] = dmar(pivot[k],a,t[k],md);            t[k] = dmar(pivot[k],a,t[k],md);
                 }      }
         }    }
         for ( j = n-1; j >= 0; j-- ) {    for ( j = n-1; j >= 0; j-- ) {
                 pivot = mat[j];      pivot = mat[j];
                 for ( i = j-1; i >= 0; i-- ) {      for ( i = j-1; i >= 0; i-- ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
                                         t[k] = dmar(pivot[k],a,t[k],md);            t[k] = dmar(pivot[k],a,t[k],md);
                 }      }
         }    }
         return 0;    return 0;
 }  }
   
 void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)  void Psolve_by_lu_gfmmat(NODE arg,VECT *rp)
 {  {
         GFMMAT lu;    GFMMAT lu;
         Q *perm,*rhs,*v;    Q *perm,*rhs,*v;
         int n,i;    int n,i;
         unsigned int md;    unsigned int md;
         unsigned int *b,*sol;    unsigned int *b,*sol;
         VECT r;    VECT r;
   
         lu = (GFMMAT)ARG0(arg);    lu = (GFMMAT)ARG0(arg);
         perm = (Q *)BDY((VECT)ARG1(arg));    perm = (Q *)BDY((VECT)ARG1(arg));
         rhs = (Q *)BDY((VECT)ARG2(arg));    rhs = (Q *)BDY((VECT)ARG2(arg));
         md = (unsigned int)QTOS((Q)ARG3(arg));    md = (unsigned int)QTOS((Q)ARG3(arg));
         n = lu->col;    n = lu->col;
         b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));    b = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
         sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));    sol = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 b[i] = QTOS(rhs[QTOS(perm[i])]);      b[i] = QTOS(rhs[QTOS(perm[i])]);
         solve_by_lu_gfmmat(lu,md,b,sol);    solve_by_lu_gfmmat(lu,md,b,sol);
         MKVECT(r,n);    MKVECT(r,n);
         for ( i = 0, v = (Q *)r->body; i < n; i++ )    for ( i = 0, v = (Q *)r->body; i < n; i++ )
                         UTOQ(sol[i],v[i]);        UTOQ(sol[i],v[i]);
         *rp = r;    *rp = r;
 }  }
   
 void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,  void solve_by_lu_gfmmat(GFMMAT lu,unsigned int md,
         unsigned int *b,unsigned int *x)    unsigned int *b,unsigned int *x)
 {  {
         int n;    int n;
         unsigned int **a;    unsigned int **a;
         unsigned int *y;    unsigned int *y;
         int i,j;    int i,j;
         unsigned int t,m;    unsigned int t,m;
   
         n = lu->col;    n = lu->col;
         a = lu->body;    a = lu->body;
         y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));    y = (unsigned int *)MALLOC_ATOMIC(n*sizeof(int));
         /* solve Ly=b */    /* solve Ly=b */
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 for ( t = b[i], j = 0; j < i; j++ )      for ( t = b[i], j = 0; j < i; j++ )
                         if ( a[i][j] ) {        if ( a[i][j] ) {
                                 m = md - a[i][j];          m = md - a[i][j];
                                 DMAR(m,y[j],t,md,t)          DMAR(m,y[j],t,md,t)
                         }        }
                 y[i] = t;      y[i] = t;
         }    }
         /* solve Ux=y */    /* solve Ux=y */
         for ( i = n-1; i >= 0; i-- ) {    for ( i = n-1; i >= 0; i-- ) {
                 for ( t = y[i], j =i+1; j < n; j++ )      for ( t = y[i], j =i+1; j < n; j++ )
                         if ( a[i][j] ) {        if ( a[i][j] ) {
                                 m = md - a[i][j];          m = md - a[i][j];
                                 DMAR(m,x[j],t,md,t)          DMAR(m,x[j],t,md,t)
                         }        }
                 /* a[i][i] = 1/U[i][i] */      /* a[i][i] = 1/U[i][i] */
                 DMAR(t,a[i][i],0,md,x[i])      DMAR(t,a[i][i],0,md,x[i])
         }    }
 }  }
   
 void Plu_mat(NODE arg,LIST *rp)  void Plu_mat(NODE arg,LIST *rp)
 {  {
         MAT m,lu;    MAT m,lu;
         Q dn;    Q dn;
         Q *v;    Q *v;
         int n,i;    int n,i;
         int *iperm;    int *iperm;
         VECT perm;    VECT perm;
         NODE n0;    NODE n0;
   
         asir_assert(ARG0(arg),O_MAT,"lu_mat");    asir_assert(ARG0(arg),O_MAT,"lu_mat");
         m = (MAT)ARG0(arg);    m = (MAT)ARG0(arg);
         n = m->row;    n = m->row;
         MKMAT(lu,n,n);    MKMAT(lu,n,n);
         lu_dec_cr(m,lu,&dn,&iperm);    lu_dec_cr(m,lu,&dn,&iperm);
         MKVECT(perm,n);    MKVECT(perm,n);
         for ( i = 0, v = (Q *)perm->body; i < n; i++ )    for ( i = 0, v = (Q *)perm->body; i < n; i++ )
                 STOQ(iperm[i],v[i]);      STOQ(iperm[i],v[i]);
         n0 = mknode(3,lu,dn,perm);    n0 = mknode(3,lu,dn,perm);
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void Plu_gfmmat(NODE arg,LIST *rp)  void Plu_gfmmat(NODE arg,LIST *rp)
 {  {
         MAT m;    MAT m;
         GFMMAT mm;    GFMMAT mm;
         unsigned int md;    unsigned int md;
         int i,row,col,status;    int i,row,col,status;
         int *iperm;    int *iperm;
         Q *v;    Q *v;
         VECT perm;    VECT perm;
         NODE n0;    NODE n0;
   
         asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");    asir_assert(ARG0(arg),O_MAT,"lu_gfmmat");
         asir_assert(ARG1(arg),O_N,"lu_gfmmat");    asir_assert(ARG1(arg),O_N,"lu_gfmmat");
         m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
         mat_to_gfmmat(m,md,&mm);    mat_to_gfmmat(m,md,&mm);
         row = m->row;    row = m->row;
         col = m->col;    col = m->col;
         iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));    iperm = (int *)MALLOC_ATOMIC(row*sizeof(int));
         status = lu_gfmmat(mm,md,iperm);    status = lu_gfmmat(mm,md,iperm);
         if ( !status )    if ( !status )
                 n0 = 0;      n0 = 0;
         else {    else {
                 MKVECT(perm,row);      MKVECT(perm,row);
                 for ( i = 0, v = (Q *)perm->body; i < row; i++ )      for ( i = 0, v = (Q *)perm->body; i < row; i++ )
                         STOQ(iperm[i],v[i]);        STOQ(iperm[i],v[i]);
                 n0 = mknode(2,mm,perm);      n0 = mknode(2,mm,perm);
         }    }
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)  void Pmat_to_gfmmat(NODE arg,GFMMAT *rp)
 {  {
         MAT m;    MAT m;
         unsigned int md;    unsigned int md;
   
         asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");    asir_assert(ARG0(arg),O_MAT,"mat_to_gfmmat");
         asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");    asir_assert(ARG1(arg),O_N,"mat_to_gfmmat");
         m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = (unsigned int)QTOS((Q)ARG1(arg));
         mat_to_gfmmat(m,md,rp);    mat_to_gfmmat(m,md,rp);
 }  }
   
 void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)  void mat_to_gfmmat(MAT m,unsigned int md,GFMMAT *rp)
 {  {
         unsigned int **wmat;    unsigned int **wmat;
         unsigned int t;    unsigned int t;
         Q **mat;    Q **mat;
         Q q;    Q q;
         int i,j,row,col;    int i,j,row,col;
   
         row = m->row; col = m->col; mat = (Q **)m->body;    row = m->row; col = m->col; mat = (Q **)m->body;
         wmat = (unsigned int **)almat(row,col);    wmat = (unsigned int **)almat(row,col);
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 bzero((char *)wmat[i],col*sizeof(unsigned int));      bzero((char *)wmat[i],col*sizeof(unsigned int));
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = mat[i][j] ) {        if ( q = mat[i][j] ) {
                                 t = (unsigned int)rem(NM(q),md);          t = (unsigned int)rem(NM(q),md);
                                 if ( SGN(q) < 0 )          if ( SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         }        }
         }    }
         TOGFMMAT(row,col,wmat,*rp);    TOGFMMAT(row,col,wmat,*rp);
 }  }
   
 void Pgeninvm_swap(NODE arg,LIST *rp)  void Pgeninvm_swap(NODE arg,LIST *rp)
 {  {
         MAT m;    MAT m;
         pointer **mat;    pointer **mat;
         Q **tmat;    Q **tmat;
         Q *tvect;    Q *tvect;
         Q q;    Q q;
         unsigned int **wmat,**invmat;    unsigned int **wmat,**invmat;
         int *index;    int *index;
         unsigned int t,md;    unsigned int t,md;
         int i,j,row,col,status;    int i,j,row,col,status;
         MAT mat1;    MAT mat1;
         VECT vect1;    VECT vect1;
         NODE node1,node2;    NODE node1,node2;
   
         asir_assert(ARG0(arg),O_MAT,"geninvm_swap");    asir_assert(ARG0(arg),O_MAT,"geninvm_swap");
         asir_assert(ARG1(arg),O_N,"geninvm_swap");    asir_assert(ARG1(arg),O_N,"geninvm_swap");
         m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));    m = (MAT)ARG0(arg); md = QTOS((Q)ARG1(arg));
         row = m->row; col = m->col; mat = m->body;    row = m->row; col = m->col; mat = m->body;
         wmat = (unsigned int **)almat(row,col+row);    wmat = (unsigned int **)almat(row,col+row);
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 bzero((char *)wmat[i],(col+row)*sizeof(int));      bzero((char *)wmat[i],(col+row)*sizeof(int));
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (Q)mat[i][j] ) {        if ( q = (Q)mat[i][j] ) {
                                 t = (unsigned int)rem(NM(q),md);          t = (unsigned int)rem(NM(q),md);
                                 if ( SGN(q) < 0 )          if ( SGN(q) < 0 )
                                         t = (md - t) % md;            t = (md - t) % md;
                                 wmat[i][j] = t;          wmat[i][j] = t;
                         }        }
                 wmat[i][col+i] = 1;      wmat[i][col+i] = 1;
         }    }
         status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);    status = gauss_elim_geninv_mod_swap(wmat,row,col,md,&invmat,&index);
         if ( status > 0 )    if ( status > 0 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 MKMAT(mat1,col,col);      MKMAT(mat1,col,col);
                 for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )      for ( i = 0, tmat = (Q **)mat1->body; i < col; i++ )
                         for ( j = 0; j < col; j++ )        for ( j = 0; j < col; j++ )
                                 UTOQ(invmat[i][j],tmat[i][j]);          UTOQ(invmat[i][j],tmat[i][j]);
                 MKVECT(vect1,row);      MKVECT(vect1,row);
                 for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )      for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
                         STOQ(index[i],tvect[i]);        STOQ(index[i],tvect[i]);
                 MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);       MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
         }    }
 }  }
   
 int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md,  int gauss_elim_geninv_mod_swap(unsigned int **mat,int row,int col,unsigned int md,
     unsigned int ***invmatp,int **indexp)      unsigned int ***invmatp,int **indexp)
 {  {
         int i,j,k,inv,a,n,m;    int i,j,k,inv,a,n,m;
         unsigned int *t,*pivot,*s;    unsigned int *t,*pivot,*s;
         int *index;    int *index;
         unsigned int **invmat;    unsigned int **invmat;
   
         n = col; m = row+col;    n = col; m = row+col;
         *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));    *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 index[i] = i;      index[i] = i;
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 for ( i = j; i < row && !mat[i][j]; i++ );      for ( i = j; i < row && !mat[i][j]; i++ );
                 if ( i == row ) {      if ( i == row ) {
                         *indexp = 0; *invmatp = 0; return 1;        *indexp = 0; *invmatp = 0; return 1;
                 }      }
                 if ( i != j ) {      if ( i != j ) {
                         t = mat[i]; mat[i] = mat[j]; mat[j] = t;        t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                         k = index[i]; index[i] = index[j]; index[j] = k;        k = index[i]; index[i] = index[j]; index[j] = k;
                 }      }
                 pivot = mat[j];      pivot = mat[j];
                 inv = (unsigned int)invm(pivot[j],md);      inv = (unsigned int)invm(pivot[j],md);
                 for ( k = j; k < m; k++ )      for ( k = j; k < m; k++ )
                         if ( pivot[k] )        if ( pivot[k] )
                                 pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);          pivot[k] = (unsigned int)dmar(pivot[k],inv,0,md);
                 for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
                                         if ( pivot[k] )            if ( pivot[k] )
                                                 t[k] = dmar(pivot[k],a,t[k],md);              t[k] = dmar(pivot[k],a,t[k],md);
                 }      }
         }    }
         for ( j = n-1; j >= 0; j-- ) {    for ( j = n-1; j >= 0; j-- ) {
                 pivot = mat[j];      pivot = mat[j];
                 for ( i = j-1; i >= 0; i-- ) {      for ( i = j-1; i >= 0; i-- ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 for ( k = j, a = md - a; k < m; k++ )          for ( k = j, a = md - a; k < m; k++ )
                                         if ( pivot[k] )            if ( pivot[k] )
                                                 t[k] = dmar(pivot[k],a,t[k],md);              t[k] = dmar(pivot[k],a,t[k],md);
                 }      }
         }    }
         *invmatp = invmat = (unsigned int **)almat(col,col);    *invmatp = invmat = (unsigned int **)almat(col,col);
         for ( i = 0; i < col; i++ )    for ( i = 0; i < col; i++ )
                 for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )      for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
                         s[j] = t[col+index[j]];        s[j] = t[col+index[j]];
         return 0;    return 0;
 }  }
   
 void Pgeninv_sf_swap(NODE arg,LIST *rp)  void Pgeninv_sf_swap(NODE arg,LIST *rp)
 {  {
         MAT m;    MAT m;
         GFS **mat,**tmat;    GFS **mat,**tmat;
         Q *tvect;    Q *tvect;
         GFS q;    GFS q;
         int **wmat,**invmat;    int **wmat,**invmat;
         int *index;    int *index;
         unsigned int t;    unsigned int t;
         int i,j,row,col,status;    int i,j,row,col,status;
         MAT mat1;    MAT mat1;
         VECT vect1;    VECT vect1;
         NODE node1,node2;    NODE node1,node2;
   
         asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");    asir_assert(ARG0(arg),O_MAT,"geninv_sf_swap");
         m = (MAT)ARG0(arg);    m = (MAT)ARG0(arg);
         row = m->row; col = m->col; mat = (GFS **)m->body;    row = m->row; col = m->col; mat = (GFS **)m->body;
         wmat = (int **)almat(row,col+row);    wmat = (int **)almat(row,col+row);
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 bzero((char *)wmat[i],(col+row)*sizeof(int));      bzero((char *)wmat[i],(col+row)*sizeof(int));
                 for ( j = 0; j < col; j++ )      for ( j = 0; j < col; j++ )
                         if ( q = (GFS)mat[i][j] )        if ( q = (GFS)mat[i][j] )
                                 wmat[i][j] = FTOIF(CONT(q));          wmat[i][j] = FTOIF(CONT(q));
                 wmat[i][col+i] = _onesf();      wmat[i][col+i] = _onesf();
         }    }
         status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);    status = gauss_elim_geninv_sf_swap(wmat,row,col,&invmat,&index);
         if ( status > 0 )    if ( status > 0 )
                 *rp = 0;      *rp = 0;
         else {    else {
                 MKMAT(mat1,col,col);      MKMAT(mat1,col,col);
                 for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )      for ( i = 0, tmat = (GFS **)mat1->body; i < col; i++ )
                         for ( j = 0; j < col; j++ )        for ( j = 0; j < col; j++ )
                                 if ( t = invmat[i][j] ) {          if ( t = invmat[i][j] ) {
                                         MKGFS(IFTOF(t),tmat[i][j]);            MKGFS(IFTOF(t),tmat[i][j]);
                                 }          }
                 MKVECT(vect1,row);      MKVECT(vect1,row);
                 for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )      for ( i = 0, tvect = (Q *)vect1->body; i < row; i++ )
                         STOQ(index[i],tvect[i]);        STOQ(index[i],tvect[i]);
                 MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);       MKNODE(node2,vect1,0); MKNODE(node1,mat1,node2); MKLIST(*rp,node1);
         }    }
 }  }
   
 int gauss_elim_geninv_sf_swap(int **mat,int row,int col,  int gauss_elim_geninv_sf_swap(int **mat,int row,int col,
         int ***invmatp,int **indexp)    int ***invmatp,int **indexp)
 {  {
         int i,j,k,inv,a,n,m,u;    int i,j,k,inv,a,n,m,u;
         int *t,*pivot,*s;    int *t,*pivot,*s;
         int *index;    int *index;
         int **invmat;    int **invmat;
   
         n = col; m = row+col;    n = col; m = row+col;
         *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));    *indexp = index = (int *)MALLOC_ATOMIC(row*sizeof(int));
         for ( i = 0; i < row; i++ )    for ( i = 0; i < row; i++ )
                 index[i] = i;      index[i] = i;
         for ( j = 0; j < n; j++ ) {    for ( j = 0; j < n; j++ ) {
                 for ( i = j; i < row && !mat[i][j]; i++ );      for ( i = j; i < row && !mat[i][j]; i++ );
                 if ( i == row ) {      if ( i == row ) {
                         *indexp = 0; *invmatp = 0; return 1;        *indexp = 0; *invmatp = 0; return 1;
                 }      }
                 if ( i != j ) {      if ( i != j ) {
                         t = mat[i]; mat[i] = mat[j]; mat[j] = t;        t = mat[i]; mat[i] = mat[j]; mat[j] = t;
                         k = index[i]; index[i] = index[j]; index[j] = k;        k = index[i]; index[i] = index[j]; index[j] = k;
                 }      }
                 pivot = mat[j];      pivot = mat[j];
                 inv = _invsf(pivot[j]);      inv = _invsf(pivot[j]);
                 for ( k = j; k < m; k++ )      for ( k = j; k < m; k++ )
                         if ( pivot[k] )        if ( pivot[k] )
                                 pivot[k] = _mulsf(pivot[k],inv);          pivot[k] = _mulsf(pivot[k],inv);
                 for ( i = j+1; i < row; i++ ) {      for ( i = j+1; i < row; i++ ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 for ( k = j, a = _chsgnsf(a); k < m; k++ )          for ( k = j, a = _chsgnsf(a); k < m; k++ )
                                         if ( pivot[k] ) {            if ( pivot[k] ) {
                                                 u = _mulsf(pivot[k],a);              u = _mulsf(pivot[k],a);
                                                 t[k] = _addsf(u,t[k]);              t[k] = _addsf(u,t[k]);
                                         }            }
                 }      }
         }    }
         for ( j = n-1; j >= 0; j-- ) {    for ( j = n-1; j >= 0; j-- ) {
                 pivot = mat[j];      pivot = mat[j];
                 for ( i = j-1; i >= 0; i-- ) {      for ( i = j-1; i >= 0; i-- ) {
                         t = mat[i];        t = mat[i];
                         if ( a = t[j] )        if ( a = t[j] )
                                 for ( k = j, a = _chsgnsf(a); k < m; k++ )          for ( k = j, a = _chsgnsf(a); k < m; k++ )
                                         if ( pivot[k] ) {            if ( pivot[k] ) {
                                                 u = _mulsf(pivot[k],a);              u = _mulsf(pivot[k],a);
                                                 t[k] = _addsf(u,t[k]);              t[k] = _addsf(u,t[k]);
                                         }            }
                 }      }
         }    }
         *invmatp = invmat = (int **)almat(col,col);    *invmatp = invmat = (int **)almat(col,col);
         for ( i = 0; i < col; i++ )    for ( i = 0; i < col; i++ )
                 for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )      for ( j = 0, s = invmat[i], t = mat[i]; j < col; j++ )
                         s[j] = t[col+index[j]];        s[j] = t[col+index[j]];
         return 0;    return 0;
 }  }
   
 void _addn(N,N,N);  void _addn(N,N,N);
Line 3309  void _muln(N,N,N);
Line 3309  void _muln(N,N,N);
   
 void inner_product_int(Q *a,Q *b,int n,Q *r)  void inner_product_int(Q *a,Q *b,int n,Q *r)
 {  {
         int la,lb,i;    int la,lb,i;
         int sgn,sgn1;    int sgn,sgn1;
         N wm,wma,sum,t;    N wm,wma,sum,t;
   
         for ( la = lb = 0, i = 0; i < n; i++ ) {    for ( la = lb = 0, i = 0; i < n; i++ ) {
                 if ( a[i] )      if ( a[i] )
                         if ( DN(a[i]) )        if ( DN(a[i]) )
                                 error("inner_product_int : invalid argument");          error("inner_product_int : invalid argument");
                         else        else
                                 la = MAX(PL(NM(a[i])),la);          la = MAX(PL(NM(a[i])),la);
                 if ( b[i] )      if ( b[i] )
                         if ( DN(b[i]) )        if ( DN(b[i]) )
                                 error("inner_product_int : invalid argument");          error("inner_product_int : invalid argument");
                         else        else
                                 lb = MAX(PL(NM(b[i])),lb);          lb = MAX(PL(NM(b[i])),lb);
         }    }
         sgn = 0;    sgn = 0;
         sum= NALLOC(la+lb+2);    sum= NALLOC(la+lb+2);
         bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));    bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
         wm = NALLOC(la+lb+2);    wm = NALLOC(la+lb+2);
         wma = NALLOC(la+lb+2);    wma = NALLOC(la+lb+2);
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 if ( !a[i] || !b[i] )      if ( !a[i] || !b[i] )
                         continue;        continue;
                 _muln(NM(a[i]),NM(b[i]),wm);      _muln(NM(a[i]),NM(b[i]),wm);
                 sgn1 = SGN(a[i])*SGN(b[i]);      sgn1 = SGN(a[i])*SGN(b[i]);
                 if ( !sgn ) {      if ( !sgn ) {
                         sgn = sgn1;        sgn = sgn1;
                         t = wm; wm = sum; sum = t;        t = wm; wm = sum; sum = t;
                 } else if ( sgn == sgn1 ) {      } else if ( sgn == sgn1 ) {
                         _addn(sum,wm,wma);        _addn(sum,wm,wma);
                         if ( !PL(wma) )        if ( !PL(wma) )
                                 sgn = 0;          sgn = 0;
                         t = wma; wma = sum; sum = t;        t = wma; wma = sum; sum = t;
                 } else {      } else {
                         /* sgn*sum+sgn1*wm = sgn*(sum-wm) */        /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
                         sgn *= _subn(sum,wm,wma);        sgn *= _subn(sum,wm,wma);
                         t = wma; wma = sum; sum = t;        t = wma; wma = sum; sum = t;
                 }      }
         }    }
         GCFREE(wm);    GCFREE(wm);
         GCFREE(wma);    GCFREE(wma);
         if ( !sgn ) {    if ( !sgn ) {
                 GCFREE(sum);      GCFREE(sum);
                 *r = 0;      *r = 0;
         } else    } else
                 NTOQ(sum,sgn,*r);      NTOQ(sum,sgn,*r);
 }  }
   
 /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */  /* (k,l) element of a*b where a: .x n matrix, b: n x . integer matrix */
   
 void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)  void inner_product_mat_int_mod(Q **a,int **b,int n,int k,int l,Q *r)
 {  {
         int la,lb,i;    int la,lb,i;
         int sgn,sgn1;    int sgn,sgn1;
         N wm,wma,sum,t;    N wm,wma,sum,t;
         Q aki;    Q aki;
         int bil,bilsgn;    int bil,bilsgn;
         struct oN tn;    struct oN tn;
   
         for ( la = 0, i = 0; i < n; i++ ) {    for ( la = 0, i = 0; i < n; i++ ) {
                 if ( aki = a[k][i] )      if ( aki = a[k][i] )
                         if ( DN(aki) )        if ( DN(aki) )
                                 error("inner_product_int : invalid argument");          error("inner_product_int : invalid argument");
                         else        else
                                 la = MAX(PL(NM(aki)),la);          la = MAX(PL(NM(aki)),la);
         }    }
         lb = 1;    lb = 1;
         sgn = 0;    sgn = 0;
         sum= NALLOC(la+lb+2);    sum= NALLOC(la+lb+2);
         bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));    bzero((char *)sum,(la+lb+3)*sizeof(unsigned int));
         wm = NALLOC(la+lb+2);    wm = NALLOC(la+lb+2);
         wma = NALLOC(la+lb+2);    wma = NALLOC(la+lb+2);
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 if ( !(aki = a[k][i]) || !(bil = b[i][l]) )      if ( !(aki = a[k][i]) || !(bil = b[i][l]) )
                         continue;        continue;
                 tn.p = 1;      tn.p = 1;
                 if ( bil > 0 ) {      if ( bil > 0 ) {
                         tn.b[0] = bil; bilsgn = 1;        tn.b[0] = bil; bilsgn = 1;
                 } else {      } else {
                         tn.b[0] = -bil; bilsgn = -1;        tn.b[0] = -bil; bilsgn = -1;
                 }      }
                 _muln(NM(aki),&tn,wm);      _muln(NM(aki),&tn,wm);
                 sgn1 = SGN(aki)*bilsgn;      sgn1 = SGN(aki)*bilsgn;
                 if ( !sgn ) {      if ( !sgn ) {
                         sgn = sgn1;        sgn = sgn1;
                         t = wm; wm = sum; sum = t;        t = wm; wm = sum; sum = t;
                 } else if ( sgn == sgn1 ) {      } else if ( sgn == sgn1 ) {
                         _addn(sum,wm,wma);        _addn(sum,wm,wma);
                         if ( !PL(wma) )        if ( !PL(wma) )
                                 sgn = 0;          sgn = 0;
                         t = wma; wma = sum; sum = t;        t = wma; wma = sum; sum = t;
                 } else {      } else {
                         /* sgn*sum+sgn1*wm = sgn*(sum-wm) */        /* sgn*sum+sgn1*wm = sgn*(sum-wm) */
                         sgn *= _subn(sum,wm,wma);        sgn *= _subn(sum,wm,wma);
                         t = wma; wma = sum; sum = t;        t = wma; wma = sum; sum = t;
                 }      }
         }    }
         GCFREE(wm);    GCFREE(wm);
         GCFREE(wma);    GCFREE(wma);
         if ( !sgn ) {    if ( !sgn ) {
                 GCFREE(sum);      GCFREE(sum);
                 *r = 0;      *r = 0;
         } else    } else
                 NTOQ(sum,sgn,*r);      NTOQ(sum,sgn,*r);
 }  }
   
 void Pmul_mat_vect_int(NODE arg,VECT *rp)  void Pmul_mat_vect_int(NODE arg,VECT *rp)
 {  {
         MAT mat;    MAT mat;
         VECT vect,r;    VECT vect,r;
         int row,col,i;    int row,col,i;
   
         mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
         vect = (VECT)ARG1(arg);    vect = (VECT)ARG1(arg);
         row = mat->row;    row = mat->row;
         col = mat->col;    col = mat->col;
         MKVECT(r,row);    MKVECT(r,row);
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);      inner_product_int((Q *)mat->body[i],(Q *)vect->body,col,(Q *)&r->body[i]);
         }    }
         *rp = r;    *rp = r;
 }  }
   
 void Pnbpoly_up2(NODE arg,GF2N *rp)  void Pnbpoly_up2(NODE arg,GF2N *rp)
 {  {
         int m,type,ret;    int m,type,ret;
         UP2 r;    UP2 r;
   
         m = QTOS((Q)ARG0(arg));    m = QTOS((Q)ARG0(arg));
         type = QTOS((Q)ARG1(arg));    type = QTOS((Q)ARG1(arg));
         ret = generate_ONB_polynomial(&r,m,type);    ret = generate_ONB_polynomial(&r,m,type);
         if ( ret == 0 )    if ( ret == 0 )
                 MKGF2N(r,*rp);      MKGF2N(r,*rp);
         else    else
                 *rp = 0;      *rp = 0;
 }  }
   
 void Px962_irredpoly_up2(NODE arg,GF2N *rp)  void Px962_irredpoly_up2(NODE arg,GF2N *rp)
 {  {
         int m,ret,w;    int m,ret,w;
         GF2N prev;    GF2N prev;
         UP2 r;    UP2 r;
   
         m = QTOS((Q)ARG0(arg));    m = QTOS((Q)ARG0(arg));
         prev = (GF2N)ARG1(arg);    prev = (GF2N)ARG1(arg);
         if ( !prev ) {    if ( !prev ) {
                 w = (m>>5)+1; NEWUP2(r,w); r->w = 0;      w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                 bzero((char *)r->b,w*sizeof(unsigned int));      bzero((char *)r->b,w*sizeof(unsigned int));
         } else {    } else {
                 r = prev->body;      r = prev->body;
                 if ( degup2(r) != m ) {      if ( degup2(r) != m ) {
                         w = (m>>5)+1; NEWUP2(r,w); r->w = 0;        w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                         bzero((char *)r->b,w*sizeof(unsigned int));        bzero((char *)r->b,w*sizeof(unsigned int));
                 }      }
         }    }
         ret = _generate_irreducible_polynomial(r,m);    ret = _generate_irreducible_polynomial(r,m);
         if ( ret == 0 )    if ( ret == 0 )
                 MKGF2N(r,*rp);      MKGF2N(r,*rp);
         else    else
                 *rp = 0;      *rp = 0;
 }  }
   
 void Pirredpoly_up2(NODE arg,GF2N *rp)  void Pirredpoly_up2(NODE arg,GF2N *rp)
 {  {
         int m,ret,w;    int m,ret,w;
         GF2N prev;    GF2N prev;
         UP2 r;    UP2 r;
   
         m = QTOS((Q)ARG0(arg));    m = QTOS((Q)ARG0(arg));
         prev = (GF2N)ARG1(arg);    prev = (GF2N)ARG1(arg);
         if ( !prev ) {    if ( !prev ) {
                 w = (m>>5)+1; NEWUP2(r,w); r->w = 0;      w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                 bzero((char *)r->b,w*sizeof(unsigned int));      bzero((char *)r->b,w*sizeof(unsigned int));
         } else {    } else {
                 r = prev->body;      r = prev->body;
                 if ( degup2(r) != m ) {      if ( degup2(r) != m ) {
                         w = (m>>5)+1; NEWUP2(r,w); r->w = 0;        w = (m>>5)+1; NEWUP2(r,w); r->w = 0;
                         bzero((char *)r->b,w*sizeof(unsigned int));        bzero((char *)r->b,w*sizeof(unsigned int));
                 }      }
         }    }
         ret = _generate_good_irreducible_polynomial(r,m);    ret = _generate_good_irreducible_polynomial(r,m);
         if ( ret == 0 )    if ( ret == 0 )
                 MKGF2N(r,*rp);      MKGF2N(r,*rp);
         else    else
                 *rp = 0;      *rp = 0;
 }  }
   
 void Pmat_swap_row_destructive(NODE arg, MAT *m)  void Pmat_swap_row_destructive(NODE arg, MAT *m)
 {  {
         int i1,i2;    int i1,i2;
         pointer *t;    pointer *t;
         MAT mat;    MAT mat;
   
         asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");    asir_assert(ARG0(arg),O_MAT,"mat_swap_row_destructive");
         asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");    asir_assert(ARG1(arg),O_N,"mat_swap_row_destructive");
         asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");    asir_assert(ARG2(arg),O_N,"mat_swap_row_destructive");
         mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
         i1 = QTOS((Q)ARG1(arg));    i1 = QTOS((Q)ARG1(arg));
         i2 = QTOS((Q)ARG2(arg));    i2 = QTOS((Q)ARG2(arg));
         if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )    if ( i1 < 0 || i2 < 0 || i1 >= mat->row || i2 >= mat->row )
                 error("mat_swap_row_destructive : Out of range");      error("mat_swap_row_destructive : Out of range");
         t = mat->body[i1];    t = mat->body[i1];
         mat->body[i1] = mat->body[i2];    mat->body[i1] = mat->body[i2];
         mat->body[i2] = t;    mat->body[i2] = t;
         *m = mat;    *m = mat;
 }  }
   
 void Pmat_swap_col_destructive(NODE arg, MAT *m)  void Pmat_swap_col_destructive(NODE arg, MAT *m)
 {  {
         int j1,j2,i,n;    int j1,j2,i,n;
         pointer *mi;    pointer *mi;
         pointer t;    pointer t;
         MAT mat;    MAT mat;
   
         asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");    asir_assert(ARG0(arg),O_MAT,"mat_swap_col_destructive");
         asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");    asir_assert(ARG1(arg),O_N,"mat_swap_col_destructive");
         asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");    asir_assert(ARG2(arg),O_N,"mat_swap_col_destructive");
         mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
         j1 = QTOS((Q)ARG1(arg));    j1 = QTOS((Q)ARG1(arg));
         j2 = QTOS((Q)ARG2(arg));    j2 = QTOS((Q)ARG2(arg));
         if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )    if ( j1 < 0 || j2 < 0 || j1 >= mat->col || j2 >= mat->col )
                 error("mat_swap_col_destructive : Out of range");      error("mat_swap_col_destructive : Out of range");
         n = mat->row;    n = mat->row;
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 mi = mat->body[i];      mi = mat->body[i];
                 t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;      t = mi[j1]; mi[j1] = mi[j2]; mi[j2] = t;
         }    }
         *m = mat;    *m = mat;
 }  }
 /*  /*
  * f = type 'type' normal polynomial of degree m if exists   * f = type 'type' normal polynomial of degree m if exists
Line 3550  void Pmat_swap_col_destructive(NODE arg, MAT *m)
Line 3550  void Pmat_swap_col_destructive(NODE arg, MAT *m)
   
 int generate_ONB_polynomial(UP2 *rp,int m,int type)  int generate_ONB_polynomial(UP2 *rp,int m,int type)
 {  {
         int i,r;    int i,r;
         int w;    int w;
         UP2 f,f0,f1,f2,t;    UP2 f,f0,f1,f2,t;
   
         w = (m>>5)+1;    w = (m>>5)+1;
         switch ( type ) {    switch ( type ) {
                 case 1:      case 1:
                         if ( !TypeT_NB_check(m,1) ) return 1;        if ( !TypeT_NB_check(m,1) ) return 1;
                         NEWUP2(f,w); *rp = f; f->w = w;        NEWUP2(f,w); *rp = f; f->w = w;
                         /* set all the bits */        /* set all the bits */
                         for ( i = 0; i < w; i++ )        for ( i = 0; i < w; i++ )
                                 f->b[i] = 0xffffffff;          f->b[i] = 0xffffffff;
                         /* mask the top word if necessary */        /* mask the top word if necessary */
                         if ( r = (m+1)&31 )        if ( r = (m+1)&31 )
                                 f->b[w-1] &= (1<<r)-1;          f->b[w-1] &= (1<<r)-1;
                         return 0;        return 0;
                         break;        break;
                 case 2:      case 2:
                         if ( !TypeT_NB_check(m,2) ) return 1;        if ( !TypeT_NB_check(m,2) ) return 1;
                         NEWUP2(f,w); *rp = f;        NEWUP2(f,w); *rp = f;
                         W_NEWUP2(f0,w);        W_NEWUP2(f0,w);
                         W_NEWUP2(f1,w);        W_NEWUP2(f1,w);
                         W_NEWUP2(f2,w);        W_NEWUP2(f2,w);
   
                         /* recursion for genrating Type II normal polynomial */        /* recursion for genrating Type II normal polynomial */
   
                         /* f0 = 1, f1 = t+1 */        /* f0 = 1, f1 = t+1 */
                         f0->w = 1; f0->b[0] = 1;        f0->w = 1; f0->b[0] = 1;
                         f1->w = 1; f1->b[0] = 3;        f1->w = 1; f1->b[0] = 3;
                         for ( i = 2; i <= m; i++ ) {        for ( i = 2; i <= m; i++ ) {
                                 /* f2 = t*f1+f0 */          /* f2 = t*f1+f0 */
                                 _bshiftup2(f1,-1,f2);          _bshiftup2(f1,-1,f2);
                                 _addup2_destructive(f2,f0);          _addup2_destructive(f2,f0);
                                 /* cyclic change of the variables */          /* cyclic change of the variables */
                                 t = f0; f0 = f1; f1 = f2; f2 = t;          t = f0; f0 = f1; f1 = f2; f2 = t;
                         }        }
                         _copyup2(f1,f);        _copyup2(f1,f);
                         return 0;        return 0;
                         break;        break;
                 default:      default:
                         return -1;        return -1;
                         break;        break;
                 }      }
 }  }
   
 /*  /*
Line 3603  int generate_ONB_polynomial(UP2 *rp,int m,int type)
Line 3603  int generate_ONB_polynomial(UP2 *rp,int m,int type)
   
 int _generate_irreducible_polynomial(UP2 f,int d)  int _generate_irreducible_polynomial(UP2 f,int d)
 {  {
         int ret,i,j,k,nz,i0,j0,k0;    int ret,i,j,k,nz,i0,j0,k0;
         int w;    int w;
         unsigned int *fd;    unsigned int *fd;
   
         /*    /*
          * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.     * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
          * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.     * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
          * otherwise i0,j0,k0 is set to 0.     * otherwise i0,j0,k0 is set to 0.
          */     */
   
         fd = f->b;    fd = f->b;
         w = (d>>5)+1;    w = (d>>5)+1;
         if ( f->w && (d==degup2(f)) ) {    if ( f->w && (d==degup2(f)) ) {
                 for ( nz = 0, i = d; i >= 0; i-- )      for ( nz = 0, i = d; i >= 0; i-- )
                         if ( fd[i>>5]&(1<<(i&31)) ) nz++;        if ( fd[i>>5]&(1<<(i&31)) ) nz++;
                 switch ( nz ) {      switch ( nz ) {
                         case 3:        case 3:
                                 for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );          for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                                 /* reset i0-th bit */          /* reset i0-th bit */
                                 fd[i0>>5] &= ~(1<<(i0&31));          fd[i0>>5] &= ~(1<<(i0&31));
                                 j0 = k0 = 0;          j0 = k0 = 0;
                                 break;          break;
                         case 5:        case 5:
                                 for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );          for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                                 /* reset i0-th bit */          /* reset i0-th bit */
                                 fd[i0>>5] &= ~(1<<(i0&31));          fd[i0>>5] &= ~(1<<(i0&31));
                                 for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );          for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
                                 /* reset j0-th bit */          /* reset j0-th bit */
                                 fd[j0>>5] &= ~(1<<(j0&31));          fd[j0>>5] &= ~(1<<(j0&31));
                                 for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );          for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
                                 /* reset k0-th bit */          /* reset k0-th bit */
                                 fd[k0>>5] &= ~(1<<(k0&31));          fd[k0>>5] &= ~(1<<(k0&31));
                                 break;          break;
                         default:        default:
                                 f->w = 0; break;          f->w = 0; break;
                 }      }
         } else    } else
                 f->w = 0;      f->w = 0;
   
         if ( !f->w ) {    if ( !f->w ) {
                 fd = f->b;      fd = f->b;
                 f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));      f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
                 i0 = j0 = k0 = 0;      i0 = j0 = k0 = 0;
         }    }
         /* if j0 > 0 then f is already a pentanomial */    /* if j0 > 0 then f is already a pentanomial */
         if ( j0 > 0 ) goto PENTA;    if ( j0 > 0 ) goto PENTA;
   
         /* searching for an irreducible trinomial */    /* searching for an irreducible trinomial */
   
         for ( i = 1; 2*i <= d; i++ ) {    for ( i = 1; 2*i <= d; i++ ) {
                 /* skip the polynomials 'before' f */      /* skip the polynomials 'before' f */
                 if ( i < i0 ) continue;      if ( i < i0 ) continue;
                 if ( i == i0 ) { i0 = 0; continue; }      if ( i == i0 ) { i0 = 0; continue; }
                 /* set i-th bit */      /* set i-th bit */
                 fd[i>>5] |= (1<<(i&31));      fd[i>>5] |= (1<<(i&31));
                 ret = irredcheck_dddup2(f);      ret = irredcheck_dddup2(f);
                 if ( ret == 1 ) return 0;      if ( ret == 1 ) return 0;
                 /* reset i-th bit */      /* reset i-th bit */
                 fd[i>>5] &= ~(1<<(i&31));      fd[i>>5] &= ~(1<<(i&31));
         }    }
   
         /* searching for an irreducible pentanomial */    /* searching for an irreducible pentanomial */
 PENTA:  PENTA:
         for ( i = 1; i < d; i++ ) {    for ( i = 1; i < d; i++ ) {
                 /* skip the polynomials 'before' f */      /* skip the polynomials 'before' f */
                 if ( i < i0 ) continue;      if ( i < i0 ) continue;
                 if ( i == i0 ) i0 = 0;      if ( i == i0 ) i0 = 0;
                 /* set i-th bit */      /* set i-th bit */
                 fd[i>>5] |= (1<<(i&31));      fd[i>>5] |= (1<<(i&31));
                 for ( j = i+1; j < d; j++ ) {      for ( j = i+1; j < d; j++ ) {
                         /* skip the polynomials 'before' f */        /* skip the polynomials 'before' f */
                         if ( j < j0 ) continue;        if ( j < j0 ) continue;
                         if ( j == j0 ) j0 = 0;        if ( j == j0 ) j0 = 0;
                         /* set j-th bit */        /* set j-th bit */
                         fd[j>>5] |= (1<<(j&31));        fd[j>>5] |= (1<<(j&31));
                         for ( k = j+1; k < d; k++ ) {        for ( k = j+1; k < d; k++ ) {
                                 /* skip the polynomials 'before' f */          /* skip the polynomials 'before' f */
                                 if ( k < k0 ) continue;          if ( k < k0 ) continue;
                                 else if ( k == k0 ) { k0 = 0; continue; }          else if ( k == k0 ) { k0 = 0; continue; }
                                 /* set k-th bit */          /* set k-th bit */
                                 fd[k>>5] |= (1<<(k&31));          fd[k>>5] |= (1<<(k&31));
                                 ret = irredcheck_dddup2(f);          ret = irredcheck_dddup2(f);
                                 if ( ret == 1 ) return 0;          if ( ret == 1 ) return 0;
                                 /* reset k-th bit */          /* reset k-th bit */
                                 fd[k>>5] &= ~(1<<(k&31));          fd[k>>5] &= ~(1<<(k&31));
                         }        }
                         /* reset j-th bit */        /* reset j-th bit */
                         fd[j>>5] &= ~(1<<(j&31));        fd[j>>5] &= ~(1<<(j&31));
                 }      }
                 /* reset i-th bit */      /* reset i-th bit */
                 fd[i>>5] &= ~(1<<(i&31));      fd[i>>5] &= ~(1<<(i&31));
         }    }
         /* exhausted */    /* exhausted */
         return 1;    return 1;
 }  }
   
 /*  /*
Line 3716  PENTA:
Line 3716  PENTA:
   
 int _generate_good_irreducible_polynomial(UP2 f,int d)  int _generate_good_irreducible_polynomial(UP2 f,int d)
 {  {
         int ret,i,j,k,nz,i0,j0,k0;    int ret,i,j,k,nz,i0,j0,k0;
         int w;    int w;
         unsigned int *fd;    unsigned int *fd;
   
         /*    /*
          * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.     * if f = x^d+x^i+1 then i0 <- i, j0 <- 0, k0 <-0.
          * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.     * if f = x^d+x^k+x^j+x^i+1 (k>j>i) then i0 <- i, j0 <- j, k0 <-k.
          * otherwise i0,j0,k0 is set to 0.     * otherwise i0,j0,k0 is set to 0.
          */     */
   
         fd = f->b;    fd = f->b;
         w = (d>>5)+1;    w = (d>>5)+1;
         if ( f->w && (d==degup2(f)) ) {    if ( f->w && (d==degup2(f)) ) {
                 for ( nz = 0, i = d; i >= 0; i-- )      for ( nz = 0, i = d; i >= 0; i-- )
                         if ( fd[i>>5]&(1<<(i&31)) ) nz++;        if ( fd[i>>5]&(1<<(i&31)) ) nz++;
                 switch ( nz ) {      switch ( nz ) {
                         case 3:        case 3:
                                 for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );          for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                                 /* reset i0-th bit */          /* reset i0-th bit */
                                 fd[i0>>5] &= ~(1<<(i0&31));          fd[i0>>5] &= ~(1<<(i0&31));
                                 j0 = k0 = 0;          j0 = k0 = 0;
                                 break;          break;
                         case 5:        case 5:
                                 for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );          for ( i0 = 1; !(fd[i0>>5]&(1<<(i0&31))) ; i0++ );
                                 /* reset i0-th bit */          /* reset i0-th bit */
                                 fd[i0>>5] &= ~(1<<(i0&31));          fd[i0>>5] &= ~(1<<(i0&31));
                                 for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );          for ( j0 = i0+1; !(fd[j0>>5]&(1<<(j0&31))) ; j0++ );
                                 /* reset j0-th bit */          /* reset j0-th bit */
                                 fd[j0>>5] &= ~(1<<(j0&31));          fd[j0>>5] &= ~(1<<(j0&31));
                                 for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );          for ( k0 = j0+1; !(fd[k0>>5]&(1<<(k0&31))) ; k0++ );
                                 /* reset k0-th bit */          /* reset k0-th bit */
                                 fd[k0>>5] &= ~(1<<(k0&31));          fd[k0>>5] &= ~(1<<(k0&31));
                                 break;          break;
                         default:        default:
                                 f->w = 0; break;          f->w = 0; break;
                 }      }
         } else    } else
                 f->w = 0;      f->w = 0;
   
         if ( !f->w ) {    if ( !f->w ) {
                 fd = f->b;      fd = f->b;
                 f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));      f->w = w; fd[0] |= 1; fd[d>>5] |= (1<<(d&31));
                 i0 = j0 = k0 = 0;      i0 = j0 = k0 = 0;
         }    }
         /* if j0 > 0 then f is already a pentanomial */    /* if j0 > 0 then f is already a pentanomial */
         if ( j0 > 0 ) goto PENTA;    if ( j0 > 0 ) goto PENTA;
   
         /* searching for an irreducible trinomial */    /* searching for an irreducible trinomial */
   
         for ( i = 1; 2*i <= d; i++ ) {    for ( i = 1; 2*i <= d; i++ ) {
                 /* skip the polynomials 'before' f */      /* skip the polynomials 'before' f */
                 if ( i < i0 ) continue;      if ( i < i0 ) continue;
                 if ( i == i0 ) { i0 = 0; continue; }      if ( i == i0 ) { i0 = 0; continue; }
                 /* set i-th bit */      /* set i-th bit */
                 fd[i>>5] |= (1<<(i&31));      fd[i>>5] |= (1<<(i&31));
                 ret = irredcheck_dddup2(f);      ret = irredcheck_dddup2(f);
                 if ( ret == 1 ) return 0;      if ( ret == 1 ) return 0;
                 /* reset i-th bit */      /* reset i-th bit */
                 fd[i>>5] &= ~(1<<(i&31));      fd[i>>5] &= ~(1<<(i&31));
         }    }
   
         /* searching for an irreducible pentanomial */    /* searching for an irreducible pentanomial */
 PENTA:  PENTA:
         for ( i = 3; i < d; i++ ) {    for ( i = 3; i < d; i++ ) {
                 /* skip the polynomials 'before' f */      /* skip the polynomials 'before' f */
                 if ( i < i0 ) continue;      if ( i < i0 ) continue;
                 if ( i == i0 ) i0 = 0;      if ( i == i0 ) i0 = 0;
                 /* set i-th bit */      /* set i-th bit */
                 fd[i>>5] |= (1<<(i&31));      fd[i>>5] |= (1<<(i&31));
                 for ( j = 2; j < i; j++ ) {       for ( j = 2; j < i; j++ ) {
                         /* skip the polynomials 'before' f */        /* skip the polynomials 'before' f */
                         if ( j < j0 ) continue;        if ( j < j0 ) continue;
                         if ( j == j0 ) j0 = 0;        if ( j == j0 ) j0 = 0;
                         /* set j-th bit */        /* set j-th bit */
                         fd[j>>5] |= (1<<(j&31));        fd[j>>5] |= (1<<(j&31));
                         for ( k = 1; k < j; k++ ) {         for ( k = 1; k < j; k++ ) {
                                 /* skip the polynomials 'before' f */          /* skip the polynomials 'before' f */
                                 if ( k < k0 ) continue;          if ( k < k0 ) continue;
                                 else if ( k == k0 ) { k0 = 0; continue; }          else if ( k == k0 ) { k0 = 0; continue; }
                                 /* set k-th bit */          /* set k-th bit */
                                 fd[k>>5] |= (1<<(k&31));          fd[k>>5] |= (1<<(k&31));
                                 ret = irredcheck_dddup2(f);          ret = irredcheck_dddup2(f);
                                 if ( ret == 1 ) return 0;          if ( ret == 1 ) return 0;
                                 /* reset k-th bit */          /* reset k-th bit */
                                 fd[k>>5] &= ~(1<<(k&31));          fd[k>>5] &= ~(1<<(k&31));
                         }        }
                         /* reset j-th bit */        /* reset j-th bit */
                         fd[j>>5] &= ~(1<<(j&31));        fd[j>>5] &= ~(1<<(j&31));
                 }      }
                 /* reset i-th bit */      /* reset i-th bit */
                 fd[i>>5] &= ~(1<<(i&31));      fd[i>>5] &= ~(1<<(i&31));
         }    }
         /* exhausted */    /* exhausted */
         return 1;    return 1;
 }  }
   
 void printqmat(Q **mat,int row,int col)  void printqmat(Q **mat,int row,int col)
 {  {
         int i,j;    int i,j;
   
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 for ( j = 0; j < col; j++ ) {      for ( j = 0; j < col; j++ ) {
                         printnum((Num)mat[i][j]); printf(" ");        printnum((Num)mat[i][j]); printf(" ");
                 }      }
                 printf("\n");      printf("\n");
         }    }
 }  }
   
 void printimat(int **mat,int row,int col)  void printimat(int **mat,int row,int col)
 {  {
         int i,j;    int i,j;
   
         for ( i = 0; i < row; i++ ) {    for ( i = 0; i < row; i++ ) {
                 for ( j = 0; j < col; j++ ) {      for ( j = 0; j < col; j++ ) {
                         printf("%d ",mat[i][j]);        printf("%d ",mat[i][j]);
                 }      }
                 printf("\n");      printf("\n");
         }    }
 }  }
   
 void Pnd_det(NODE arg,P *rp)  void Pnd_det(NODE arg,P *rp)
 {  {
         if ( argc(arg) == 1 )    if ( argc(arg) == 1 )
                 nd_det(0,ARG0(arg),rp);      nd_det(0,ARG0(arg),rp);
         else    else
                 nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);      nd_det(QTOS((Q)ARG1(arg)),ARG0(arg),rp);
 }  }
   
 void Pmat_col(NODE arg,VECT *rp)  void Pmat_col(NODE arg,VECT *rp)
 {  {
         int i,j,n;    int i,j,n;
         MAT mat;    MAT mat;
         VECT vect;    VECT vect;
   
         asir_assert(ARG0(arg),O_MAT,"mat_col");    asir_assert(ARG0(arg),O_MAT,"mat_col");
         asir_assert(ARG1(arg),O_N,"mat_col");    asir_assert(ARG1(arg),O_N,"mat_col");
         mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
         j = QTOS((Q)ARG1(arg));    j = QTOS((Q)ARG1(arg));
         if ( j < 0 || j >= mat->col) {    if ( j < 0 || j >= mat->col) {
                 error("mat_col : Out of range");      error("mat_col : Out of range");
         }    }
         n = mat->row;    n = mat->row;
         MKVECT(vect,n);    MKVECT(vect,n);
         for(i=0; i<n; i++) {    for(i=0; i<n; i++) {
                 BDY(vect)[i] = BDY(mat)[i][j];      BDY(vect)[i] = BDY(mat)[i][j];
         }    }
         *rp = vect;    *rp = vect;
 }  }
   
 NODE triangleq(NODE e)  NODE triangleq(NODE e)
Legend:
Removed from v.1.75  
changed lines
  Added in v.1.76
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>