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

Diff for /OpenXM_contrib2/asir2000/builtin/gf.c between version 1.9 and 1.16

version 1.9, 2001/06/28 08:57:20 version 1.16, 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/gf.c,v 1.8 2001/06/25 10:01:27 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/gf.c,v 1.15 2002/03/15 02:52:09 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
   
 struct resf_dlist {  struct resf_dlist {
         int ib,id;    int ib,id;
 };  };
   
 int resf_degtest(int,int *,int,struct resf_dlist *);  int resf_degtest(int,int *,int,struct resf_dlist *);
 void uhensel(P,NODE,int,int,NODE *);  void uhensel(P,NODE,int,int,NODE *);
   void uhensel_incremental(P,NODE,int,int,int,NODE *);
 void resf_hensel(int,P,int,P *,ML *);  void resf_hensel(int,P,int,P *,ML *);
 void resf_dtest(P,ML,int,int *,int *,DCP *);  void resf_dtest(P,ML,int,int *,int *,DCP *);
 void resf_dtest_special(P,ML,int,int *,int *,DCP *);  void resf_dtest_special(P,ML,int,int *,int *,DCP *);
Line 67  void nullspace_lm(LM **,int,int *);
Line 68  void nullspace_lm(LM **,int,int *);
 void nullspace_gf2n(GF2N **,int,int *);  void nullspace_gf2n(GF2N **,int,int *);
 void nullspace_gfpn(GFPN **,int,int *);  void nullspace_gfpn(GFPN **,int,int *);
 void nullspace_gfs(GFS **,int,int *);  void nullspace_gfs(GFS **,int,int *);
   void nullspace_gfsn(GFSN **,int,int *);
 void null_to_sol(int **,int *,int,int,UM *);  void null_to_sol(int **,int *,int,int,UM *);
   
 void showgfmat(UM **,int);  void showgfmat(UM **,int);
Line 74  void pwr_mod(P,P,V,P,int,N,P *);
Line 76  void pwr_mod(P,P,V,P,int,N,P *);
 void rem_mod(P,P,V,P,int,P *);  void rem_mod(P,P,V,P,int,P *);
   
 void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();  void Pnullspace(),Pgcda_mod(),Pftest(),Presfmain(),Ppwr_mod(),Puhensel();
   void Puhensel_incremental();
 void Psfuhensel();  void Psfuhensel();
   
 void Pnullspace_ff();  void Pnullspace_ff();
Line 89  void sfuhensel(P,NODE,GFS,int,NODE *);
Line 92  void sfuhensel(P,NODE,GFS,int,NODE *);
 extern int current_ff;  extern int current_ff;
   
 struct ftab gf_tab[] = {  struct ftab gf_tab[] = {
         {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},    {"solve_linear_equation_gf2n",Psolve_linear_equation_gf2n,1},
         {"nullspace",Pnullspace,3},    {"nullspace",Pnullspace,3},
         {"nullspace_ff",Pnullspace_ff,1},    {"nullspace_ff",Pnullspace_ff,1},
 /*      {"gcda_mod",Pgcda_mod,5}, */  /*  {"gcda_mod",Pgcda_mod,5}, */
         {"ftest",Pftest,2},    {"ftest",Pftest,2},
         {"resfmain",Presfmain,4},    {"resfmain",Presfmain,4},
         {"pwr_mod",Ppwr_mod,6},    {"pwr_mod",Ppwr_mod,6},
         {"uhensel",Puhensel,4},    {"uhensel",Puhensel,4},
         {"sfuhensel",Psfuhensel,4},    {"uhensel_incremental",Puhensel_incremental,5},
         {0,0,0},    {"sfuhensel",Psfuhensel,4},
     {0,0,0},
 };  };
   
 int resf_degtest(k,in,nb,dlist)  int resf_degtest(k,in,nb,dlist)
Line 107  int *in;
Line 111  int *in;
 int nb;  int nb;
 struct resf_dlist *dlist;  struct resf_dlist *dlist;
 {  {
         int i,d0,d;    int i,d0,d;
         int dl_i;    int dl_i;
         struct resf_dlist *t;    struct resf_dlist *t;
   
         if ( k < nb )    if ( k < nb )
                 return 0;      return 0;
         if ( nb == 1 )    if ( nb == 1 )
                 return 1;      return 1;
         d0 = 0; d = 0; dl_i = 0;    d0 = 0; d = 0; dl_i = 0;
         for ( i = 0; i < k; i++ ) {    for ( i = 0; i < k; i++ ) {
                 t = &dlist[in[i]];      t = &dlist[in[i]];
                 if ( t->ib > dl_i + 1 )      if ( t->ib > dl_i + 1 )
                         return 0;        return 0;
                 else if ( t->ib == dl_i )      else if ( t->ib == dl_i )
                         d += t->id;        d += t->id;
                 else if ( !d || (dl_i && d0 != d) )      else if ( !d || (dl_i && d0 != d) )
                         return 0;        return 0;
                 else {      else {
                         d0 = d; dl_i++; d = t->id;        d0 = d; dl_i++; d = t->id;
                 }      }
         }    }
         if ( dl_i != nb - 1 || d0 != d )    if ( dl_i != nb - 1 || d0 != d )
                 return 0;      return 0;
         else    else
                 return 1;      return 1;
 }  }
   
 void Puhensel(arg,rp)  void Puhensel(arg,rp)
 NODE arg;  NODE arg;
 LIST *rp;  LIST *rp;
 {  {
         P f;    P f;
         NODE mfl,r;    NODE mfl,r;
         int mod,bound;    int mod,bound;
   
         f = (P)ARG0(arg);    f = (P)ARG0(arg);
         mfl = BDY((LIST)ARG1(arg));    mfl = BDY((LIST)ARG1(arg));
         mod = QTOS((Q)ARG2(arg));    mod = QTOS((Q)ARG2(arg));
         bound = QTOS((Q)ARG3(arg));    bound = QTOS((Q)ARG3(arg));
         uhensel(f,mfl,mod,bound,&r);    uhensel(f,mfl,mod,bound,&r);
         MKLIST(*rp,r);    MKLIST(*rp,r);
 }  }
   
   void Puhensel_incremental(arg,rp)
   NODE arg;
   LIST *rp;
   {
     P f;
     NODE mfl,r;
     int mod,bound,start;
   
     f = (P)ARG0(arg);
     mfl = BDY((LIST)ARG1(arg));
     mod = QTOS((Q)ARG2(arg));
     start = QTOS((Q)ARG3(arg));
     bound = QTOS((Q)ARG4(arg));
     uhensel_incremental(f,mfl,mod,start,bound,&r);
     MKLIST(*rp,r);
   }
   
 void uhensel(f,mfl,mod,bound,rp)  void uhensel(f,mfl,mod,bound,rp)
 P f;  P f;
 NODE mfl;  NODE mfl;
 int mod,bound;  int mod,bound;
 NODE *rp;  NODE *rp;
 {  {
         ML blist,clist,rlist;    ML blist,clist,rlist;
         LUM fl;    LUM fl;
         int nf,i;    int nf,i;
         P s;    P s;
         V v;    V v;
         NODE t,top;    NODE t,top;
   
         nf = length(mfl);    nf = length(mfl);
         blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;    blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
         for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {    for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
                 blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));      blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
                 ptoum(mod,(P)BDY(t),blist->c[i]);      ptoum(mod,(P)BDY(t),blist->c[i]);
         }    }
         gcdgen(f,blist,&clist);    gcdgen(f,blist,&clist);
         blist->bound = clist->bound = bound;    blist->bound = clist->bound = bound;
         W_LUMALLOC((int)UDEG(f),bound,fl);    W_LUMALLOC((int)UDEG(f),bound,fl);
         ptolum(mod,bound,f,fl);    ptolum(mod,bound,f,fl);
         henmain(fl,blist,clist,&rlist);    henmain(fl,blist,clist,&rlist);
         v = VR(f);    v = VR(f);
         for ( i = nf-1, top = 0; i >= 0; i-- ) {    for ( i = nf-1, top = 0; i >= 0; i-- ) {
                 lumtop(v,mod,bound,rlist->c[i],&s);      lumtop(v,mod,bound,rlist->c[i],&s);
                 MKNODE(t,s,top); top = t;      MKNODE(t,s,top); top = t;
         }    }
         *rp = top;    *rp = top;
 }  }
   
   void uhensel_incremental(f,mfl,mod,start,bound,rp)
   P f;
   NODE mfl;
   int mod,start,bound;
   NODE *rp;
   {
     ML blist,clist,rlist;
     LUM fl;
     LUM *lblist;
     int nf,i,j,k;
     int **p;
     P s;
     V v;
     NODE t,top;
   
     nf = length(mfl);
     blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
     lblist = (LUM *)MALLOC(nf*sizeof(LUM));
     for ( i = 0, t = mfl; i < nf; i++, t = NEXT(t) ) {
       blist->c[i] = (pointer)UMALLOC(UDEG((P)BDY(t)));
       ptoum(mod,(P)BDY(t),blist->c[i]);
       W_LUMALLOC((int)UDEG((P)BDY(t)),bound,lblist[i]);
       ptolum(mod,start,(P)BDY(t),lblist[i]);
       p = lblist[i]->c;
       for ( j = DEG(lblist[i]); j >= 0; j-- )
         for ( k = start; k < bound; k++ )
           p[j][k] = 0;
     }
     gcdgen(f,blist,&clist);
     clist->bound = bound;
     W_LUMALLOC((int)UDEG(f),bound,fl);
     ptolum(mod,bound,f,fl);
     henmain_incremental(fl,lblist,clist,nf,mod,start,bound);
     v = VR(f);
     for ( i = nf-1, top = 0; i >= 0; i-- ) {
       lumtop_unsigned(v,mod,bound,lblist[i],&s);
       MKNODE(t,s,top); top = t;
     }
     *rp = top;
   }
   
 void Psfuhensel(arg,rp)  void Psfuhensel(arg,rp)
 NODE arg;  NODE arg;
 LIST *rp;  LIST *rp;
 {  {
         P f;    P f;
         int bound;    int bound;
         NODE r,mfl;    NODE r,mfl;
         GFS ev;    GFS ev;
   
         f = (P)ARG0(arg);    f = (P)ARG0(arg);
         mfl = BDY((LIST)ARG1(arg));    mfl = BDY((LIST)ARG1(arg));
         ev = (GFS)ARG2(arg);    ev = (GFS)ARG2(arg);
         bound = QTOS((Q)ARG3(arg));    bound = QTOS((Q)ARG3(arg));
         sfuhensel(f,mfl,ev,bound,&r);    sfuhensel(f,mfl,ev,bound,&r);
         MKLIST(*rp,r);    MKLIST(*rp,r);
 }  }
   
 void sfuhensel(f,mfl,ev,bound,rp)  void sfuhensel(f,mfl,ev,bound,rp)
