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