version 1.3, 2009/10/13 01:28:00 |
version 1.19, 2011/03/30 05:07:01 |
|
|
|
/* $OpenXM: OpenXM/src/asir-contrib/testing/noro/ndbf.rr,v 1.18 2011/01/16 08:46:10 noro Exp $ */ |
/* requires 'primdec' */ |
/* requires 'primdec' */ |
|
|
#define TMP_H hhhhhhhh |
#define TMP_H hhhhhhhh |
#define TMP_S ssssssss |
#define TMP_S ssssssss |
#define TMP_DS dssssssss |
#define TMP_DS dssssssss |
#define TMP_T t |
#define TMP_T tttttttt |
#define TMP_DT dt |
#define TMP_DT dtttttttt |
#define TMP_Y1 yyyyyyyy1 |
#define TMP_Y1 yyyyyyyy1 |
#define TMP_DY1 dyyyyyyyy1 |
#define TMP_DY1 dyyyyyyyy1 |
#define TMP_Y2 yyyyyyyy2 |
#define TMP_Y2 yyyyyyyy2 |
|
|
|
|
if (!module_definedp("gr")) load("gr")$ else{ }$ |
if (!module_definedp("gr")) load("gr")$ else{ }$ |
if (!module_definedp("primdec")) load("primdec")$ else{ }$ |
if (!module_definedp("primdec")) load("primdec")$ else{ }$ |
|
if (!module_definedp("newsyz")) load("noro_module_syz.rr")$ else{ }$ |
/* Empty for now. It will be used in a future. */ |
/* Empty for now. It will be used in a future. */ |
|
|
/* toplevel */ |
/* toplevel */ |
|
|
/* bfunction */ |
/* bfunction */ |
|
|
localf bfunction, in_ww, in_ww_main, ann, ann_n$ |
localf bfunction, in_ww, in_ww_main, ann, ann_n$ |
localf ann0, psi, ww_weight, compare_first, generic_bfct$ |
localf ann0, ann_fa, psi, ww_weight, compare_first, generic_bfct$ |
localf generic_bfct_1, initial_part, bfct, indicial1, bfct_via_gbfct$ |
localf generic_bfct_1, initial_part, bfct, indicial1, bfct_via_gbfct$ |
localf bfct_via_gbfct_weight, bfct_via_gbfct_weight_1, bfct_via_gbfct_weight_2$ |
localf bfct_via_gbfct_weight, bfct_via_gbfct_weight_1, bfct_via_gbfct_weight_2$ |
localf weyl_minipolym, weyl_minipoly, weyl_nf, weyl_nf_quo_check$ |
localf weyl_minipolym, weyl_minipoly, weyl_nf, weyl_nf_quo_check$ |
Line 29 localf weyl_nf_quo, weyl_nf_mod, b_subst, v_factorial, |
|
Line 31 localf weyl_nf_quo, weyl_nf_mod, b_subst, v_factorial, |
|
localf replace_vars_f, replace_vars_v, replace_var$ |
localf replace_vars_f, replace_vars_v, replace_var$ |
localf action_on_gfs, action_on_gfs_1$ |
localf action_on_gfs, action_on_gfs_1$ |
localf nd_gb_candidate$ |
localf nd_gb_candidate$ |
|
localf in_gb_oaku, homogenize_oaku$ |
|
|
/* stratification */ |
/* stratification */ |
|
|
Line 42 localf ideal_intersection$ |
|
Line 45 localf ideal_intersection$ |
|
|
|
def bfunction(F) |
def bfunction(F) |
{ |
{ |
|
if ( member(s,vars(F)) ) |
|
error("ann : the variable 's' is reserved."); |
|
/* F -> F/Fcont */ |
|
F1 = ptozp(F); Fcont = sdiv(F,F1); F = F1; |
|
|
if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; |
if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; |
if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; |
if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; |
if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; |
if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; |
|
if ( type(Op=getopt(op)) == -1 ) Op = 0; |
L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord); |
L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord); |
Indata = L[0]; AllData = L[1]; VData = L[2]; |
Indata = L[0]; AllData = L[1]; VData = L[2]; |
GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4]; |
GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4]; |
|
|
dp_set_weight(W); |
dp_set_weight(W); |
B = weyl_minipoly(GIN,VDV,0,WVDV); |
B = weyl_minipoly(GIN,VDV,0,WVDV); |
dp_set_weight(0); |
dp_set_weight(0); |
return subst(B,s,-s-1); |
if ( !Op ) return subst(B,s,-s-1); |
|
|
|
V0 = VData[0]; DV0 = VData[1]; T = VData[2]; DT = VData[3]; |
|
BPT = weyl_subst(B,T*DT,VDV); |
|
|
|
/* computation using G0,GIN0,VDV0 */ |
|
G0 = AllData[0]; GIN0 = AllData[1]; VDV0 = AllData[2]; WtV0 = AllData[5]; |
|
dp_set_weight(WtV0); dp_ord(0); |
|
PS = map(dp_ptod,GIN0,VDV0); Len = length(PS); |
|
for ( I = Len-1, Ind = []; I >= 0; I-- ) Ind = cons(I,Ind); |
|
/* QR = [D,M,Coef] */ |
|
Ax = 1; |
|
AxBPT = dp_ptod(Ax*BPT,VDV0); |
|
QR = weyl_nf_quo(Ind,AxBPT,1,PS); |
|
if ( !weyl_nf_quo_check(AxBPT,PS,QR) ) |
|
error("bfunction : invalid quotient"); |
|
if ( QR[0] ) error("bfunction : invalid quotient"); |
|
Den = QR[1]; Coef = QR[2]; |
|
for ( I = 0, R = Den*AxBPT; I < Len; I++ ) |
|
R -= dp_weyl_mul(Coef[I],dp_ptod(G0[I],VDV0)); |
|
R = dp_dtop(R,VDV0); |
|
CR = conv_tdt(R,F,V0,DV0,T,DT); |
|
|
|
dp_set_weight(0); |
|
Cont = cont(CR); CR /= Cont; |
|
Cont *= dn(Fcont); Den *= nm(Fcont); |
|
Gcd = igcd(Den,Cont); |
|
return [subst(B,s,-s-1),(Cont*CR)/(Den*Ax)]; |
} |
} |
|
|
/* |
/* |
|
|
{ |
{ |
if ( member(s,vars(F)) ) |
if ( member(s,vars(F)) ) |
error("ann : the variable 's' is reserved."); |
error("ann : the variable 's' is reserved."); |
|
if ( type(Vord=getopt(vord)) == -1 ) Vord = 0; |
|
F = ptozp(F); |
V = vars(F); |
V = vars(F); |
|
if ( Vord ) { |
|
Param = setminus(V,Vord); |
|
V = Vord; |
|
} |
N = length(V); |
N = length(V); |
D = newvect(N); |
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++ ) |
Wt1 = vector(N); |
D[I] = [deg(F,V[I]),V[I]]; |
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); |
qsort(D,compare_first); |
for ( V = [], I = N-1; I >= 0; I-- ) |
for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V); |
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-- ) |
for ( I = N-1, DV = []; I >= 0; I-- ) |
DV = cons(strtov("d"+rtostr(V[I])),DV); |
DV = cons(strtov("d"+rtostr(V[I])),DV); |
|
|
B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B); |
B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B); |
} |
} |
|
|
/* homogenized (heuristics) */ |
Tdeg = w_tdeg(F,V,Wt); |
dp_nelim(2); |
/* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */ |
G0 = nd_weyl_gr(B,append(W,DW),0,[[0,2],[0,length(W)*2-2]]); |
/* 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 = []; |
G1 = []; |
for ( T = G0; T != []; T = cdr(T) ) { |
for ( T = G0; T != []; T = cdr(T) ) { |
E = car(T); VL = vars(E); |
E = car(T); VL = vars(E); |
|
|
return G3; |
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)]; |
|
} |
|
|
|
/* homogenization w.r.t. (-W,W)-weight */ |
|
/* VDV = [x1,...,xn,dx1,...,dxn] */ |
|
/* homogenize F w.r.t. (W,-W,1) for (x,dx,y) */ |
|
|
|
def homogenize_oaku(F,VDV,W,Y) |
|
{ |
|
N = length(VDV); |
|
if ( N%2 ) error("invalid variable list"); |
|
N2 = N/2; |
|
if ( length(W) != N2 ) error("inconsistent weight vector"); |
|
W0 = dp_set_weight(); |
|
Wt = append(W,append(vtol(-ltov(W)),[1])); |
|
dp_set_weight(Wt); |
|
H = homogenize(F,VDV,Y); |
|
dp_set_weight(W0); |
|
if ( type(Vars=getopt(vars)) != -1 && Vars ) { |
|
DY = strtov("d"+rtostr(Y)); |
|
for ( I = 0, T = VDV, V = []; I < N2; I++, T = cdr(T) ) |
|
V = cons(car(T),V); |
|
T = cons(Y,append(T,[DY])); |
|
for ( ; V != []; V = cdr(V) ) T = cons(car(V),T); |
|
return [H,T]; |
|
} else return H; |
|
} |
|
|
/* F = [F0,F1,...] */ |
/* F = [F0,F1,...] */ |
|
|
def ann_n(F) |
def ann_n(F) |
|
|
VA = append(U,append(W,V)); |
VA = append(U,append(W,V)); |
DVA = append(DU,append(DW,DV)); |
DVA = append(DU,append(DW,DV)); |
VDV = append(VA,DVA); |
VDV = append(VA,DVA); |
|
#if 0 |
G0 = nd_weyl_gr(B,VDV,0,[[0,2*L],[0,length(VDV)-2*L]]); |
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 = []; |
G1 = []; |
for ( T = G0; T != []; T = cdr(T) ) { |
for ( T = G0; T != []; T = cdr(T) ) { |
E = car(T); VL = vars(E); |
E = car(T); VL = vars(E); |
|
|
return [Min,map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds)]; |
return [Min,map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds)]; |
} |
} |
|
|
|
/* |
|
* For a polynomial F and a scalar A, |
|
* compute generators of Ann(F^A). |
|
*/ |
|
|
|
def ann_fa(F,A) |
|
{ |
|
if ( type(Syz=getopt(syz)) == -1 ) Syz = 0; |
|
|
|
F = subst(F,s,TMP_S); |
|
Ann = ann(F); |
|
Bf = bfunction(F); |
|
|
|
FList = cdr(fctr(Bf)); |
|
for ( T = FList, Min = 0; T != []; T = cdr(T) ) { |
|
LF = car(car(T)); |
|
Root = -coef(LF,0)/coef(LF,1); |
|
if ( dn(Root) == 1 && Root < Min ) |
|
Min = Root; |
|
} |
|
D = A-Min; |
|
if ( dn(D) != 1 || D <= 0 ) |
|
return map(ptozp,map(subst,Ann,s,A,TMP_S,s,TMP_DS,ds)); |
|
|
|
V = vars(F); |
|
for ( I = length(V)-1, DV = []; I >= 0; I-- ) |
|
DV = cons(strtov("d"+rtostr(V[I])),DV); |
|
VDV = append(V,DV); |
|
R = map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds); |
|
F = ptozp(F); |
|
|
|
if ( Syz ) { |
|
/* syzygy method */ |
|
S = newsyz.module_syz(cons(F^D,R),VDV,1,0|weyl=1); |
|
B = car(S); |
|
for ( R = []; B != []; B = cdr(B) ) |
|
if ( H = car(car(B)) ) |
|
R = cons(ptozp(H),R); |
|
} else { |
|
/* colon method */ |
|
for ( I = 0; I < D; I++ ) |
|
R = weyl_ideal_quotient(R,F,VDV); |
|
} |
|
return R; |
|
} |
|
|
def psi0(F,T,DT) |
def psi0(F,T,DT) |
{ |
{ |
D = dp_ptod(F,[T,DT]); |
D = dp_ptod(F,[T,DT]); |
Line 1086 def replace_var(V,X,Y) |
|
Line 1293 def replace_var(V,X,Y) |
|
|
|
def action_on_gfs(P,V,GFS) |
def action_on_gfs(P,V,GFS) |
{ |
{ |
|
for ( T = V, DV = []; T != []; T = cdr(T) ) |
|
DV = cons(strtov("d"+rtostr(car(T))),DV); |
|
V = append(append(V,[s]),reverse(cons(ds,DV))); |
DP = dp_ptod(P,V); |
DP = dp_ptod(P,V); |
N = length(V)/2; |
N = length(V)/2; |
for ( I = N-1, V0 = []; I >= 0; I-- ) |
for ( I = N-1, V0 = []; I >= 0; I-- ) |
Line 1227 def weyl_ideal_quotient(B,F,VDV) |
|
Line 1437 def weyl_ideal_quotient(B,F,VDV) |
|
O1 = [[0,1],[0,N+1]]; |
O1 = [[0,1],[0,N+1]]; |
GJ = nd_weyl_gr(J,VDV1,0,O1); |
GJ = nd_weyl_gr(J,VDV1,0,O1); |
R = elimination(GJ,VDV); |
R = elimination(GJ,VDV); |
return R; |
return map(weyl_divide_by_right,R,F,VDV,0); |
} |
} |
|
|
def bf_strat(F) |
def bf_strat(F) |
{ |
{ |
|
if ( member(s,vars(F)) ) |
|
error("ann : the variable 's' is reserved."); |
dp_ord(0); |
dp_ord(0); |
T0 = time(); |
T0 = time(); |
if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; |
if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; |
if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; |
if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; |
if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; |
if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; |
|
if ( type(ElimIdeal=getopt(elimideal)) == -1 ) ElimIdeal = 0; |
L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord); |
L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord); |
T1 = time(); |
T1 = time(); |
print(["in_ww",(T1[0]+T1[1])-(T0[0]+T0[1])]); |
print(["in_ww",(T1[0]+T1[1])-(T0[0]+T0[1])]); |
Line 1275 def bf_strat(F) |
|
Line 1488 def bf_strat(F) |
|
} |
} |
|
|
L2 = bf_strat_stage2(L); |
L2 = bf_strat_stage2(L); |
|
if ( ElimIdeal ) return L2; |
S = bf_strat_stage3(L2); |
S = bf_strat_stage3(L2); |
R = []; |
R = []; |
for ( T = S; T != []; T = cdr(T) ) { |
for ( T = S; T != []; T = cdr(T) ) { |
Line 1354 def bf_strat_stage3(L) |
|
Line 1568 def bf_strat_stage3(L) |
|
} |
} |
dp_set_weight(0); |
dp_set_weight(0); |
if ( K < J ) { |
if ( K < J ) { |
for ( L = 0, T = []; L < NDI; L++ ) |
for ( L = 0, T = []; L < NDI; L++ ) { |
|
#if 0 |
|
NewId = DK[L][1]; |
|
#else |
|
NewId = ideal_intersection(DK[L][1],DI[L][1],V0,0); |
|
#endif |
T = cons([[DK[L][0][0]*DI[L][0][0],DK[L][0][1]], |
T = cons([[DK[L][0][0]*DI[L][0][0],DK[L][0][1]], |
DK[L][1],DK[L][2]],T); |
NewId,DK[L][2]],T); |
|
} |
Data[K] = reverse(T); |
Data[K] = reverse(T); |
} else |
} else |
Data[J++] = DI; |
Data[J++] = DI; |
Line 1379 def bf_strat_stage3(L) |
|
Line 1599 def bf_strat_stage3(L) |
|
|
|
def bf_local(F,P) |
def bf_local(F,P) |
{ |
{ |
|
if ( member(s,vars(F)) ) |
|
error("ann : the variable 's' is reserved."); |
|
/* F -> F/Fcont */ |
|
F1 = ptozp(F); Fcont = sdiv(F,F1); F = F1; |
if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; |
if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0; |
if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; |
if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0; |
if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; |
if ( type(Wt=getopt(weight)) == -1 ) Wt = 0; |
Line 1443 def bf_local(F,P) |
|
Line 1667 def bf_local(F,P) |
|
CR = conv_tdt(R,F,V0,DV0,T,DT); |
CR = conv_tdt(R,F,V0,DV0,T,DT); |
|
|
dp_set_weight(0); |
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 */ |
/* t^(l+k)*dt^l (k>l) -> (s-k)(s-k-1)...(s-(k+l-1))t^k */ |