[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.7 and 1.8

version 1.7, 2004/12/10 02:45:02 version 1.8, 2005/01/23 14:03:47
Line 1 
Line 1 
 /*  /*
  * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.6 2004/12/07 15:15:52 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/engine/dalg.c,v 1.7 2004/12/10 02:45:02 noro Exp $
 */  */
   
 #include "ca.h"  #include "ca.h"
Line 8 
Line 8 
 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);  void simpdalg(DAlg da,DAlg *r);
 void invdalg(DAlg a,DAlg *c);  int invdalg(DAlg a,DAlg *c);
 void rmcontdalg(DAlg a, DAlg *c);  void rmcontdalg(DAlg a, DAlg *c);
 void algtodalg(Alg a,DAlg *r);  void algtodalg(Alg a,DAlg *r);
 void dalgtoalg(DAlg da,Alg *r);  void dalgtoalg(DAlg da,Alg *r);
Line 62  void setfield_dalg(NODE alist)
Line 62  void setfield_dalg(NODE alist)
         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(t);                  mb[dim-i-1] = (DP)BDY(t);
 }  }
   
 void qtodalg(Q q,DAlg *r)  void qtodalg(Q q,DAlg *r)
Line 463  void muldalg(DAlg a,DAlg b,DAlg *c)
Line 463  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,t;          DAlg inv,t;
           int ret;
   
         if ( !current_numberfield )          if ( !current_numberfield )
                 error("divdalg : current_numberfield is not set");                  error("divdalg : current_numberfield is not set");
Line 472  void divdalg(DAlg a,DAlg b,DAlg *c)
Line 473  void divdalg(DAlg a,DAlg b,DAlg *c)
                 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;
                 invdalg(b,&inv);                  ret = invdalg(b,&inv);
                   if ( !ret ) {
                           error("divdalg : the denominator is not invertible");
                   }
                 muldalg(a,inv,c);                  muldalg(a,inv,c);
         }          }
 }  }
Line 496  void rmcontdalg(DAlg a, DAlg *r)
Line 500  void rmcontdalg(DAlg a, DAlg *r)
         }          }
 }  }
   
 void invdalg(DAlg a,DAlg *c)  int invdalg(DAlg a,DAlg *c)
 {  {
         NumberField nf;          NumberField nf;
         int dim,n,i,j;          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;
Line 510  void invdalg(DAlg a,DAlg *c)
Line 514  void invdalg(DAlg a,DAlg *c)
         Q **mat,**solmat;          Q **mat,**solmat;
         MP mp0,mp;          MP mp0,mp;
         int *rinfo,*cinfo;          int *rinfo,*cinfo;
           int rank,nparam;
           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;
Line 530  void invdalg(DAlg a,DAlg *c)
Line 536  void invdalg(DAlg a,DAlg *c)
         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 = dim-1; i >= 0; i-- ) {          for ( i = 0; i < dim; i++ ) {
                 m = mb[i];                  m = mb[i];
                 for ( j = i+1; j < dim; j++ )                  for ( j = i-1; j >= 0; j-- )
                         if ( dp_redble(m,mb[j]) )                          if ( dp_redble(m,mb[j]) )
                                 break;                                  break;
                 if ( j < dim ) {                  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);
Line 551  void invdalg(DAlg a,DAlg *c)
Line 557  void invdalg(DAlg a,DAlg *c)
         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[dim-1][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 = 0, mp = BDY(simp[j]->nm); mp && i < dim; 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);
         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);
         solmat = (Q **)BDY(sol);          if ( cinfo[0] == dim ) {
         for ( i = 0, mp0 = 0; i < dim; i++ )                  /* the input is invertible */
                 if ( solmat[i][0] ) {                  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;
   
           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;
           }
           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 0;
           } else {
                   /* the input is not invertible */
                   nparam = (dim+1)-rank;
                   /* the index 'dim' should not be in cinfo[] */
                   solmat = (Q **)BDY(sol);
                   for ( k = 0; k < nparam; k++ )
                           if ( cinfo[k] == dim )
                                   error("invdalg : cannot happen");
                   nd0 = 0;
                   for ( k = 0; k < nparam; k++ ) {
                           m = mb[cinfo[k]];
                           for ( ndt = nd0; ndt; ndt = NEXT(ndt) ) {
                                   if ( dp_redble(m,(DP)BDY(ndt)) ) break;
                           }
                           /* skip a redundunt basis element */
                           if ( ndt ) continue;
                           /* construct a new basis element */
                           mp0 = 0;
                         NEXTMP(mp0,mp);                          NEXTMP(mp0,mp);
                         mp->c = (P)solmat[i][0];                          chsgnq(dnsol,&dn1); mp->c = (P)dn1;
                         mp->dl = BDY(mb[i])->dl;                          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(mp) = 0; MKDP(n,mp0,u);                  NEXT(nd) = 0;
         mulq(dnsol,nmc,&dn1);                  return nd0;
         MKDAlg(u,dn1,r);          }
         rmcontdalg(r,c);  
 }  }
   
 void chsgndalg(DAlg a,DAlg *c)  void chsgndalg(DAlg a,DAlg *c)
Line 597  void pwrdalg(DAlg a,Q e,DAlg *c)
Line 730  void pwrdalg(DAlg a,Q e,DAlg *c)
         Q q;          Q q;
         N en,qn;          N en,qn;
         int r;          int r;
           int ret;
   
         if ( !(nf=current_numberfield) )          if ( !(nf=current_numberfield) )
                 error("pwrdalg : current_numberfield is not set");                  error("pwrdalg : current_numberfield is not set");
Line 610  void pwrdalg(DAlg a,Q e,DAlg *c)
Line 744  void pwrdalg(DAlg a,Q e,DAlg *c)
                 *c = a;                  *c = a;
         else {          else {
                 if ( SGN(e) < 0 ) {                  if ( SGN(e) < 0 ) {
                         invdalg(a,&t); a = t;                          ret = invdalg(a,&t);
                           if ( !ret )
                                   error("pwrdalg : the denominator is not invertible");
                           a = t;
                 }                  }
                 en = NM(e);                  en = NM(e);
                 y = nf->one;                  y = nf->one;
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.8
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>