![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to ndbf.rr CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / testing / noro|
Return to ndbf.rr CVS log | Up to [local] / OpenXM / src / asir-contrib / testing / noro |
1.1 noro 1: /* requires 'primdec' */
2:
3: #define TMP_H hhhhhhhh
4: #define TMP_S ssssssss
5: #define TMP_DS dssssssss
6: #define TMP_T t
7: #define TMP_DT dt
8: #define TMP_Y1 yyyyyyyy1
9: #define TMP_DY1 dyyyyyyyy1
10: #define TMP_Y2 yyyyyyyy2
11: #define TMP_DY2 dyyyyyyyy2
12:
13: if (!module_definedp("gr")) load("gr")$ else{ }$
14: if (!module_definedp("primdec")) load("primdec")$ else{ }$
15: /* Empty for now. It will be used in a future. */
16:
17: /* toplevel */
18:
19: module ndbf$
20:
21: /* bfunction */
22:
23: localf bfunction, in_ww, in_ww_main, ann, ann_n$
24: localf ann0, psi, ww_weight, compare_first, generic_bfct$
25: localf generic_bfct_1, initial_part, bfct, indicial1, bfct_via_gbfct$
26: localf bfct_via_gbfct_weight, bfct_via_gbfct_weight_1, bfct_via_gbfct_weight_2$
27: localf weyl_minipolym, weyl_minipoly, weyl_nf, weyl_nf_quo_check$
28: localf weyl_nf_quo, weyl_nf_mod, b_subst, v_factorial, w_tdeg$
29: localf replace_vars_f, replace_vars_v, replace_var$
30: localf action_on_gfs, action_on_gfs_1$
1.2 noro 31: localf nd_gb_candidate$
1.1 noro 32:
33: /* stratification */
34:
35: localf weyl_subst, bfactor, gen_a, gen_d$
36: localf gen_dm, elimination, weyl_ideal_quotient, psi0$
37: localf bf_strat, bf_strat_stage2, bf_strat_stage3, bf_local$
38: localf conv_tdt, merge_tower, stratify_bf, tdt_homogenize$
39: localf sing, tower_in_p, subst_vars, compute_exponent$
40: localf zero_inclusion, weyl_divide_by_right, elim_mat, int_cons$
41: localf ideal_intersection$
42:
43: def bfunction(F)
44: {
45: if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
46: if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
47: if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
48: L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
49: Indata = L[0]; AllData = L[1]; VData = L[2];
50: GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4];
51: W = Indata[4];
52: dp_set_weight(W);
53: B = weyl_minipoly(GIN,VDV,0,WVDV);
54: dp_set_weight(0);
55: return subst(B,s,-s-1);
56: }
57:
58: /*
59: returns [InData,AllData,VData]
60: InData = [GIN,VDV,V,DV,WtV]
61: AllData = [G0,GIN0,VDV0,W,WVDV,WtV0]
62: VData = [V0,DV0,T,DT]
63: GIN0 = ini_(-W,W)(G0)
64: WVDV = W[0]*V[0]*DV[0]+...
65: */
66:
67: def in_ww(F)
68: {
1.2 noro 69: F = ptozp(F);
1.1 noro 70: V = vars(F);
71: N = length(V);
72: D = newvect(N);
73: Wt = getopt(weight);
74: Vord = getopt(vord);
75: if ( type(Wt) != 4 ) {
76: if ( type(Vord) != 4 ) {
77: for ( I = 0; I < N; I++ )
78: D[I] = [deg(F,V[I]),V[I]];
79: qsort(D,compare_first);
80: for ( V = [], I = 0; I < N; I++ )
81: V = cons(D[I][1],V);
82: V = reverse(V);
83: } else
84: V = Vord;
85: for ( I = 0, Wt = []; I < N; I++ )
86: Wt = cons(1,Wt);
87: } else {
88: Wt1 = vector(N);
89: if ( type(Vord) != 4 ) {
90: for ( I = 0, F1 =F; I < N; I++ ) {
91: VI = Wt[2*I]; WI = Wt[2*I+1];
92: for ( J = 0; J < N; J++ )
93: if ( VI == V[J] ) break;
94: F1 = subst(F1,VI,VI^WI);
95: }
96: for ( I = 0; I < N; I++ )
97: D[I] = [deg(F1,V[I]),V[I]];
98: qsort(D,compare_first);
99: for ( V = [], I = 0; I < N; I++ )
100: V = cons(D[I][1],V);
101: V = reverse(V);
102: } else
103: V = Vord;
104: for ( I = 0; I < N; I++ ) {
105: VI = Wt[2*I]; WI = Wt[2*I+1];
106: for ( J = 0; J < N; J++ )
107: if ( VI == V[J] ) break;
108: Wt1[J] = WI;
109: }
110: Wt = vtol(Wt1);
111: }
112: Tdeg = w_tdeg(F,V,Wt);
113: /* weight for [t,x1,...,xn,dt,dx1,...,dxn,h] */
114: WtV = newvect(2*(N+1)+1);
115: WtV[0] = Tdeg; WtV[N+1] = 1; WtV[2*(N+1)] = 1;
116: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
117: for ( I = 1; I <= N; I++ ) {
118: WtV[I] = Wt[I-1]; WtV[N+1+I] = Tdeg-Wt[I-1]+1;
119: }
120: for ( I = N-1, DV = []; I >= 0; I-- )
121: DV = cons(strtov("d"+rtostr(V[I])),DV);
122:
123: B = [TMP_T-F];
124: for ( I = 0; I < N; I++ ) {
125: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
126: }
127: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
128: W = newvect(N+1); W[0] = 1;
129: VW1 = [V1,DV1,WtV,W];
130:
131: /* weight for [x1,...,xn,t,dx1,...,dxn,dt,h] */
132: WtV = newvect(2*(N+1)+1);
133: WtV[N] = Tdeg; WtV[2*N+1] = 1; WtV[2*(N+1)] = 1;
134: for ( I = 0; I < N; I++ ) {
135: WtV[I] = Wt[I]; WtV[N+1+I] = Tdeg-Wt[I]+1;
136: }
137: for ( I = N-1, DV = []; I >= 0; I-- )
138: DV = cons(strtov("d"+rtostr(V[I])),DV);
139:
140: B = [TMP_T-F];
141: for ( I = 0; I < N; I++ ) {
142: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
143: }
144: V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_DT]);
145: W = newvect(N+1); W[N] = 1;
146: VW2 = [V1,DV1,WtV,W];
147:
148: Heu = getopt(heuristic);
149: if ( type(Heu) != -1 && Heu )
150: L = in_ww_main(B,VW1,VW2);
151: else
152: L = in_ww_main(B,VW1,0);
153: return append(L,[[V,DV,TMP_T,TMP_DT,Wt]]);
154: }
155:
156: /*
157: returns [InData,AllData]
158: InData = [GIN,VDV,V,DV,WtV]
159: AllData = [G0,GIN0,VDV0,W,WVDV,WtV0]
160: GIN0 = ini_(-W,W)(G0)
161: WVDV = W[0]*V[0]*DV[0]+...
162: */
163:
164: def in_ww_main(F,VW1,VW2)
165: {
166: V = VW1[0]; DV = VW1[1]; WtV = VW1[2]; W = VW1[3];
167: dp_set_weight(WtV);
168:
169: N = length(V);
170: N2 = N*2;
171:
172: /* If W is a list, convert it to a vector */
173: if ( type(W) == 4 )
174: W = newvect(length(W),W);
175: dp_weyl_set_weight(W);
176:
177: /* create a term order M in D<x,d> (DRL) */
178: M = newmat(N2,N2);
179: for ( J = 0; J < N2; J++ )
180: M[0][J] = 1;
181: for ( I = 1; I < N2; I++ )
182: M[I][N2-I] = -1;
183:
184: VDV = append(V,DV);
185:
186: /* create a non-term order MW in D<x,d> */
187: MW = newmat(N2+1,N2);
188: for ( J = 0; J < N; J++ )
189: MW[0][J] = -W[J];
190: for ( ; J < N2; J++ )
191: MW[0][J] = W[J-N];
192: for ( I = 1; I <= N2; I++ )
193: for ( J = 0; J < N2; J++ )
194: MW[I][J] = M[I-1][J];
195:
196: /* create a homogenized term order MWH in D<x,d,h> */
197: MWH = newmat(N2+2,N2+1);
198: for ( J = 0; J <= N2; J++ )
199: MWH[0][J] = 1;
200: for ( I = 1; I <= N2+1; I++ )
201: for ( J = 0; J < N2; J++ )
202: MWH[I][J] = MW[I-1][J];
203:
204: /* homogenize F */
205: VDVH = append(VDV,[TMP_H]);
206: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
207:
1.2 noro 208: /*
209: * FH is a GB w.r.t. any term order s.t. LT(FH)=[t,dx1,...,dxn]
210: * Compute a groebner basis of FH w.r.t. MWH by modular change of
211: * ordering.
212: * Since F is Z-coef, LC(FH)=[1,...,1] and we can use any prime p
213: * for trace algorithm.
214: */
1.1 noro 215: /* dp_gr_flags(["Top",1,"NoRA",1]); */
1.2 noro 216: for ( I = 0, HC=[]; I <= N; I++ ) HC = cons(1,HC);
217: GH = nd_gb_candidate(FH,VDVH,11,0,HC,1);
1.1 noro 218: /* dp_gr_flags(["Top",0,"NoRA",0]); */
219:
220: /* dehomigenize GH */
221: G = map(subst,GH,TMP_H,1);
222:
223: /* G is a groebner basis w.r.t. a non term order MW */
224: /* take the initial part w.r.t. (-W,W) */
225: GIN = map(initial_part,G,VDV,MW,W);
226:
227: /* GIN is a groebner basis w.r.t. a term order M */
228: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
229:
230: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
231: for ( I = 0, T = 0; I < N; I++ )
232: T += W[I]*V[I]*DV[I];
233:
234: AllData = [G,GIN,VDV,W,T,WtV];
235: if ( VW2 ) {
1.2 noro 236: /* take LC(GIN) w.r.t. DRL */
237: dp_set_weight(WtV); dp_ord(0);
238: HC = map(dp_hc,map(dp_ptod,GIN,VDV));
1.1 noro 239: V2 = VW2[0]; DV2 = VW2[1]; WtV2 = VW2[2];
240: VDV2 = append(V2,DV2);
241: dp_set_weight(WtV2);
1.2 noro 242: GIN2 = nd_gb_candidate(GIN,VDV2,0,0,HC,1);
1.1 noro 243: InData = [GIN2,VDV2,V2,DV2,WtV2];
244: } else {
245: if ( 0 ) {
246: dp_set_weight(WtV);
247: GIN1 = nd_weyl_gr_postproc(GIN,VDV,0,0,0);
248: InData = [GIN1,VDV,V,DV,WtV];
249: } else
250: InData = [GIN,VDV,V,DV,WtV];
251: }
252:
253: /* B = weyl_minipoly(GIN2,VDV2,0,T); */ /* M represents DRL order */
254: WtV = dp_set_weight();
255: dp_set_weight(0);
256:
257: return [InData,AllData];
258: }
259:
260: /* annihilating ideal of F^s */
261:
262: def ann(F)
263: {
264: if ( member(s,vars(F)) )
265: error("ann : the variable 's' is reserved.");
1.4 noro 266: F = ptozp(F);
1.1 noro 267: V = vars(F);
268: N = length(V);
269: D = newvect(N);
1.4 noro 270: if ( type(Wt=getopt(weight)) == -1 )
271: for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt);
1.1 noro 272:
1.4 noro 273: Wt1 = vector(N);
274: for ( I = 0, F1 =F; I < N; I++ ) {
275: VI = Wt[2*I]; WI = Wt[2*I+1];
276: for ( J = 0; J < N; J++ )
277: if ( VI == V[J] ) break;
278: F1 = subst(F1,VI,VI^WI);
279: }
280: for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]];
1.1 noro 281: qsort(D,compare_first);
1.4 noro 282: for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V);
283: V = reverse(V);
284: for ( I = 0; I < N; I++ ) {
285: VI = Wt[2*I]; WI = Wt[2*I+1];
286: for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break;
287: Wt1[J] = WI;
288: }
289: Wt = vtol(Wt1);
1.1 noro 290:
291: for ( I = N-1, DV = []; I >= 0; I-- )
292: DV = cons(strtov("d"+rtostr(V[I])),DV);
293:
294: W = append([TMP_Y1,TMP_Y2,TMP_T],V);
295: DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV);
296:
297: B = [1-TMP_Y1*TMP_Y2,TMP_T-TMP_Y1*F];
298: for ( I = 0; I < N; I++ ) {
299: B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
300: }
301:
1.4 noro 302: Tdeg = w_tdeg(F,V,Wt);
303: /* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */
304: /* weight for [y1,y2,t, x1,...,xn, dy1,dy2, dt,dx1,...,dxn, h] */
305: /* 0 1 2 3 N3-1 N3 N3+1 N3+2 2*N3 */
1.7 ! noro 306: /* 1 1 D+1 w1 wn 1 1 1 D D 1 */
1.4 noro 307: N3 = N+3;
308: WtV = newvect(2*N3+1);
309: WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1;
1.7 ! noro 310: for ( I = 3; I < N3; I++ ) WtV[I] = Wt[I-3];
! 311: for ( ; I <= N3+2; I++ ) WtV[I] = 1;
1.4 noro 312: for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg;
313: WtV[2*N3] = 1;
314:
315: /* B is already a GB => modular change of ordering can be applied */
316: /* any prime is available => HC=[1] */
317: dp_set_weight(WtV);
318: G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1);
319: dp_set_weight(0);
1.1 noro 320: G1 = [];
321: for ( T = G0; T != []; T = cdr(T) ) {
322: E = car(T); VL = vars(E);
323: if ( !member(TMP_Y1,VL) && !member(TMP_Y2,VL) )
324: G1 = cons(E,G1);
325: }
326: G2 = map(psi,G1,TMP_T,TMP_DT);
327: G3 = map(subst,G2,TMP_T,-1-s);
328: return G3;
329: }
330:
1.7 ! noro 331: def in_gb_oaku(F)
! 332: {
! 333: if ( member(s,vars(F)) )
! 334: error("ann : the variable 's' is reserved.");
! 335: F = ptozp(F);
! 336: V = vars(F);
! 337: N = length(V);
! 338: D = newvect(N);
! 339: if ( type(Wt=getopt(weight)) == -1 )
! 340: for ( I = N-1, Wt = []; I >= 0; I-- ) Wt = append([V[I],1],Wt);
! 341:
! 342: Wt1 = vector(N);
! 343: for ( I = 0, F1 =F; I < N; I++ ) {
! 344: VI = Wt[2*I]; WI = Wt[2*I+1];
! 345: for ( J = 0; J < N; J++ )
! 346: if ( VI == V[J] ) break;
! 347: F1 = subst(F1,VI,VI^WI);
! 348: }
! 349: for ( I = 0; I < N; I++ ) D[I] = [deg(F1,V[I]),V[I]];
! 350: qsort(D,compare_first);
! 351: for ( V = [], I = 0; I < N; I++ ) V = cons(D[I][1],V);
! 352: V = reverse(V);
! 353: for ( I = 0; I < N; I++ ) {
! 354: VI = Wt[2*I]; WI = Wt[2*I+1];
! 355: for ( J = 0; J < N; J++ ) if ( VI == V[J] ) break;
! 356: Wt1[J] = WI;
! 357: }
! 358: Wt = vtol(Wt1);
! 359:
! 360: for ( I = N-1, DV = []; I >= 0; I-- )
! 361: DV = cons(strtov("d"+rtostr(V[I])),DV);
! 362:
! 363: W = append([TMP_Y1,TMP_Y2,TMP_T],V);
! 364: DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV);
! 365:
! 366: B = [TMP_T-TMP_Y1*F];
! 367: for ( I = 0; I < N; I++ ) {
! 368: B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
! 369: }
! 370:
! 371: Tdeg = w_tdeg(F,V,Wt);
! 372: /* y1*y2-1, t-y1*f, dx1+y1*df/dx1*dt ... */
! 373: /* weight for [y1,y2,t, x1,...,xn, dy1,dy2, dt,dx1,...,dxn, h] */
! 374: /* 0 1 2 3 N3-1 N3 N3+1 N3+2 2*N3 */
! 375: /* 1 1 D+1 1 1 1 1 1 D D 1 */
! 376: N3 = N+3;
! 377: WtV = newvect(2*N3+1);
! 378: WtV[0] = WtV[1] = 1; WtV[2] = Tdeg+1;
! 379: for ( I = 3; I <= N3+2; I++ ) WtV[I] = 1;
! 380: for ( ; I < 2*N3; I++ ) WtV[I] = Tdeg;
! 381: WtV[2*N3] = 1;
! 382:
! 383: /* B is already a GB => modular change of ordering can be applied */
! 384: /* any prime is available => HC=[1] */
! 385: dp_set_weight(WtV);
! 386: G0 = nd_gb_candidate(B,append(W,DW),[[0,2],[0,length(W)*2-2]],0,[1],1);
! 387: dp_set_weight(0);
! 388: G1 = map(subst,G0,TMP_Y1,1);
! 389: return [G1,append(V,DV)];
! 390: }
! 391:
1.1 noro 392: /* F = [F0,F1,...] */
393:
394: def ann_n(F)
395: {
396: L = length(F);
397: V = vars(F);
398: N = length(V);
399: D = newvect(N);
400:
401: for ( I = N-1, DV = []; I >= 0; I-- )
402: DV = cons(strtov("d"+rtostr(V[I])),DV);
403: W = []; DW = [];
404: for ( I = L-1; I >= 0; I-- ) {
405: SI = rtostr(I);
406: W = cons(strtov("t"+SI),W);
407: DW = cons(strtov("dt"+SI),DW);
408: }
409: U = []; DU = [];
410: for ( I = L-1; I >= 0; I-- ) {
411: SI = rtostr(I);
412: U = append([strtov("u"+SI),strtov("v"+SI)],U);
413: DU = append([strtov("du"+SI),strtov("dv"+SI)],DU);
414: }
415:
416: B = [];
417: for ( I = 0; I < N; I++ ) {
418: T = DV[I];
419: for ( J = 0; J < L; J++ )
420: T += U[2*J]*diff(F[J],V[I])*DW[J];
421: B = cons(T,B);
422: }
423: for ( I = 0; I < L; I++ )
424: B = append([W[I]-U[2*I]*F[I],1-U[2*I]*U[2*I+1]],B);
425:
426: VA = append(U,append(W,V));
427: DVA = append(DU,append(DW,DV));
428: VDV = append(VA,DVA);
1.4 noro 429: #if 0
1.1 noro 430: G0 = nd_weyl_gr(B,VDV,0,[[0,2*L],[0,length(VDV)-2*L]]);
1.4 noro 431: #else
432: G0 = nd_gb_candidate(B,VDV,[[0,2*L],[0,length(VDV)-2*L]],0,[1],1);
433: #endif
1.1 noro 434: G1 = [];
435: for ( T = G0; T != []; T = cdr(T) ) {
436: E = car(T); VL = vars(E);
437: for ( TV = U; TV != []; TV = cdr(TV) )
438: if ( member(car(TV),VL) ) break;
439: if ( TV == [] )
440: G1 = cons(E,G1);
441: }
442: G2 = G1;
443: for ( I = 0; I < L; I++ ) {
444: G2 = map(psi,G2,W[I],DW[I]);
445: G2 = map(subst,G2,W[I],-1-strtov("s"+rtostr(I)));
446: }
447: return G2;
448: }
449:
450: /*
451: * compute J_f|s=r, where r = the minimal integral root of global b_f(s)
452: * ann0(F) returns [MinRoot,Ideal]
453: */
454:
455: def ann0(F)
456: {
457: F = subst(F,s,TMP_S);
458: Ann = ann(F);
459: Bf = bfunction(F);
460:
461: FList = cdr(fctr(Bf));
462: for ( T = FList, Min = 0; T != []; T = cdr(T) ) {
463: LF = car(car(T));
464: Root = -coef(LF,0)/coef(LF,1);
465: if ( dn(Root) == 1 && Root < Min )
466: Min = Root;
467: }
468: return [Min,map(subst,Ann,s,Min,TMP_S,s,TMP_DS,ds)];
469: }
470:
471: def psi0(F,T,DT)
472: {
473: D = dp_ptod(F,[T,DT]);
474: Wmax = ww_weight(D);
475: D1 = dp_rest(D);
476: for ( ; D1; D1 = dp_rest(D1) )
477: if ( ww_weight(D1) > Wmax )
478: Wmax = ww_weight(D1);
479: for ( D1 = D, Dmax = 0; D1; D1 = dp_rest(D1) )
480: if ( ww_weight(D1) == Wmax )
481: Dmax += dp_hm(D1);
482: if ( Wmax >= 0 )
483: Dmax = dp_weyl_mul(<<Wmax,0>>,Dmax);
484: else
485: Dmax = dp_weyl_mul(<<0,-Wmax>>,Dmax);
486: Rmax = dp_dtop(Dmax,[T,DT]);
487: return Rmax;
488: }
489:
490: def psi(F,T,DT)
491: {
492: Rmax = psi0(F,T,DT);
493: R = b_subst(subst(Rmax,DT,1),T);
494: return R;
495: }
496:
497: def ww_weight(D)
498: {
499: V = dp_etov(D);
500: return V[1]-V[0];
501: }
502:
503: def compare_first(A,B)
504: {
505: A0 = car(A);
506: B0 = car(B);
507: if ( A0 > B0 )
508: return 1;
509: else if ( A0 < B0 )
510: return -1;
511: else
512: return 0;
513: }
514:
515: /* generic b-function w.r.t. weight vector W */
516:
517: def generic_bfct(F,V,DV,W)
518: {
519: N = length(V);
520: N2 = N*2;
521:
522: /* If W is a list, convert it to a vector */
523: if ( type(W) == 4 )
524: W = newvect(length(W),W);
525: dp_weyl_set_weight(W);
526:
527: /* create a term order M in D<x,d> (DRL) */
528: M = newmat(N2,N2);
529: for ( J = 0; J < N2; J++ )
530: M[0][J] = 1;
531: for ( I = 1; I < N2; I++ )
532: M[I][N2-I] = -1;
533:
534: VDV = append(V,DV);
535:
536: /* create a non-term order MW in D<x,d> */
537: MW = newmat(N2+1,N2);
538: for ( J = 0; J < N; J++ )
539: MW[0][J] = -W[J];
540: for ( ; J < N2; J++ )
541: MW[0][J] = W[J-N];
542: for ( I = 1; I <= N2; I++ )
543: for ( J = 0; J < N2; J++ )
544: MW[I][J] = M[I-1][J];
545:
546: /* create a homogenized term order MWH in D<x,d,h> */
547: MWH = newmat(N2+2,N2+1);
548: for ( J = 0; J <= N2; J++ )
549: MWH[0][J] = 1;
550: for ( I = 1; I <= N2+1; I++ )
551: for ( J = 0; J < N2; J++ )
552: MWH[I][J] = MW[I-1][J];
553:
554: /* homogenize F */
555: VDVH = append(VDV,[TMP_H]);
556: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
557:
558: /* compute a groebner basis of FH w.r.t. MWH */
559: dp_gr_flags(["Top",1,"NoRA",1]);
560: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
561: dp_gr_flags(["Top",0,"NoRA",0]);
562:
563: /* dehomigenize GH */
564: G = map(subst,GH,TMP_H,1);
565:
566: /* G is a groebner basis w.r.t. a non term order MW */
567: /* take the initial part w.r.t. (-W,W) */
568: GIN = map(initial_part,G,VDV,MW,W);
569:
570: /* GIN is a groebner basis w.r.t. a term order M */
571: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
572:
573: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
574: for ( I = 0, T = 0; I < N; I++ )
575: T += W[I]*V[I]*DV[I];
576: B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
577: return B;
578: }
579:
580: /* all term reduction + interreduce */
581: def generic_bfct_1(F,V,DV,W)
582: {
583: N = length(V);
584: N2 = N*2;
585:
586: /* If W is a list, convert it to a vector */
587: if ( type(W) == 4 )
588: W = newvect(length(W),W);
589: dp_weyl_set_weight(W);
590:
591: /* create a term order M in D<x,d> (DRL) */
592: M = newmat(N2,N2);
593: for ( J = 0; J < N2; J++ )
594: M[0][J] = 1;
595: for ( I = 1; I < N2; I++ )
596: M[I][N2-I] = -1;
597:
598: VDV = append(V,DV);
599:
600: /* create a non-term order MW in D<x,d> */
601: MW = newmat(N2+1,N2);
602: for ( J = 0; J < N; J++ )
603: MW[0][J] = -W[J];
604: for ( ; J < N2; J++ )
605: MW[0][J] = W[J-N];
606: for ( I = 1; I <= N2; I++ )
607: for ( J = 0; J < N2; J++ )
608: MW[I][J] = M[I-1][J];
609:
610: /* create a homogenized term order MWH in D<x,d,h> */
611: MWH = newmat(N2+2,N2+1);
612: for ( J = 0; J <= N2; J++ )
613: MWH[0][J] = 1;
614: for ( I = 1; I <= N2+1; I++ )
615: for ( J = 0; J < N2; J++ )
616: MWH[I][J] = MW[I-1][J];
617:
618: /* homogenize F */
619: VDVH = append(VDV,[TMP_H]);
620: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
621:
622: /* compute a groebner basis of FH w.r.t. MWH */
623: /* dp_gr_flags(["Top",1,"NoRA",1]); */
624: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
625: /* dp_gr_flags(["Top",0,"NoRA",0]); */
626:
627: /* dehomigenize GH */
628: G = map(subst,GH,TMP_H,1);
629:
630: /* G is a groebner basis w.r.t. a non term order MW */
631: /* take the initial part w.r.t. (-W,W) */
632: GIN = map(initial_part,G,VDV,MW,W);
633:
634: /* GIN is a groebner basis w.r.t. a term order M */
635: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
636:
637: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
638: for ( I = 0, T = 0; I < N; I++ )
639: T += W[I]*V[I]*DV[I];
640: B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
641: return B;
642: }
643:
644: def initial_part(F,V,MW,W)
645: {
646: N2 = length(V);
647: N = N2/2;
648: dp_ord(MW);
649: DF = dp_ptod(F,V);
650: R = dp_hm(DF);
651: DF = dp_rest(DF);
652:
653: E = dp_etov(R);
654: for ( I = 0, TW = 0; I < N; I++ )
655: TW += W[I]*(-E[I]+E[N+I]);
656: RW = TW;
657:
658: for ( ; DF; DF = dp_rest(DF) ) {
659: E = dp_etov(DF);
660: for ( I = 0, TW = 0; I < N; I++ )
661: TW += W[I]*(-E[I]+E[N+I]);
662: if ( TW == RW )
663: R += dp_hm(DF);
664: else if ( TW < RW )
665: break;
666: else
667: error("initial_part : cannot happen");
668: }
669: return dp_dtop(R,V);
670:
671: }
672:
673: /* b-function of F ? */
674:
675: def bfct(F)
676: {
677: /* XXX */
678: F = replace_vars_f(F);
679:
680: V = vars(F);
681: N = length(V);
682: D = newvect(N);
683:
684: for ( I = 0; I < N; I++ )
685: D[I] = [deg(F,V[I]),V[I]];
686: qsort(D,compare_first);
687: for ( V = [], I = 0; I < N; I++ )
688: V = cons(D[I][1],V);
689: for ( I = N-1, DV = []; I >= 0; I-- )
690: DV = cons(strtov("d"+rtostr(V[I])),DV);
691: V1 = cons(s,V); DV1 = cons(ds,DV);
692:
693: G0 = indicial1(F,reverse(V));
694: G1 = dp_weyl_gr_main(G0,append(V1,DV1),0,1,0);
695: Minipoly = weyl_minipoly(G1,append(V1,DV1),0,s);
696: return Minipoly;
697: }
698:
699: /* called from bfct() only */
700:
701: def indicial1(F,V)
702: {
703: W = append([y1,t],V);
704: N = length(V);
705: B = [t-y1*F];
706: for ( I = N-1, DV = []; I >= 0; I-- )
707: DV = cons(strtov("d"+rtostr(V[I])),DV);
708: DW = append([dy1,dt],DV);
709: for ( I = 0; I < N; I++ ) {
710: B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
711: }
712: dp_nelim(1);
713:
714: /* homogenized (heuristics) */
715: G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
716: G1 = map(subst,G0,y1,1);
717: G2 = map(psi,G1,t,dt);
718: G3 = map(subst,G2,t,-s-1);
719: return G3;
720: }
721:
722: /* b-function computation via generic_bfct() (experimental) */
723:
724: def bfct_via_gbfct(F)
725: {
726: V = vars(F);
727: N = length(V);
728: D = newvect(N);
729:
730: for ( I = 0; I < N; I++ )
731: D[I] = [deg(F,V[I]),V[I]];
732: qsort(D,compare_first);
733: for ( V = [], I = 0; I < N; I++ )
734: V = cons(D[I][1],V);
735: V = reverse(V);
736: for ( I = N-1, DV = []; I >= 0; I-- )
737: DV = cons(strtov("d"+rtostr(V[I])),DV);
738:
739: B = [TMP_T-F];
740: for ( I = 0; I < N; I++ ) {
741: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
742: }
743: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
744: W = newvect(N+1);
745: W[0] = 1;
746: R = generic_bfct(B,V1,DV1,W);
747:
748: return subst(R,s,-s-1);
749: }
750:
751: /* use an order s.t. [t,x,y,z,...,dt,dx,dy,dz,...,h] */
752:
753: def bfct_via_gbfct_weight(F,V)
754: {
755: N = length(V);
756: D = newvect(N);
757: Wt = getopt(weight);
758: if ( type(Wt) != 4 ) {
759: for ( I = 0, Wt = []; I < N; I++ )
760: Wt = cons(1,Wt);
761: }
762: Tdeg = w_tdeg(F,V,Wt);
763: WtV = newvect(2*(N+1)+1);
764: WtV[0] = Tdeg;
765: WtV[N+1] = 1;
766: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
767: for ( I = 1; I <= N; I++ ) {
768: WtV[I] = Wt[I-1];
769: WtV[N+1+I] = Tdeg-Wt[I-1]+1;
770: }
771: WtV[2*(N+1)] = 1;
772: dp_set_weight(WtV);
773: for ( I = N-1, DV = []; I >= 0; I-- )
774: DV = cons(strtov("d"+rtostr(V[I])),DV);
775:
776: B = [TMP_T-F];
777: for ( I = 0; I < N; I++ ) {
778: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
779: }
780: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
781: W = newvect(N+1);
782: W[0] = 1;
783: R = generic_bfct_1(B,V1,DV1,W);
784: dp_set_weight(0);
785: return subst(R,s,-s-1);
786: }
787:
788: /* use an order s.t. [x,y,z,...,t,dx,dy,dz,...,dt,h] */
789:
790: def bfct_via_gbfct_weight_1(F,V)
791: {
792: N = length(V);
793: D = newvect(N);
794: Wt = getopt(weight);
795: if ( type(Wt) != 4 ) {
796: for ( I = 0, Wt = []; I < N; I++ )
797: Wt = cons(1,Wt);
798: }
799: Tdeg = w_tdeg(F,V,Wt);
800: WtV = newvect(2*(N+1)+1);
801: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
802: for ( I = 0; I < N; I++ ) {
803: WtV[I] = Wt[I];
804: WtV[N+1+I] = Tdeg-Wt[I]+1;
805: }
806: WtV[N] = Tdeg;
807: WtV[2*N+1] = 1;
808: WtV[2*(N+1)] = 1;
809: dp_set_weight(WtV);
810: for ( I = N-1, DV = []; I >= 0; I-- )
811: DV = cons(strtov("d"+rtostr(V[I])),DV);
812:
813: B = [TMP_T-F];
814: for ( I = 0; I < N; I++ ) {
815: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
816: }
817: V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_DT]);
818: W = newvect(N+1);
819: W[N] = 1;
820: R = generic_bfct_1(B,V1,DV1,W);
821: dp_set_weight(0);
822: return subst(R,s,-s-1);
823: }
824:
825: def bfct_via_gbfct_weight_2(F,V)
826: {
827: N = length(V);
828: D = newvect(N);
829: Wt = getopt(weight);
830: if ( type(Wt) != 4 ) {
831: for ( I = 0, Wt = []; I < N; I++ )
832: Wt = cons(1,Wt);
833: }
834: Tdeg = w_tdeg(F,V,Wt);
835:
836: /* a weight for the first GB computation */
837: /* [t,x1,...,xn,dt,dx1,...,dxn,h] */
838: WtV = newvect(2*(N+1)+1);
839: WtV[0] = Tdeg;
840: WtV[N+1] = 1;
841: WtV[2*(N+1)] = 1;
842: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
843: for ( I = 1; I <= N; I++ ) {
844: WtV[I] = Wt[I-1];
845: WtV[N+1+I] = Tdeg-Wt[I-1]+1;
846: }
847: dp_set_weight(WtV);
848:
849: /* a weight for the second GB computation */
850: /* [x1,...,xn,t,dx1,...,dxn,dt,h] */
851: WtV2 = newvect(2*(N+1)+1);
852: WtV2[N] = Tdeg;
853: WtV2[2*N+1] = 1;
854: WtV2[2*(N+1)] = 1;
855: for ( I = 0; I < N; I++ ) {
856: WtV2[I] = Wt[I];
857: WtV2[N+1+I] = Tdeg-Wt[I]+1;
858: }
859:
860: for ( I = N-1, DV = []; I >= 0; I-- )
861: DV = cons(strtov("d"+rtostr(V[I])),DV);
862:
863: B = [TMP_T-F];
864: for ( I = 0; I < N; I++ ) {
865: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
866: }
867: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
868: V2 = append(V,[TMP_T]); DV2 = append(DV,[TMP_DT]);
869: W = newvect(N+1,[1]);
870: dp_weyl_set_weight(W);
871:
872: VDV = append(V1,DV1);
873: N1 = length(V1);
874: N2 = N1*2;
875:
876: /* create a non-term order MW in D<x,d> */
877: MW = newmat(N2+1,N2);
878: for ( J = 0; J < N1; J++ ) {
879: MW[0][J] = -W[J]; MW[0][N1+J] = W[J];
880: }
881: for ( J = 0; J < N2; J++ ) MW[1][J] = 1;
882: for ( I = 2; I <= N2; I++ ) MW[I][N2-I+1] = -1;
883:
884: /* homogenize F */
885: VDVH = append(VDV,[TMP_H]);
886: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,B,VDV)),VDVH);
887:
888: /* compute a groebner basis of FH w.r.t. MWH */
889: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
890:
891: /* dehomigenize GH */
892: G = map(subst,GH,TMP_H,1);
893:
894: /* G is a groebner basis w.r.t. a non term order MW */
895: /* take the initial part w.r.t. (-W,W) */
896: GIN = map(initial_part,G,VDV,MW,W);
897:
898: /* GIN is a groebner basis w.r.t. a term order M */
899: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
900:
901: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
902: for ( I = 0, T = 0; I < N1; I++ )
903: T += W[I]*V1[I]*DV1[I];
904:
905: /* change of ordering from VDV to VDV2 */
906: VDV2 = append(V2,DV2);
907: dp_set_weight(WtV2);
908: for ( Pind = 0; ; Pind++ ) {
909: Prime = lprime(Pind);
910: GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-Prime,0);
911: if ( GIN2 ) break;
912: }
913:
914: R = weyl_minipoly(GIN2,VDV2,0,T); /* M represents DRL order */
915: dp_set_weight(0);
916: return subst(R,s,-s-1);
917: }
918:
919: /* minimal polynomial of s; modular computation */
920:
921: def weyl_minipolym(G,V,O,M,V0)
922: {
923: N = length(V);
924: Len = length(G);
925: dp_ord(O);
926: setmod(M);
927: PS = newvect(Len);
928: PS0 = newvect(Len);
929:
930: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
931: PS0[I] = dp_ptod(car(T),V);
932: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
933: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
934:
935: for ( I = Len - 1, GI = []; I >= 0; I-- )
936: GI = cons(I,GI);
937:
938: U = dp_mod(dp_ptod(V0,V),M,[]);
939: U = dp_weyl_nf_mod(GI,U,PS,1,M);
940:
941: T = dp_mod(<<0>>,M,[]);
942: TT = dp_mod(dp_ptod(1,V),M,[]);
943: G = H = [[TT,T]];
944:
945: for ( I = 1; ; I++ ) {
946: if ( dp_gr_print() )
947: print(".",2);
948: T = dp_mod(<<I>>,M,[]);
949:
950: TT = dp_weyl_nf_mod(GI,dp_weyl_mul_mod(TT,U,M),PS,1,M);
951: H = cons([TT,T],H);
952: L = dp_lnf_mod([TT,T],G,M);
953: if ( !L[0] ) {
954: if ( dp_gr_print() )
955: print("");
956: return dp_dtop(L[1],[t]); /* XXX */
957: } else
958: G = insert(G,L);
959: }
960: }
961:
962: /* minimal polynomial of s over Q */
963:
964: def weyl_minipoly(G0,V0,O0,P)
965: {
966: HM = hmlist(G0,V0,O0);
967:
968: N = length(V0);
969: Len = length(G0);
970: dp_ord(O0);
971: PS = newvect(Len);
972: for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
973: PS[I] = dp_ptod(car(T),V0);
974: for ( I = Len - 1, GI = []; I >= 0; I-- )
975: GI = cons(I,GI);
976: PSM = newvect(Len);
977: DP = dp_ptod(P,V0);
978:
979: for ( Pind = 0; ; Pind++ ) {
980: Prime = lprime(Pind);
981: if ( !valid_modulus(HM,Prime) )
982: continue;
983: setmod(Prime);
984: for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
985: PSM[I] = dp_mod(dp_ptod(car(T),V0),Prime,[]);
986:
987: NFP = weyl_nf(GI,DP,1,PS);
988: NFPM = dp_mod(NFP[0],Prime,[])/ptomp(NFP[1],Prime);
989:
990: NF = [[dp_ptod(1,V0),1]];
991: LCM = 1;
992:
993: TM = dp_mod(<<0>>,Prime,[]);
994: TTM = dp_mod(dp_ptod(1,V0),Prime,[]);
995: GM = NFM = [[TTM,TM]];
996:
997: for ( D = 1; ; D++ ) {
998: if ( dp_gr_print() )
999: print(".",2);
1000: NFPrev = car(NF);
1001: NFJ = weyl_nf(GI,
1002: dp_weyl_mul(NFP[0],NFPrev[0]),NFP[1]*NFPrev[1],PS);
1003: NFJ = remove_cont(NFJ);
1004: NF = cons(NFJ,NF);
1005: LCM = ilcm(LCM,NFJ[1]);
1006:
1007: /* modular computation */
1008: TM = dp_mod(<<D>>,Prime,[]);
1009: TTM = dp_mod(NFJ[0],Prime,[])/ptomp(NFJ[1],Prime);
1010: NFM = cons([TTM,TM],NFM);
1011: LM = dp_lnf_mod([TTM,TM],GM,Prime);
1012: if ( !LM[0] )
1013: break;
1014: else
1015: GM = insert(GM,LM);
1016: }
1017:
1018: if ( dp_gr_print() )
1019: print("");
1020: U = NF[0][0]*idiv(LCM,NF[0][1]);
1021: Coef = [];
1022: for ( J = D-1; J >= 0; J-- ) {
1023: Coef = cons(strtov("u"+rtostr(J)),Coef);
1024: U += car(Coef)*NF[D-J][0]*idiv(LCM,NF[D-J][1]);
1025: }
1026:
1027: for ( UU = U, Eq = []; UU; UU = dp_rest(UU) )
1028: Eq = cons(dp_hc(UU),Eq);
1029: M = etom([Eq,Coef]);
1030: B = henleq(M,Prime);
1031: if ( dp_gr_print() )
1032: print("");
1033: if ( B ) {
1034: R = 0;
1035: for ( I = 0; I < D; I++ )
1036: R += B[0][I]*s^I;
1037: R += B[1]*s^D;
1038: return R;
1039: }
1040: }
1041: }
1042:
1043: def weyl_nf(B,G,M,PS)
1044: {
1045: for ( D = 0; G; ) {
1046: for ( U = 0, L = B; L != []; L = cdr(L) ) {
1047: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
1048: GCD = igcd(dp_hc(G),dp_hc(R));
1049: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
1050: U = CG*G-dp_weyl_mul(CR*dp_subd(G,R),R);
1051: if ( !U )
1052: return [D,M];
1053: D *= CG; M *= CG;
1054: break;
1055: }
1056: }
1057: if ( U )
1058: G = U;
1059: else {
1060: D += dp_hm(G); G = dp_rest(G);
1061: }
1062: }
1063: return [D,M];
1064: }
1065:
1066: def weyl_nf_quo_check(G,PS,R)
1067: {
1068: D = R[0]; M = R[1]; Coef = R[2];
1069: Len = length(PS);
1070: T = 0;
1071: for ( I = 0; I < Len; I++ )
1072: T += dp_weyl_mul(Coef[I],PS[I]);
1073: return (M*G-T)==D;
1074: }
1075:
1076: def weyl_nf_quo(B,G,M,PS)
1077: {
1078: Len = length(PS);
1079: Coef = vector(Len);
1080: for ( D = 0; G; ) {
1081: for ( U = 0, L = B; L != []; L = cdr(L) ) {
1082: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
1083: GCD = igcd(dp_hc(G),dp_hc(R));
1084: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
1085: for ( I = 0; I < Len; I++ ) Coef[I] *= CG;
1086: Q = CR*dp_subd(G,R);
1087: Coef[car(L)] += Q;
1088: U = CG*G-dp_weyl_mul(Q,R);
1089: D *= CG; M *= CG;
1090: if ( !U )
1091: return [D,M,Coef];
1092: break;
1093: }
1094: }
1095: if ( U )
1096: G = U;
1097: else {
1098: D += dp_hm(G); G = dp_rest(G);
1099: }
1100: }
1101: return [D,M,Coef];
1102: }
1103:
1104: def weyl_nf_mod(B,G,PS,Mod)
1105: {
1106: for ( D = 0; G; ) {
1107: for ( U = 0, L = B; L != []; L = cdr(L) ) {
1108: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
1109: CR = dp_hc(G)/dp_hc(R);
1110: U = G-dp_weyl_mul_mod(CR*dp_mod(dp_subd(G,R),Mod,[]),R,Mod);
1111: if ( !U )
1112: return D;
1113: break;
1114: }
1115: }
1116: if ( U )
1117: G = U;
1118: else {
1119: D += dp_hm(G); G = dp_rest(G);
1120: }
1121: }
1122: return D;
1123: }
1124:
1125: def b_subst(F,V)
1126: {
1127: D = deg(F,V);
1128: C = newvect(D+1);
1129: for ( I = D; I >= 0; I-- )
1130: C[I] = coef(F,I,V);
1131: for ( I = 0, R = 0; I <= D; I++ )
1132: if ( C[I] )
1133: R += C[I]*v_factorial(V,I);
1134: return R;
1135: }
1136:
1137: def v_factorial(V,N)
1138: {
1139: for ( J = N-1, R = 1; J >= 0; J-- )
1140: R *= V-J;
1141: return R;
1142: }
1143:
1144: def w_tdeg(F,V,W)
1145: {
1146: dp_set_weight(newvect(length(W),W));
1147: T = dp_ptod(F,V);
1148: for ( R = 0; T; T = cdr(T) ) {
1149: D = dp_td(T);
1150: if ( D > R ) R = D;
1151: }
1152: return R;
1153: }
1154:
1155: def replace_vars_f(F)
1156: {
1157: return subst(F,s,TMP_S,t,TMP_T,y1,TMP_Y1,y2,TMP_Y2);
1158: }
1159:
1160: def replace_vars_v(V)
1161: {
1162: V = replace_var(V,s,TMP_S);
1163: V = replace_var(V,t,TMP_T);
1164: V = replace_var(V,y1,TMP_Y1);
1165: V = replace_var(V,y2,TMP_Y2);
1166: return V;
1167: }
1168:
1169: def replace_var(V,X,Y)
1170: {
1171: V = reverse(V);
1172: for ( R = []; V != []; V = cdr(V) ) {
1173: Z = car(V);
1174: if ( Z == X )
1175: R = cons(Y,R);
1176: else
1177: R = cons(Z,R);
1178: }
1179: return R;
1180: }
1181:
1182: def action_on_gfs(P,V,GFS)
1183: {
1184: DP = dp_ptod(P,V);
1185: N = length(V)/2;
1186: for ( I = N-1, V0 = []; I >= 0; I-- )
1187: V0 = cons(V[I],V0);
1188: R = [];
1189: for ( T = DP; T; T = dp_rest(T) )
1190: R = cons(action_on_gfs_1(dp_hm(T),N,V0,GFS),R);
1191: D = coef(car(R)[2],0);
1192: for ( T = cdr(R); T != []; T = cdr(T) ) {
1193: Di = coef(car(T)[2],0);
1194: if ( Di < D )
1195: D = Di;
1196: }
1197: F = GFS[1];
1198: for ( T = R, G = 0; T != []; T = cdr(T) )
1199: G += car(T)[0]*F^(car(T)[2]-(s+D));
1200: while ( 1 ) {
1201: G1 = tdiv(G,F);
1202: if ( G1 ) {
1203: G = G1;
1204: D++;
1205: } else
1206: return [G,F,s+D];
1207: }
1208: }
1209:
1210: def action_on_gfs_1(M,N,V,GFS)
1211: {
1212: G = GFS[0];
1213: F = GFS[1];
1214: S = GFS[2];
1215: C = dp_hc(M);
1216: E = dp_etov(M);
1217: for ( I = 0; I < N; I++ ) {
1218: VI = V[I];
1219: C *= VI^E[I];
1220: DFVI = diff(F,VI);
1221: for ( J = 0, EI = E[I+N]; J < EI; J++, S-- )
1222: G = diff(G,VI)*F+S*G*DFVI;
1223: }
1224: return [C*G,F,S];
1225: }
1226:
1227: /* stratification */
1228:
1229: def weyl_subst(F,P,V)
1230: {
1231: VF = var(F);
1232: D = deg(F,VF);
1233: P = dp_ptod(P,V);
1234: One = dp_ptod(1,V);
1235: for ( R = 0, I = D; I >= 0; I-- )
1236: R = dp_weyl_mul(R,P)+coef(F,I,VF)*One;
1237: return dp_dtop(R,V);
1238: }
1239:
1240: def bfactor(F)
1241: {
1242: L=length(F);
1243: for(I=0,B=1;I<L;I++)B*=F[I][0]^F[I][1];
1244: return fctr(B);
1245: }
1246:
1247: def gen_a(K)
1248: {
1249: D = x^(K+1);
1250: W = [];
1251: for ( I = 1; I <= K; I++ ) {
1252: D += (V=strtov("u"+rtostr(K-I+1)))*x^(K-I);
1253: W = cons(I+1,cons(V,W));
1254: }
1255: F = res(x,D,diff(D,x));
1256: return [D,F,reverse(W)];
1257: }
1258:
1259: def gen_d(K)
1260: {
1261: D = x^2*y+y^(K-1)+u1+u2*x+u3*x^2;
1262: W = reverse([u1,2*K-2,u2,K,u3,2]);
1263: U = [u3,u2,u1];
1264: for ( I = 4; I <= K; I++ ) {
1265: D += (V=strtov("u"+rtostr(I)))*y^(I-3);
1266: W = cons((2*K-2)-2*(I-3),cons(V,W));
1267: U = cons(V,U);
1268: }
1269: B = [D,diff(D,x),diff(D,y)];
1270: G = nd_gr_trace(B,append([x,y],U),1,1,0);
1271: G = nd_gr_trace(G,append([x,y],U),1,-1,[[0,2],[0,K]]);
1272: E = elimination(G,U);
1273: F = E[0];
1274: return [D,F,reverse(W)];
1275: }
1276:
1277: def gen_dm(K)
1278: {
1279: D = x^2*y-y^(K-1)+u1+u2*x+u3*x^2;
1280: W = reverse([u1,2*K-2,u2,K,u3,2]);
1281: U = [u3,u2,u1];
1282: for ( I = 4; I <= K; I++ ) {
1283: D += (V=strtov("u"+rtostr(I)))*y^(I-3);
1284: W = cons((2*K-2)-2*(I-3),cons(V,W));
1285: U = cons(V,U);
1286: }
1287: B = [D,diff(D,x),diff(D,y)];
1288: G = nd_gr_trace(B,append([x,y],U),1,1,0);
1289: G = nd_gr_trace(G,append([x,y],U),1,-1,[[0,2],[0,K]]);
1290: E = elimination(G,U);
1291: F = E[0];
1292: return [D,F,reverse(W)];
1293: }
1294:
1295: def elimination(G,V)
1296: {
1297: ANS=[];
1298: NG=length(G);
1299:
1300: for (I=NG-1;I>=0;I--)
1301: {
1302: VSet=vars(G[I]);
1303: DIFF=setminus(VSet,V);
1304:
1305: if ( DIFF ==[] )
1306: {
1307: ANS=cons(G[I],ANS);
1308: }
1309: }
1310: return ANS;
1311: }
1312:
1313: def weyl_ideal_quotient(B,F,VDV)
1314: {
1315: T = ttttt; DT = dttttt;
1316: J = cons((1-T)*F,vtol(ltov(B)*T));
1317: N = length(VDV); N1 = N/2;
1318: for ( I = N1-1, V1 = []; I >= 0; I-- )
1319: V1 = cons(VDV[I],V1);
1320: for ( I = 0, VDV1 = VDV; I < N1; I++ ) VDV1 = cdr(VDV1);
1321: VDV1 = append(cons(T,V1),cons(DT,VDV1));
1322: O1 = [[0,1],[0,N+1]];
1323: GJ = nd_weyl_gr(J,VDV1,0,O1);
1324: R = elimination(GJ,VDV);
1325: return R;
1326: }
1327:
1328: def bf_strat(F)
1329: {
1330: dp_ord(0);
1331: T0 = time();
1332: if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
1333: if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
1334: if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
1335: L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
1336: T1 = time();
1337: print(["in_ww",(T1[0]+T1[1])-(T0[0]+T0[1])]);
1338:
1339: /* shortcuts */
1340: V = vars(F);
1341: dp_set_weight(0);
1342: dp_ord(0);
1343: Sing = sing(F,V);
1344: if ( Sing[0] == 1 || Sing[0] == -1 ) {
1345: return [[[F],[1],[[s+1,1]]],[[0],[F],[]]];
1346: } else if ( zero_dim(Sing,V,0) ) {
1347: N = length(V);
1348: P0 = [];
1349: for ( I = 0; I < N; I++ ) {
1350: M = minipoly(Sing,V,0,V[I],TMP_S);
1351: MF = cdr(fctr(M));
1352: if ( length(MF) == 1 && deg(MF[0][0],TMP_S)==1 ) {
1353: P0 = cons(subst(MF[0][0],TMP_S,V[I]),P0);
1354: } else break;
1355: }
1356: if ( I == N ) {
1357: Indata = L[0]; AllData = L[1]; VData = L[2];
1358: GIN = Indata[0]; VDV = Indata[1]; WVDV = AllData[4];
1359: W = Indata[4];
1360: dp_set_weight(W);
1361: B = weyl_minipoly(GIN,VDV,0,WVDV);
1362: B = subst(B,s,-s-1);
1363: dp_set_weight(0);
1364: return [
1365: [P0,[1],cdr(fctr(B))],
1366: [[F],P0,[[s+1,1]]],
1367: [[0],[F],[]]
1368: ];
1369: }
1370: }
1371:
1372: L2 = bf_strat_stage2(L);
1373: S = bf_strat_stage3(L2);
1374: R = [];
1375: for ( T = S; T != []; T = cdr(T) ) {
1376: Str = car(T);
1377: B = Str[2];
1378: B1 = [];
1379: for ( U = B; U != []; U = cdr(U) )
1380: B1 = cons([subst(car(U)[0],s,-s-1),car(U)[1]],B1);
1381: B1 = reverse(B1);
1382: R = cons([Str[0],Str[1],B1],R);
1383: }
1384: return reverse(R);
1385: }
1386:
1387: /* returns [QQ,V,V0,B,BF,W] */
1388: /* QQ : ideal in C[x,s] (s=tdt), V=[x1,..,xn,t], V0 = [x1,..,xn] */
1389: /* B : global b-function, BF : factor list of B, W : weight */
1390:
1391: def bf_strat_stage2(L)
1392: {
1393: T0 = time();
1394: InData = L[0]; VData = L[2];
1395: G1 = InData[0]; VDV = InData[1]; W = InData[4]; W0 = VData[4];
1396: N = length(VDV); N1 = N/2;
1397: V = InData[2]; DV = InData[3];
1398: T = VData[2]; DT = VData[3];
1399: V0 = VData[0]; DVR = VData[1];
1400: dp_set_weight(W);
1401: for ( I = 0; DVR != []; I++, DVR = cdr(DVR) ) {
1402: DVRV = cons(DT,append(cdr(DVR),V));
1403: M = elim_mat(VDV,DVRV);
1404: for ( K = 0; K < N; K++ )
1405: M[1][K] = W[K];
1.2 noro 1406: dp_ord(0); D1 = map(dp_ptod,G1,VDV);
1407: H1 = map(dp_ht,D1); HC1 = map(dp_hc,D1);
1.1 noro 1408: dp_ord(M); H2 = map(dp_ht,map(dp_ptod,G1,VDV));
1409: if ( H1 == H2 )
1410: G2 = G1;
1411: else
1.2 noro 1412: G2 = nd_gb_candidate(G1,VDV,M,0,HC1,1);
1.1 noro 1413: G1 = elimination(G2,DVRV);
1414: }
1415: T1 = time();
1416: B = weyl_minipoly(G1,VDV,0,T*DT);
1417: T2 = time();
1418: BF = cdr(fctr(B));
1419:
1420: dp_set_weight(0);
1421: G1 = map(psi0,G1,T,DT);
1422: QQ = map(subst,map(b_subst,map(subst,G1,DT,1),T),T,var(B));
1423: if ( type(getopt(ideal)) != -1 ) return [QQ,V];
1424: print(["elim",(T1[0]+T1[1])-(T0[0]+T0[1])]);
1425: print(["globalb",(T2[0]+T2[1])-(T1[0]+T1[1])]);
1426: return [QQ,V,V0,B,BF,W0,DV];
1427: }
1428:
1429: def bf_strat_stage3(L)
1430: {
1431: T0 = time();
1432: QQ = L[0]; V0 = L[2]; B = L[3]; BF = L[4]; W0 = L[5];
1433: NF = length(BF);
1434: Data = vector(NF);
1.2 noro 1435: W1 = W0? cons(1,append(W0,[1])) : 0;
1.1 noro 1436: for ( I = J = 0; I < NF; I++ ) {
1437: DI = tower_in_p(QQ,B,BF[I],V0,W0);
1438: NDI = length(DI);
1.2 noro 1439: dp_set_weight(W1);
1.1 noro 1440: for ( K = 0; K < J; K++ ) {
1441: if ( length(DK=Data[K]) == NDI ) {
1442: for ( L = 0; L < NDI; L++ ) {
1443: CL = DI[L][1]; CK = DK[L][1];
1444: if ( !zero_inclusion(CL,CK,V0)
1445: || !zero_inclusion(CK,CL,V0) ) break;
1446: }
1447: if ( L == NDI ) break;
1448: }
1449: }
1.2 noro 1450: dp_set_weight(0);
1.1 noro 1451: if ( K < J ) {
1452: for ( L = 0, T = []; L < NDI; L++ )
1453: T = cons([[DK[L][0][0]*DI[L][0][0],DK[L][0][1]],
1454: DK[L][1],DK[L][2]],T);
1455: Data[K] = reverse(T);
1456: } else
1457: Data[J++] = DI;
1458: }
1459: Data1 = vector(J);
1460: for ( I = 0; I < J; I++ )
1461: Data1[I] = Data[I];
1462: T1 = time();
1.2 noro 1463: Str = stratify_bf(Data1,V0,W0);
1.1 noro 1464: T2 = time();
1465: print(["tower",(T1[0]+T1[1])-(T0[0]+T0[1])]);
1466: print(["strat",(T2[0]+T2[1])-(T1[0]+T1[1])]);
1467: return Str;
1468: }
1469:
1470: /*
1471: InData = [GIN,VDV,V,DV,WtV]
1472: AllData = [G0,GIN0,VDV0,W,WVDV,WtV0]
1473: */
1474:
1475: def bf_local(F,P)
1476: {
1.6 noro 1477: /* F -> F/Fcont */
1478: F1 = ptozp(F); Fcont = sdiv(F,F1); F = F1;
1.1 noro 1479: if ( type(Heu=getopt(heuristic)) == -1 ) Heu = 0;
1480: if ( type(Vord=getopt(vord)) == -1 || type(Vord) != 4 ) Vord = 0;
1481: if ( type(Wt=getopt(weight)) == -1 ) Wt = 0;
1.3 noro 1482: if ( type(Op=getopt(op)) == -1 ) Op = 0;
1.1 noro 1483: L = in_ww(F|weight=Wt,heuristic=Heu,vord=Vord);
1484: InData = L[0]; AllData = L[1]; VData = L[2];
1485: G = InData[0]; VDV = InData[1];
1486: V = InData[2]; DV = InData[3];
1487:
1488: V0 = VData[0]; DV0 = VData[1]; T = VData[2]; DT = VData[3]; W0 = VData[4];
1489:
1490: L2 = bf_strat_stage2(L);
1491:
1492: /* L2 = [QQ,V,V0,B,BF,W] */
1493: /* QQ : ideal in C[x,s] (s=tdt), V=[t,x1,..,xn], V0 = [x1,..,xn] */
1494: /* B : global b-function, BF : factor list of B, W : weight */
1495:
1496: QQ = L2[0]; B = L2[3]; BF = L2[4]; W = L2[5];
1497:
1498: NF = length(BF);
1499: BP = [];
1500: dp_set_weight(0);
1501: for ( I = J = 0; I < NF; I++ ) {
1502: List = compute_exponent(QQ,B,BF[I],P,V0,W0);
1503: DI = List[0]; QQI = List[1];
1504: if ( DI )
1505: BP = cons([BF[I][0],DI],BP);
1506: if ( I == 0 )
1507: Id = QQI;
1508: else
1509: Id = ideal_intersection(Id,QQI,V0,0);
1510: }
1511: for ( List = Id; List != []; List = cdr(List) )
1512: if ( subst_vars(car(List),P) )
1513: break;
1514: if ( List == [] ) error("bf_local : inconsitent intersection");
1515: Ax = car(List);
1.3 noro 1516: LB = [];
1517: for ( BPT = 1, List = BP; List != []; List = cdr(List) ) {
1.1 noro 1518: BPT *= car(List)[0]^car(List)[1];
1.3 noro 1519: LB = cons([subst(car(List)[0],s,-s-1),car(List)[1]],LB);
1520: }
1521: LB = reverse(LB);
1522: if ( !Op ) return LB;
1523:
1.1 noro 1524: BPT = weyl_subst(BPT,T*DT,VDV);
1525:
1526: /* computation using G0,GIN0,VDV0 */
1527: G0 = AllData[0]; GIN0 = AllData[1]; VDV0 = AllData[2]; WtV0 = AllData[5];
1528: dp_set_weight(WtV0); dp_ord(0);
1529: PS = map(dp_ptod,GIN0,VDV0); Len = length(PS);
1530: for ( I = Len-1, Ind = []; I >= 0; I-- ) Ind = cons(I,Ind);
1531: /* QR = [D,M,Coef] */
1532: AxBPT = dp_ptod(Ax*BPT,VDV0);
1533: QR = weyl_nf_quo(Ind,AxBPT,1,PS);
1534: if ( !weyl_nf_quo_check(AxBPT,PS,QR) ) error("bf_local : invalid quotient");
1535: if ( QR[0] ) error("bf_local : invalid quotient");
1536: Den = QR[1]; Coef = QR[2];
1537: for ( I = 0, R = Den*AxBPT; I < Len; I++ )
1538: R -= dp_weyl_mul(Coef[I],dp_ptod(G0[I],VDV0));
1539: R = dp_dtop(R,VDV0);
1540: CR = conv_tdt(R,F,V0,DV0,T,DT);
1541:
1542: dp_set_weight(0);
1.6 noro 1543: Cont = cont(CR); CR /= Cont;
1544: Cont *= dn(Fcont); Den *= nm(Fcont);
1.5 noro 1545: Gcd = igcd(Den,Cont);
1.6 noro 1546: return [LB,(Den/Gcd)*Ax,(Cont/Gcd)*CR];
1.1 noro 1547: }
1548:
1549: /* t^(l+k)*dt^l (k>l) -> (s-k)(s-k-1)...(s-(k+l-1))t^k */
1550: def conv_tdt(P,F,V0,DV0,T,DT)
1551: {
1552: DP = dp_ptod(P,[T,DT]);
1553: VDV = append(cons(T,V0),cons(DT,DV0));
1554: R = 0;
1555: DF = dp_ptod(F,VDV);
1556: for ( ; DP; DP = dp_rest(DP) ) {
1557: C = dp_hc(DP);
1558: E = dp_etov(dp_ht(DP));
1559: L = E[1]; K = E[0]-E[1];
1560: S = 1;
1561: for ( I = 0; I < L; I++ )
1562: S *= ((-T-1)-K-I);
1563: U = dp_ptod(C*S,VDV);
1564: for ( I = 1; I < K; I++ )
1565: U = dp_weyl_mul(U,DF);
1566: R += dp_dtop(U,VDV);
1567: }
1568: return subst(R,T,s);
1569: }
1570:
1.2 noro 1571: /* W1=[W,1], W2=[1,W,1] */
1572:
1573: def merge_tower(Str,Tower,V,W1,W2)
1.1 noro 1574: {
1575: Prev = car(Tower); T = cdr(Tower);
1576: Str1 = [];
1577: for ( ; T != []; T = cdr(T) ) {
1578: Cur = car(T);
1579: Str1 = cons([Cur[1],Prev[1],[Prev[0]]],Str1);
1580: Prev = Cur;
1581: }
1582: Str1 = cons([[0],Prev[1],[]],Str1);
1583: Str1 = reverse(Str1);
1584: if ( Str == [] ) return Str1;
1585: U = [];
1586: for ( T = Str; T != []; T = cdr(T) ) {
1587: for ( S = Str1; S != []; S = cdr(S) ) {
1.2 noro 1588: Int = int_cons(T[0],S[0],V,W1,W2);
1.1 noro 1589: if ( Int[0] != [1] )
1590: U = cons(append(Int,[append(T[0][2],S[0][2])]),U);
1591: }
1592: }
1593: U = reverse(U);
1594: return U;
1595: }
1596:
1.2 noro 1597: def stratify_bf(Data,V,W)
1.1 noro 1598: {
1599: N = length(Data);
1600: R = [];
1.2 noro 1601: if ( W ) {
1602: W1 = append(W,[1]);
1603: W2 = cons(1,W1);
1604: } else
1605: W1 = W2 = 0;
1.1 noro 1606: for ( I = 0; I < N; I++ )
1.2 noro 1607: R = merge_tower(R,Data[I],V,W1,W2);
1.1 noro 1608: return R;
1609: }
1610:
1611: def tdt_homogenize(F,L)
1612: {
1613: TY1 = L[0]; T = TY1[0]; Y1 = TY1[1];
1614: TY2 = L[1]; DT = TY2[0]; Y2 = TY2[1];
1615: DF = dp_ptod(F,[T,DT]);
1616: for ( R = 0; DF; DF = dp_rest(DF) ) {
1617: M = dp_hm(DF);
1618: E = dp_etov(M);
1619: W = E[1]-E[0];
1620: if ( W > 0 ) R += Y1^W*dp_dtop(M,[T,DT]);
1621: else R += Y2^W*dp_dtop(M,[T,DT]);
1622: }
1623: return R;
1624: }
1625:
1626: def sing(F,V)
1627: {
1628: R = [F];
1629: for ( T = V; T != []; T = cdr(T) )
1630: R = cons(diff(F,car(T)),R);
1631: G = nd_gr_trace(R,V,1,1,0);
1632: return G;
1633: }
1634:
1635: def tower_in_p(B,F,FD,V,W)
1636: {
1637: TT = ttttt;
1638: N = length(V); S = var(F); SV = cons(S,V); V1 = cons(TT,SV);
1639: Wt = append(append([1,1],W),[1]);
1640: dp_set_weight(Wt);
1641:
1642: F1 = FD[0]; D = FD[1];
1643: O1 = [[0,1],[0,N+1]]; O2 = [[0,1],[0,N]];
1644:
1645: TF = sdiv(F,F1^D);
1646:
1647: T = nd_gr_trace(cons((1-TT)*TF,vtol(TT*ltov(B))),V1,1,1,O1);
1648: T = elimination(T,SV);
1649: Q = map(sdiv,T,TF);
1650: dp_set_weight(cdr(Wt));
1651: QQ = elimination(nd_gr(Q,SV,0,O2),V);
1652: E = [[[F1,0],QQ,PD]];
1653:
1654: for ( I = D-1; I >= 0; I-- ) {
1655: dp_set_weight(Wt);
1656: T = nd_gr_trace(cons((1-TT)*F1,vtol(TT*ltov(Q))),V1,1,1,O1);
1657: T = elimination(T,SV);
1658: Q = map(sdiv,T,F1);
1659: dp_set_weight(cdr(Wt));
1660: QQ = elimination(nd_gr(Q,SV,0,O2),V);
1661: E = cons([[F1,D-I],QQ,PD],E);
1662: }
1663: dp_set_weight(0);
1664: return E;
1665: }
1666:
1667: def subst_vars(F,P)
1668: {
1669: for ( ; P != []; P = cdr(cdr(P)) )
1670: F = subst(F,P[0],P[1]);
1671: return F;
1672: }
1673:
1674: def compute_exponent(B,F,FD,P,V,W)
1675: {
1676: TT = ttttt;
1677: N = length(V); S = var(F); SV = cons(S,V); V1 = cons(TT,SV);
1678: F1 = FD[0]; D = FD[1];
1679:
1680: Wt = append(append([1,1],W),[1]);
1681: dp_set_weight(Wt);
1682: O1 = [[0,1],[0,N+1]]; O2 = [[0,1],[0,N]];
1683:
1684: TF = sdiv(F,F1^D);
1685:
1686: T = nd_gr_trace(cons((1-TT)*TF,vtol(TT*ltov(B))),V1,0,1,O1);
1687: T = elimination(T,SV);
1688: Q = map(sdiv,T,TF);
1689: QQ = elimination(nd_gr(Q,SV,0,O2),V);
1690:
1691: for ( I = 0; I < D; I++ ) {
1692: for ( T = QQ; T != []; T = cdr(T) )
1693: if ( subst_vars(car(T),P) ) {
1694: dp_set_weight(0);
1695: return [I,QQ];
1696: }
1697: T = nd_gr_trace(cons((1-TT)*F1,vtol(TT*ltov(Q))),V1,0,1,O1);
1698: T = elimination(T,SV);
1699: Q = map(sdiv,T,F1);
1700: QQ = elimination(nd_gr(Q,SV,0,O2),V);
1701: }
1702: dp_set_weight(0);
1703: return [D,QQ];
1704: }
1705:
1706: /* V(B) subset V(A) ? */
1707:
1708: def zero_inclusion(A,B,V)
1709: {
1710: NV = ttttt;
1711: for ( T = A; T != []; T = cdr(T) ) {
1712: F = car(T);
1.2 noro 1713: G = nd_gr_trace(cons(NV*F-1,B),cons(NV,V),1,1,0);
1.1 noro 1714: if ( type(car(G)) != 1 ) return 0;
1715: }
1716: return 1;
1717: }
1718:
1719: def weyl_divide_by_right(G,H,V,O)
1720: {
1721: dp_ord(O); G = dp_ptod(G,V); H = dp_ptod(H,V);
1722: CH = dp_hc(H);
1723: for ( Q = 0; G; ) {
1724: if ( !dp_redble(G,H) ) return 0;
1725: CG = dp_hc(G);
1726: Q1 = CG/CH*dp_subd(G,H);
1727: G -= dp_weyl_mul(Q1,H);
1728: Q += Q1;
1729: }
1730: return dp_dtop(Q,V);
1731: }
1732:
1733: def elim_mat(V,W)
1734: {
1735: N = length(V);
1736: M = matrix(N+1,N);
1737: for ( J = 0; J < N; J++ ) if ( !member(V[J],W) ) M[0][J] = 1;
1738: for ( J = 0; J < N; J++ ) M[1][J] = 1;
1739: for ( I = 2; I <= N; I++ ) M[I][N-I+1] = -1;
1740: return M;
1741: }
1742:
1743: /* (P-Q)cap(R-S)=(P cap Q^c)cap(R cap S^c)=(P cap R)cap(Q cup S)^c
1744: =(P cap R)-(Q cup S)
1745: */
1746:
1.2 noro 1747: def int_cons(A,B,V,W1,W2)
1.1 noro 1748: {
1749: AZ = A[0]; AN = A[1];
1750: BZ = B[0]; BN = B[1];
1.2 noro 1751: if ( W1 ) dp_set_weight(W1);
1.1 noro 1752: CZ = nd_gr_trace(append(AZ,BZ),V,1,1,0);
1.2 noro 1753: if ( W2 ) dp_set_weight(W2);
1.1 noro 1754: CN = ideal_intersection(AN,BN,V,0);
1.2 noro 1755: ZI = zero_inclusion(CN,CZ,V);
1756: dp_set_weight(0);
1757: if ( ZI )
1.1 noro 1758: return [[1],[]];
1759: else
1760: return [CZ,CN];
1761: }
1762:
1763: def ideal_intersection(A,B,V,Ord)
1764: {
1765: T = ttttt;
1766: G = nd_gr_trace(append(vtol(ltov(A)*T),vtol(ltov(B)*(1-T))),
1767: cons(T,V),1,1,[[0,1],[Ord,length(V)]]);
1768: return elimination(G,V);
1769: }
1.2 noro 1770:
1771: def nd_gb_candidate(G,V,Ord,Homo,HC,Weyl)
1772: {
1773: Ind = 0;
1774: N = length(HC);
1775: while ( 1 ) {
1776: P = lprime(Ind++);
1777: for ( I = 0; I < N; I++ )
1778: if ( !(HC[I]%P) ) break;
1779: if ( I < N ) continue;
1780: if ( Weyl )
1781: G = nd_weyl_gr_trace(G,V,Homo,-P,Ord);
1782: else
1783: G = nd_gr_trace(G,V,Homo,-P,Ord);
1784: if ( G ) return G;
1785: }
1786: }
1787:
1.1 noro 1788: endmodule $
1789: end$
1790:
1791: