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