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