Line 206  GFS ev;
Line 268  GFS ev;
 int bound;  int bound;
 NODE *rp;  NODE *rp;
 {  {
         BM fl;    BM fl;
         BM *r;    BM *r;
         VL vl,nvl;    VL vl,nvl;
         int i,fn,dx,dy;    int i,fn,dx,dy,d;
         NODE t,top;    NODE t,top;
         UM fm,hm,q;    UM fm,hm,q;
         UM *gm;    UM *gm;
         V x,y;    V x,y;
         P g,s,u;    P g,s,u;
   
         clctv(CO,f,&vl);    clctv(CO,f,&vl);
         if ( !vl || !vl->next || vl->next->next )    if ( !vl || !vl->next || vl->next->next )
                 error("sfuhensel : f must be a bivariate poly");      error("sfuhensel : f must be a bivariate poly");
   
         for ( i = 0, t = mfl; t; i++, t = NEXT(t) );    for ( i = 0, t = mfl; t; i++, t = NEXT(t) );
         fn = i;    fn = i;
   
         gm = (UM *)MALLOC(fn*sizeof(UM));    gm = (UM *)MALLOC(fn*sizeof(UM));
   
         /* XXX : more severe check is necessary */    /* XXX : more severe check is necessary */
         x = VR((P)BDY(mfl));    x = VR((P)BDY(mfl));
         y = vl->v == x ? vl->next->v : vl->v;    y = vl->v == x ? vl->next->v : vl->v;
   
         for ( i = 0, t = mfl; i < fn; i++, t = NEXT(t) ) {    for ( i = 0, t = mfl, d = 0; i < fn; i++, t = NEXT(t) ) {
                 gm[i] = (pointer)UMALLOC(getdeg(x,(P)BDY(t)));      gm[i] = (pointer)UMALLOC(getdeg(x,(P)BDY(t)));
                 ptosfum((P)BDY(t),gm[i]);      ptosfum((P)BDY(t),gm[i]);
         }      d += DEG(gm[i]);
     }
   
         /* reorder f if necessary */    /* reorder f if necessary */
         if ( vl->v != x ) {    if ( vl->v != x ) {
                 reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&g);      reordvar(vl,x,&nvl); reorderp(nvl,vl,f,&g);
                 vl = nvl; f = g;      vl = nvl; f = g;
         }    }
         dx = getdeg(x,f);    dx = getdeg(x,f);
         dy = getdeg(y,f);    if ( dx != d )
         if ( bound < dy+1 ) bound = dy+1;      error("sfuhensel : product of factors has incompatible degree");
         fl = BMALLOC(dx,bound);  
         ptosfbm(bound,f,fl);  
         shiftsfbm(bound,fl,FTOIF(CONT(ev)));  
   
         /* fm = fl mod y */    dy = getdeg(y,f);
         fm = W_UMALLOC(dx);    dy = MAX(dy,bound);
         cpyum(COEF(fl)[0],fm);    fl = BMALLOC(dx,dy);
         hm = W_UMALLOC(dx);    ptosfbm(dy,f,fl);
     if ( ev ) shiftsfbm(fl,FTOIF(CONT(ev)));
   
         q = W_UMALLOC(dx);    /* fm = fl mod y */
         r = (BM *)MLALLOC(fn*sizeof(BM));    fm = W_UMALLOC(dx);
         for ( i = 0; i < fn-1; i++ ) {    cpyum(COEF(fl)[0],fm);
                 /* fl = gm[i]*hm mod y */    hm = W_UMALLOC(dx);
                 divsfum(fm,gm[i],hm);  
                 /* fl is replaced by the cofactor of gk mod y^bound */  
                 /* r[i] = gk */  
                 sfhenmain2(fl,gm[i],hm,bound,r+i);  
                 cpyum(hm,fm);  
         }  
         /* finally, fl must be the lift of gm[fn-1] */  
         r[i] = fl;  
   
         for ( i = fn-1, top = 0; i >= 0; i-- ) {    q = W_UMALLOC(dx);
                 sfbmtop(bound,r[i],x,y,&s);    r = (BM *)MLALLOC(fn*sizeof(BM));
                 reorderp(CO,vl,s,&u);    for ( i = 0; i < fn-1; i++ ) {
                 MKNODE(t,u,top); top = t;      /* fl = gm[i]*hm mod y */
         }      divsfum(fm,gm[i],hm);
         *rp = top;      /* fl is replaced by the cofactor of gk mod y^bound */
       /* r[i] = gk */
       sfhenmain2(fl,gm[i],hm,bound,r+i);
       cpyum(hm,fm);
     }
     /* finally, fl must be the lift of gm[fn-1] */
     r[i] = fl;
   
     for ( i = fn-1, top = 0; i >= 0; i-- ) {
       sfbmtop(r[i],x,y,&s);
       reorderp(CO,vl,s,&u);
       MKNODE(t,u,top); top = t;
     }
     *rp = top;
 }  }
   
 void Presfmain(arg,rp)  void Presfmain(arg,rp)
 NODE arg;  NODE arg;
 LIST *rp;  LIST *rp;
 {  {
         P f;    P f;
         int mod,n,nb,i,j,k;    int mod,n,nb,i,j,k;
         int *nf,*md;    int *nf,*md;
         P *mf;    P *mf;
         NODE mfl,mdl,t,s,u;    NODE mfl,mdl,t,s,u;
         ML list;    ML list;
         DCP dc;    DCP dc;
         int sflag;    int sflag;
   
         f = (P)ARG0(arg);    f = (P)ARG0(arg);
         mod = QTOS((Q)ARG1(arg));    mod = QTOS((Q)ARG1(arg));
         mfl = BDY((LIST)ARG2(arg));    mfl = BDY((LIST)ARG2(arg));
         mdl = BDY((LIST)ARG3(arg));    mdl = BDY((LIST)ARG3(arg));
         for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )    for ( n = nb = 0, t = mfl; t; nb++, t = NEXT(t) )
                 for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );      for ( s = BDY((LIST)BDY(t)); s; n++, s = NEXT(s) );
         if ( n == nb ) {    if ( n == nb ) {
                 /* f must be irreducible! */      /* f must be irreducible! */
                 NEWDC(dc);      NEWDC(dc);
                 dc->c = f; dc->d = ONE;      dc->c = f; dc->d = ONE;
         } else {    } else {
                 nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));      nf = W_ALLOC(nb); md = W_ALLOC(nb); mf = (P *)ALLOCA(n*sizeof(P));
                 for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;      for ( k = i = 0, t = mfl, u = mdl, sflag = 1; t;
                         i++, t = NEXT(t), u = NEXT(u) ) {        i++, t = NEXT(t), u = NEXT(u) ) {
                         if ( sflag && length(BDY((LIST)BDY(t))) != 2 )        if ( sflag && length(BDY((LIST)BDY(t))) != 2 )
                                 sflag = 0;          sflag = 0;
                         for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )        for ( j = 0, s = BDY((LIST)BDY(t)); s; j++, s = NEXT(s) )
                                 mf[k++] = (P)BDY(s);          mf[k++] = (P)BDY(s);
                         nf[i] = j; md[i] = QTOS((Q)BDY(u));        nf[i] = j; md[i] = QTOS((Q)BDY(u));
                 }      }
                 resf_hensel(mod,f,n,mf,&list);      resf_hensel(mod,f,n,mf,&list);
                 if ( sflag )      if ( sflag )
                         resf_dtest_special(f,list,nb,nf,md,&dc);        resf_dtest_special(f,list,nb,nf,md,&dc);
                 else      else
                         resf_dtest(f,list,nb,nf,md,&dc);        resf_dtest(f,list,nb,nf,md,&dc);
         }    }
         dcptolist(dc,rp);    dcptolist(dc,rp);
 }  }
   
 void resf_hensel(mod,f,nf,mfl,listp)  void resf_hensel(mod,f,nf,mfl,listp)
