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