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