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

Diff for /OpenXM_contrib2/asir2000/builtin/algnum.c between version 1.2 and 1.9

version 1.2, 2000/08/21 08:31:18 version 1.9, 2004/12/06 09:29:34
Line 23 
Line 23 
  * shall be made on your publication or presentation in any form of the   * shall be made on your publication or presentation in any form of the
  * results obtained by use of the SOFTWARE.   * results obtained by use of the SOFTWARE.
  * (4) In the event that you modify the SOFTWARE, you shall notify FLL by   * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
  * e-mail at risa-admin@flab.fujitsu.co.jp of the detailed specification   * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
  * for such modification or the source code of the modified part of the   * for such modification or the source code of the modified part of the
  * SOFTWARE.   * SOFTWARE.
  *   *
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.1.1.1 1999/12/03 07:39:07 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/algnum.c,v 1.8 2004/12/06 01:15:18 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
   
 void Pdefpoly(), Pnewalg(), Pmainalg(), Palgtorat(), Prattoalg(), Pgetalg();  void Pdefpoly(), Pnewalg(), Pmainalg(), Palgtorat(), Prattoalg(), Pgetalg();
 void Palg(), Palgv(), Pgetalgtree();  void Palg(), Palgv(), Pgetalgtree();
   void Pinvalg_le();
   void Pset_field(),Palgtodalg(),Pdalgtoalg();
   
 #if defined(THINK_C)  
 void mkalg(P,Alg *);  void mkalg(P,Alg *);
 int cmpalgp(P,P);  int cmpalgp(P,P);
 void algptop(P,P *);  void algptop(P,P *);
 void algtorat(Num,Obj *);  void algtorat(Num,Obj *);
 void rattoalg(Obj,Alg *);  void rattoalg(Obj,Alg *);
 void ptoalgp(P,P *);  void ptoalgp(P,P *);
 #else  void clctalg(P,VL *);
 void mkalg();  void get_algtree(Obj f,VL *r);
 int cmpalgp();  
 void algptop();  
 void algtorat();  
 void rattoalg();  
 void ptoalgp();  
 void clctalg();  
 #endif  
   
 struct ftab alg_tab[] = {  struct ftab alg_tab[] = {
           {"set_field",Pset_field,1},
           {"algtodalg",Palgtodalg,1},
           {"dalgtoalg",Pdalgtoalg,1},
           {"invalg_le",Pinvalg_le,1},
         {"defpoly",Pdefpoly,1},          {"defpoly",Pdefpoly,1},
         {"newalg",Pnewalg,1},          {"newalg",Pnewalg,1},
         {"mainalg",Pmainalg,1},          {"mainalg",Pmainalg,1},
Line 85  struct ftab alg_tab[] = {
Line 83  struct ftab alg_tab[] = {
   
 static int UCN,ACNT;  static int UCN,ACNT;
   
   void Pset_field(NODE arg,Q *rp)
   {
           setfield_dalg(BDY((LIST)ARG0(arg)));
           *rp = 0;
   }
   
   void Palgtodalg(NODE arg,DAlg *rp)
   {
           algtodalg((Alg)ARG0(arg),rp);
   }
   
   void Pdalgtoalg(NODE arg,Alg *rp)
   {
           dalgtoalg((DAlg)ARG0(arg),rp);
   }
   
 void Pnewalg(arg,rp)  void Pnewalg(arg,rp)
 NODE arg;  NODE arg;
 Alg *rp;  Alg *rp;
Line 253  LIST *rp;
Line 267  LIST *rp;
         Alg b;          Alg b;
         NODE n0,n;          NODE n0,n;
   
   #if 0
         if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )          if ( !(a = (Num)ARG0(arg)) || NID(a) <= N_R )
                 vl = 0;                  vl = 0;
         else {          else {
                 t = BDY((Alg)a);                  t = BDY((Alg)a);
                 switch ( OID(t) ) {                  switch ( OID(t) ) {
                         case O_P:                          case O_P:
                                 clctalg(t,&vl); break;                                  clctalg((P)t,&vl); break;
                         case O_R:                          case O_R:
                                 clctalg(NM((R)t),&vl1);                                  clctalg(NM((R)t),&vl1);
                                 clctalg(DN((R)t),&vl2);                                  clctalg(DN((R)t),&vl2);
Line 268  LIST *rp;
Line 283  LIST *rp;
                                 vl = 0; break;                                  vl = 0; break;
                 }                  }
         }          }
   #else
           get_algtree((Obj)ARG0(arg),&vl);
   #endif
         for ( n0 = 0; vl; vl = NEXT(vl) ) {          for ( n0 = 0; vl; vl = NEXT(vl) ) {
                 NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;                  NEXTNODE(n0,n); MKV(vl->v,p); MKAlg(p,b); BDY(n) = (pointer)b;
         }          }
Line 419  P p,*r;
Line 437  P p,*r;
                 }                  }
                 NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);                  NEXT(dcr) = 0; MKP((V)(VR(p)->priv),dcr0,*r);
         }          }
   }
   
   void invalg_le(Alg a,LIST *r);
   
   void Pinvalg_le(NODE arg,LIST *r)
   {
           invalg_le((Alg)ARG0(arg),r);
   }
   
   typedef struct oMono_nf {
           DP mono;
           DP nf;
           Q dn;
   } *Mono_nf;
   
   void invalg_le(Alg a,LIST *r)
   {
           Alg inv;
           MAT mobj,sol;
           int *rinfo,*cinfo;
           P p,dn,dn1,ap;
           VL vl,tvl;
           Q c1,c2,c3,cont,c,two,iq,dn0,mul,dnsol;
           int i,j,n,len,k;
           MP mp,mp0;
           DP dp,nm,nm1,m,d,u,u1;
           NODE b,b1,hlist,mblist,t,s,rev0,rev,hist;
           DP *ps;
           struct order_spec *spec;
           Mono_nf h,h1;
           N nq,nr,nl,ng;
           Q **mat,**solmat;
           Q *w;
           int *wi;
   
           ap = (P)BDY(a);
           asir_assert(ap,O_P,"invalg_le");
   
           /* collecting algebraic numbers */
           clctalg(ap,&vl);
   
           /* setup */
           ptozp(ap,1,&c,&p);
           STOQ(2,two); create_order_spec(0,(Obj)two,&spec); initd(spec);
           for ( n = 0, tvl = vl; tvl; tvl = NEXT(tvl), n++ );
           ps = (DP *)ALLOCA(n*sizeof(DP));
   
           /* conversion to DP */
           for ( i = 0, tvl = vl; i < n; i++, tvl = NEXT(tvl) ) {
                   ptod(ALG,vl,tvl->v->attr,&ps[i]);
           }
           ptod(ALG,vl,p,&dp);
           /* index list */
           for ( b = 0, i = 0; i < n; i++ ) {
                   STOQ(i,iq); MKNODE(b1,(pointer)iq,b); b = b1;
           }
           /* simplification */
           dp_true_nf(b,dp,ps,1,&nm,&dn);
   
           /* construction of NF table */
   
           /* stdmono: <<0,...,0>> < ... < max */
           for ( hlist = 0, i = 0; i < n; i++ ) {
                   MKNODE(b1,(pointer)ps[i],hlist); hlist = b1;
           }
           dp_mbase(hlist,&rev0);
           for ( mblist = 0, rev = rev0; rev; rev = NEXT(rev) ) {
                   MKNODE(b1,BDY(rev),mblist); mblist = b1;
           }
           dn0 = ONE;
           for ( hist = 0, t = mblist; t; t = NEXT(t) ) {
                   /* searching a predecessor */
                   for ( m = (DP)BDY(t), s = hist; s; s = NEXT(s) ) {
                           h = (Mono_nf)BDY(s);
                           if ( dp_redble(m,h->mono) )
                                   break;
                   }
                   h1 = (Mono_nf)ALLOCA(sizeof(struct oMono_nf));
                   if ( s ) {
                           dp_subd(m,h->mono,&d);
                           muld(CO,d,h->nf,&u);
                           dp_true_nf(b,u,ps,1,&nm1,&dn1);
                           mulq(h->dn,(Q)dn1,&h1->dn);
                   } else {
                           muld(CO,m,nm,&u);
                           dp_true_nf(b,u,ps,1,&nm1,&dn1);
                           h1->dn = (Q)dn1;
                   }
                   h1->mono = m;
                   h1->nf = nm1;
                   MKNODE(b1,(pointer)h1,hist); hist = b1;
   
                   /* dn0 = LCM(dn0,h1->dn) */
                   gcdn(NM(dn0),NM(h1->dn),&ng); divn(NM(dn0),ng,&nq,&nr);
                   muln(nq,NM(h1->dn),&nl); NTOQ(nl,1,dn0);
           }
           /* create a matrix */
           len = length(mblist);
           MKMAT(mobj,len,len+1);
           mat = (Q **)BDY(mobj);
           mat[len-1][len] = dn0;
           for ( j = 0, t = hist; j < len; j++, t = NEXT(t) ) {
                   h = (Mono_nf)BDY(t);
                   nm1 = h->nf;
                   divq((Q)dn0,h->dn,&mul);
                   for ( i = 0, rev = rev0, mp = BDY(nm1); mp && i < len; i++, rev = NEXT(rev) )
                           if ( dl_equal(n,BDY((DP)BDY(rev))->dl,mp->dl) ) {
                                   mulq(mul,(Q)mp->c,&mat[i][j]);
                                   mp = NEXT(mp);
                           }
           }
   #if 0
           w = (Q *)ALLOCA((len+1)*sizeof(Q));
           wi = (int *)ALLOCA((len+1)*sizeof(int));
           for ( i = 0; i < len; i++ ) {
                   for ( j = 0, k = 0; j <= len; j++ )
                           if ( mat[i][j] ) {
                                   w[k] = mat[i][j];
                                   wi[k] = j;
                                   k++;
                           }
                   removecont_array(w,k);
                   for ( j = 0; j < k; j++ )
                           mat[i][wi[j]] = w[j];
           }
   #endif
           generic_gauss_elim_hensel(mobj,&sol,&dnsol,&rinfo,&cinfo);
           solmat = (Q **)BDY(sol);
           for ( i = 0, t = rev0, mp0 = 0; i < len; i++, t = NEXT(t) )
                   if ( solmat[i][0] ) {
                           NEXTMP(mp0,mp);
                           mp->c = (P)solmat[i][0];
                           mp->dl = BDY((DP)BDY(t))->dl;
                   }
           NEXT(mp) = 0; MKDP(n,mp0,u);
           dp_ptozp(u,&u1);
           divq((Q)BDY(u)->c,(Q)BDY(u1)->c,&cont);
           dtop(ALG,vl,u1,&ap);
           MKAlg(ap,inv);
           mulq(dnsol,(Q)dn,&c1);
           mulq(c1,c,&c2);
           divq(c2,cont,&c3);
           b = mknode(2,inv,c3);
           MKLIST(*r,b);
   }
   
   void get_algtree(Obj f,VL *r)
   {
           VL vl1,vl2,vl3;
           Obj t;
           DCP dc;
           NODE b;
           pointer *a;
           pointer **m;
           int len,row,col,i,j,l;
   
           if ( !f ) *r = 0;
           else
                   switch ( OID(f) ) {
                           case O_N:
                                   if ( NID((Num)f) != N_A ) *r = 0;
                                   else  {
                                           t = BDY((Alg)f);
                                           switch ( OID(t) ) {
                                                   case O_P:
                                                           clctalg((P)t,r); break;
                                                   case O_R:
                                                           clctalg(NM((R)t),&vl1);
                                                           clctalg(DN((R)t),&vl2);
                                                           mergev(ALG,vl1,vl2,r); break;
                                                   default:
                                                           *r = 0; break;
                                           }
                                   }
                                   break;
                           case O_P:
                                   vl1 = 0;
                                   for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
                                           get_algtree((Obj)COEF(dc),&vl2);
                                           mergev(ALG,vl1,vl2,&vl3);
                                           vl1 = vl3;
                                   }
                                   *r = vl1;
                                   break;
                           case O_R:
                                   get_algtree((Obj)NM((R)f),&vl1);
                                   get_algtree((Obj)DN((R)f),&vl2);
                                   mergev(ALG,vl1,vl2,r);
                                   break;
                           case O_LIST:
                                   vl1 = 0;
                                   for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
                                           get_algtree((Obj)BDY(b),&vl2);
                                           mergev(ALG,vl1,vl2,&vl3);
                                           vl1 = vl3;
                                   }
                                   *r = vl1;
                                   break;
                           case O_VECT:
                                   vl1 = 0;
                                   l = ((VECT)f)->len;
                                   a = BDY((VECT)f);
                                   for ( i = 0; i < l; i++ ) {
                                           get_algtree((Obj)a[i],&vl2);
                                           mergev(ALG,vl1,vl2,&vl3);
                                           vl1 = vl3;
                                   }
                                   *r = vl1;
                                   break;
                           case O_MAT:
                                   vl1 = 0;
                                   row = ((MAT)f)->row; col = ((MAT)f)->col;
                                   m = BDY((MAT)f);
                                   for ( i = 0; i < row; i++ )
                                           for ( j = 0; j < col; j++ ) {
                                                   get_algtree((Obj)m[i][j],&vl2);
                                                   mergev(ALG,vl1,vl2,&vl3);
                                                   vl1 = vl3;
                                           }
                                   *r = vl1;
                                   break;
                           default:
                                   *r = 0;
                                   break;
                   }
   }
   
   void algobjtorat(Obj f,Obj *r)
   {
           Obj t;
           DCP dc,dcr,dcr0;
           P p,nm,dn;
           R rat;
           NODE b,s,s0;
           VECT v;
           MAT mat;
           LIST list;
           pointer *a;
           pointer **m;
           int len,row,col,i,j,l;
   
           if ( !f ) *r = 0;
           else
                   switch ( OID(f) ) {
                           case O_N:
                                   algtorat((Num)f,r);
                                   break;
                           case O_P:
                                   dcr0 = 0;
                                   for ( dc = DC((P)f); dc; dc = NEXT(dc) ) {
                                           NEXTDC(dcr0,dcr);
                                           algobjtorat((Obj)COEF(dc),&t);
                                           COEF(dcr) = (P)t;
                                           DEG(dcr) = DEG(dc);
                                   }
                                   NEXT(dcr) = 0; MKP(VR((P)f),dcr0,p); *r = (Obj)p;
                                   break;
                           case O_R:
                                   algobjtorat((Obj)NM((R)f),&t); nm = (P)t;
                                   algobjtorat((Obj)DN((R)f),&t); dn = (P)t;
                                   MKRAT(nm,dn,0,rat); *r = (Obj)rat;
                                   break;
                           case O_LIST:
                                   s0 = 0;
                                   for ( b = BDY((LIST)f); b; b = NEXT(b) ) {
                                           NEXTNODE(s0,s);
                                           algobjtorat((Obj)BDY(b),&t);
                                           BDY(s) = (pointer)t;
                                   }
                                   NEXT(s) = 0;
                                   MKLIST(list,s0);
                                   *r = (Obj)list;
                                   break;
                           case O_VECT:
                                   l = ((VECT)f)->len;
                                   a = BDY((VECT)f);
                                   MKVECT(v,l);
                                   for ( i = 0; i < l; i++ ) {
                                           algobjtorat((Obj)a[i],&t);
                                           BDY(v)[i] = (pointer)t;
                                   }
                                   *r = (Obj)v;
                                   break;
                           case O_MAT:
                                   row = ((MAT)f)->row; col = ((MAT)f)->col;
                                   m = BDY((MAT)f);
                                   MKMAT(mat,row,col);
                                   for ( i = 0; i < row; i++ )
                                           for ( j = 0; j < col; j++ ) {
                                                   algobjtorat((Obj)m[i][j],&t);
                                                   BDY(mat)[i][j] = (pointer)t;
                                           }
                                   *r = (Obj)mat;
                                   break;
                           default:
                                   *r = f;
                                   break;
                   }
 }  }
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.9
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>