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