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

Diff for /OpenXM_contrib2/asir2000/engine/dalg.c between version 1.10 and 1.18

version 1.10, 2005/08/02 07:16:42 version 1.18, 2018/03/29 01:32:51
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.9 2005/07/11 00:24:02 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.17 2017/08/31 02:36:21 noro Exp $
 */  */
   
 #include "ca.h"  #include "ca.h"
Line 15  void dalgtoalg(DAlg da,Alg *r);
Line 15  void dalgtoalg(DAlg da,Alg *r);
   
 NumberField get_numberfield()  NumberField get_numberfield()
 {  {
         return current_numberfield;    return current_numberfield;
 }  }
   
 void setfield_dalg(NODE alist)  void setfield_dalg(NODE alist)
 {  {
         NumberField nf;    NumberField nf;
         VL vl,vl1,vl2;    VL vl,vl1,vl2;
         int n,i,dim;    int n,i,dim;
         Alg *gen;    Alg *gen;
         P *defpoly;    P *defpoly;
         P p;    P p;
         Q c,iq,two;    Q c,iq,two;
         DP *ps,*mb;    DP *ps,*mb;
         DP one;    DP one;
         NODE t,b,b1,b2,hlist,mblist;    NODE t,b,b1,b2,hlist,mblist;
         struct order_spec *current_spec;    struct order_spec *current_spec;
   
         nf = (NumberField)MALLOC(sizeof(struct oNumberField));    nf = (NumberField)MALLOC(sizeof(struct oNumberField));
         current_numberfield = nf;    current_numberfield = nf;
         vl = 0;    vl = 0;
         for ( t = alist; t; t = NEXT(t) ) {    for ( t = alist; t; t = NEXT(t) ) {
                 clctalg(BDY((Alg)BDY(t)),&vl1);      clctalg((P)BDY((Alg)BDY(t)),&vl1);
                 mergev(ALG,vl,vl1,&vl2); vl = vl2;      mergev(ALG,vl,vl1,&vl2); vl = vl2;
         }    }
         for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ );    for ( n = 0, vl1 = vl; vl1; vl1 = NEXT(vl1), n++ );
         nf->n = n;    nf->n = n;
         nf->vl = vl;    nf->vl = vl;
         nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P));    nf->defpoly = defpoly = (P *)MALLOC(n*sizeof(P));
         nf->ps = ps = (DP *)MALLOC(n*sizeof(DP));    nf->ps = ps = (DP *)MALLOC(n*sizeof(DP));
         current_spec = dp_current_spec;    current_spec = dp_current_spec;
         STOQ(2,two);    STOQ(2,two);
         create_order_spec(0,(Obj)two,&nf->spec);    create_order_spec(0,(Obj)two,&nf->spec);
         initd(nf->spec);    initd(nf->spec);
         for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) {    for ( b = hlist = 0, i = 0, vl1 = vl; i < n; vl1 = NEXT(vl1), i++ ) {
                 ptozp(vl1->v->attr,1,&c,&defpoly[i]);      ptozp(vl1->v->attr,1,&c,&defpoly[i]);
                 ptod(ALG,vl,defpoly[i],&ps[i]);      ptod(ALG,vl,defpoly[i],&ps[i]);
                 STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;      STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
                 MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;      MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
         }    }
         ptod(ALG,vl,(P)ONE,&one);    ptod(ALG,vl,(P)ONE,&one);
         MKDAlg(one,ONE,nf->one);    MKDAlg(one,ONE,nf->one);
         nf->ind = b;    nf->ind = b;
         dp_mbase(hlist,&mblist);    dp_mbase(hlist,&mblist);
         initd(current_spec);    initd(current_spec);
         nf->dim = dim = length(mblist);    nf->dim = dim = length(mblist);
         nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));    nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
         for ( i = 0, t = mblist; t; t = NEXT(t), i++ )    for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
                 mb[dim-i-1] = (DP)BDY(t);      mb[dim-i-1] = (DP)BDY(t);
 }  }
   
 void setfield_gb(NODE gb,VL vl,struct order_spec *spec)  void setfield_gb(NODE gb,VL vl,struct order_spec *spec)
 {  {
         NumberField nf;    NumberField nf;
         VL vl1,vl2;    VL vl1,vl2;
         int n,i,dim;    int n,i,dim;
         Alg *gen;    Alg *gen;
         P *defpoly;    P *defpoly;
         P p;    P p;
         Q c,iq,two;    Q c,iq,two;
         DP *ps,*mb;    DP *ps,*mb;
         DP one;    DP one;
         NODE t,b,b1,b2,hlist,mblist;    NODE t,b,b1,b2,hlist,mblist;
         struct order_spec *current_spec;    struct order_spec *current_spec;
   
         nf = (NumberField)MALLOC(sizeof(struct oNumberField));    nf = (NumberField)MALLOC(sizeof(struct oNumberField));
         current_numberfield = nf;    current_numberfield = nf;
         for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ );    for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ );
         nf->n = n;    nf->n = n;
         nf->psn = length(gb);    nf->psn = length(gb);
         nf->vl = vl;    nf->vl = vl;
         nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P));    nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P));
         nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP));    nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP));
         current_spec = dp_current_spec;    current_spec = dp_current_spec;
         nf->spec = spec;    nf->spec = spec;
         initd(nf->spec);    initd(nf->spec);
         for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) {    for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) {
                 ptozp((P)BDY(t),1,&c,&defpoly[i]);      ptozp((P)BDY(t),1,&c,&defpoly[i]);
                 ptod(CO,vl,defpoly[i],&ps[i]);      ptod(CO,vl,defpoly[i],&ps[i]);
                 STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;      STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
                 MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;      MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
         }    }
         ptod(ALG,vl,(P)ONE,&one);    ptod(ALG,vl,(P)ONE,&one);
         MKDAlg(one,ONE,nf->one);    MKDAlg(one,ONE,nf->one);
         nf->ind = b;    nf->ind = b;
         dp_mbase(hlist,&mblist);    dp_mbase(hlist,&mblist);
         initd(current_spec);    initd(current_spec);
         nf->dim = dim = length(mblist);    nf->dim = dim = length(mblist);
         nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));    nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
         for ( i = 0, t = mblist; t; t = NEXT(t), i++ )    for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
                 mb[dim-i-1] = (DP)BDY(t);      mb[dim-i-1] = (DP)BDY(t);
 }  }
   
 void qtodalg(Q q,DAlg *r)  void qtodalg(Q q,DAlg *r)
 {  {
         NumberField nf;    NumberField nf;
         Q t;    Q t;
         DP nm;    DP nm;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("qtodalg : current_numberfield is not set");      error("qtodalg : current_numberfield is not set");
         if ( !q )    if ( !q )
                 *r = 0;      *r = 0;
         else if ( NID(q) == N_DA )    else if ( NID(q) == N_DA )
                 *r = (DAlg)q;      *r = (DAlg)q;
         else if ( NID(q) == N_Q ) {    else if ( NID(q) == N_Q ) {
                 if ( INT(q) ) {      if ( INT(q) ) {
                         muldc(CO,nf->one->nm,(P)q,&nm);        muldc(CO,nf->one->nm,(Obj)q,&nm);
                         MKDAlg(nm,ONE,*r);        MKDAlg(nm,ONE,*r);
                 } else {      } else {
                         NTOQ(NM(q),SGN(q),t);        NTOQ(NM(q),SGN(q),t);
                         muldc(CO,nf->one->nm,(P)t,&nm);        muldc(CO,nf->one->nm,(Obj)t,&nm);
                         NTOQ(DN(q),1,t);        NTOQ(DN(q),1,t);
                         MKDAlg(nm,t,*r);        MKDAlg(nm,t,*r);
                 }      }
         } else    } else
                 error("qtodalg : invalid argument");      error("qtodalg : invalid argument");
 }  }
   
 void obj_algtodalg(Obj obj,Obj *r)  void obj_algtodalg(Obj obj,Obj *r)
 {  {
         DAlg d;    DAlg d;
         DCP dc,dcr0,dcr;    DCP dc,dcr0,dcr;
         P c,p;    P c,p;
         Obj t;    Obj t;
         Obj nm,dn;    Obj nm,dn;
         NODE b,s,s0;    NODE b,s,s0;
         R rat;    R rat;
         VECT v;    VECT v;
         MAT mat;    MAT mat;
         LIST list;    LIST list;
         pointer *a;    pointer *a;
         pointer **m;    pointer **m;
         int len,row,col,i,j,l;    int len,row,col,i,j,l;
   
         if ( !obj ) {    if ( !obj ) {
                 *r = 0;      *r = 0;
                 return;      return;
         }    }
         switch ( OID(obj) ) {    switch ( OID(obj) ) {
                 case O_N:      case O_N:
                         algtodalg((Alg)obj,&d); *r = (Obj)d;        algtodalg((Alg)obj,&d); *r = (Obj)d;
                         break;        break;
                 case O_P:      case O_P:
                         for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {        for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
                                 obj_algtodalg((Obj)COEF(dc),&t);          obj_algtodalg((Obj)COEF(dc),&t);
                                 if ( t ) {          if ( t ) {
                                         NEXTDC(dcr0,dcr);            NEXTDC(dcr0,dcr);
                                         COEF(dcr) = (P)t;            COEF(dcr) = (P)t;
                                         DEG(dcr) = DEG(dc);            DEG(dcr) = DEG(dc);
                                 }          }
                         }        }
                         if ( dcr0 ) {        if ( dcr0 ) {
                                 MKP(VR((P)obj),dcr0,p);          MKP(VR((P)obj),dcr0,p);
                                 *r = (Obj)p;          *r = (Obj)p;
                         } else        } else
                                 *r = 0;          *r = 0;
                         break;        break;
                 case O_R:      case O_R:
                         obj_algtodalg((Obj)NM((R)obj),&nm);        obj_algtodalg((Obj)NM((R)obj),&nm);
                         obj_algtodalg((Obj)DN((R)obj),&dn);        obj_algtodalg((Obj)DN((R)obj),&dn);
                         if ( !dn )        if ( !dn )
                                 error("obj_algtodalg : division by 0");          error("obj_algtodalg : division by 0");
                         if ( !nm )        if ( !nm )
                                 *r = 0;          *r = 0;
                         else {        else {
                                 MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;          MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
                         }        }
                         break;        break;
                 case O_LIST:      case O_LIST:
                         s0 = 0;        s0 = 0;
                         for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {        for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
                                 NEXTNODE(s0,s);          NEXTNODE(s0,s);
                                 obj_algtodalg((Obj)BDY(b),&t);          obj_algtodalg((Obj)BDY(b),&t);
                                 BDY(s) = (pointer)t;          BDY(s) = (pointer)t;
                         }        }
                         NEXT(s) = 0;        NEXT(s) = 0;
                         MKLIST(list,s0);        MKLIST(list,s0);
                         *r = (Obj)list;        *r = (Obj)list;
                         break;        break;
                 case O_VECT:      case O_VECT:
                         l = ((VECT)obj)->len;        l = ((VECT)obj)->len;
                         a = BDY((VECT)obj);        a = BDY((VECT)obj);
                         MKVECT(v,l);        MKVECT(v,l);
                         for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
                                 obj_algtodalg((Obj)a[i],&t);          obj_algtodalg((Obj)a[i],&t);
                                 BDY(v)[i] = (pointer)t;          BDY(v)[i] = (pointer)t;
                         }        }
                         *r = (Obj)v;        *r = (Obj)v;
                         break;        break;
                 case O_MAT:      case O_MAT:
                         row = ((MAT)obj)->row; col = ((MAT)obj)->col;        row = ((MAT)obj)->row; col = ((MAT)obj)->col;
                         m = BDY((MAT)obj);        m = BDY((MAT)obj);
                         MKMAT(mat,row,col);        MKMAT(mat,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++ ) {
                                         obj_algtodalg((Obj)m[i][j],&t);            obj_algtodalg((Obj)m[i][j],&t);
                                         BDY(mat)[i][j] = (pointer)t;            BDY(mat)[i][j] = (pointer)t;
                                 }          }
                         *r = (Obj)mat;        *r = (Obj)mat;
                         break;        break;
                 default:      default:
                         *r = obj;        *r = obj;
                         break;        break;
         }    }
 }  }
   
 void obj_dalgtoalg(Obj obj,Obj *r)  void obj_dalgtoalg(Obj obj,Obj *r)
 {  {
         Alg d;    Alg d;
         DCP dc,dcr0,dcr;    DCP dc,dcr0,dcr;
         P c,p;    P c,p;
         Obj t;    Obj t;
         Obj nm,dn;    Obj nm,dn;
         NODE b,s,s0;    NODE b,s,s0;
         R rat;    R rat;
         VECT v;    VECT v;
         MAT mat;    MAT mat;
         LIST list;    LIST list;
         pointer *a;    pointer *a;
         pointer **m;    pointer **m;
         int len,row,col,i,j,l;    int len,row,col,i,j,l;
   
         if ( !obj ) {    if ( !obj ) {
                 *r = 0;      *r = 0;
                 return;      return;
         }    }
         switch ( OID(obj) ) {    switch ( OID(obj) ) {
                 case O_N:      case O_N:
                         dalgtoalg((DAlg)obj,&d); *r = (Obj)d;        dalgtoalg((DAlg)obj,&d); *r = (Obj)d;
                         break;        break;
                 case O_P:      case O_P:
                         for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {        for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
                                 obj_dalgtoalg((Obj)COEF(dc),&t);          obj_dalgtoalg((Obj)COEF(dc),&t);
                                 if ( t ) {          if ( t ) {
                                         NEXTDC(dcr0,dcr);            NEXTDC(dcr0,dcr);
                                         COEF(dcr) = (P)t;            COEF(dcr) = (P)t;
                                         DEG(dcr) = DEG(dc);            DEG(dcr) = DEG(dc);
                                 }          }
                         }        }
                         if ( dcr0 ) {        if ( dcr0 ) {
                                 MKP(VR((P)obj),dcr0,p);          MKP(VR((P)obj),dcr0,p);
                                 *r = (Obj)p;          *r = (Obj)p;
                         } else        } else
                                 *r = 0;          *r = 0;
                         break;        break;
                 case O_R:      case O_R:
                         obj_dalgtoalg((Obj)NM((R)obj),&nm);        obj_dalgtoalg((Obj)NM((R)obj),&nm);
                         obj_dalgtoalg((Obj)DN((R)obj),&dn);        obj_dalgtoalg((Obj)DN((R)obj),&dn);
                         if ( !dn )        if ( !dn )
                                 error("obj_dalgtoalg : division by 0");          error("obj_dalgtoalg : division by 0");
                         if ( !nm )        if ( !nm )
                                 *r = 0;          *r = 0;
                         else {        else {
                                 MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;          MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
                         }        }
                         break;        break;
                 case O_LIST:      case O_LIST:
                         s0 = 0;        s0 = 0;
                         for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {        for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
                                 NEXTNODE(s0,s);          NEXTNODE(s0,s);
                                 obj_dalgtoalg((Obj)BDY(b),&t);          obj_dalgtoalg((Obj)BDY(b),&t);
                                 BDY(s) = (pointer)t;          BDY(s) = (pointer)t;
                         }        }
                         NEXT(s) = 0;        NEXT(s) = 0;
                         MKLIST(list,s0);        MKLIST(list,s0);
                         *r = (Obj)list;        *r = (Obj)list;
                         break;        break;
                 case O_VECT:      case O_VECT:
                         l = ((VECT)obj)->len;        l = ((VECT)obj)->len;
                         a = BDY((VECT)obj);        a = BDY((VECT)obj);
                         MKVECT(v,l);        MKVECT(v,l);
                         for ( i = 0; i < l; i++ ) {        for ( i = 0; i < l; i++ ) {
                                 obj_dalgtoalg((Obj)a[i],&t);          obj_dalgtoalg((Obj)a[i],&t);
                                 BDY(v)[i] = (pointer)t;          BDY(v)[i] = (pointer)t;
                         }        }
                         *r = (Obj)v;        *r = (Obj)v;
                         break;        break;
                 case O_MAT:      case O_MAT:
                         row = ((MAT)obj)->row; col = ((MAT)obj)->col;        row = ((MAT)obj)->row; col = ((MAT)obj)->col;
                         m = BDY((MAT)obj);        m = BDY((MAT)obj);
                         MKMAT(mat,row,col);        MKMAT(mat,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++ ) {
                                         obj_dalgtoalg((Obj)m[i][j],&t);            obj_dalgtoalg((Obj)m[i][j],&t);
                                         BDY(mat)[i][j] = (pointer)t;            BDY(mat)[i][j] = (pointer)t;
                                 }          }
                         *r = (Obj)mat;        *r = (Obj)mat;
                         break;        break;
                 default:      default:
                         *r = obj;        *r = obj;
                         break;        break;
         }    }
 }  }
   
 void algtodalg(Alg a,DAlg *r)  void algtodalg(Alg a,DAlg *r)
 {  {
         P ap,p,p1;    P ap,p,p1;
         Q c,c1,d1,dn,nm;    Q c,c1,d1,dn,nm;
         DP dp;    DP dp;
         DAlg da;    DAlg da;
         NumberField nf;    NumberField nf;
         struct order_spec *current_spec;    struct order_spec *current_spec;
         VL vl,tvl,svl;    VL vl,tvl,svl;
         V v;    V v;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("algtodalg : current_numberfield is not set");      error("algtodalg : current_numberfield is not set");
         if ( !a ) {    if ( !a ) {
                 *r = 0;      *r = 0;
                 return;      return;
         }    }
         switch (NID((Num)a) ) {    switch (NID((Num)a) ) {
                 case N_Q:      case N_Q:
                         c = (Q)a;        c = (Q)a;
                         if ( INT(c) ) {        if ( INT(c) ) {
                                 muldc(CO,nf->one->nm,(P)c,&dp);          muldc(CO,nf->one->nm,(Obj)c,&dp);
                                 MKDAlg(dp,ONE,*r);          MKDAlg(dp,ONE,*r);
                         } else {        } else {
                                 NTOQ(NM(c),SGN(c),c1);          NTOQ(NM(c),SGN(c),c1);
                                 NTOQ(DN(c),1,d1);          NTOQ(DN(c),1,d1);
                                 muldc(CO,nf->one->nm,(P)c1,&dp);          muldc(CO,nf->one->nm,(Obj)c1,&dp);
                                 MKDAlg(dp,d1,*r);          MKDAlg(dp,d1,*r);
                         }        }
                         break;        break;
                 case N_A:      case N_A:
                         ap = (P)BDY(a);        ap = (P)BDY(a);
                         ptozp(ap,1,&c,&p);        ptozp(ap,1,&c,&p);
                         if ( INT(c) ) {        if ( INT(c) ) {
                                 p = ap;          p = ap;
                                 dn = ONE;          dn = ONE;
                         } else {        } else {
                                 NTOQ(NM(c),SGN(c),nm);          NTOQ(NM(c),SGN(c),nm);
                                 NTOQ(DN(c),1,dn);          NTOQ(DN(c),1,dn);
                                 mulpq(p,(P)nm,&p1); p = p1;          mulpq(p,(P)nm,&p1); p = p1;
                         }        }
                         current_spec = dp_current_spec; initd(nf->spec);        current_spec = dp_current_spec; initd(nf->spec);
                         get_vars(p,&vl);        get_vars((Obj)p,&vl);
                         for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {        for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
                                 v = tvl->v;          v = tvl->v;
                                 for ( svl = nf->vl; svl; svl = NEXT(svl) )          for ( svl = nf->vl; svl; svl = NEXT(svl) )
                                         if ( v == svl->v )            if ( v == svl->v )
                                                 break;              break;
                                 if ( !svl )          if ( !svl )
                                         error("algtodalg : incompatible numberfield");            error("algtodalg : incompatible numberfield");
                         }        }
                         ptod(ALG,nf->vl,p,&dp);        ptod(ALG,nf->vl,p,&dp);
                         MKDAlg(dp,dn,da);        MKDAlg(dp,dn,da);
                         simpdalg(da,r);        simpdalg(da,r);
                         break;        break;
                 default:      default:
                         error("algtodalg : invalid argument");        error("algtodalg : invalid argument");
                         break;        break;
         }    }
 }  }
   
 void dalgtoalg(DAlg da,Alg *r)  void dalgtoalg(DAlg da,Alg *r)
 {  {
         NumberField nf;    NumberField nf;
         P p,p1;    P p,p1;
         Q inv;    Q inv;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("dalgtoalg : current_numberfield is not set");      error("dalgtoalg : current_numberfield is not set");
         dtop(ALG,nf->vl,da->nm,&p);    if ( !da ) *r = 0;
         invq(da->dn,&inv);    else {
         mulpq(p,(P)inv,&p1);      dtop(ALG,nf->vl,da->nm,(Obj *)&p);
         MKAlg(p1,*r);      invq(da->dn,&inv);
       mulpq(p,(P)inv,&p1);
       MKAlg(p1,*r);
     }
 }  }
   
 void simpdalg(DAlg da,DAlg *r)  void simpdalg(DAlg da,DAlg *r)
 {  {
         NumberField nf;    NumberField nf;
         DP nm;    DP nm;
         DAlg d;    DAlg d;
         Q dn,dn1;    Q dn,dn1;
         struct order_spec *current_spec;    struct order_spec *current_spec;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("simpdalg : current_numberfield is not set");      error("simpdalg : current_numberfield is not set");
         if ( !da ) {    if ( !da ) {
                 *r = 0;      *r = 0;
                 return;      return;
         }    }
         current_spec = dp_current_spec; initd(nf->spec);    current_spec = dp_current_spec; initd(nf->spec);
         dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);    dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,(P *)&dn);
         if ( !nm ) *r = 0;    if ( !nm ) *r = 0;
         else {    else {
                 initd(current_spec);      initd(current_spec);
                 mulq(da->dn,dn,&dn1);      mulq(da->dn,dn,&dn1);
                 MKDAlg(nm,dn1,d);      MKDAlg(nm,dn1,d);
                 rmcontdalg(d,r);      rmcontdalg(d,r);
         }    }
 }  }
   
 void adddalg(DAlg a,DAlg b,DAlg *c)  void adddalg(DAlg a,DAlg b,DAlg *c)
 {  {
         NumberField nf;    NumberField nf;
         Q dna,dnb,a1,b1,dn,g;    Q dna,dnb,a1,b1,dn,g;
         N an,bn,gn;    N an,bn,gn;
         DAlg t;    DAlg t;
         DP ta,tb,nm;    DP ta,tb,nm;
         struct order_spec *current_spec;    struct order_spec *current_spec;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("adddalg : current_numberfield is not set");      error("adddalg : current_numberfield is not set");
         if ( !a )    if ( !a )
                 *c = b;      *c = b;
         else if ( !b )    else if ( !b )
                 *c = a;      *c = a;
         else {    else {
                 qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;      qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                 dna = a->dn;      dna = a->dn;
                 dnb = b->dn;      dnb = b->dn;
                 gcdn(NM(dna),NM(dnb),&gn);      gcdn(NM(dna),NM(dnb),&gn);
                 divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);      divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
                 NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);      NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
                 /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */      /* nma/dna+nmb/dnb = (nma*b1+nmb*a1)/(dna*b1) */
                 muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);      muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb);
                 current_spec = dp_current_spec; initd(nf->spec);      current_spec = dp_current_spec; initd(nf->spec);
                 addd(CO,ta,tb,&nm);      addd(CO,ta,tb,&nm);
                 initd(current_spec);      initd(current_spec);
                 if ( !nm )      if ( !nm )
                         *c = 0;        *c = 0;
                 else {      else {
                         mulq(dna,b1,&dn);        mulq(dna,b1,&dn);
                         MKDAlg(nm,dn,*c);        MKDAlg(nm,dn,*c);
                 }      }
         }    }
 }  }
   
 void subdalg(DAlg a,DAlg b,DAlg *c)  void subdalg(DAlg a,DAlg b,DAlg *c)
 {  {
         NumberField nf;    NumberField nf;
         Q dna,dnb,a1,b1,dn,g;    Q dna,dnb,a1,b1,dn,g;
         N an,bn,gn;    N an,bn,gn;
         DP ta,tb,nm;    DP ta,tb,nm;
         DAlg t;    DAlg t;
         struct order_spec *current_spec;    struct order_spec *current_spec;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("subdalg : current_numberfield is not set");      error("subdalg : current_numberfield is not set");
         if ( !a )    if ( !a )
                 *c = b;      *c = b;
         else if ( !b )    else if ( !b )
                 *c = a;      *c = a;
         else {    else {
                 qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;      qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                 dna = a->dn;      dna = a->dn;
                 dnb = b->dn;      dnb = b->dn;
                 gcdn(NM(dna),NM(dnb),&gn);      gcdn(NM(dna),NM(dnb),&gn);
                 divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);      divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
                 NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);      NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),b1);
                 /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */      /* nma/dna-nmb/dnb = (nma*b1-nmb*a1)/(dna*b1) */
                 muldc(CO,a->nm,(P)b1,&ta); muldc(CO,b->nm,(P)a1,&tb);      muldc(CO,a->nm,(Obj)b1,&ta); muldc(CO,b->nm,(Obj)a1,&tb);
                 current_spec = dp_current_spec; initd(nf->spec);      current_spec = dp_current_spec; initd(nf->spec);
                 subd(CO,ta,tb,&nm);      subd(CO,ta,tb,&nm);
                 initd(current_spec);      initd(current_spec);
                 if ( !nm )      if ( !nm )
                         *c = 0;        *c = 0;
                 else {      else {
                         mulq(dna,b1,&dn);        mulq(dna,b1,&dn);
                         MKDAlg(nm,dn,*c);        MKDAlg(nm,dn,*c);
                 }      }
         }    }
 }  }
   
 void muldalg(DAlg a,DAlg b,DAlg *c)  void muldalg(DAlg a,DAlg b,DAlg *c)
 {  {
         NumberField nf;    NumberField nf;
         DP nm;    DP nm;
         Q dn;    Q dn;
         DAlg t;    DAlg t;
         struct order_spec *current_spec;    struct order_spec *current_spec;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("muldalg : current_numberfield is not set");      error("muldalg : current_numberfield is not set");
         if ( !a || !b )    if ( !a || !b )
                 *c = 0;      *c = 0;
         else {    else {
                 qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;      qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                 current_spec = dp_current_spec; initd(nf->spec);      current_spec = dp_current_spec; initd(nf->spec);
                 muld(CO,a->nm,b->nm,&nm);      muld(CO,a->nm,b->nm,&nm);
                 initd(current_spec);      initd(current_spec);
                 mulq(a->dn,b->dn,&dn);      mulq(a->dn,b->dn,&dn);
                 MKDAlg(nm,dn,t);      MKDAlg(nm,dn,t);
                 simpdalg(t,c);      simpdalg(t,c);
         }    }
 }  }
   
   
 void divdalg(DAlg a,DAlg b,DAlg *c)  void divdalg(DAlg a,DAlg b,DAlg *c)
 {  {
         DAlg inv,t;    DAlg inv,t;
         int ret;    int ret;
   
         if ( !current_numberfield )    if ( !current_numberfield )
                 error("divdalg : current_numberfield is not set");      error("divdalg : current_numberfield is not set");
         if ( !b )    if ( !b )
                 error("divdalg : division by 0");      error("divdalg : division by 0");
         if ( !a )    if ( !a )
                 c = 0;      c = 0;
         else {    else {
                 qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;      qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                 ret = invdalg(b,&inv);      ret = invdalg(b,&inv);
                 if ( !ret ) {      if ( !ret ) {
                         error("divdalg : the denominator is not invertible");        error("divdalg : the denominator is not invertible");
                 }      }
                 muldalg(a,inv,c);      muldalg(a,inv,c);
         }    }
 }  }
   
 void rmcontdalg(DAlg a, DAlg *r)  void rmcontdalg(DAlg a, DAlg *r)
 {  {
         DP u,u1;    DP u,u1;
         Q cont,c,d;    Q cont,c,d;
         N gn,cn,dn;    N gn,cn,dn;
   
         if ( !a )    if ( !a )
                 *r = a;      *r = a;
         else {    else {
                 dp_ptozp(a->nm,&u);      dp_ptozp(a->nm,&u);
                 divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);      divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
                 gcdn(NM(cont),NM(a->dn),&gn);      gcdn(NM(cont),NM(a->dn),&gn);
                 divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);      divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
                 divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);      divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
                 muldc(CO,u,(P)c,&u1);      muldc(CO,u,(Obj)c,&u1);
                 MKDAlg(u1,d,*r);      MKDAlg(u1,d,*r);
         }    }
 }  }
   
 int invdalg(DAlg a,DAlg *c)  int invdalg(DAlg a,DAlg *c)
 {  {
         NumberField nf;    NumberField nf;
         int dim,n,i,j,k,l;    int dim,n,i,j,k,l;
         DP *mb;    DP *mb;
         DP m,d,u;    DP m,d,u;
         N ln,gn,qn;    N ln,gn,qn;
         DAlg *simp;    DAlg *simp;
         DAlg t,a0,r;    DAlg t,a0,r;
         Q dn,dnsol,mul,nmc,dn1;    Q dn,dnsol,mul,nmc,dn1;
         MAT mobj,sol;    MAT mobj,sol;
         Q **mat,**solmat;    Q **mat,**solmat;
         MP mp0,mp;    MP mp0,mp;
         int *rinfo,*cinfo;    int *rinfo,*cinfo;
         int rank,nparam;    int rank,nparam;
         NODE nd0,nd,ndt;    NODE nd0,nd,ndt;
         struct order_spec *current_spec;    struct order_spec *current_spec;
         struct oEGT eg0,eg1;    struct oEGT eg0,eg1;
         extern struct oEGT eg_le;    extern struct oEGT eg_le;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("invdalg : current_numberfield is not set");      error("invdalg : current_numberfield is not set");
         if ( !a )    if ( !a )
                 error("invdalg : division by 0");      error("invdalg : division by 0");
         else if ( NID(a) == N_Q ) {    else if ( NID(a) == N_Q ) {
                 invq((Q)a,&dn); *c = (DAlg)dn;      invq((Q)a,&dn); *c = (DAlg)dn;
                 return;      return 1;
         }    }
         dim = nf->dim;    dim = nf->dim;
         mb = nf->mb;    mb = nf->mb;
         n = nf->n;    n = nf->n;
         ln = ONEN;    ln = ONEN;
         dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);    dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
         MKDAlg(u,ONE,a0);    MKDAlg(u,ONE,a0);
         simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));    simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
         current_spec = dp_current_spec; initd(nf->spec);    current_spec = dp_current_spec; initd(nf->spec);
         for ( i = 0; i < dim; i++ ) {    for ( i = 0; i < dim; i++ ) {
                 m = mb[i];      m = mb[i];
                 for ( j = i-1; j >= 0; j-- )      for ( j = i-1; j >= 0; j-- )
                         if ( dp_redble(m,mb[j]) )        if ( dp_redble(m,mb[j]) )
                                 break;          break;
                 if ( j >= 0 ) {      if ( j >= 0 ) {
                         dp_subd(m,mb[j],&d);        dp_subd(m,mb[j],&d);
                         muld(CO,d,simp[j]->nm,&u);        muld(CO,d,simp[j]->nm,&u);
                         MKDAlg(u,simp[j]->dn,t);        MKDAlg(u,simp[j]->dn,t);
                         simpdalg(t,&simp[i]);        simpdalg(t,&simp[i]);
                 } else {      } else {
                         MKDAlg(m,ONE,t);        MKDAlg(m,ONE,t);
                         muldalg(t,a0,&simp[i]);        muldalg(t,a0,&simp[i]);
                 }      }
                 gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);      gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
                 muln(NM(simp[i]->dn),qn,&ln);      muln(NM(simp[i]->dn),qn,&ln);
         }    }
         initd(current_spec);    initd(current_spec);
         NTOQ(ln,1,dn);    NTOQ(ln,1,dn);
         MKMAT(mobj,dim,dim+1);    MKMAT(mobj,dim,dim+1);
         mat = (Q **)BDY(mobj);    mat = (Q **)BDY(mobj);
         mulq(dn,a->dn,&mat[0][dim]);    mulq(dn,a->dn,&mat[0][dim]);
         for ( j = 0; j < dim; j++ ) {    for ( j = 0; j < dim; j++ ) {
                 divq(dn,simp[j]->dn,&mul);      divq(dn,simp[j]->dn,&mul);
                 for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )      for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
                         if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {        if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
                                 mulq(mul,(Q)mp->c,&mat[i][j]);          mulq(mul,(Q)mp->c,&mat[i][j]);
                                 mp = NEXT(mp);          mp = NEXT(mp);
                         }        }
         }    }
         get_eg(&eg0);    get_eg(&eg0);
         rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);    rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
         get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);    get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
         if ( cinfo[0] == dim ) {    if ( cinfo[0] == dim ) {
                 /* the input is invertible */      /* the input is invertible */
                 solmat = (Q **)BDY(sol);      solmat = (Q **)BDY(sol);
                 for ( i = dim-1, mp0 = 0; i >= 0; i-- )      for ( i = dim-1, mp0 = 0; i >= 0; i-- )
                         if ( solmat[i][0] ) {        if ( solmat[i][0] ) {
                                 NEXTMP(mp0,mp);          NEXTMP(mp0,mp);
                                 mp->c = (P)solmat[i][0];          mp->c = (Obj)solmat[i][0];
                                 mp->dl = BDY(mb[i])->dl;          mp->dl = BDY(mb[i])->dl;
                         }        }
                 NEXT(mp) = 0; MKDP(n,mp0,u);      NEXT(mp) = 0; MKDP(n,mp0,u);
                 mulq(dnsol,nmc,&dn1);      mulq(dnsol,nmc,&dn1);
                 MKDAlg(u,dn1,r);      MKDAlg(u,dn1,r);
                 rmcontdalg(r,c);      rmcontdalg(r,c);
                 return 1;      return 1;
         } else    } else
                 return 0;      return 0;
 }  }
   
 NODE inv_or_split_dalg(DAlg a,DAlg *c)  NODE inv_or_split_dalg(DAlg a,DAlg *c)
 {  {
         NumberField nf;    NumberField nf;
         int dim,n,i,j,k,l;    int dim,n,i,j,k,l;
         DP *mb;    DP *mb;
         DP m,d,u;    DP m,d,u;
         N ln,gn,qn;    N ln,gn,qn;
         DAlg *simp;    DAlg *simp;
         DAlg t,a0,r;    DAlg t,a0,r;
         Q dn,dnsol,mul,nmc,dn1;    Q dn,dnsol,mul,nmc,dn1;
         MAT mobj,sol;    MAT mobj,sol;
         Q **mat,**solmat;    Q **mat,**solmat;
         MP mp0,mp;    MP mp0,mp;
         int *rinfo,*cinfo;    int *rinfo,*cinfo;
         int rank,nparam;    int rank,nparam;
         NODE nd0,nd,ndt;    NODE nd0,nd,ndt;
         struct order_spec *current_spec;    struct order_spec *current_spec;
         struct oEGT eg0,eg1;    struct oEGT eg0,eg1;
         extern struct oEGT eg_le;    extern struct oEGT eg_le;
     extern int DP_Print;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("invdalg : current_numberfield is not set");      error("invdalg : current_numberfield is not set");
         if ( !a )    if ( !a )
                 error("invdalg : division by 0");      error("invdalg : division by 0");
         else if ( NID(a) == N_Q ) {    else if ( NID(a) == N_Q ) {
                 invq((Q)a,&dn); *c = (DAlg)dn;      invq((Q)a,&dn); *c = (DAlg)dn;
                 return;      return 0;
         }    }
         dim = nf->dim;    dim = nf->dim;
         mb = nf->mb;    mb = nf->mb;
         n = nf->n;    n = nf->n;
         ln = ONEN;    ln = ONEN;
         dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);    dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
         MKDAlg(u,ONE,a0);    MKDAlg(u,ONE,a0);
         simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));    simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
         current_spec = dp_current_spec; initd(nf->spec);    current_spec = dp_current_spec; initd(nf->spec);
         for ( i = 0; i < dim; i++ ) {    for ( i = 0; i < dim; i++ ) {
                 m = mb[i];      if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
                 for ( j = i-1; j >= 0; j-- )      m = mb[i];
                         if ( dp_redble(m,mb[j]) )      for ( j = i-1; j >= 0; j-- )
                                 break;        if ( dp_redble(m,mb[j]) )
                 if ( j >= 0 ) {          break;
                         dp_subd(m,mb[j],&d);      if ( j >= 0 ) {
                         muld(CO,d,simp[j]->nm,&u);        dp_subd(m,mb[j],&d);
                         MKDAlg(u,simp[j]->dn,t);        if ( simp[j] ) {
                         simpdalg(t,&simp[i]);          muld(CO,d,simp[j]->nm,&u);
                 } else {          MKDAlg(u,simp[j]->dn,t);
                         MKDAlg(m,ONE,t);          simpdalg(t,&simp[i]);
                         muldalg(t,a0,&simp[i]);        } else
                 }          simp[i] = 0;
                 gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);      } else {
                 muln(NM(simp[i]->dn),qn,&ln);        MKDAlg(m,ONE,t);
         }        muldalg(t,a0,&simp[i]);
         initd(current_spec);      }
         NTOQ(ln,1,dn);      if ( simp[i] ) {
         MKMAT(mobj,dim,dim+1);        gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
         mat = (Q **)BDY(mobj);        muln(NM(simp[i]->dn),qn,&ln);
         mulq(dn,a->dn,&mat[0][dim]);      }
         for ( j = 0; j < dim; j++ ) {    }
                 divq(dn,simp[j]->dn,&mul);    initd(current_spec);
                 for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )    NTOQ(ln,1,dn);
                         if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {    MKMAT(mobj,dim,dim+1);
                                 mulq(mul,(Q)mp->c,&mat[i][j]);    mat = (Q **)BDY(mobj);
                                 mp = NEXT(mp);    mulq(dn,a->dn,&mat[0][dim]);
                         }    for ( j = 0; j < dim; j++ ) {
         }      if ( simp[j] ) {
         get_eg(&eg0);        divq(dn,simp[j]->dn,&mul);
         rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);        for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
         get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);          if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
         if ( cinfo[0] == dim ) {            mulq(mul,(Q)mp->c,&mat[i][j]);
                 /* the input is invertible */            mp = NEXT(mp);
                 solmat = (Q **)BDY(sol);          }
                 for ( i = dim-1, mp0 = 0; i >= 0; i-- )      }
                         if ( solmat[i][0] ) {    }
                                 NEXTMP(mp0,mp);    get_eg(&eg0);
                                 mp->c = (P)solmat[i][0];    rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
                                 mp->dl = BDY(mb[i])->dl;    get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
                         }    if ( cinfo[0] == dim ) {
                 NEXT(mp) = 0; MKDP(n,mp0,u);      /* the input is invertible */
                 mulq(dnsol,nmc,&dn1);      solmat = (Q **)BDY(sol);
                 MKDAlg(u,dn1,r);      for ( i = dim-1, mp0 = 0; i >= 0; i-- )
                 rmcontdalg(r,c);        if ( solmat[i][0] ) {
                 return 0;          NEXTMP(mp0,mp);
         } else {          mp->c = (Obj)solmat[i][0];
                 /* the input is not invertible */          mp->dl = BDY(mb[i])->dl;
                 nparam = (dim+1)-rank;        }
                 /* the index 'dim' should not be in cinfo[] */      NEXT(mp) = 0; MKDP(n,mp0,u);
                 solmat = (Q **)BDY(sol);      mulq(dnsol,nmc,&dn1);
                 for ( k = 0; k < nparam; k++ )      MKDAlg(u,dn1,r);
                         if ( cinfo[k] == dim )      rmcontdalg(r,c);
                                 error("invdalg : cannot happen");      return 0;
                 nd0 = 0;    } else {
                 for ( k = 0; k < nparam; k++ ) {      /* the input is not invertible */
                         m = mb[cinfo[k]];      nparam = sol->col;
                         for ( ndt = nd0; ndt; ndt = NEXT(ndt) ) {      solmat = (Q **)BDY(sol);
                                 if ( dp_redble(m,(DP)BDY(ndt)) ) break;      nd0 = 0;
                         }      for ( k = 0; k < nparam; k++ ) {
                         /* skip a redundunt basis element */        /* construct a new basis element */
                         if ( ndt ) continue;        m = mb[cinfo[k]];
                         /* construct a new basis element */        mp0 = 0;
                         mp0 = 0;        NEXTMP(mp0,mp);
                         NEXTMP(mp0,mp);        chsgnq(dnsol,&dn1); mp->c = (Obj)dn1;
                         chsgnq(dnsol,&dn1); mp->c = (P)dn1;        mp->dl = BDY(m)->dl;
                         mp->dl = BDY(m)->dl;        /* skip the last parameter */
                         /* skip the last parameter */        for ( l = rank-2; l >= 0; l-- ) {
                         for ( l = rank-2; l >= 0; l-- ) {          if ( solmat[l][k] ) {
                                 if ( solmat[l][k] ) {            NEXTMP(mp0,mp);
                                         NEXTMP(mp0,mp);            mp->c = (Obj)solmat[l][k];
                                         mp->c = (P)solmat[l][k];            mp->dl = BDY(mb[rinfo[l]])->dl;
                                         mp->dl = BDY(mb[rinfo[l]])->dl;          }
                                 }        }
                         }        NEXT(mp) = 0; MKDP(n,mp0,u);
                         NEXT(mp) = 0; MKDP(n,mp0,u);        NEXTNODE(nd0,nd);
                         NEXTNODE(nd0,nd);        BDY(nd) = (pointer)u;
                         BDY(nd) = (pointer)u;        NEXT(nd) = 0;
                         NEXT(nd) = 0;      }
                 }      NEXT(nd) = 0;
                 NEXT(nd) = 0;      return nd0;
                 return nd0;    }
         }  
 }  }
   
   NODE dp_inv_or_split(NODE gb,DP f,struct order_spec *spec, DP *inv)
   {
     int dim,n,i,j,k,l,nv;
     DP *mb,*ps;
     DP m,d,u,nm;
     N ln,gn,qn;
     DAlg *simp;
     DAlg a0,r;
     Q dn,dnsol,mul,nmc,dn1,iq;
     MAT mobj,sol;
     Q **mat,**solmat;
     MP mp0,mp;
     int *rinfo,*cinfo;
     int rank,nparam;
     NODE nd0,nd,ndt,ind,indt,t,mblist;
     struct oEGT eg0,eg1;
     extern struct oEGT eg_le;
     extern int DP_Print;
     initd(spec);
     dp_ptozp(f,&u); f = u;
   
     n = length(gb);
     ps = (DP *)MALLOC(n*sizeof(DP));
     for ( ind = 0, i = 0, t = gb; i < n; i++, t = NEXT(t) ) {
       ps[i] = (DP)BDY(t);
       NEXTNODE(ind,indt);
       STOQ(i,iq); BDY(indt) = iq;
     }
     if ( ind ) NEXT(indt) = 0;
     dp_true_nf(ind,f,ps,1,&nm,(P *)&dn);
     if ( !nm ) error("dp_inv_or_split : input is 0");
     f = nm;
   
     dp_mbase(gb,&mblist);
     dim = length(mblist);
     mb = (DP *)MALLOC(dim*sizeof(DP));
     for ( i = 0, t = mblist; i < dim; i++, t = NEXT(t) )
       mb[dim-i-1] = (DP)BDY(t);
     nv = mb[0]->nv;
     ln = ONEN;
     simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
     for ( i = 0; i < dim; i++ ) {
       if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
       m = mb[i];
       for ( j = i-1; j >= 0; j-- )
         if ( dp_redble(m,mb[j]) )
           break;
       if ( j >= 0 ) {
         dp_subd(m,mb[j],&d);
         if ( simp[j] ) {
           muld(CO,d,simp[j]->nm,&u);
           dp_true_nf(ind,u,ps,1,&nm,(P *)&dn);
           mulq(simp[j]->dn,dn,&dn1);
           MKDAlg(nm,dn1,simp[i]);
         } else
           simp[i] = 0;
       } else {
         dp_true_nf(ind,f,ps,1,&nm,(P *)&dn);
         MKDAlg(nm,dn,simp[i]);
       }
       if ( simp[i] ) {
         gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
         muln(NM(simp[i]->dn),qn,&ln);
       }
     }
     NTOQ(ln,1,dn);
     MKMAT(mobj,dim,dim+1);
     mat = (Q **)BDY(mobj);
     mat[0][dim] = dn;
     for ( j = 0; j < dim; j++ ) {
       if ( simp[j] ) {
         divq(dn,simp[j]->dn,&mul);
         for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
           if ( dl_equal(nv,BDY(mb[i])->dl,mp->dl) ) {
             mulq(mul,(Q)mp->c,&mat[i][j]);
             mp = NEXT(mp);
           }
       }
     }
     get_eg(&eg0);
     rank = generic_gauss_elim_hensel_dalg(mobj,mb,&sol,&dnsol,&rinfo,&cinfo);
     get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
     if ( cinfo[0] == dim ) {
       /* the input is invertible */
       solmat = (Q **)BDY(sol);
       for ( i = dim-1, mp0 = 0; i >= 0; i-- )
         if ( solmat[i][0] ) {
           NEXTMP(mp0,mp);
           mp->c = (Obj)solmat[i][0];
           mp->dl = BDY(mb[i])->dl;
         }
       NEXT(mp) = 0; MKDP(nv,mp0,*inv);
       return 0;
     } else {
       /* the input is not invertible */
       nparam = sol->col;
       solmat = (Q **)BDY(sol);
       nd0 = 0;
       for ( k = 0; k < nparam; k++ ) {
         /* construct a new basis element */
         m = mb[cinfo[k]];
         mp0 = 0;
         NEXTMP(mp0,mp);
         chsgnq(dnsol,&dn1); mp->c = (Obj)dn1;
         mp->dl = BDY(m)->dl;
         /* skip the last parameter */
         for ( l = rank-2; l >= 0; l-- ) {
           if ( solmat[l][k] ) {
             NEXTMP(mp0,mp);
             mp->c = (Obj)solmat[l][k];
             mp->dl = BDY(mb[rinfo[l]])->dl;
           }
         }
         NEXT(mp) = 0; MKDP(nv,mp0,u);
         NEXTNODE(nd0,nd);
         BDY(nd) = (pointer)u;
         NEXT(nd) = 0;
       }
       NEXT(nd) = 0;
       return nd0;
     }
   }
   
 void chsgndalg(DAlg a,DAlg *c)  void chsgndalg(DAlg a,DAlg *c)
 {  {
         DP nm;    DP nm;
         Q t;    Q t;
   
         if ( !a ) *c = 0;    if ( !a ) *c = 0;
         else if ( NID(a) == N_Q ) {    else if ( NID(a) == N_Q ) {
                 chsgnq((Q)a,&t); *c = (DAlg)t;      chsgnq((Q)a,&t); *c = (DAlg)t;
         } else {    } else {
                 chsgnd(a->nm,&nm);      chsgnd(a->nm,&nm);
                 MKDAlg(nm,a->dn,*c);      MKDAlg(nm,a->dn,*c);
         }    }
 }  }
   
 void pwrdalg(DAlg a,Q e,DAlg *c)  void pwrdalg(DAlg a,Q e,DAlg *c)
 {  {
         NumberField nf;    NumberField nf;
         DAlg t,z,y;    DAlg t,z,y;
         Q q;    Q q;
         N en,qn;    N en,qn;
         int r;    int r;
         int ret;    int ret;
   
         if ( !(nf=current_numberfield) )    if ( !(nf=current_numberfield) )
                 error("pwrdalg : current_numberfield is not set");      error("pwrdalg : current_numberfield is not set");
         if ( !a )    if ( !a )
                 *c = !e ? (DAlg)ONE : 0;      *c = !e ? (DAlg)ONE : 0;
         else if ( NID(a) == N_Q ) {    else if ( NID(a) == N_Q ) {
                 pwrq((Q)a,e,&q); *c = (DAlg)q;      pwrq((Q)a,e,&q); *c = (DAlg)q;
         } else if ( !e )    } else if ( !e )
                 *c = nf->one;      *c = nf->one;
         else if ( UNIQ(e) )    else if ( UNIQ(e) )
                 *c = a;      *c = a;
         else {    else {
                 if ( SGN(e) < 0 ) {      if ( SGN(e) < 0 ) {
                         ret = invdalg(a,&t);        ret = invdalg(a,&t);
                         if ( !ret )        if ( !ret )
                                 error("pwrdalg : the denominator is not invertible");          error("pwrdalg : the denominator is not invertible");
                         a = t;        a = t;
                 }      }
                 en = NM(e);      en = NM(e);
                 y = nf->one;      y = nf->one;
                 z = a;      z = a;
                 while ( 1 ) {      while ( 1 ) {
                         r = divin(en,2,&qn); en = qn;        r = divin(en,2,&qn); en = qn;
                         if ( r ) {        if ( r ) {
                                 muldalg(z,y,&t); y = t;          muldalg(z,y,&t); y = t;
                                 if ( !en ) {          if ( !en ) {
                                         *c = y;            *c = y;
                                         return;            return;
                                 }          }
                         }        }
                         muldalg(z,z,&t); z = t;        muldalg(z,z,&t); z = t;
                 }      }
         }    }
 }  }
   
 int cmpdalg(DAlg a,DAlg b)  int cmpdalg(DAlg a,DAlg b)
 {  {
         DAlg c;    DAlg c;
   
         subdalg(a,b,&c);    subdalg(a,b,&c);
         if ( !c ) return 0;    if ( !c ) return 0;
         else    else
                 return SGN((Q)BDY(c->nm)->c);      return SGN((Q)BDY(c->nm)->c);
 }  }
   
 /* convert da to a univariate poly; return the position of variable */  /* convert da to a univariate poly; return the position of variable */
   
 int dalgtoup(DAlg da,P *up,Q *dn)  int dalgtoup(DAlg da,P *up,Q *dn)
 {  {
         int nv,i,hi,current_d;    int nv,i,hi,current_d;
         DCP dc0,dc;    DCP dc0,dc;
         MP h,mp0,mp,t;    MP h,mp0,mp,t;
         DL hd,d;    DL hd,d;
         DP c;    DP c;
         DAlg cc;    DAlg cc;
         P v;    P v;
   
         nv = da->nm->nv;    nv = da->nm->nv;
         h = BDY(da->nm);    h = BDY(da->nm);
         *dn = da->dn;    *dn = da->dn;
         hd = h->dl;    hd = h->dl;
         for ( i = 0; i < nv; i++ )    for ( i = 0; i < nv; i++ )
                 if ( hd->d[i] ) break;      if ( hd->d[i] ) break;
         hi = i;    hi = i;
         current_d = hd->d[i];    current_d = hd->d[i];
         dc0 = 0;    dc0 = 0;
         mp0 = 0;    mp0 = 0;
         for ( t = h; t; t = NEXT(t) ) {    for ( t = h; t; t = NEXT(t) ) {
                 NEWDL(d,nv);      NEWDL(d,nv);
                 for ( i = 0; i <= hi; i++ ) d->d[i] = 0;      for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
                 for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];      for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
                 d->td = t->dl->td - t->dl->d[hi];      d->td = t->dl->td - t->dl->d[hi];
                 if ( t->dl->d[hi] != current_d ) {      if ( t->dl->d[hi] != current_d ) {
                         NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);        NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
                         NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;        NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
                         current_d = t->dl->d[hi];        current_d = t->dl->d[hi];
                         mp0 = 0;        mp0 = 0;
                 }      }
                 NEXTMP(mp0,mp);      NEXTMP(mp0,mp);
                 mp->c = t->c; mp->dl = d;      mp->c = t->c; mp->dl = d;
         }    }
         NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);    NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
         NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;    NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
         NEXT(dc) = 0;    NEXT(dc) = 0;
         makevar("x",&v);    makevar("x",&v);
         MKP(VR(v),dc0,*up);    MKP(VR(v),dc0,*up);
         return hi;    return hi;
 }  }
   
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.18
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>