Line 321  int nf;
Line 387  int nf;
 P *mfl;  P *mfl;
 ML *listp;  ML *listp;
 {  {
         register int i,j;    register int i,j;
         int q,n,bound,inv,lc;    int q,n,bound,inv,lc;
         int *p;    int *p;
         int **pp;    int **pp;
         ML blist,clist,bqlist,cqlist,rlist;    ML blist,clist,bqlist,cqlist,rlist;
         UM *b;    UM *b;
         LUM fl,tl;    LUM fl,tl;
         LUM *l;    LUM *l;
         UM w;    UM w;
   
         w = W_UMALLOC(UDEG(f));    w = W_UMALLOC(UDEG(f));
         blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;    blist = MLALLOC(nf); blist->n = nf; blist->mod = mod;
   
         /* c[0] must have lc(f) */    /* c[0] must have lc(f) */
         blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));    blist->c[0] = (pointer)UMALLOC(UDEG(mfl[0]));
         ptoum(mod,mfl[0],w);    ptoum(mod,mfl[0],w);
         inv = invm(w->c[UDEG(mfl[0])],mod);    inv = invm(w->c[UDEG(mfl[0])],mod);
         lc = rem(NM((Q)LC(f)),mod);    lc = rem(NM((Q)LC(f)),mod);
         if ( SGN((Q)LC(f)) < 0 )    if ( SGN((Q)LC(f)) < 0 )
                 lc = (mod-lc)%mod;      lc = (mod-lc)%mod;
         lc = dmar(inv,lc,0,mod);    lc = dmar(inv,lc,0,mod);
         if ( lc == 1 )    if ( lc == 1 )
                 copyum(w,blist->c[0]);      copyum(w,blist->c[0]);
         else    else
                 mulsum(mod,w,lc,blist->c[0]);      mulsum(mod,w,lc,blist->c[0]);
   
         /* c[i] (i=1,...,nf-1) must be monic */    /* c[i] (i=1,...,nf-1) must be monic */
         for ( i = 1; i < nf; i++ ) {    for ( i = 1; i < nf; i++ ) {
                 blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));      blist->c[i] = (pointer)UMALLOC(UDEG(mfl[i]));
                 ptoum(mod,mfl[i],w);      ptoum(mod,mfl[i],w);
                 inv = invm(w->c[UDEG(mfl[i])],mod);      inv = invm(w->c[UDEG(mfl[i])],mod);
                 if ( inv == 1 )      if ( inv == 1 )
                         copyum(w,blist->c[i]);        copyum(w,blist->c[i]);
                 else      else
                         mulsum(mod,w,inv,blist->c[i]);        mulsum(mod,w,inv,blist->c[i]);
         }    }
   
         gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);    gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
         n = bqlist->n; q = bqlist->mod;    n = bqlist->n; q = bqlist->mod;
         bqlist->bound = cqlist->bound = bound = mignotte(q,f);    bqlist->bound = cqlist->bound = bound = mignotte(q,f);
         if ( bound == 1 ) {    if ( bound == 1 ) {
                 *listp = rlist = MLALLOC(n);      *listp = rlist = MLALLOC(n);
                 rlist->n = n; rlist->mod = q; rlist->bound = bound;      rlist->n = n; rlist->mod = q; rlist->bound = bound;
                 for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {      for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
                         tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);        tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
                         for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )        for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
                                         pp[j][0] = p[j];            pp[j][0] = p[j];
                 }      }
         } else {    } else {
                 W_LUMALLOC((int)UDEG(f),bound,fl);      W_LUMALLOC((int)UDEG(f),bound,fl);
                 ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);      ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
         }    }
 }  }
   
 void resf_dtest(f,list,nb,nfl,mdl,dcp)  void resf_dtest(f,list,nb,nfl,mdl,dcp)
Line 382  int nb;
Line 448  int nb;
 int *nfl,*mdl;  int *nfl,*mdl;
 DCP *dcp;  DCP *dcp;
 {  {
         int n,np,bound,q;    int n,np,bound,q;
         int i,j,k;    int i,j,k;
         int *win;    int *win;
         P g,factor,cofactor;    P g,factor,cofactor;
         Q csum,csumt;    Q csum,csumt;
         DCP dcf,dcf0;    DCP dcf,dcf0;
         LUM *c;    LUM *c;
         ML wlist;    ML wlist;
         struct resf_dlist *dlist;    struct resf_dlist *dlist;
         int z;    int z;
   
         n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;    n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
         win = W_ALLOC(np+1);    win = W_ALLOC(np+1);
         ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;    ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
         wlist = W_MLALLOC(np); wlist->n = list->n;    wlist = W_MLALLOC(np); wlist->n = list->n;
         wlist->mod = list->mod; wlist->bound = list->bound;    wlist->mod = list->mod; wlist->bound = list->bound;
         c = (LUM *)COEF(wlist);    c = (LUM *)COEF(wlist);
         bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));    bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
         dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));    dlist = (struct resf_dlist *)ALLOCA(np*sizeof(struct resf_dlist));
         for ( i = 0, j = 0; j < nb; j++ )    for ( i = 0, j = 0; j < nb; j++ )
                 for ( k = 0; k < nfl[j]; k++, i++ ) {      for ( k = 0; k < nfl[j]; k++, i++ ) {
                         dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];        dlist[i].ib = j; dlist[i].id = DEG(c[i])/mdl[j];
                 }      }
         k = nb;    k = nb;
         for ( i = 0; i < nb; i++ )    for ( i = 0; i < nb; i++ )
                 win[i] = i+1;      win[i] = i+1;
         for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {    for ( g = f, dcf = dcf0 = 0, --np, z = 0; ; ) {
 #if 0  #if 0
                 if ( !(z++ % 10000) )      if ( !(z++ % 10000) )
                         fputc('.',stderr);        fputc('.',stderr);
 #endif  #endif
                 if ( resf_degtest(k,win,nb,dlist) &&      if ( resf_degtest(k,win,nb,dlist) &&
                         dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {        dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
                         NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;        NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
                         g = cofactor;        g = cofactor;
                         ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;        ucsump(g,&csum); mulq(csum,(Q)COEF(DC(g)),&csumt); csum = csumt;
                         for ( i = 0; i < k - 1; i++ )        for ( i = 0; i < k - 1; i++ )
                                 for ( j = win[i] + 1; j < win[i + 1]; j++ ) {          for ( j = win[i] + 1; j < win[i + 1]; j++ ) {
                                         c[j-i-1] = c[j];            c[j-i-1] = c[j];
                                         dlist[j-i-1] = dlist[j];            dlist[j-i-1] = dlist[j];
                                 }          }
                         for ( j = win[k-1] + 1; j <= np; j++ ) {        for ( j = win[k-1] + 1; j <= np; j++ ) {
                                 c[j-k] = c[j];          c[j-k] = c[j];
                                 dlist[j-k] = dlist[j];          dlist[j-k] = dlist[j];
                         }        }
                         if ( ( np -= k ) < k )        if ( ( np -= k ) < k )
                                 break;          break;
                         if ( np - win[0] + 1 < k )        if ( np - win[0] + 1 < k )
                                 if ( ++k > np )          if ( ++k > np )
                                         break;            break;
                                 else          else
                                         for ( i = 0; i < k; i++ )            for ( i = 0; i < k; i++ )
                                                 win[i] = i + 1;              win[i] = i + 1;
                         else        else
                                 for ( i = 1; i < k; i++ )          for ( i = 1; i < k; i++ )
                                         win[i] = win[0] + i;            win[i] = win[0] + i;
                 } else if ( !ncombi(1,np,k,win) )      } else if ( !ncombi(1,np,k,win) )
                         if ( k == np )        if ( k == np )
                                 break;          break;
                         else        else
                                 for ( i = 0, ++k; i < k; i++ )          for ( i = 0, ++k; i < k; i++ )
                                         win[i] = i + 1;            win[i] = i + 1;
         }    }
         NEXTDC(dcf0,dcf); COEF(dcf) = g;    NEXTDC(dcf0,dcf); COEF(dcf) = g;
         DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;    DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
 }  }
   
 void resf_dtest_special(f,list,nb,nfl,mdl,dcp)  void resf_dtest_special(f,list,nb,nfl,mdl,dcp)
Line 456  int nb;
Line 522  int nb;
 int *nfl,*mdl;  int *nfl,*mdl;
 DCP *dcp;  DCP *dcp;
 {  {
         int n,np,bound,q;    int n,np,bound,q;
         int i,j;    int i,j;
         int *win;    int *win;
         P g,factor,cofactor;    P g,factor,cofactor;
         Q csum,csumt;    Q csum,csumt;
         DCP dcf,dcf0;    DCP dcf,dcf0;
         LUM *c;    LUM *c;
         ML wlist;    ML wlist;
         int max;    int max;
   
         n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;    n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
         win = W_ALLOC(np+1);    win = W_ALLOC(np+1);
         ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;    ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
         wlist = W_MLALLOC(np); wlist->n = list->n;    wlist = W_MLALLOC(np); wlist->n = list->n;
         wlist->mod = list->mod; wlist->bound = list->bound;    wlist->mod = list->mod; wlist->bound = list->bound;
         c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));    c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
         max = 1<<nb;    max = 1<<nb;
         for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {    for ( g = f, dcf = dcf0 = 0, i = 0; i < max; i++ ) {
 #if 0  #if 0
                 if ( !(i % 1000) )      if ( !(i % 1000) )
                         fprintf(stderr,"i=%d\n",i);        fprintf(stderr,"i=%d\n",i);
 #endif  #endif
                 for ( j = 0; j < nb; j++ )      for ( j = 0; j < nb; j++ )
                         win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;        win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
                 if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {      if ( dtestmain(g,csum,wlist,nb,win,&factor,&cofactor) ) {
 #if 0  #if 0
                         fprintf(stderr,"success : i=%d\n",i);        fprintf(stderr,"success : i=%d\n",i);
 #endif  #endif
                         NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;        NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = factor;
                         NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;        NEXTDC(dcf0,dcf); DEG(dcf) = ONE; COEF(dcf) = cofactor;
                         NEXT(dcf) = 0;*dcp = dcf0;        NEXT(dcf) = 0;*dcp = dcf0;
                         return;        return;
                 }      }
         }    }
         NEXTDC(dcf0,dcf); COEF(dcf) = g;    NEXTDC(dcf0,dcf); COEF(dcf) = g;
         DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;    DEG(dcf) = ONE; NEXT(dcf) = 0;*dcp = dcf0;
 }  }
   
 #if 0  #if 0
