=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/engine/nd.c,v retrieving revision 1.114 retrieving revision 1.118 diff -u -p -r1.114 -r1.118 --- OpenXM_contrib2/asir2000/engine/nd.c 2004/10/25 04:19:50 1.114 +++ OpenXM_contrib2/asir2000/engine/nd.c 2004/12/04 09:39:27 1.118 @@ -1,4 +1,4 @@ -/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.113 2004/10/06 12:09:19 noro Exp $ */ +/* $OpenXM: OpenXM_contrib2/asir2000/engine/nd.c,v 1.117 2004/12/03 08:57:30 noro Exp $ */ #include "nd.h" @@ -8,6 +8,7 @@ NM _nm_free_list; ND _nd_free_list; ND_pairs _ndp_free_list; +static int nd_nalg; #if 0 static int ndv_alloc; #endif @@ -43,6 +44,7 @@ static int nd_demand; UINT *nd_det_compute_bound(NDV **dm,int n,int j); void nd_det_reconstruct(NDV **dm,int n,int j,NDV d); +int nd_monic(int m,ND *p); void nd_free_private_storage() { @@ -1580,7 +1582,7 @@ NODE nd_gb(int m,int ishomo,int checkonly) NODE r,g,t; ND_pairs d; ND_pairs l; - ND h,nf,s,head; + ND h,nf,s,head,nf1; NDV nfv; Q q,num,den; union oNDC dn; @@ -1619,6 +1621,10 @@ again: if ( checkonly ) return 0; if ( DP_Print ) { printf("+"); fflush(stdout); } nd_removecont(m,nf); + if ( !m && nd_nalg ) { + nd_monic(0,&nf); + nd_removecont(m,nf); + } nfv = ndtondv(m,nf); nd_free(nf); nh = ndv_newps(m,nfv,0); d = update_pairs(d,g,nh); @@ -1685,6 +1691,9 @@ void do_diagonalize_trace(int sugar,int m) } } +static struct oEGT eg_invdalg; +struct oEGT eg_le; + NODE nd_gb_trace(int m,int ishomo) { int i,nh,sugar,stat; @@ -1695,7 +1704,11 @@ NODE nd_gb_trace(int m,int ishomo) NDV nfv,nfqv; Q q,den,num; union oNDC dn; + struct oEGT eg_monic,egm0,egm1; + init_eg(&eg_monic); + init_eg(&eg_invdalg); + init_eg(&eg_le); g = 0; d = 0; for ( i = 0; i < nd_psn; i++ ) { d = update_pairs(d,g,i); @@ -1743,8 +1756,17 @@ again: if ( !rem(NM(HCQ(nfq)),m) ) return 0; if ( DP_Print ) { printf("+"); fflush(stdout); } - nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf); - nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq); + if ( nd_nalg ) { + /* m|DN(HC(nf)^(-1)) => failure */ + get_eg(&egm0); + if ( !nd_monic(m,&nfq) ) return 0; + get_eg(&egm1); add_eg(&eg_monic,&egm0,&egm1); + nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq); + nfv = ndv_dup(0,nfqv); ndv_mod(m,nfv); nd_free(nf); + } else { + nd_removecont(0,nfq); nfqv = ndtondv(0,nfq); nd_free(nfq); + nd_removecont(m,nf); nfv = ndtondv(m,nf); nd_free(nf); + } nh = ndv_newps(0,nfv,nfqv); d = update_pairs(d,g,nh); g = update_base(g,nh); @@ -1762,6 +1784,11 @@ again: else for ( t = g; t; t = NEXT(t) ) BDY(t) = (pointer)nd_ps_trace[(int)BDY(t)]; + if ( nd_nalg ) { + print_eg("monic",&eg_monic); + print_eg("invdalg",&eg_invdalg); + print_eg("le",&eg_le); + } return g; } @@ -2622,6 +2649,7 @@ void nd_mul_c(int mod,ND p,int mul) int c,c1; if ( !p ) return; + if ( mul == 1 ) return; if ( mod == -1 ) for ( m = BDY(p); m; m = NEXT(m) ) CM(m) = _mulsf(CM(m),mul); @@ -2637,6 +2665,7 @@ void nd_mul_c_q(ND p,Q mul) Q c; if ( !p ) return; + if ( UNIQ(mul) ) return; for ( m = BDY(p); m; m = NEXT(m) ) { mulq(CQ(m),mul,&c); CQ(m) = c; } @@ -3715,6 +3744,10 @@ void nd_nf_p(P f,LIST g,LIST v,int m,struct order_spec int stat,nvar,max,e; union oNDC dn; + if ( !f ) { + *rp = 0; + return; + } pltovl(v,&vv); for ( nvar = 0, tv = vv; tv; tv = NEXT(tv), nvar++ ); @@ -5089,4 +5122,131 @@ UINT *nd_det_compute_bound(NDV **dm,int n,int j) r = (UINT *)ALLOCA(nd_wpd*sizeof(UINT)); for ( k = 0; k < nd_wpd; k++ ) r[k] = d0[k]; return r; +} + +DL nd_separate_d(UINT *d,UINT *trans) +{ + int n,ntrans,td,i,e; + DL a; + + n = nd_nvar; ntrans = n-nd_nalg; + ndl_zero(trans); + td = 0; + for ( i = 0; i < ntrans; i++ ) { + e = GET_EXP(d,i); + PUT_EXP(trans,i,e); + td += MUL_WEIGHT(e,i); + } + TD(trans) = td; + if ( nd_blockmask) ndl_weight_mask(trans); + NEWDL(a,nd_nalg); + td = 0; + for ( ; i < n; i++ ) { + e = GET_EXP(d,i); + a->d[i-ntrans] = e; + td += e; + } + a->td = td; + return a; +} + +NumberField get_numberfield(); + +int nd_monic(int mod,ND *p) +{ + UINT *trans,*t; + DL alg; + MP mp0,mp; + NM m,m0,m1,ma0,ma,mb,mr0,mr; + ND r; + DL dl; + DP nm; + NDV ndv; + DAlg inv,cd; + ND s,c; + Q l,mul; + N ln; + int n,ntrans,i,e,td,is_lc,len; + NumberField nf; + struct oEGT eg0,eg1; + + if ( !(nf = get_numberfield()) ) + error("nd_monic : current_numberfield is not set"); + + n = nd_nvar; ntrans = n-nd_nalg; + /* Q coef -> DAlg coef */ + NEWNM(ma0); ma = ma0; + m = BDY(*p); + is_lc = 1; + while ( 1 ) { + NEWMP(mp0); mp = mp0; + mp->c = (P)CQ(m); + mp->dl = nd_separate_d(DL(m),DL(ma)); + NEWNM(mb); + for ( m = NEXT(m); m; m = NEXT(m) ) { + alg = nd_separate_d(DL(m),DL(mb)); + if ( !ndl_equal(DL(ma),DL(mb)) ) + break; + NEXTMP(mp0,mp); mp->c = (P)CQ(m); mp->dl = alg; + } + NEXT(mp) = 0; + MKDP(nd_nalg,mp0,nm); + MKDAlg(nm,ONE,cd); + if ( is_lc == 1 ) { + /* if the lc is a rational number, we have nothing to do */ + if ( !mp0->dl->td ) + return 1; + + get_eg(&eg0); + invdalg(cd,&inv); + get_eg(&eg1); add_eg(&eg_invdalg,&eg0,&eg1); + /* check the validity of inv */ + if ( mod && !rem(NM(inv->dn),mod) ) + return 0; + CA(ma) = nf->one; + is_lc = 0; + ln = ONEN; + } else { + muldalg(cd,inv,&CA(ma)); + lcmn(ln,NM(CA(ma)->dn),&ln); + } + if ( m ) { + NEXT(ma) = mb; ma = mb; + } else { + NEXT(ma) = 0; + break; + } + } + /* l = lcm(denoms) */ + NTOQ(ln,1,l); + for ( mr0 = 0, m = ma0; m; m = NEXT(m) ) { + divq(l,CA(m)->dn,&mul); + for ( mp = BDY(CA(m)->nm); mp; mp = NEXT(mp) ) { + NEXTNM(mr0,mr); + mulq((Q)mp->c,mul,&CQ(mr)); + dl = mp->dl; + td = TD(DL(m)); + ndl_copy(DL(m),DL(mr)); + for ( i = ntrans; i < n; i++ ) { + e = dl->d[i-ntrans]; + PUT_EXP(DL(mr),i,e); + td += MUL_WEIGHT(e,i); + } + TD(DL(mr)) = td; + if ( nd_blockmask) ndl_weight_mask(DL(mr)); + } + } + NEXT(mr) = 0; + for ( len = 0, mr = mr0; mr; mr = NEXT(mr), len++ ); + MKND(NV(*p),mr0,len,r); + /* XXX */ + SG(r) = SG(*p); + nd_free(*p); + *p = r; + return 1; +} + +void nd_set_nalg(int nalg) +{ + nd_nalg = nalg; }