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