Line 499  void Pftest(arg,rp)
Line 565  void Pftest(arg,rp)
 NODE arg;  NODE arg;
 P *rp;  P *rp;
 {  {
         ML list;    ML list;
         DCP dc;    DCP dc;
         P p;    P p;
         P *mfl;    P *mfl;
   
         p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);    p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
         hensel_special(4,1,p,mfl,&list);    hensel_special(4,1,p,mfl,&list);
         dtest_special(p,list,rp);    dtest_special(p,list,rp);
 }  }
   
 void dtest_special(f,list,pr)  void dtest_special(f,list,pr)
Line 514  P f;
Line 580  P f;
 ML list;  ML list;
 P *pr;  P *pr;
 {  {
         int n,np,bound,q;    int n,np,bound,q;
         int i,j,k;    int i,j,k;
         int *win;    int *win;
         P g,factor,cofactor;    P g,factor,cofactor;
         Q csum,csumt;    Q csum,csumt;
         DCP dc,dcf,dcf0;    DCP dc,dcf,dcf0;
         LUM *c;    LUM *c;
         ML wlist;    ML wlist;
   
         n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;    n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
         win = W_ALLOC(np+1);    win = W_ALLOC(np+1);
         ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;    ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
         wlist = W_MLALLOC(np); wlist->n = list->n;    wlist = W_MLALLOC(np); wlist->n = list->n;
         wlist->mod = list->mod; wlist->bound = list->bound;    wlist->mod = list->mod; wlist->bound = list->bound;
         c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));    c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
         for ( g = f, i = 130000; i < (1<<20); i++ ) {    for ( g = f, i = 130000; i < (1<<20); i++ ) {
 #if 0  #if 0
                 if ( !(i % 1000) )      if ( !(i % 1000) )
                         fprintf(stderr,"i=%d\n",i);        fprintf(stderr,"i=%d\n",i);
 #endif  #endif
                 for ( j = 0; j < 20; j++ )      for ( j = 0; j < 20; j++ )
                         win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;        win[j] = (i&(1<<j)) ? 2*j+1 : 2*j;
                 if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {      if ( dtestmain(g,csum,wlist,20,win,&factor,&cofactor) ) {
 #if 0  #if 0
                         fprintf(stderr,"success : i=%d\n",i);        fprintf(stderr,"success : i=%d\n",i);
 #endif  #endif
                         *pr = factor; return;        *pr = factor; return;
                 }      }
         }    }
 }  }
   
 void hensel_special(index,count,f,mfl,listp)  void hensel_special(index,count,f,mfl,listp)
Line 551  P f;
Line 617  P f;
 P *mfl;  P *mfl;
 ML *listp;  ML *listp;
 {  {
         register int i,j;    register int i,j;
         int q,n,t,d,r,u,br,tmp,bound;    int q,n,t,d,r,u,br,tmp,bound;
         int *c,*p,*m,*w;    int *c,*p,*m,*w;
         int **pp;    int **pp;
         DCP dc;    DCP dc;
         ML blist,clist,bqlist,cqlist,rlist;    ML blist,clist,bqlist,cqlist,rlist;
         UM *b;    UM *b;
         LUM fl,tl;    LUM fl,tl;
         LUM *l;    LUM *l;
   
         blist = MLALLOC(40); blist->n = 40; blist->mod = 11;    blist = MLALLOC(40); blist->n = 40; blist->mod = 11;
         for ( i = 0; i < 40; i++ ) {    for ( i = 0; i < 40; i++ ) {
                 blist->c[i] = (pointer)UMALLOC(6);      blist->c[i] = (pointer)UMALLOC(6);
                 ptoum(11,mfl[i],blist->c[i]);      ptoum(11,mfl[i],blist->c[i]);
         }    }
         gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);    gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
         n = bqlist->n; q = bqlist->mod;    n = bqlist->n; q = bqlist->mod;
         bqlist->bound = cqlist->bound = bound = mignotte(q,f);    bqlist->bound = cqlist->bound = bound = mignotte(q,f);
         if ( bound == 1 ) {    if ( bound == 1 ) {
                 *listp = rlist = MLALLOC(n);      *listp = rlist = MLALLOC(n);
                 rlist->n = n; rlist->mod = q; rlist->bound = bound;      rlist->n = n; rlist->mod = q; rlist->bound = bound;
                 for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {      for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
                         tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);        tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
                         for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )        for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
                                         pp[j][0] = p[j];            pp[j][0] = p[j];
                 }      }
         } else {    } else {
                 W_LUMALLOC(UDEG(f),bound,fl);      W_LUMALLOC(UDEG(f),bound,fl);
                 ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);      ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
         }    }
 }  }
 #endif  #endif
   
Line 589  void Pftest(arg,rp)
Line 655  void Pftest(arg,rp)
 NODE arg;  NODE arg;
 P *rp;  P *rp;
 {  {
         ML list;    ML list;
         DCP dc;    DCP dc;
         P p;    P p;
         P *mfl;    P *mfl;
   
         p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);    p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
         hensel_special(2,1,p,mfl,&list);    hensel_special(2,1,p,mfl,&list);
         dtest_special(p,list,rp);    dtest_special(p,list,rp);
 }  }
   
 void dtest_special(f,list,pr)  void dtest_special(f,list,pr)
Line 604  P f;
Line 670  P f;
 ML list;  ML list;
 P *pr;  P *pr;
 {  {
         int n,np,bound,q;    int n,np,bound,q;
         int i,j,k,t,b0;    int i,j,k,t,b0;
         int *win;    int *win;
         P g,factor,cofactor;    P g,factor,cofactor;
         Q csum,csumt;    Q csum,csumt;
         DCP dc,dcf,dcf0;    DCP dc,dcf,dcf0;
         LUM *c;    LUM *c;
         ML wlist;    ML wlist;
         static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};    static int nbits[16] = {0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
   
         n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;    n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
         win = W_ALLOC(np+1);    win = W_ALLOC(np+1);
         ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;    ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
         wlist = W_MLALLOC(np); wlist->n = list->n;    wlist = W_MLALLOC(np); wlist->n = list->n;
         wlist->mod = list->mod; wlist->bound = list->bound;    wlist->mod = list->mod; wlist->bound = list->bound;
         c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));    c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
         for ( g = f, i = 0; i < (1<<23); i++ ) {    for ( g = f, i = 0; i < (1<<23); i++ ) {
 #if 0  #if 0
                 if ( !(i % 1000) )      if ( !(i % 1000) )
                 fprintf(stderr,"i=%d\n",i);      fprintf(stderr,"i=%d\n",i);
 #endif  #endif
                 t = i>>20; b0 = nbits[t];      t = i>>20; b0 = nbits[t];
                 if ( !b0 )      if ( !b0 )
                         continue;        continue;
                 for ( j = 1; j < 6; j++ ) {      for ( j = 1; j < 6; j++ ) {
                         t = (i>>(20-4*j))&0xf;        t = (i>>(20-4*j))&0xf;
                         if ( nbits[t] != b0 )        if ( nbits[t] != b0 )
                                 break;          break;
                 }      }
                 if ( j != 6 )      if ( j != 6 )
                         continue;        continue;
                 for ( j = k = 0; j < 24; j++ )      for ( j = k = 0; j < 24; j++ )
                         if ( i & (1<<(23-j)) )        if ( i & (1<<(23-j)) )
                                 win[k++] = j;          win[k++] = j;
                 if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {      if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
 #if 0  #if 0
                         fprintf(stderr,"success : i=%d\n",i);        fprintf(stderr,"success : i=%d\n",i);
 #endif  #endif
                         *pr = factor; return;        *pr = factor; return;
                 }      }
         }    }
         *pr = f;    *pr = f;
 }  }
   
 void hensel_special(index,count,f,mfl,listp)  void hensel_special(index,count,f,mfl,listp)
Line 654  P f;
Line 720  P f;
 P *mfl;  P *mfl;
 ML *listp;  ML *listp;
 {  {
         register int i,j;    register int i,j;
         int q,n,t,d,r,u,br,tmp,bound;    int q,n,t,d,r,u,br,tmp,bound;
         int *c,*p,*m,*w;    int *c,*p,*m,*w;
         int **pp;    int **pp;
         DCP dc;    DCP dc;
         ML blist,clist,bqlist,cqlist,rlist;    ML blist,clist,bqlist,cqlist,rlist;
         UM *b;    UM *b;
         LUM fl,tl;    LUM fl,tl;
         LUM *l;    LUM *l;
   
         blist = MLALLOC(24); blist->n = 24; blist->mod = 5;    blist = MLALLOC(24); blist->n = 24; blist->mod = 5;
         for ( i = 0; i < 24; i++ ) {    for ( i = 0; i < 24; i++ ) {
                 blist->c[i] = (pointer)UMALLOC(7);      blist->c[i] = (pointer)UMALLOC(7);
                 ptoum(5,mfl[i],blist->c[i]);      ptoum(5,mfl[i],blist->c[i]);
         }    }
         gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);    gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
         n = bqlist->n; q = bqlist->mod;    n = bqlist->n; q = bqlist->mod;
         bqlist->bound = cqlist->bound = bound = mignotte(q,f);    bqlist->bound = cqlist->bound = bound = mignotte(q,f);
         if ( bound == 1 ) {    if ( bound == 1 ) {
                 *listp = rlist = MLALLOC(n);      *listp = rlist = MLALLOC(n);
                 rlist->n = n; rlist->mod = q; rlist->bound = bound;      rlist->n = n; rlist->mod = q; rlist->bound = bound;
                 for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {      for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
                         tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);        tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
                         for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )        for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
                                         pp[j][0] = p[j];            pp[j][0] = p[j];
                 }      }
         } else {    } else {
                 W_LUMALLOC(UDEG(f),bound,fl);      W_LUMALLOC(UDEG(f),bound,fl);
                 ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);      ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
         }    }
 }  }
 #endif  #endif
   
