[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.1 and 1.16

version 1.1, 2004/12/02 08:23:25 version 1.16, 2013/11/05 02:55:03
Line 1 
Line 1 
 /*  /*
  * $OpenXM$   * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.15 2007/02/13 07:12:54 noro Exp $
 */  */
   
 #include "ca.h"  #include "ca.h"
 #include "base.h"  #include "base.h"
   
 typedef struct oNumberField {  
         int n;  
         int dim;  
         VL vl;  
         P *defpoly;  
         DP *mb;  
         DP *ps;  
         NODE ind;  
         struct order_spec *spec;  
 } *NumberField;  
   
 typedef struct oDAlg {  
         short id;  
         char nid;  
         char pad;  
         DP nm;  
         Q dn;  
 } *DAlg;  
   
 #define N_DA 11  
 #define NEWDAlg(r) ((r)=(DAlg)MALLOC(sizeof(struct oDAlg)),OID(r)=O_N,NID(r)=N_DA)  
 #define MKDAlg(dp,dn,r) (NEWDAlg(r),(r)->nm = (dp),(r)->dn=(dn))  
   
 static NumberField current_numberfield;  static NumberField current_numberfield;
 extern struct order_spec *dp_current_spec;  extern struct order_spec *dp_current_spec;
   void simpdalg(DAlg da,DAlg *r);
   int invdalg(DAlg a,DAlg *c);
   void rmcontdalg(DAlg a, DAlg *c);
   void algtodalg(Alg a,DAlg *r);
   void dalgtoalg(DAlg da,Alg *r);
   
   NumberField get_numberfield()
   {
           return current_numberfield;
   }
   
 void setfield_dalg(NODE alist)  void setfield_dalg(NODE alist)
 {  {
         NumberField nf;          NumberField nf;
Line 41  void setfield_dalg(NODE alist)
Line 28  void setfield_dalg(NODE alist)
         P p;          P p;
         Q c,iq,two;          Q c,iq,two;
         DP *ps,*mb;          DP *ps,*mb;
           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;
   
Line 64  void setfield_dalg(NODE alist)
Line 52  void setfield_dalg(NODE alist)
                 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);
           MKDAlg(one,ONE,nf->one);
           nf->ind = b;
           dp_mbase(hlist,&mblist);
         initd(current_spec);          initd(current_spec);
           nf->dim = dim = length(mblist);
           nf->mb = mb = (DP *)MALLOC(dim*sizeof(DP));
           for ( i = 0, t = mblist; t; t = NEXT(t), i++ )
                   mb[dim-i-1] = (DP)BDY(t);
   }
   
   void setfield_gb(NODE gb,VL vl,struct order_spec *spec)
   {
           NumberField nf;
           VL vl1,vl2;
           int n,i,dim;
           Alg *gen;
           P *defpoly;
           P p;
           Q c,iq,two;
           DP *ps,*mb;
           DP one;
           NODE t,b,b1,b2,hlist,mblist;
           struct order_spec *current_spec;
   
           nf = (NumberField)MALLOC(sizeof(struct oNumberField));
           current_numberfield = nf;
           for ( vl1 = vl, n = 0; vl1; vl1 = NEXT(vl1), n++ );
           nf->n = n;
           nf->psn = length(gb);
           nf->vl = vl;
           nf->defpoly = defpoly = (P *)MALLOC(nf->psn*sizeof(P));
           nf->ps = ps = (DP *)MALLOC(nf->psn*sizeof(DP));
           current_spec = dp_current_spec;
           nf->spec = spec;
           initd(nf->spec);
           for ( b = hlist = 0, i = 0, t = gb; i < nf->psn; t = NEXT(t), i++ ) {
                   ptozp((P)BDY(t),1,&c,&defpoly[i]);
                   ptod(CO,vl,defpoly[i],&ps[i]);
                   STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
                   MKNODE(b2,(pointer)ps[i],hlist); hlist = b2;
           }
           ptod(ALG,vl,(P)ONE,&one);
           MKDAlg(one,ONE,nf->one);
         nf->ind = b;          nf->ind = b;
         dp_base(hlist,&mblist);          dp_mbase(hlist,&mblist);
           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[i] = (DP)BDY(mblist);                  mb[dim-i-1] = (DP)BDY(t);
 }  }
   
   void qtodalg(Q q,DAlg *r)
   {
           NumberField nf;
           Q t;
           DP nm;
   
           if ( !(nf=current_numberfield) )
                   error("qtodalg : current_numberfield is not set");
           if ( !q )
                   *r = 0;
           else if ( NID(q) == N_DA )
                   *r = (DAlg)q;
           else if ( NID(q) == N_Q ) {
                   if ( INT(q) ) {
                           muldc(CO,nf->one->nm,(P)q,&nm);
                           MKDAlg(nm,ONE,*r);
                   } else {
                           NTOQ(NM(q),SGN(q),t);
                           muldc(CO,nf->one->nm,(P)t,&nm);
                           NTOQ(DN(q),1,t);
                           MKDAlg(nm,t,*r);
                   }
           } else
                   error("qtodalg : invalid argument");
   }
   
   void obj_algtodalg(Obj obj,Obj *r)
   {
           DAlg d;
           DCP dc,dcr0,dcr;
           P c,p;
           Obj t;
           Obj nm,dn;
           NODE b,s,s0;
           R rat;
           VECT v;
           MAT mat;
           LIST list;
           pointer *a;
           pointer **m;
           int len,row,col,i,j,l;
   
           if ( !obj ) {
                   *r = 0;
                   return;
           }
           switch ( OID(obj) ) {
                   case O_N:
                           algtodalg((Alg)obj,&d); *r = (Obj)d;
                           break;
                   case O_P:
                           for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
                                   obj_algtodalg((Obj)COEF(dc),&t);
                                   if ( t ) {
                                           NEXTDC(dcr0,dcr);
                                           COEF(dcr) = (P)t;
                                           DEG(dcr) = DEG(dc);
                                   }
                           }
                           if ( dcr0 ) {
                                   MKP(VR((P)obj),dcr0,p);
                                   *r = (Obj)p;
                           } else
                                   *r = 0;
                           break;
                   case O_R:
                           obj_algtodalg((Obj)NM((R)obj),&nm);
                           obj_algtodalg((Obj)DN((R)obj),&dn);
                           if ( !dn )
                                   error("obj_algtodalg : division by 0");
                           if ( !nm )
                                   *r = 0;
                           else {
                                   MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
                           }
                           break;
                   case O_LIST:
                           s0 = 0;
                           for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
                                   NEXTNODE(s0,s);
                                   obj_algtodalg((Obj)BDY(b),&t);
                                   BDY(s) = (pointer)t;
                           }
                           NEXT(s) = 0;
                           MKLIST(list,s0);
                           *r = (Obj)list;
                           break;
                   case O_VECT:
                           l = ((VECT)obj)->len;
                           a = BDY((VECT)obj);
                           MKVECT(v,l);
                           for ( i = 0; i < l; i++ ) {
                                   obj_algtodalg((Obj)a[i],&t);
                                   BDY(v)[i] = (pointer)t;
                           }
                           *r = (Obj)v;
                           break;
                   case O_MAT:
                           row = ((MAT)obj)->row; col = ((MAT)obj)->col;
                           m = BDY((MAT)obj);
                           MKMAT(mat,row,col);
                           for ( i = 0; i < row; i++ )
                                   for ( j = 0; j < col; j++ ) {
                                           obj_algtodalg((Obj)m[i][j],&t);
                                           BDY(mat)[i][j] = (pointer)t;
                                   }
                           *r = (Obj)mat;
                           break;
                   default:
                           *r = obj;
                           break;
           }
   }
   
   void obj_dalgtoalg(Obj obj,Obj *r)
   {
           Alg d;
           DCP dc,dcr0,dcr;
           P c,p;
           Obj t;
           Obj nm,dn;
           NODE b,s,s0;
           R rat;
           VECT v;
           MAT mat;
           LIST list;
           pointer *a;
           pointer **m;
           int len,row,col,i,j,l;
   
           if ( !obj ) {
                   *r = 0;
                   return;
           }
           switch ( OID(obj) ) {
                   case O_N:
                           dalgtoalg((DAlg)obj,&d); *r = (Obj)d;
                           break;
                   case O_P:
                           for ( dcr0 = 0, dc = DC((P)obj); dc; dc = NEXT(dc) ) {
                                   obj_dalgtoalg((Obj)COEF(dc),&t);
                                   if ( t ) {
                                           NEXTDC(dcr0,dcr);
                                           COEF(dcr) = (P)t;
                                           DEG(dcr) = DEG(dc);
                                   }
                           }
                           if ( dcr0 ) {
                                   MKP(VR((P)obj),dcr0,p);
                                   *r = (Obj)p;
                           } else
                                   *r = 0;
                           break;
                   case O_R:
                           obj_dalgtoalg((Obj)NM((R)obj),&nm);
                           obj_dalgtoalg((Obj)DN((R)obj),&dn);
                           if ( !dn )
                                   error("obj_dalgtoalg : division by 0");
                           if ( !nm )
                                   *r = 0;
                           else {
                                   MKRAT((P)nm,(P)dn,0,rat); *r = (Obj)rat;
                           }
                           break;
                   case O_LIST:
                           s0 = 0;
                           for ( b = BDY((LIST)obj); b; b = NEXT(b) ) {
                                   NEXTNODE(s0,s);
                                   obj_dalgtoalg((Obj)BDY(b),&t);
                                   BDY(s) = (pointer)t;
                           }
                           NEXT(s) = 0;
                           MKLIST(list,s0);
                           *r = (Obj)list;
                           break;
                   case O_VECT:
                           l = ((VECT)obj)->len;
                           a = BDY((VECT)obj);
                           MKVECT(v,l);
                           for ( i = 0; i < l; i++ ) {
                                   obj_dalgtoalg((Obj)a[i],&t);
                                   BDY(v)[i] = (pointer)t;
                           }
                           *r = (Obj)v;
                           break;
                   case O_MAT:
                           row = ((MAT)obj)->row; col = ((MAT)obj)->col;
                           m = BDY((MAT)obj);
                           MKMAT(mat,row,col);
                           for ( i = 0; i < row; i++ )
                                   for ( j = 0; j < col; j++ ) {
                                           obj_dalgtoalg((Obj)m[i][j],&t);
                                           BDY(mat)[i][j] = (pointer)t;
                                   }
                           *r = (Obj)mat;
                           break;
                   default:
                           *r = obj;
                           break;
           }
   }
   
 void algtodalg(Alg a,DAlg *r)  void algtodalg(Alg a,DAlg *r)
 {  {
         P ap,p,p1;          P ap,p,p1;
         Q c,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;
           V v;
   
         if ( !(nf=current_numberfield) )          if ( !(nf=current_numberfield) )
                 error("algtodalg : current_numberfield is not set");                  error("algtodalg : current_numberfield is not set");
         ap = (P)BDY(a);          if ( !a ) {
         ptozp(ap,1,&c,&p);                  *r = 0;
         if ( INT(c) ) {                  return;
                 p = ap;  
                 dn = ONE;  
         } else {  
                 NTOQ(NM(c),SGN(c),nm);  
                 NTOQ(DN(c),1,dn);  
                 mulpq(p,(P)nm,&p1); p = p1;  
         }          }
         current_spec = dp_current_spec;          switch (NID((Num)a) ) {
         initd(nf->spec);                  case N_Q:
         ptod(ALG,nf->vl,p,&dp);                          c = (Q)a;
         MKDAlg(dp,dn,da);                          if ( INT(c) ) {
         *r = da;                                  muldc(CO,nf->one->nm,(P)c,&dp);
                                   MKDAlg(dp,ONE,*r);
                           } else {
                                   NTOQ(NM(c),SGN(c),c1);
                                   NTOQ(DN(c),1,d1);
                                   muldc(CO,nf->one->nm,(P)c1,&dp);
                                   MKDAlg(dp,d1,*r);
                           }
                           break;
                   case N_A:
                           ap = (P)BDY(a);
                           ptozp(ap,1,&c,&p);
                           if ( INT(c) ) {
                                   p = ap;
                                   dn = ONE;
                           } else {
                                   NTOQ(NM(c),SGN(c),nm);
                                   NTOQ(DN(c),1,dn);
                                   mulpq(p,(P)nm,&p1); p = p1;
                           }
                           current_spec = dp_current_spec; initd(nf->spec);
                           get_vars(p,&vl);
                           for ( tvl = vl; tvl; tvl = NEXT(tvl) ) {
                                   v = tvl->v;
                                   for ( svl = nf->vl; svl; svl = NEXT(svl) )
                                           if ( v == svl->v )
                                                   break;
                                   if ( !svl )
                                           error("algtodalg : incompatible numberfield");
                           }
                           ptod(ALG,nf->vl,p,&dp);
                           MKDAlg(dp,dn,da);
                           simpdalg(da,r);
                           break;
                   default:
                           error("algtodalg : invalid argument");
                           break;
           }
 }  }
   
 void dalgtoalg(DAlg da,Num *a)  void dalgtoalg(DAlg da,Alg *r)
 {  {
         if ( !current_numberfield )          NumberField nf;
                 error("algtodalg : current_numberfield is not set");          P p,p1;
           Q inv;
   
           if ( !(nf=current_numberfield) )
                   error("dalgtoalg : current_numberfield is not set");
           if ( !da ) *r = 0;
           else {
                   dtop(ALG,nf->vl,da->nm,&p);
                   invq(da->dn,&inv);
                   mulpq(p,(P)inv,&p1);
                   MKAlg(p1,*r);
           }
 }  }
   
 void simpdalg(DAlg da,DAlg *r)  void simpdalg(DAlg da,DAlg *r)
 {  {
         if ( !current_numberfield )          NumberField nf;
                 error("algtodalg : current_numberfield is not set");          DP nm;
           DAlg d;
           Q dn,dn1;
           struct order_spec *current_spec;
   
           if ( !(nf=current_numberfield) )
                   error("simpdalg : current_numberfield is not set");
           if ( !da ) {
                   *r = 0;
                   return;
           }
           current_spec = dp_current_spec; initd(nf->spec);
           dp_true_nf(nf->ind,da->nm,nf->ps,1,&nm,&dn);
           if ( !nm ) *r = 0;
           else {
                   initd(current_spec);
                   mulq(da->dn,dn,&dn1);
                   MKDAlg(nm,dn1,d);
                   rmcontdalg(d,r);
           }
 }  }
   
 void adddalg(DAlg a,DAlg b,DAlg *c)  void adddalg(DAlg a,DAlg b,DAlg *c)
 {  {
         if ( !current_numberfield )          NumberField nf;
                 error("algtodalg : current_numberfield is not set");          Q dna,dnb,a1,b1,dn,g;
           N an,bn,gn;
           DAlg t;
           DP ta,tb,nm;
           struct order_spec *current_spec;
   
           if ( !(nf=current_numberfield) )
                   error("adddalg : current_numberfield is not set");
           if ( !a )
                   *c = b;
           else if ( !b )
                   *c = a;
           else {
                   qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                   dna = a->dn;
                   dnb = b->dn;
                   gcdn(NM(dna),NM(dnb),&gn);
                   divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
                   NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),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);
                   current_spec = dp_current_spec; initd(nf->spec);
                   addd(CO,ta,tb,&nm);
                   initd(current_spec);
                   if ( !nm )
                           *c = 0;
                   else {
                           mulq(dna,b1,&dn);
                           MKDAlg(nm,dn,*c);
                   }
           }
 }  }
   
 void subdalg(DAlg a,DAlg b,DAlg *c)  void subdalg(DAlg a,DAlg b,DAlg *c)
 {  {
         if ( !current_numberfield )          NumberField nf;
                 error("algtodalg : current_numberfield is not set");          Q dna,dnb,a1,b1,dn,g;
           N an,bn,gn;
           DP ta,tb,nm;
           DAlg t;
           struct order_spec *current_spec;
   
           if ( !(nf=current_numberfield) )
                   error("subdalg : current_numberfield is not set");
           if ( !a )
                   *c = b;
           else if ( !b )
                   *c = a;
           else {
                   qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                   dna = a->dn;
                   dnb = b->dn;
                   gcdn(NM(dna),NM(dnb),&gn);
                   divsn(NM(dna),gn,&an); divsn(NM(dnb),gn,&bn);
                   NTOQ(an,SGN(dna),a1); NTOQ(bn,SGN(dnb),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);
                   current_spec = dp_current_spec; initd(nf->spec);
                   subd(CO,ta,tb,&nm);
                   initd(current_spec);
                   if ( !nm )
                           *c = 0;
                   else {
                           mulq(dna,b1,&dn);
                           MKDAlg(nm,dn,*c);
                   }
           }
 }  }
   
   void muldalg(DAlg a,DAlg b,DAlg *c)
   {
           NumberField nf;
           DP nm;
           Q dn;
           DAlg t;
           struct order_spec *current_spec;
   
           if ( !(nf=current_numberfield) )
                   error("muldalg : current_numberfield is not set");
           if ( !a || !b )
                   *c = 0;
           else {
                   qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                   current_spec = dp_current_spec; initd(nf->spec);
                   muld(CO,a->nm,b->nm,&nm);
                   initd(current_spec);
                   mulq(a->dn,b->dn,&dn);
                   MKDAlg(nm,dn,t);
                   simpdalg(t,c);
           }
   }
   
   
 void divdalg(DAlg a,DAlg b,DAlg *c)  void divdalg(DAlg a,DAlg b,DAlg *c)
 {  {
           DAlg inv,t;
           int ret;
   
         if ( !current_numberfield )          if ( !current_numberfield )
                 error("algtodalg : current_numberfield is not set");                  error("divdalg : current_numberfield is not set");
           if ( !b )
                   error("divdalg : division by 0");
           if ( !a )
                   c = 0;
           else {
                   qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                   ret = invdalg(b,&inv);
                   if ( !ret ) {
                           error("divdalg : the denominator is not invertible");
                   }
                   muldalg(a,inv,c);
           }
 }  }
   
 void invdalg(DAlg a,DAlg *c)  void rmcontdalg(DAlg a, DAlg *r)
 {  {
         if ( !current_numberfield )          DP u,u1;
                 error("algtodalg : current_numberfield is not set");          Q cont,c,d;
           N gn,cn,dn;
   
           if ( !a )
                   *r = a;
           else {
                   dp_ptozp(a->nm,&u);
                   divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&cont);
                   gcdn(NM(cont),NM(a->dn),&gn);
                   divsn(NM(cont),gn,&cn); NTOQ(cn,SGN(cont),c);
                   divsn(NM(a->dn),gn,&dn); NTOQ(dn,SGN(a->dn),d);
                   muldc(CO,u,(P)c,&u1);
                   MKDAlg(u1,d,*r);
           }
 }  }
   
   int invdalg(DAlg a,DAlg *c)
   {
           NumberField nf;
           int dim,n,i,j,k,l;
           DP *mb;
           DP m,d,u;
           N ln,gn,qn;
           DAlg *simp;
           DAlg t,a0,r;
           Q dn,dnsol,mul,nmc,dn1;
           MAT mobj,sol;
           Q **mat,**solmat;
           MP mp0,mp;
           int *rinfo,*cinfo;
           int rank,nparam;
           NODE nd0,nd,ndt;
           struct order_spec *current_spec;
           struct oEGT eg0,eg1;
           extern struct oEGT eg_le;
   
           if ( !(nf=current_numberfield) )
                   error("invdalg : current_numberfield is not set");
           if ( !a )
                   error("invdalg : division by 0");
           else if ( NID(a) == N_Q ) {
                   invq((Q)a,&dn); *c = (DAlg)dn;
                   return 1;
           }
           dim = nf->dim;
           mb = nf->mb;
           n = nf->n;
           ln = ONEN;
           dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
           MKDAlg(u,ONE,a0);
           simp = (DAlg *)ALLOCA(dim*sizeof(DAlg));
           current_spec = dp_current_spec; initd(nf->spec);
           for ( i = 0; i < dim; i++ ) {
                   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);
                           muld(CO,d,simp[j]->nm,&u);
                           MKDAlg(u,simp[j]->dn,t);
                           simpdalg(t,&simp[i]);
                   } else {
                           MKDAlg(m,ONE,t);
                           muldalg(t,a0,&simp[i]);
                   }
                   gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
                   muln(NM(simp[i]->dn),qn,&ln);
           }
           initd(current_spec);
           NTOQ(ln,1,dn);
           MKMAT(mobj,dim,dim+1);
           mat = (Q **)BDY(mobj);
           mulq(dn,a->dn,&mat[0][dim]);
           for ( j = 0; j < dim; j++ ) {
                   divq(dn,simp[j]->dn,&mul);
                   for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
                           if ( dl_equal(n,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(mobj,&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 = (P)solmat[i][0];
                                   mp->dl = BDY(mb[i])->dl;
                           }
                   NEXT(mp) = 0; MKDP(n,mp0,u);
                   mulq(dnsol,nmc,&dn1);
                   MKDAlg(u,dn1,r);
                   rmcontdalg(r,c);
                   return 1;
           } else
                   return 0;
   }
   
   NODE inv_or_split_dalg(DAlg a,DAlg *c)
   {
           NumberField nf;
           int dim,n,i,j,k,l;
           DP *mb;
           DP m,d,u;
           N ln,gn,qn;
           DAlg *simp;
           DAlg t,a0,r;
           Q dn,dnsol,mul,nmc,dn1;
           MAT mobj,sol;
           Q **mat,**solmat;
           MP mp0,mp;
           int *rinfo,*cinfo;
           int rank,nparam;
           NODE nd0,nd,ndt;
           struct order_spec *current_spec;
           struct oEGT eg0,eg1;
           extern struct oEGT eg_le;
           extern int DP_Print;
   
           if ( !(nf=current_numberfield) )
                   error("invdalg : current_numberfield is not set");
           if ( !a )
                   error("invdalg : division by 0");
           else if ( NID(a) == N_Q ) {
                   invq((Q)a,&dn); *c = (DAlg)dn;
                   return 1;
           }
           dim = nf->dim;
           mb = nf->mb;
           n = nf->n;
           ln = ONEN;
           dp_ptozp(a->nm,&u); divq((Q)BDY(a->nm)->c,(Q)BDY(u)->c,&nmc);
           MKDAlg(u,ONE,a0);
           simp = (DAlg *)MALLOC(dim*sizeof(DAlg));
           current_spec = dp_current_spec; initd(nf->spec);
           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);
                                   MKDAlg(u,simp[j]->dn,t);
                                   simpdalg(t,&simp[i]);
                           } else
                                   simp[i] = 0;
                   } else {
                           MKDAlg(m,ONE,t);
                           muldalg(t,a0,&simp[i]);
                   }
                   if ( simp[i] ) {
                           gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
                           muln(NM(simp[i]->dn),qn,&ln);
                   }
           }
           initd(current_spec);
           NTOQ(ln,1,dn);
           MKMAT(mobj,dim,dim+1);
           mat = (Q **)BDY(mobj);
           mulq(dn,a->dn,&mat[0][dim]);
           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(n,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 = (P)solmat[i][0];
                                   mp->dl = BDY(mb[i])->dl;
                           }
                   NEXT(mp) = 0; MKDP(n,mp0,u);
                   mulq(dnsol,nmc,&dn1);
                   MKDAlg(u,dn1,r);
                   rmcontdalg(r,c);
                   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 = (P)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 = (P)solmat[l][k];
                                           mp->dl = BDY(mb[rinfo[l]])->dl;
                                   }
                           }
                           NEXT(mp) = 0; MKDP(n,mp0,u);
                           NEXTNODE(nd0,nd);
                           BDY(nd) = (pointer)u;
                           NEXT(nd) = 0;
                   }
                   NEXT(nd) = 0;
                   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,&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,&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,&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 = (P)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 = (P)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 = (P)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;
           Q t;
   
           if ( !a ) *c = 0;
           else if ( NID(a) == N_Q ) {
                   chsgnq((Q)a,&t); *c = (DAlg)t;
           } else {
                   chsgnd(a->nm,&nm);
                   MKDAlg(nm,a->dn,*c);
           }
 }  }
   
 void pwrdalg(DAlg a,Q b,DAlg *c)  void pwrdalg(DAlg a,Q e,DAlg *c)
 {  {
         if ( !current_numberfield )          NumberField nf;
                 error("algtodalg : current_numberfield is not set");          DAlg t,z,y;
           Q q;
           N en,qn;
           int r;
           int ret;
   
           if ( !(nf=current_numberfield) )
                   error("pwrdalg : current_numberfield is not set");
           if ( !a )
                   *c = !e ? (DAlg)ONE : 0;
           else if ( NID(a) == N_Q ) {
                   pwrq((Q)a,e,&q); *c = (DAlg)q;
           } else if ( !e )
                   *c = nf->one;
           else if ( UNIQ(e) )
                   *c = a;
           else {
                   if ( SGN(e) < 0 ) {
                           ret = invdalg(a,&t);
                           if ( !ret )
                                   error("pwrdalg : the denominator is not invertible");
                           a = t;
                   }
                   en = NM(e);
                   y = nf->one;
                   z = a;
                   while ( 1 ) {
                           r = divin(en,2,&qn); en = qn;
                           if ( r ) {
                                   muldalg(z,y,&t); y = t;
                                   if ( !en ) {
                                           *c = y;
                                           return;
                                   }
                           }
                           muldalg(z,z,&t); z = t;
                   }
           }
 }  }
   
 int cmpgdalg(DAlg a,DAlg b)  int cmpdalg(DAlg a,DAlg b)
 {  {
         if ( !current_numberfield )          DAlg c;
                 error("algtodalg : current_numberfield is not set");  
           subdalg(a,b,&c);
           if ( !c ) return 0;
           else
                   return SGN((Q)BDY(c->nm)->c);
 }  }
   
   /* convert da to a univariate poly; return the position of variable */
   
   int dalgtoup(DAlg da,P *up,Q *dn)
   {
           int nv,i,hi,current_d;
           DCP dc0,dc;
           MP h,mp0,mp,t;
           DL hd,d;
           DP c;
           DAlg cc;
           P v;
   
           nv = da->nm->nv;
           h = BDY(da->nm);
           *dn = da->dn;
           hd = h->dl;
           for ( i = 0; i < nv; i++ )
                   if ( hd->d[i] ) break;
           hi = i;
           current_d = hd->d[i];
           dc0 = 0;
           mp0 = 0;
           for ( t = h; t; t = NEXT(t) ) {
                   NEWDL(d,nv);
                   for ( i = 0; i <= hi; i++ ) d->d[i] = 0;
                   for ( ; i < nv; i++ ) d->d[i] = t->dl->d[i];
                   d->td = t->dl->td - t->dl->d[hi];
                   if ( t->dl->d[hi] != current_d ) {
                           NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
                           NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
                           current_d = t->dl->d[hi];
                           mp0 = 0;
                   }
                   NEXTMP(mp0,mp);
                   mp->c = t->c; mp->dl = d;
           }
           NEXT(mp) = 0; MKDP(nv,mp0,c); MKDAlg(c,ONE,cc);
           NEXTDC(dc0,dc); STOQ(current_d,DEG(dc)); COEF(dc) = (P)cc;
           NEXT(dc) = 0;
           makevar("x",&v);
           MKP(VR(v),dc0,*up);
           return hi;
   }
   
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.16
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>