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