Line 691  void Pftest(arg,rp)
Line 757  void Pftest(arg,rp)
 NODE arg;  NODE arg;
 P *rp;  P *rp;
 {  {
         ML list;    ML list;
         P p;    P p;
         P *mfl;    P *mfl;
   
         p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);    p = (P)ARG0(arg); mfl = (P *)(((VECT)ARG1(arg))->body);
         hensel_special(5,1,p,mfl,&list);    hensel_special(5,1,p,mfl,&list);
         dtest_special(p,list,rp);    dtest_special(p,list,rp);
 }  }
   
 int nbits(a)  int nbits(a)
 int a;  int a;
 {  {
         int i,s;    int i,s;
   
         for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )    for ( i = 0, s = 0; a && (i < 20); i++, a >>= 1 )
                 if ( a & 1 ) s++;      if ( a & 1 ) s++;
         return s;    return s;
 }  }
   
 void dtest_special(f,list,pr)  void dtest_special(f,list,pr)
Line 715  P f;
Line 781  P f;
 ML list;  ML list;
 P *pr;  P *pr;
 {  {
         int n,np,bound,q;    int n,np,bound,q;
         int i,j,k,b0;    int i,j,k,b0;
         int *win;    int *win;
         P g,factor,cofactor;    P g,factor,cofactor;
         Q csum,csumt;    Q csum,csumt;
         LUM *c;    LUM *c;
         ML wlist;    ML wlist;
   
         n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;    n = UDEG(f); np = list->n; bound = list->bound; q = list->mod;
         win = W_ALLOC(np+1);    win = W_ALLOC(np+1);
         ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;    ucsump(f,&csum); mulq(csum,(Q)COEF(DC(f)),&csumt); csum = csumt;
         wlist = W_MLALLOC(np); wlist->n = list->n;    wlist = W_MLALLOC(np); wlist->n = list->n;
         wlist->mod = list->mod; wlist->bound = list->bound;    wlist->mod = list->mod; wlist->bound = list->bound;
         c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));    c = (LUM *)COEF(wlist); bcopy((char *)COEF(list),(char *)c,(int)(sizeof(LUM)*np));
         for ( g = f, i = 0; i < (1<<19); i++ ) {    for ( g = f, i = 0; i < (1<<19); i++ ) {
 #if 0  #if 0
                 if ( !(i % 10000) )      if ( !(i % 10000) )
                 fprintf(stderr,"i=%d\n",i);      fprintf(stderr,"i=%d\n",i);
 #endif  #endif
                 b0 = nbits(i>>10);      b0 = nbits(i>>10);
                 if ( !b0 || (nbits(i&0x3ff) != b0) )      if ( !b0 || (nbits(i&0x3ff) != b0) )
                         continue;        continue;
                 for ( j = k = 0; j < 20; j++ )      for ( j = k = 0; j < 20; j++ )
                         if ( i & (1<<(19-j)) )        if ( i & (1<<(19-j)) )
                                 win[k++] = j;          win[k++] = j;
                 if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {      if ( dtestmain(g,csum,wlist,k,win,&factor,&cofactor) ) {
 #if 0  #if 0
                         fprintf(stderr,"success : i=%d\n",i);        fprintf(stderr,"success : i=%d\n",i);
 #endif  #endif
                         *pr = factor; return;        *pr = factor; return;
                 }      }
         }    }
         *pr = f;    *pr = f;
 }  }
   
 void hensel_special(index,count,f,mfl,listp)  void hensel_special(index,count,f,mfl,listp)
Line 756  P f;
Line 822  P f;
 P *mfl;  P *mfl;
 ML *listp;  ML *listp;
 {  {
         register int i,j;    register int i,j;
         int q,n,bound;    int q,n,bound;
         int *p;    int *p;
         int **pp;    int **pp;
         ML blist,clist,bqlist,cqlist,rlist;    ML blist,clist,bqlist,cqlist,rlist;
         UM *b;    UM *b;
         LUM fl,tl;    LUM fl,tl;
         LUM *l;    LUM *l;
   
         blist = MLALLOC(20); blist->n = 20; blist->mod = 11;    blist = MLALLOC(20); blist->n = 20; blist->mod = 11;
         for ( i = 0; i < 20; i++ ) {    for ( i = 0; i < 20; i++ ) {
                 blist->c[i] = (pointer)UMALLOC(10);      blist->c[i] = (pointer)UMALLOC(10);
                 ptoum(11,mfl[i],blist->c[i]);      ptoum(11,mfl[i],blist->c[i]);
         }    }
         gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);    gcdgen(f,blist,&clist); henprep(f,blist,clist,&bqlist,&cqlist);
         n = bqlist->n; q = bqlist->mod;    n = bqlist->n; q = bqlist->mod;
         bqlist->bound = cqlist->bound = bound = mignotte(q,f);    bqlist->bound = cqlist->bound = bound = mignotte(q,f);
         if ( bound == 1 ) {    if ( bound == 1 ) {
                 *listp = rlist = MLALLOC(n);      *listp = rlist = MLALLOC(n);
                 rlist->n = n; rlist->mod = q; rlist->bound = bound;      rlist->n = n; rlist->mod = q; rlist->bound = bound;
                 for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {      for ( i = 0, b = (UM *)bqlist->c, l = (LUM *)rlist->c; i < n; i++ ) {
                         tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);        tl = LUMALLOC(DEG(b[i]),1); l[i] = tl; p = COEF(b[i]);
                         for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )        for ( j = 0, pp = COEF(tl); j <= DEG(tl); j++ )
                                         pp[j][0] = p[j];            pp[j][0] = p[j];
                 }      }
         } else {    } else {
                 W_LUMALLOC((int)UDEG(f),bound,fl);      W_LUMALLOC((int)UDEG(f),bound,fl);
                 ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);      ptolum(q,bound,f,fl); henmain(fl,bqlist,cqlist,listp);
         }    }
 }  }
   
 void Pnullspace(arg,rp)  void Pnullspace(arg,rp)
 NODE arg;  NODE arg;
 LIST *rp;  LIST *rp;
 {  {
         int i,j,n,mod;    int i,j,n,mod;
         MAT mat,r;    MAT mat,r;
         VECT u;    VECT u;
         V v;    V v;
         P p,z;    P p,z;
         Q q;    Q q;
         UM **w;    UM **w;
         UM mp;    UM mp;
         P *t;    P *t;
         UM *s;    UM *s;
         int *ind;    int *ind;
         NODE n0,n1;    NODE n0,n1;
   
         mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
         p = (P)ARG1(arg);    p = (P)ARG1(arg);
         v = VR(p);    v = VR(p);
         mod = QTOS((Q)ARG2(arg));    mod = QTOS((Q)ARG2(arg));
         n = mat->row;    n = mat->row;
         w = (UM **)almat_pointer(n,n);    w = (UM **)almat_pointer(n,n);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {      for ( j = 0, t = (P *)mat->body[i], s = w[i]; j < n; j++ ) {
                         ptomp(mod,t[j],&z);        ptomp(mod,t[j],&z);
                         s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);        s[j] = W_UMALLOC((z&&!NUM(z))?UDEG(z):0);
                         mptoum(z,s[j]);        mptoum(z,s[j]);
                 }      }
         mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);    mp = W_UMALLOC(UDEG(p)); ptoum(mod,p,mp);
         ind = (int *)ALLOCA(n*sizeof(int));    ind = (int *)ALLOCA(n*sizeof(int));
         nullspace(w,mp,mod,n,ind);    nullspace(w,mp,mod,n,ind);
         MKMAT(r,n,n);    MKMAT(r,n,n);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )      for ( j = 0, t = (P *)r->body[i], s = w[i]; j < n; j++ )
                         umtop(v,s[j],&t[j]);        umtop(v,s[j],&t[j]);
         MKVECT(u,n);    MKVECT(u,n);
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 STOQ(ind[i],q); u->body[i] = (pointer)q;      STOQ(ind[i],q); u->body[i] = (pointer)q;
         }    }
         MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);    MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
 }  }
   
 void nullspace(mat,p,mod,n,ind)  void nullspace(mat,p,mod,n,ind)
