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