[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.16

version 1.10, 2005/08/02 07:16:42 version 1.16, 2013/11/05 02:55:03
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.15 2007/02/13 07:12:54 noro Exp $
 */  */
   
 #include "ca.h"  #include "ca.h"
Line 378  void dalgtoalg(DAlg da,Alg *r)
Line 378  void dalgtoalg(DAlg da,Alg *r)
   
         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,&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)
Line 571  int invdalg(DAlg a,DAlg *c)
Line 574  int invdalg(DAlg a,DAlg *c)
                 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;
Line 651  NODE inv_or_split_dalg(DAlg a,DAlg *c)
Line 654  NODE inv_or_split_dalg(DAlg a,DAlg *c)
         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");
Line 658  NODE inv_or_split_dalg(DAlg a,DAlg *c)
Line 662  NODE inv_or_split_dalg(DAlg a,DAlg *c)
                 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;
Line 666  NODE inv_or_split_dalg(DAlg a,DAlg *c)
Line 670  NODE inv_or_split_dalg(DAlg a,DAlg *c)
         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++ ) {
                   if ( DP_Print ) { fprintf(asir_out,"."); fflush(asir_out); }
                 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);                          if ( simp[j] ) {
                         MKDAlg(u,simp[j]->dn,t);                                  muld(CO,d,simp[j]->nm,&u);
                         simpdalg(t,&simp[i]);                                  MKDAlg(u,simp[j]->dn,t);
                                   simpdalg(t,&simp[i]);
                           } else
                                   simp[i] = 0;
                 } 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);                  if ( simp[i] ) {
                 muln(NM(simp[i]->dn),qn,&ln);                          gcdn(NM(simp[i]->dn),ln,&gn); divsn(ln,gn,&qn);
                           muln(NM(simp[i]->dn),qn,&ln);
                   }
         }          }
         initd(current_spec);          initd(current_spec);
         NTOQ(ln,1,dn);          NTOQ(ln,1,dn);
Line 691  NODE inv_or_split_dalg(DAlg a,DAlg *c)
Line 701  NODE inv_or_split_dalg(DAlg a,DAlg *c)
         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);                  if ( simp[j] ) {
                 for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )                          divq(dn,simp[j]->dn,&mul);
                         if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {                          for ( i = dim-1, mp = BDY(simp[j]->nm); mp && i >= 0; i-- )
                                 mulq(mul,(Q)mp->c,&mat[i][j]);                                  if ( dl_equal(n,BDY(mb[i])->dl,mp->dl) ) {
                                 mp = NEXT(mp);                                          mulq(mul,(Q)mp->c,&mat[i][j]);
                         }                                          mp = NEXT(mp);
                                   }
                   }
         }          }
         get_eg(&eg0);          get_eg(&eg0);
         rank = generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);          rank = generic_gauss_elim_hensel_dalg(mobj,mb,&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 */
Line 717  NODE inv_or_split_dalg(DAlg a,DAlg *c)
Line 729  NODE inv_or_split_dalg(DAlg a,DAlg *c)
                 return 0;                  return 0;
         } else {          } else {
                 /* the input is not invertible */                  /* the input is not invertible */
                 nparam = (dim+1)-rank;                  nparam = sol->col;
                 /* the index 'dim' should not be in cinfo[] */  
                 solmat = (Q **)BDY(sol);                  solmat = (Q **)BDY(sol);
                 for ( k = 0; k < nparam; k++ )  
                         if ( cinfo[k] == dim )  
                                 error("invdalg : cannot happen");  
                 nd0 = 0;                  nd0 = 0;
                 for ( k = 0; k < nparam; k++ ) {                  for ( k = 0; k < nparam; k++ ) {
                           /* construct a new basis element */
                         m = mb[cinfo[k]];                          m = mb[cinfo[k]];
                         for ( ndt = nd0; ndt; ndt = NEXT(ndt) ) {                          mp0 = 0;
                                 if ( dp_redble(m,(DP)BDY(ndt)) ) break;                          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;
                                   }
                         }                          }
                         /* skip a redundunt basis element */                          NEXT(mp) = 0; MKDP(n,mp0,u);
                         if ( ndt ) continue;                          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 */                          /* construct a new basis element */
                           m = mb[cinfo[k]];
                         mp0 = 0;                          mp0 = 0;
                         NEXTMP(mp0,mp);                          NEXTMP(mp0,mp);
                         chsgnq(dnsol,&dn1); mp->c = (P)dn1;                          chsgnq(dnsol,&dn1); mp->c = (P)dn1;
Line 744  NODE inv_or_split_dalg(DAlg a,DAlg *c)
Line 870  NODE inv_or_split_dalg(DAlg a,DAlg *c)
                                         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(nv,mp0,u);
                         NEXTNODE(nd0,nd);                          NEXTNODE(nd0,nd);
                         BDY(nd) = (pointer)u;                          BDY(nd) = (pointer)u;
                         NEXT(nd) = 0;                          NEXT(nd) = 0;
Legend:
Removed from v.1.10  
changed lines
  Added in v.1.16
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>