[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.3 and 1.6

version 1.3, 2004/12/02 13:48:43 version 1.6, 2004/12/07 15:15:52
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.2 2004/12/02 08:39:54 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.5 2004/12/04 09:39:27 noro Exp $
 */  */
   
 #include "ca.h"  #include "ca.h"
Line 10  extern struct order_spec *dp_current_spec;
Line 10  extern struct order_spec *dp_current_spec;
 void simpdalg(DAlg da,DAlg *r);  void simpdalg(DAlg da,DAlg *r);
 void invdalg(DAlg a,DAlg *c);  void invdalg(DAlg a,DAlg *c);
 void rmcontdalg(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 58  void setfield_dalg(NODE alist)
Line 65  void setfield_dalg(NODE alist)
                 mb[i] = (DP)BDY(t);                  mb[i] = (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;
Line 139  void simpdalg(DAlg da,DAlg *r)
Line 348  void simpdalg(DAlg da,DAlg *r)
         DP nm;          DP nm;
         DAlg d;          DAlg d;
         Q dn,dn1;          Q dn,dn1;
           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");
Line 146  void simpdalg(DAlg da,DAlg *r)
Line 356  void simpdalg(DAlg da,DAlg *r)
                 *r = 0;                  *r = 0;
                 return;                  return;
         }          }
           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,&dn);
           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);
Line 157  void adddalg(DAlg a,DAlg b,DAlg *c)
Line 369  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;
         DP ta,tb,nm;          DP ta,tb,nm;
         struct order_spec *current_spec;          struct order_spec *current_spec;
   
Line 167  void adddalg(DAlg a,DAlg b,DAlg *c)
Line 380  void adddalg(DAlg a,DAlg b,DAlg *c)
         else if ( !b )          else if ( !b )
                 *c = a;                  *c = a;
         else {          else {
                   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);
Line 192  void subdalg(DAlg a,DAlg b,DAlg *c)
Line 406  void subdalg(DAlg a,DAlg b,DAlg *c)
         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;
         struct order_spec *current_spec;          struct order_spec *current_spec;
   
         if ( !(nf=current_numberfield) )          if ( !(nf=current_numberfield) )
Line 201  void subdalg(DAlg a,DAlg b,DAlg *c)
Line 416  void subdalg(DAlg a,DAlg b,DAlg *c)
         else if ( !b )          else if ( !b )
                 *c = a;                  *c = a;
         else {          else {
                   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);
Line 233  void muldalg(DAlg a,DAlg b,DAlg *c)
Line 449  void muldalg(DAlg a,DAlg b,DAlg *c)
         if ( !a || !b )          if ( !a || !b )
                 *c = 0;                  *c = 0;
         else {          else {
                   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);
Line 245  void muldalg(DAlg a,DAlg b,DAlg *c)
Line 462  void muldalg(DAlg a,DAlg b,DAlg *c)
   
 void divdalg(DAlg a,DAlg b,DAlg *c)  void divdalg(DAlg a,DAlg b,DAlg *c)
 {  {
         DAlg inv;          DAlg inv,t;
   
         if ( !current_numberfield )          if ( !current_numberfield )
                 error("divdalg : current_numberfield is not set");                  error("divdalg : current_numberfield is not set");
Line 254  void divdalg(DAlg a,DAlg b,DAlg *c)
Line 471  void divdalg(DAlg a,DAlg b,DAlg *c)
         if ( !a )          if ( !a )
                 c = 0;                  c = 0;
         else {          else {
                   qtodalg((Q)a,&t); a = t; qtodalg((Q)b,&t); b = t;
                 invdalg(b,&inv);                  invdalg(b,&inv);
                 muldalg(a,inv,c);                  muldalg(a,inv,c);
         }          }
Line 293  void invdalg(DAlg a,DAlg *c)
Line 511  void invdalg(DAlg a,DAlg *c)
         MP mp0,mp;          MP mp0,mp;
         int *rinfo,*cinfo;          int *rinfo,*cinfo;
         struct order_spec *current_spec;          struct order_spec *current_spec;
           struct oEGT eg0,eg1;
           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 ) {
                   invq((Q)a,&dn); *c = (DAlg)dn;
                   return;
           }
         dim = nf->dim;          dim = nf->dim;
         mb = nf->mb;          mb = nf->mb;
         n = nf->n;          n = nf->n;
Line 335  void invdalg(DAlg a,DAlg *c)
Line 559  void invdalg(DAlg a,DAlg *c)
                                 mp = NEXT(mp);                                  mp = NEXT(mp);
                         }                          }
         }          }
           get_eg(&eg0);
         generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);          generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
           get_eg(&eg1); add_eg(&eg_le,&eg0,&eg1);
         solmat = (Q **)BDY(sol);          solmat = (Q **)BDY(sol);
         for ( i = 0, mp0 = 0; i < dim; i++ )          for ( i = 0, mp0 = 0; i < dim; i++ )
                 if ( solmat[i][0] ) {                  if ( solmat[i][0] ) {
Line 351  void invdalg(DAlg a,DAlg *c)
Line 577  void invdalg(DAlg a,DAlg *c)
 void chsgndalg(DAlg a,DAlg *c)  void chsgndalg(DAlg a,DAlg *c)
 {  {
         DP nm;          DP nm;
           Q t;
   
         if ( !a ) *c = 0;          if ( !a ) *c = 0;
         else {          else if ( NID(a) == N_Q ) {
                   chsgnq((Q)a,&t); *c = (DAlg)t;
           } else {
                 chsgnd(a->nm,&nm);                  chsgnd(a->nm,&nm);
                 MKDAlg(nm,a->dn,*c);                  MKDAlg(nm,a->dn,*c);
         }          }
Line 363  void pwrdalg(DAlg a,Q e,DAlg *c)
Line 592  void pwrdalg(DAlg a,Q e,DAlg *c)
 {  {
         NumberField nf;          NumberField nf;
         DAlg t,z,y;          DAlg t,z,y;
           Q q;
         N en,qn;          N en,qn;
         int r;          int r;
   
         if ( !(nf=current_numberfield) )          if ( !(nf=current_numberfield) )
                 error("pwrdalg : current_numberfield is not set");                  error("pwrdalg : current_numberfield is not set");
         if ( !e )          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;                  *c = nf->one;
         else if ( UNIQ(e) )          else if ( UNIQ(e) )
                 *c = a;                  *c = a;
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.6
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>