Line 836  UM p;
Line 902  UM p;
 int mod,n;  int mod,n;
 int *ind;  int *ind;
 {  {
         int i,j,l,s,d;    int i,j,l,s,d;
         UM q,w,w1,h,inv;    UM q,w,w1,h,inv;
         UM *t,*u;    UM *t,*u;
   
         d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);    d = DEG(p); inv = W_UMALLOC(d); q = W_UMALLOC(2*d);
         w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);    w = W_UMALLOC(2*d); w1 = W_UMALLOC(2*d); h = W_UMALLOC(d);
         bzero(ind,n*sizeof(int));    bzero(ind,n*sizeof(int));
         ind[0] = 0;    ind[0] = 0;
         for ( i = j = 0; j < n; i++, j++ ) {    for ( i = j = 0; j < n; i++, j++ ) {
                 for ( ; j < n; j++ ) {      for ( ; j < n; j++ ) {
                         for ( l = i; l < n; l++ )        for ( l = i; l < n; l++ )
                                 if ( DEG(mat[l][j])>=0 )          if ( DEG(mat[l][j])>=0 )
                                         break;            break;
                         if ( l < n ) {        if ( l < n ) {
                                 t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;          t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                         } else        } else
                                 ind[j] = 1;          ind[j] = 1;
                 }      }
                 if ( j == n )      if ( j == n )
                         break;        break;
                 invum(mod,p,mat[i][j],inv);      invum(mod,p,mat[i][j],inv);
                 for ( s = j, t = mat[i]; s < n; s++ ) {      for ( s = j, t = mat[i]; s < n; s++ ) {
                         mulum(mod,t[s],inv,w);        mulum(mod,t[s],inv,w);
                         DEG(w) = divum(mod,w,p,q);        DEG(w) = divum(mod,w,p,q);
                         cpyum(w,t[s]);        cpyum(w,t[s]);
                 }      }
                 for ( l = 0; l < n; l++ ) {      for ( l = 0; l < n; l++ ) {
                         if ( l == i )        if ( l == i )
                                 continue;          continue;
                         u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);        u = mat[l]; DEG(w) = -1; subum(mod,w,u[j],h);
                         for ( s = j; s < n; s++ ) {        for ( s = j; s < n; s++ ) {
                                 mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);          mulum(mod,h,t[s],w); addum(mod,w,u[s],w1);
                                 DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);          DEG(w1) = divum(mod,w1,p,q); cpyum(w1,u[s]);
                         }        }
                 }      }
         }    }
 }  }
   
 void Pnullspace_ff(arg,rp)  void Pnullspace_ff(arg,rp)
 NODE arg;  NODE arg;
 LIST *rp;  LIST *rp;
 {  {
         int i,j,n;    int i,j,n;
         Q mod;    MAT mat,r;
         MAT mat,r;    VECT u;
         VECT u;    Q q;
         Q q;    Obj **w;
         Obj **w;    Obj *t;
         Obj *t;    Obj *s;
         Obj *s;    int *ind;
         int *ind;    NODE n0,n1;
         NODE n0,n1;  
   
         mat = (MAT)ARG0(arg);    mat = (MAT)ARG0(arg);
         n = mat->row;    n = mat->row;
         w = (Obj **)almat_pointer(n,n);    w = (Obj **)almat_pointer(n,n);
         for ( i = 0; i < n; i++ )    for ( i = 0; i < n; i++ )
                 for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )      for ( j = 0, t = (Obj *)mat->body[i], s = w[i]; j < n; j++ )
                         s[j] = t[j];        s[j] = t[j];
         ind = (int *)ALLOCA(n*sizeof(int));    ind = (int *)ALLOCA(n*sizeof(int));
         switch ( current_ff ) {    switch ( current_ff ) {
                 case FF_GFP:      case FF_GFP:
                         nullspace_lm((LM **)w,n,ind); break;        nullspace_lm((LM **)w,n,ind); break;
                 case FF_GF2N:      case FF_GF2N:
                         nullspace_gf2n((GF2N **)w,n,ind); break;        nullspace_gf2n((GF2N **)w,n,ind); break;
                 case FF_GFPN:      case FF_GFPN:
                         nullspace_gfpn((GFPN **)w,n,ind); break;        nullspace_gfpn((GFPN **)w,n,ind); break;
                 case FF_GFS:      case FF_GFS:
                         nullspace_gfs((GFS **)w,n,ind); break;        nullspace_gfs((GFS **)w,n,ind); break;
                 default:      case FF_GFSN:
                         error("nullspace_ff : current_ff is not set");        nullspace_gfsn((GFSN **)w,n,ind); break;
         }      default:
         MKMAT(r,n,n);        error("nullspace_ff : current_ff is not set");
         for ( i = 0; i < n; i++ )    }
                 for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )    MKMAT(r,n,n);
                         t[j] = s[j];    for ( i = 0; i < n; i++ )
         MKVECT(u,n);      for ( j = 0, t = (Obj *)r->body[i], s = w[i]; j < n; j++ )
         for ( i = 0; i < n; i++ ) {        t[j] = s[j];
                 STOQ(ind[i],q); u->body[i] = (pointer)q;    MKVECT(u,n);
         }    for ( i = 0; i < n; i++ ) {
         MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);      STOQ(ind[i],q); u->body[i] = (pointer)q;
     }
     MKNODE(n1,u,0); MKNODE(n0,r,n1); MKLIST(*rp,n0);
 }  }
   
 void nullspace_lm(mat,n,ind)  void nullspace_lm(mat,n,ind)
Line 924  LM **mat;
Line 991  LM **mat;
 int n;  int n;
 int *ind;  int *ind;
 {  {
         int i,j,l,s;    int i,j,l,s;
         Q q,mod;    Q q,mod;
         N lm;    N lm;
         LM w,w1,h,inv;    LM w,w1,h,inv;
         LM *t,*u;    LM *t,*u;
   
         getmod_lm(&lm); NTOQ(lm,1,mod);    getmod_lm(&lm); NTOQ(lm,1,mod);
   
         bzero(ind,n*sizeof(int));    bzero(ind,n*sizeof(int));
         ind[0] = 0;    ind[0] = 0;
         for ( i = j = 0; j < n; i++, j++ ) {    for ( i = j = 0; j < n; i++, j++ ) {
                 for ( ; j < n; j++ ) {      for ( ; j < n; j++ ) {
                         for ( l = i; l < n; l++ ) {        for ( l = i; l < n; l++ ) {
                                 simplm(mat[l][j],&w); mat[l][j] = w;          simplm(mat[l][j],&w); mat[l][j] = w;
                                 if ( mat[l][j] )          if ( mat[l][j] )
                                         break;            break;
                         }        }
                         if ( l < n ) {        if ( l < n ) {
                                 t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;          t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                         } else        } else
                                 ind[j] = 1;          ind[j] = 1;
                 }      }
                 if ( j == n )      if ( j == n )
                         break;        break;
                 NTOQ(mat[i][j]->body,1,q); invl(q,mod,(Q *)&inv);      NTOQ(mat[i][j]->body,1,q); invl(q,mod,(Q *)&inv);
                 for ( s = j, t = mat[i]; s < n; s++ ) {      for ( s = j, t = mat[i]; s < n; s++ ) {
                         mullm(t[s],inv,&w); t[s] = w;        mullm(t[s],inv,&w); t[s] = w;
                 }      }
                 for ( l = 0; l < n; l++ ) {      for ( l = 0; l < n; l++ ) {
                         if ( l == i )        if ( l == i )
                                 continue;          continue;
                         u = mat[l]; chsgnlm(u[j],&h);        u = mat[l]; chsgnlm(u[j],&h);
                         for ( s = j; s < n; s++ ) {        for ( s = j; s < n; s++ ) {
                                 mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;          mullm(h,t[s],&w); addlm(w,u[s],&w1); u[s] = w1;
                         }        }
                 }      }
         }    }
 }  }
   
 void nullspace_gf2n(mat,n,ind)  void nullspace_gf2n(mat,n,ind)
Line 968  GF2N **mat;
Line 1035  GF2N **mat;
 int n;  int n;
 int *ind;  int *ind;
 {  {
         int i,j,l,s;    int i,j,l,s;
         GF2N w,w1,h,inv;    GF2N w,w1,h,inv;
         GF2N *t,*u;    GF2N *t,*u;
         extern gf2n_lazy;    extern gf2n_lazy;
   
         bzero(ind,n*sizeof(int));    bzero(ind,n*sizeof(int));
         ind[0] = 0;    ind[0] = 0;
         for ( i = j = 0; j < n; i++, j++ ) {    for ( i = j = 0; j < n; i++, j++ ) {
                 for ( ; j < n; j++ ) {      for ( ; j < n; j++ ) {
                         for ( l = i; l < n; l++ ) {        for ( l = i; l < n; l++ ) {
                                 simpgf2n(mat[l][j],&w); mat[l][j] = w;          simpgf2n(mat[l][j],&w); mat[l][j] = w;
                                 if ( mat[l][j] )          if ( mat[l][j] )
                                         break;            break;
                         }        }
                         if ( l < n ) {        if ( l < n ) {
                                 t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;          t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                         } else        } else
                                 ind[j] = 1;          ind[j] = 1;
                 }      }
                 if ( j == n )      if ( j == n )
                         break;        break;
                 invgf2n(mat[i][j],&inv);      invgf2n(mat[i][j],&inv);
                 for ( s = j, t = mat[i]; s < n; s++ ) {      for ( s = j, t = mat[i]; s < n; s++ ) {
                         mulgf2n(t[s],inv,&w); t[s] = w;        mulgf2n(t[s],inv,&w); t[s] = w;
                 }      }
                 for ( l = 0; l < n; l++ ) {      for ( l = 0; l < n; l++ ) {
                         if ( l == i )        if ( l == i )
                                 continue;          continue;
                         u = mat[l]; h = u[j];        u = mat[l]; h = u[j];
                         for ( s = j; s < n; s++ ) {        for ( s = j; s < n; s++ ) {
                                 mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;          mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
                         }        }
                 }      }
         }    }
 }  }
   
 void nullspace_gfpn(mat,n,ind)  void nullspace_gfpn(mat,n,ind)
