=================================================================== RCS file: /home/cvs/OpenXM/src/asir-contrib/testing/noro/ndbf.rr,v retrieving revision 1.3 retrieving revision 1.8 diff -u -p -r1.3 -r1.8 --- OpenXM/src/asir-contrib/testing/noro/ndbf.rr 2009/10/13 01:28:00 1.3 +++ OpenXM/src/asir-contrib/testing/noro/ndbf.rr 2009/10/27 02:53:31 1.8 @@ -29,6 +29,7 @@ localf weyl_nf_quo, weyl_nf_mod, b_subst, v_factorial, localf replace_vars_f, replace_vars_v, replace_var$ localf action_on_gfs, action_on_gfs_1$ localf nd_gb_candidate$ +localf in_gb_oaku$ /* stratification */ @@ -263,15 +264,30 @@ def ann(F) { if ( member(s,vars(F)) ) error("ann : the variable 's' is reserved."); + F = ptozp(F); V = vars(F); N = length(V); D = newvect(N); + if ( type(Wt=getopt(weight)) == -1 ) + for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt); - for ( I = 0; I < N; I++ ) - D[I] = [deg(F,V[I]),V[I]]; + Wt1 = vector(N); + for ( I = 0, F1 =F; I < N; I++ ) { + VI = Wt[2*I]; WI = Wt[2*I+1]; + for ( J = 0; J < N; J++ ) + if ( VI == V[J] ) break; + F1 = subst(F1,VI,VI^WI); + } + for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]]; qsort(D,compare_first); - for ( V = [], I = N-1; I >= 0; I-- ) - V = cons(D[I][1],V); + for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V); + V = reverse(V); + for ( I = 0; I < N; I++ ) { + VI = Wt[2*I]; WI = Wt[2*I+1]; + for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break; + Wt1[J] = WI; + } + Wt = vtol(Wt1); for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); @@ -284,9 +300,24 @@ def ann(F) B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B); } - /* homogenized (heuristics) */ - dp_nelim(2); - G0 = nd_weyl_gr(B,append(W,DW),0,[[0,2],[0,length(W)*2-2]]); + Tdeg = w_tdeg(F,V,Wt); + /* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */ + /* weight for [y1,y2,t, x1,...,xn, dy1,dy2, dt,dx1,...,dxn, h] */ + /* 0 1 2 3 N3-1 N3 N3+1 N3+2 2*N3 */ + /* 1 1 D+1 w1 wn 1 1 1 D D 1 */ + N3 = N+3; + WtV = newvect(2*N3+1); + WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1; + for ( I = 3; I < N3; I++ ) WtV[I] = Wt[I-3]; + for ( ; I <= N3+2; I++ ) WtV[I] = 1; + for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg; + WtV[2*N3] = 1; + + /* B is already a GB => modular change of ordering can be applied */ + /* any prime is available => HC=[1] */ + dp_set_weight(WtV); + G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1); + dp_set_weight(0); G1 = []; for ( T = G0; T != []; T = cdr(T) ) { E = car(T); VL = vars(E); @@ -298,6 +329,67 @@ def ann(F) return G3; } +def in_gb_oaku(F) +{ + if ( member(s,vars(F)) ) + error("ann : the variable 's' is reserved."); + F = ptozp(F); + V = vars(F); + N = length(V); + D = newvect(N); + if ( type(Wt=getopt(weight)) == -1 ) + for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt); + + Wt1 = vector(N); + for ( I = 0, F1 =F; I < N; I++ ) { + VI = Wt[2*I]; WI = Wt[2*I+1]; + for ( J = 0; J < N; J++ ) + if ( VI == V[J] ) break; + F1 = subst(F1,VI,VI^WI); + } + for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]]; + qsort(D,compare_first); + for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V); + V = reverse(V); + for ( I = 0; I < N; I++ ) { + VI = Wt[2*I]; WI = Wt[2*I+1]; + for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break; + Wt1[J] = WI; + } + Wt = vtol(Wt1); + + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + + W = append([TMP_Y1,TMP_Y2,TMP_T],V); + DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV); + + B = [TMP_T-TMP_Y1*F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B); + } + + Tdeg = w_tdeg(F,V,Wt); + /* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */ + /* weight for [y1,y2,t, x1,...,xn, dy1,dy2, dt,dx1,...,dxn, h] */ + /* 0 1 2 3 N3-1 N3 N3+1 N3+2 2*N3 */ + /* 1 1 D+1 1 1 1 1 1 D D 1 */ + N3 = N+3; + WtV = newvect(2*N3+1); + WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1; + for ( I = 3; I <= N3+2; I++ ) WtV[I] = 1; + for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg; + WtV[2*N3] = 1; + + /* B is already a GB => modular change of ordering can be applied */ + /* any prime is available => HC=[1] */ + dp_set_weight(WtV); + G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1); + dp_set_weight(0); + G1 = map(subst,G0,TMP_Y1,1); + return [G1,append(V,DV)]; +} + /* F = [F0,F1,...] */ def ann_n(F) @@ -335,7 +427,11 @@ def ann_n(F) VA = append(U,append(W,V)); DVA = append(DU,append(DW,DV)); VDV = append(VA,DVA); +#if 0 G0 = nd_weyl_gr(B,VDV,0,[[0,2*L],[0,length(VDV)-2*L]]); +#else + G0 = nd_gb_candidate(B,VDV,[[0,2*L],[0,length(VDV)-2*L]],0,[1],1); +#endif G1 = []; for ( T = G0; T != []; T = cdr(T) ) { E = car(T); VL = vars(E); @@ -1379,6 +1475,8 @@ def bf_strat_stage3(L) def bf_local(F,P) { + /* F -> F/Fcont */ + F1 = ptozp(F); Fcont = sdiv(F,F1); F = F1; if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; @@ -1443,7 +1541,10 @@ def bf_local(F,P) CR = conv_tdt(R,F,V0,DV0,T,DT); dp_set_weight(0); - return [LB,Ax,CR]; + Cont = cont(CR); CR /= Cont; + Cont *= dn(Fcont); Den *= nm(Fcont); + Gcd = igcd(Den,Cont); + return [LB,(Den/Gcd)*Ax,(Cont/Gcd)*CR]; } /* t^(l+k)*dt^l (k>l) -> (s-k)(s-k-1)...(s-(k+l-1))t^k */