Line 1009  GFPN **mat;
Line 1076  GFPN **mat;
 int n;  int n;
 int *ind;  int *ind;
 {  {
         int i,j,l,s;    int i,j,l,s;
         GFPN w,w1,h,inv;    GFPN w,w1,h,inv;
         GFPN *t,*u;    GFPN *t,*u;
   
         bzero(ind,n*sizeof(int));    bzero(ind,n*sizeof(int));
         ind[0] = 0;    ind[0] = 0;
         for ( i = j = 0; j < n; i++, j++ ) {    for ( i = j = 0; j < n; i++, j++ ) {
                 for ( ; j < n; j++ ) {      for ( ; j < n; j++ ) {
                         for ( l = i; l < n; l++ ) {        for ( l = i; l < n; l++ ) {
                                 simpgfpn(mat[l][j],&w); mat[l][j] = w;          simpgfpn(mat[l][j],&w); mat[l][j] = w;
                                 if ( mat[l][j] )          if ( mat[l][j] )
                                         break;            break;
                         }        }
                         if ( l < n ) {        if ( l < n ) {
                                 t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;          t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                         } else        } else
                                 ind[j] = 1;          ind[j] = 1;
                 }      }
                 if ( j == n )      if ( j == n )
                         break;        break;
                 divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);      divgfpn((GFPN)ONE,(GFPN)mat[i][j],&inv);
                 for ( s = j, t = mat[i]; s < n; s++ ) {      for ( s = j, t = mat[i]; s < n; s++ ) {
                         mulgfpn(t[s],inv,&w); t[s] = w;        mulgfpn(t[s],inv,&w); t[s] = w;
                 }      }
                 for ( l = 0; l < n; l++ ) {      for ( l = 0; l < n; l++ ) {
                         if ( l == i )        if ( l == i )
                                 continue;          continue;
                         u = mat[l]; chsgngfpn(u[j],&h);        u = mat[l]; chsgngfpn(u[j],&h);
                         for ( s = j; s < n; s++ ) {        for ( s = j; s < n; s++ ) {
                                 mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;          mulgfpn(h,t[s],&w); addgfpn(w,u[s],&w1); u[s] = w1;
                         }        }
                 }      }
         }    }
 }  }
   
 void nullspace_gfs(mat,n,ind)  void nullspace_gfs(mat,n,ind)
Line 1049  GFS **mat;
Line 1116  GFS **mat;
 int n;  int n;
 int *ind;  int *ind;
 {  {
         int i,j,l,s;    int i,j,l,s;
         GFS w,w1,h,inv;    GFS w,w1,h,inv;
         GFS *t,*u;    GFS *t,*u;
         GFS one;    GFS one;
   
         bzero(ind,n*sizeof(int));    bzero(ind,n*sizeof(int));
         ind[0] = 0;    ind[0] = 0;
         mqtogfs(ONEM,&one);    mqtogfs(ONEM,&one);
   
         for ( i = j = 0; j < n; i++, j++ ) {    for ( i = j = 0; j < n; i++, j++ ) {
                 for ( ; j < n; j++ ) {      for ( ; j < n; j++ ) {
                         for ( l = i; l < n; l++ )        for ( l = i; l < n; l++ )
                                 if ( mat[l][j] )          if ( mat[l][j] )
                                         break;            break;
                         if ( l < n ) {        if ( l < n ) {
                                 t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;          t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                         } else        } else
                                 ind[j] = 1;          ind[j] = 1;
                 }      }
                 if ( j == n )      if ( j == n )
                         break;        break;
                 divgfs(one,mat[i][j],&inv);      divgfs(one,mat[i][j],&inv);
                 for ( s = j, t = mat[i]; s < n; s++ ) {      for ( s = j, t = mat[i]; s < n; s++ ) {
                         mulgfs(t[s],inv,&w); t[s] = w;        mulgfs(t[s],inv,&w); t[s] = w;
                 }      }
                 for ( l = 0; l < n; l++ ) {      for ( l = 0; l < n; l++ ) {
                         if ( l == i )        if ( l == i )
                                 continue;          continue;
                         u = mat[l];        u = mat[l];
                         chsgngfs(u[j],&h);        chsgngfs(u[j],&h);
                         for ( s = j; s < n; s++ ) {        for ( s = j; s < n; s++ ) {
                                 mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;          mulgfs(h,t[s],&w); addgfs(w,u[s],&w1); u[s] = w1;
                         }        }
                 }      }
         }    }
 }  }
   
   void nullspace_gfsn(mat,n,ind)
   GFSN **mat;
   int n;
   int *ind;
   {
     int i,j,l,s;
     GFSN w,w1,h,inv;
     GFSN *t,*u;
   
     bzero(ind,n*sizeof(int));
     ind[0] = 0;
   
     for ( i = j = 0; j < n; i++, j++ ) {
       for ( ; j < n; j++ ) {
         for ( l = i; l < n; l++ )
           if ( mat[l][j] )
             break;
         if ( l < n ) {
           t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
         } else
           ind[j] = 1;
       }
       if ( j == n )
         break;
       invgfsn(mat[i][j],&inv);
       for ( s = j, t = mat[i]; s < n; s++ ) {
         mulgfsn(t[s],inv,&w); t[s] = w;
       }
       for ( l = 0; l < n; l++ ) {
         if ( l == i )
           continue;
         u = mat[l];
         chsgngfsn(u[j],&h);
         for ( s = j; s < n; s++ ) {
           mulgfsn(h,t[s],&w); addgfsn(w,u[s],&w1); u[s] = w1;
         }
       }
     }
   }
   
 /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */  /* p = a(0)vl[0]+a(1)vl[1]+...+a(m-1)vl[m-1]+a(m) -> array = [a(0) a(1) ... a(m)] */
   
 void linear_form_to_array(p,vl,m,array)  void linear_form_to_array(p,vl,m,array)
Line 1094  VL vl;
Line 1201  VL vl;
 int m;  int m;
 Num *array;  Num *array;
 {  {
         int i;    int i;
         DCP dc;    DCP dc;
   
         bzero((char *)array,(m+1)*sizeof(Num *));    bzero((char *)array,(m+1)*sizeof(Num *));
         for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {    for ( i = 0; p && vl; vl = NEXT(vl), i++ ) {
                 if ( ID(p) == O_N )      if ( ID(p) == O_N )
                         break;        break;
                 else if ( VR(p) == vl->v ) {      else if ( VR(p) == vl->v ) {
                         dc = DC(p);        dc = DC(p);
                         array[i] = (Num)COEF(dc);        array[i] = (Num)COEF(dc);
                         dc = NEXT(dc);        dc = NEXT(dc);
                         p = dc ? COEF(dc) : 0;        p = dc ? COEF(dc) : 0;
                 }      }
         }    }
         array[m] = (Num)p;    array[m] = (Num)p;
 }  }
   
 void array_to_linear_form(array,vl,m,r)  void array_to_linear_form(array,vl,m,r)
Line 1117  VL vl;
Line 1224  VL vl;
 int m;  int m;
 P *r;  P *r;
 {  {
         P t;    P t;
         DCP dc0,dc1;    DCP dc0,dc1;
   
         if ( !m )    if ( !m )
                 *r = (P)array[0];      *r = (P)array[0];
         else {    else {
                 array_to_linear_form(array+1,NEXT(vl),m-1,&t);      array_to_linear_form(array+1,NEXT(vl),m-1,&t);
                 if ( !array[0] )      if ( !array[0] )
                         *r = t;        *r = t;
                 else {      else {
                         NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];        NEWDC(dc0); DEG(dc0) = ONE; COEF(dc0) = (P)array[0];
                         if ( !t )        if ( !t )
                                 NEXT(dc0) = 0;          NEXT(dc0) = 0;
                         else {        else {
                                 NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;          NEWDC(dc1); DEG(dc1) = 0; COEF(dc1) = t;
                                 NEXT(dc1) = 0;          NEXT(dc1) = 0;
                                 NEXT(dc0) = dc1;          NEXT(dc0) = dc1;
                         }        }
                         MKP(vl->v,dc0,*r);        MKP(vl->v,dc0,*r);
                 }      }
         }    }
 }  }
   
 void Psolve_linear_equation_gf2n(arg,rp)  void Psolve_linear_equation_gf2n(arg,rp)
 NODE arg;  NODE arg;
 LIST *rp;  LIST *rp;
 {  {
         NODE eqs,tn;    NODE eqs,tn;
         VL vars,tvl;    VL vars,tvl;
         int i,j,n,m,dim,codim;    int i,j,n,m,dim,codim;
         GF2N **w;    GF2N **w;
         int *ind;    int *ind;
         NODE n0,n1;    NODE n0,n1;
   
         get_vars(ARG0(arg),&vars);    get_vars(ARG0(arg),&vars);
         eqs = BDY((LIST)ARG0(arg));    eqs = BDY((LIST)ARG0(arg));
         for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);    for ( n = 0, tn = eqs; tn; tn = NEXT(tn), n++);
         for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);    for ( m = 0, tvl = vars; tvl; tvl = NEXT(tvl), m++);
         w = (GF2N **)almat_pointer(n,m+1);    w = (GF2N **)almat_pointer(n,m+1);
         for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )    for ( i = 0, tn = eqs; i < n; i++, tn = NEXT(tn) )
                 linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);      linear_form_to_array(BDY(tn),vars,m,(Num *)w[i]);
         ind = (int *)ALLOCA(m*sizeof(int));    ind = (int *)ALLOCA(m*sizeof(int));
         solve_linear_equation_gf2n(w,n,m,ind);    solve_linear_equation_gf2n(w,n,m,ind);
         for ( j = 0, dim = 0; j < m; j++ )    for ( j = 0, dim = 0; j < m; j++ )
                 if ( ind[j] )      if ( ind[j] )
                         dim++;        dim++;
         codim = m-dim;    codim = m-dim;
         for ( i = codim; i < n; i++ )    for ( i = codim; i < n; i++ )
                 if ( w[i][m] ) {      if ( w[i][m] ) {
                         MKLIST(*rp,0); return;        MKLIST(*rp,0); return;
                 }      }
         for ( i = 0, n0 = 0; i < codim; i++ ) {    for ( i = 0, n0 = 0; i < codim; i++ ) {
                 NEXTNODE(n0,n1);      NEXTNODE(n0,n1);
                 array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));      array_to_linear_form((Num *)w[i],vars,m,(P *)&BDY(n1));
         }    }
         if ( n0 )    if ( n0 )
                 NEXT(n1) = 0;      NEXT(n1) = 0;
         MKLIST(*rp,n0);    MKLIST(*rp,n0);
 }  }
   
 void solve_linear_equation_gf2n(mat,n,m,ind)  void solve_linear_equation_gf2n(mat,n,m,ind)
Line 1182  GF2N **mat;
Line 1289  GF2N **mat;
 int n;  int n;
 int *ind;  int *ind;
 {  {
         int i,j,l,s;    int i,j,l,s;
         GF2N w,w1,h,inv;    GF2N w,w1,h,inv;
         GF2N *t,*u;    GF2N *t,*u;
         extern gf2n_lazy;    extern gf2n_lazy;
   
         bzero(ind,m*sizeof(int));    bzero(ind,m*sizeof(int));
         ind[0] = 0;    ind[0] = 0;
         for ( i = j = 0; j < m; i++, j++ ) {    for ( i = j = 0; j < m; i++, j++ ) {
                 for ( ; j < m; j++ ) {      for ( ; j < m; j++ ) {
                         for ( l = i; l < n; l++ ) {        for ( l = i; l < n; l++ ) {
                                 simpgf2n(mat[l][j],&w); mat[l][j] = w;          simpgf2n(mat[l][j],&w); mat[l][j] = w;
                                 if ( mat[l][j] )          if ( mat[l][j] )
                                         break;            break;
                         }        }
                         if ( l < n ) {        if ( l < n ) {
                                 t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;          t = mat[i]; mat[i] = mat[l]; mat[l] = t; break;
                         } else        } else
                                 ind[j] = 1;          ind[j] = 1;
                 }      }
                 if ( j == m )      if ( j == m )
                         break;        break;
                 invgf2n(mat[i][j],&inv);      invgf2n(mat[i][j],&inv);
                 for ( s = j, t = mat[i]; s <= m; s++ ) {      for ( s = j, t = mat[i]; s <= m; s++ ) {
                         mulgf2n(t[s],inv,&w); t[s] = w;        mulgf2n(t[s],inv,&w); t[s] = w;
                 }      }
                 for ( l = 0; l < n; l++ ) {      for ( l = 0; l < n; l++ ) {
                         if ( l == i )        if ( l == i )
                                 continue;          continue;
                         u = mat[l]; h = u[j];        u = mat[l]; h = u[j];
                         for ( s = j; s <= m; s++ ) {        for ( s = j; s <= m; s++ ) {
                                 mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;          mulgf2n(h,t[s],&w); addgf2n(w,u[s],&w1); u[s] = w1;
                         }        }
                 }      }
         }    }
 }  }
   
 /*  /*
Line 1225  int *ind;
Line 1332  int *ind;
 int mod,n;  int mod,n;
 UM *r;  UM *r;
 {  {
         int i,j,k,l;    int i,j,k,l;
         int *c;    int *c;
         UM w;    UM w;
   
         for ( i = 0, l = 0; i < n; i++ ) {    for ( i = 0, l = 0; i < n; i++ ) {
                 if ( !ind[i] )      if ( !ind[i] )
                         continue;        continue;
                 w = UMALLOC(n);      w = UMALLOC(n);
                 for ( j = k = 0, c = COEF(w); j < n; j++ )      for ( j = k = 0, c = COEF(w); j < n; j++ )
                         if ( ind[j] )        if ( ind[j] )
                                 c[j] = 0;          c[j] = 0;
                         else        else
                                 c[j] = mat[k++][i];          c[j] = mat[k++][i];
                 c[i] = mod-1;      c[i] = mod-1;
                 for ( j = n; j >= 0; j-- )      for ( j = n; j >= 0; j-- )
                         if ( c[j] )        if ( c[j] )
                                 break;          break;
                 DEG(w) = j;      DEG(w) = j;
                 r[l++] = w;      r[l++] = w;
         }    }
 }  }
 */  */
   
Line 1252  void showgfmat(mat,n)
Line 1359  void showgfmat(mat,n)
 UM **mat;  UM **mat;
 int n;  int n;
 {  {
         int i,j,k;    int i,j,k;
         int *c;    int *c;
         UM p;    UM p;
   
         for ( i = 0; i < n; i++ ) {    for ( i = 0; i < n; i++ ) {
                 for ( j = 0; j < n; j++ ) {      for ( j = 0; j < n; j++ ) {
                         p = mat[i][j];        p = mat[i][j];
                         if ( DEG(p) < 0 )        if ( DEG(p) < 0 )
                                 fprintf(asir_out,"0");          fprintf(asir_out,"0");
                         else        else
                                 for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {          for ( p = mat[i][j], k = DEG(p), c = COEF(p); k >= 0; k-- ) {
                                         if ( c[k] )            if ( c[k] )
                                                 fprintf(asir_out,"+%d",c[k]);              fprintf(asir_out,"+%d",c[k]);
                                         if ( k > 1 )            if ( k > 1 )
                                                 fprintf(asir_out,"a^%d",k);              fprintf(asir_out,"a^%d",k);
                                         else if ( k == 1 )            else if ( k == 1 )
                                                 fprintf(asir_out,"a",k);              fprintf(asir_out,"a",k);
                                 }          }
                         fprintf(asir_out," ");        fprintf(asir_out," ");
                 }      }
                 fprintf(asir_out,"\n");      fprintf(asir_out,"\n");
         }    }
 }  }
   
 #if 0  #if 0
Line 1281  void Pgcda_mod(arg,rp)
Line 1388  void Pgcda_mod(arg,rp)
 NODE arg;  NODE arg;
 P *rp;  P *rp;
 {  {
         p1 = (P)ARG0(arg);    p1 = (P)ARG0(arg);
         p2 = (P)ARG1(arg);    p2 = (P)ARG1(arg);
         v = VR((P)ARG2(arg));    v = VR((P)ARG2(arg));
         d = (P)ARG3(arg);    d = (P)ARG3(arg);
         m = QTOS((Q)ARG4(arg));    m = QTOS((Q)ARG4(arg));
         reordvar(CO,v,&vl);    reordvar(CO,v,&vl);
         reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);    reorderp(vl,CO,p1,&t); ptomp(m,t,&m1);
         reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);    reorderp(vl,CO,p2,&t); ptomp(m,t,&m2);
         if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {    if ( NUM(m1) || NUM(m2) || VR(m1) != v || VR(m2) != v ) {
                 *rp = ONE; return;      *rp = ONE; return;
         }    }
         if ( deg(v,m1) >= deg(v,m2) ) {    if ( deg(v,m1) >= deg(v,m2) ) {
                 t = m1; m1 = m2; m2 = t;      t = m1; m1 = m2; m2 = t;
         }    }
         while ( 1 ) {    while ( 1 ) {
                 inva_mod(COEF(DC(m2)),d,m,&inv);      inva_mod(COEF(DC(m2)),d,m,&inv);
                 NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);      NEWDC(dc); NEXT(dc) = 0; MKP(v,dc,h);
                 d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));      d0 = deg(v,m1)-deg(v,m2); STOQ(d0,DEG(dc));
                 mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));      mulgq(m,d,inv,COEF(DC(m1)),&COEF(dc));
                 mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);      mulgp(vl,m,d,m2,h,&t); subgp(vl,m,d,m1,t,&s);
         }    }
 }  }
 #endif  #endif
   
Line 1309  void Ppwr_mod(arg,rp)
Line 1416  void Ppwr_mod(arg,rp)
 NODE arg;  NODE arg;
 P *rp;  P *rp;
 {  {
         P p,a,d,r;    P p,a,d,r;
         int m;    int m;
         Q q;    Q q;
         N n;    N n;
   
         m = QTOS((Q)ARG4(arg)); q = (Q)ARG5(arg); n = q ? NM(q) : 0;    m = QTOS((Q)ARG4(arg)); q = (Q)ARG5(arg); n = q ? NM(q) : 0;
         ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);    ptomp(m,(P)ARG0(arg),&p); ptomp(m,(P)ARG1(arg),&a);
         ptomp(m,(P)ARG3(arg),&d);    ptomp(m,(P)ARG3(arg),&d);
         pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);    pwr_mod(p,a,VR((P)ARG2(arg)),d,m,n,&r);
         mptop(r,rp);    mptop(r,rp);
 }  }
   
 void pwr_mod(p,a,v,d,m,n,rp)  void pwr_mod(p,a,v,d,m,n,rp)
Line 1328  int m;
Line 1435  int m;
 N n;  N n;
 P *rp;  P *rp;
 {  {
         int b;    int b;
         P t,s,r;    P t,s,r;
         N n1;    N n1;
   
         if ( !n )    if ( !n )
                 *rp = (P)ONEM;      *rp = (P)ONEM;
         else if ( UNIN(n) )    else if ( UNIN(n) )
                 *rp = p;      *rp = p;
         else {    else {
                 b = divin(n,2,&n1);      b = divin(n,2,&n1);
                 pwr_mod(p,a,v,d,m,n1,&t);      pwr_mod(p,a,v,d,m,n1,&t);
                 mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);      mulmp(CO,m,t,t,&s); rem_mod(s,a,v,d,m,&r);
                 if ( b ) {      if ( b ) {
                         mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);        mulmp(CO,m,r,p,&t); rem_mod(t,a,v,d,m,rp);
                 } else      } else
                         *rp = r;        *rp = r;
         }    }
 }  }
   
 void rem_mod(p,a,v,d,m,rp)  void rem_mod(p,a,v,d,m,rp)
Line 1353  V v;
Line 1460  V v;
 int m;  int m;
 P *rp;  P *rp;
 {  {
         P tmp,r1,r2;    P tmp,r1,r2;
   
         divsrmp(CO,m,p,d,&tmp,&r1);    divsrmp(CO,m,p,d,&tmp,&r1);
         divsrmp(CO,m,r1,a,&tmp,&r2);    divsrmp(CO,m,r1,a,&tmp,&r2);
         divsrmp(CO,m,r2,d,&tmp,rp);    divsrmp(CO,m,r2,d,&tmp,rp);
 }  }
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.16
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>