Annotation of OpenXM_contrib2/asir2000/lib/bfct, Revision 1.24
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.24 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.23 2003/04/20 11:59:57 noro Exp $
1.10 noro 49: */
1.1 noro 50: /* requires 'primdec' */
1.22 noro 51:
1.23 noro 52: #define TMP_S ssssssss
1.24 ! noro 53: #define TMP_DS dssssssss
! 54: #define TMP_T dtttttttt
! 55: #define TMP_DT tttttttt
1.23 noro 56: #define TMP_Y1 yyyyyyyy1
1.24 ! noro 57: #define TMP_DY1 dyyyyyyyy1
1.23 noro 58: #define TMP_Y2 yyyyyyyy2
1.24 ! noro 59: #define TMP_DY2 dyyyyyyyy2
1.23 noro 60:
1.22 noro 61: extern LIBRARY_GR_LOADED$
62: extern LIBRARY_PRIMDEC_LOADED$
63:
64: if(!LIBRARY_GR_LOADED) load("gr"); else ; LIBRARY_GR_LOADED = 1$
65: if(!LIBRARY_PRIMDEC_LOADED) load("primdec"); else ; LIBRARY_PRIMDEC_LOADED = 1$
66:
67: /* toplevel */
68:
69: def bfunction(F)
70: {
71: V = vars(F);
72: N = length(V);
73: D = newvect(N);
74:
75: for ( I = 0; I < N; I++ )
76: D[I] = [deg(F,V[I]),V[I]];
77: qsort(D,compare_first);
78: for ( V = [], I = 0; I < N; I++ )
79: V = cons(D[I][1],V);
80: return bfct_via_gbfct_weight(F,V);
81: }
1.1 noro 82:
1.6 noro 83: /* annihilating ideal of F^s */
1.1 noro 84:
85: def ann(F)
86: {
1.24 ! noro 87: if ( member(s,vars(F)) )
! 88: error("ann : the variable 's' is reserved.");
1.1 noro 89: V = vars(F);
90: N = length(V);
1.8 noro 91: D = newvect(N);
92:
93: for ( I = 0; I < N; I++ )
94: D[I] = [deg(F,V[I]),V[I]];
95: qsort(D,compare_first);
96: for ( V = [], I = N-1; I >= 0; I-- )
97: V = cons(D[I][1],V);
98:
1.1 noro 99: for ( I = N-1, DV = []; I >= 0; I-- )
100: DV = cons(strtov("d"+rtostr(V[I])),DV);
1.8 noro 101:
1.24 ! noro 102: W = append([TMP_Y1,TMP_Y2,TMP_T],V);
! 103: DW = append([TMP_DY1,TMP_DY2,TMP_DT],DV);
1.8 noro 104:
1.24 ! noro 105: B = [1-TMP_Y1*TMP_Y2,TMP_T-TMP_Y1*F];
1.1 noro 106: for ( I = 0; I < N; I++ ) {
1.24 ! noro 107: B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
1.1 noro 108: }
1.10 noro 109:
110: /* homogenized (heuristics) */
1.1 noro 111: dp_nelim(2);
1.10 noro 112: G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
1.1 noro 113: G1 = [];
114: for ( T = G0; T != []; T = cdr(T) ) {
115: E = car(T); VL = vars(E);
1.24 ! noro 116: if ( !member(TMP_Y1,VL) && !member(TMP_Y2,VL) )
1.1 noro 117: G1 = cons(E,G1);
118: }
1.24 ! noro 119: G2 = map(psi,G1,TMP_T,TMP_DT);
! 120: G3 = map(subst,G2,TMP_T,-1-s);
1.12 noro 121: return G3;
1.1 noro 122: }
123:
1.10 noro 124: /*
125: * compute J_f|s=r, where r = the minimal integral root of global b_f(s)
126: * ann0(F) returns [MinRoot,Ideal]
127: */
128:
129: def ann0(F)
130: {
131: V = vars(F);
132: N = length(V);
133: D = newvect(N);
134:
135: for ( I = 0; I < N; I++ )
136: D[I] = [deg(F,V[I]),V[I]];
137: qsort(D,compare_first);
138: for ( V = [], I = 0; I < N; I++ )
139: V = cons(D[I][1],V);
140:
141: for ( I = N-1, DV = []; I >= 0; I-- )
142: DV = cons(strtov("d"+rtostr(V[I])),DV);
143:
144: /* XXX : heuristics */
1.24 ! noro 145: W = append([TMP_Y1,TMP_Y2,TMP_T],reverse(V));
! 146: DW = append([TMP_DY1,TMP_DY2,TMP_DT],reverse(DV));
1.10 noro 147: WDW = append(W,DW);
148:
1.24 ! noro 149: B = [1-TMP_Y1*TMP_Y2,TMP_T-TMP_Y1*F];
1.10 noro 150: for ( I = 0; I < N; I++ ) {
1.24 ! noro 151: B = cons(DV[I]+TMP_Y1*diff(F,V[I])*TMP_DT,B);
1.10 noro 152: }
153:
154: /* homogenized (heuristics) */
155: dp_nelim(2);
156: G0 = dp_weyl_gr_main(B,WDW,1,0,6);
157: G1 = [];
158: for ( T = G0; T != []; T = cdr(T) ) {
159: E = car(T); VL = vars(E);
1.24 ! noro 160: if ( !member(TMP_Y1,VL) && !member(TMP_Y2,VL) )
1.10 noro 161: G1 = cons(E,G1);
162: }
1.24 ! noro 163: G2 = map(psi,G1,TMP_T,TMP_DT);
! 164: G3 = map(subst,G2,TMP_T,-1-TMP_S);
1.10 noro 165:
1.12 noro 166: /* G3 = J_f(s) */
1.10 noro 167:
1.24 ! noro 168: V1 = cons(TMP_S,V); DV1 = cons(TMP_DS,DV); V1DV1 = append(V1,DV1);
1.12 noro 169: G4 = dp_weyl_gr_main(cons(F,G3),V1DV1,0,1,0);
1.24 ! noro 170: Bf = weyl_minipoly(G4,V1DV1,0,TMP_S);
1.10 noro 171:
172: FList = cdr(fctr(Bf));
173: for ( T = FList, Min = 0; T != []; T = cdr(T) ) {
174: LF = car(car(T));
175: Root = -coef(LF,0)/coef(LF,1);
176: if ( dn(Root) == 1 && Root < Min )
177: Min = Root;
178: }
1.24 ! noro 179: return [Min,map(subst,G3,TMP_S,Min)];
1.6 noro 180: }
181:
182: def psi(F,T,DT)
183: {
184: D = dp_ptod(F,[T,DT]);
185: Wmax = weight(D);
186: D1 = dp_rest(D);
187: for ( ; D1; D1 = dp_rest(D1) )
188: if ( weight(D1) > Wmax )
189: Wmax = weight(D1);
190: for ( D1 = D, Dmax = 0; D1; D1 = dp_rest(D1) )
191: if ( weight(D1) == Wmax )
192: Dmax += dp_hm(D1);
193: if ( Wmax >= 0 )
194: Dmax = dp_weyl_mul(<<Wmax,0>>,Dmax);
195: else
196: Dmax = dp_weyl_mul(<<0,-Wmax>>,Dmax);
197: Rmax = dp_dtop(Dmax,[T,DT]);
198: R = b_subst(subst(Rmax,DT,1),T);
199: return R;
200: }
201:
202: def weight(D)
203: {
204: V = dp_etov(D);
205: return V[1]-V[0];
206: }
207:
208: def compare_first(A,B)
209: {
210: A0 = car(A);
211: B0 = car(B);
212: if ( A0 > B0 )
213: return 1;
214: else if ( A0 < B0 )
215: return -1;
216: else
217: return 0;
218: }
219:
1.13 noro 220: /* generic b-function w.r.t. weight vector W */
221:
222: def generic_bfct(F,V,DV,W)
223: {
224: N = length(V);
225: N2 = N*2;
226:
1.16 noro 227: /* If W is a list, convert it to a vector */
228: if ( type(W) == 4 )
229: W = newvect(length(W),W);
1.15 noro 230: dp_weyl_set_weight(W);
231:
1.14 noro 232: /* create a term order M in D<x,d> (DRL) */
1.13 noro 233: M = newmat(N2,N2);
234: for ( J = 0; J < N2; J++ )
235: M[0][J] = 1;
236: for ( I = 1; I < N2; I++ )
237: M[I][N2-I] = -1;
238:
239: VDV = append(V,DV);
240:
241: /* create a non-term order MW in D<x,d> */
242: MW = newmat(N2+1,N2);
243: for ( J = 0; J < N; J++ )
244: MW[0][J] = -W[J];
245: for ( ; J < N2; J++ )
246: MW[0][J] = W[J-N];
247: for ( I = 1; I <= N2; I++ )
248: for ( J = 0; J < N2; J++ )
249: MW[I][J] = M[I-1][J];
250:
251: /* create a homogenized term order MWH in D<x,d,h> */
252: MWH = newmat(N2+2,N2+1);
253: for ( J = 0; J <= N2; J++ )
254: MWH[0][J] = 1;
255: for ( I = 1; I <= N2+1; I++ )
256: for ( J = 0; J < N2; J++ )
257: MWH[I][J] = MW[I-1][J];
258:
259: /* homogenize F */
260: VDVH = append(VDV,[h]);
261: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
262:
263: /* compute a groebner basis of FH w.r.t. MWH */
1.21 noro 264: dp_gr_flags(["Top",1,"NoRA",1]);
1.15 noro 265: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
1.21 noro 266: dp_gr_flags(["Top",0,"NoRA",0]);
1.13 noro 267:
268: /* dehomigenize GH */
269: G = map(subst,GH,h,1);
270:
271: /* G is a groebner basis w.r.t. a non term order MW */
272: /* take the initial part w.r.t. (-W,W) */
273: GIN = map(initial_part,G,VDV,MW,W);
274:
275: /* GIN is a groebner basis w.r.t. a term order M */
276: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
277:
278: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
279: for ( I = 0, T = 0; I < N; I++ )
280: T += W[I]*V[I]*DV[I];
1.14 noro 281: B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
1.13 noro 282: return B;
283: }
284:
1.18 noro 285: /* all term reduction + interreduce */
286: def generic_bfct_1(F,V,DV,W)
287: {
288: N = length(V);
289: N2 = N*2;
290:
291: /* If W is a list, convert it to a vector */
292: if ( type(W) == 4 )
293: W = newvect(length(W),W);
294: dp_weyl_set_weight(W);
295:
296: /* create a term order M in D<x,d> (DRL) */
297: M = newmat(N2,N2);
298: for ( J = 0; J < N2; J++ )
299: M[0][J] = 1;
300: for ( I = 1; I < N2; I++ )
301: M[I][N2-I] = -1;
302:
303: VDV = append(V,DV);
304:
305: /* create a non-term order MW in D<x,d> */
306: MW = newmat(N2+1,N2);
307: for ( J = 0; J < N; J++ )
308: MW[0][J] = -W[J];
309: for ( ; J < N2; J++ )
310: MW[0][J] = W[J-N];
311: for ( I = 1; I <= N2; I++ )
312: for ( J = 0; J < N2; J++ )
313: MW[I][J] = M[I-1][J];
314:
315: /* create a homogenized term order MWH in D<x,d,h> */
316: MWH = newmat(N2+2,N2+1);
317: for ( J = 0; J <= N2; J++ )
318: MWH[0][J] = 1;
319: for ( I = 1; I <= N2+1; I++ )
320: for ( J = 0; J < N2; J++ )
321: MWH[I][J] = MW[I-1][J];
322:
323: /* homogenize F */
324: VDVH = append(VDV,[h]);
325: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
326:
327: /* compute a groebner basis of FH w.r.t. MWH */
328: /* dp_gr_flags(["Top",1,"NoRA",1]); */
329: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
330: /* dp_gr_flags(["Top",0,"NoRA",0]); */
331:
332: /* dehomigenize GH */
333: G = map(subst,GH,h,1);
334:
335: /* G is a groebner basis w.r.t. a non term order MW */
336: /* take the initial part w.r.t. (-W,W) */
337: GIN = map(initial_part,G,VDV,MW,W);
338:
339: /* GIN is a groebner basis w.r.t. a term order M */
340: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
341:
342: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
343: for ( I = 0, T = 0; I < N; I++ )
344: T += W[I]*V[I]*DV[I];
345: B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
346: return B;
347: }
348:
1.13 noro 349: def initial_part(F,V,MW,W)
350: {
351: N2 = length(V);
352: N = N2/2;
353: dp_ord(MW);
354: DF = dp_ptod(F,V);
355: R = dp_hm(DF);
356: DF = dp_rest(DF);
357:
358: E = dp_etov(R);
359: for ( I = 0, TW = 0; I < N; I++ )
360: TW += W[I]*(-E[I]+E[N+I]);
361: RW = TW;
362:
363: for ( ; DF; DF = dp_rest(DF) ) {
364: E = dp_etov(DF);
365: for ( I = 0, TW = 0; I < N; I++ )
366: TW += W[I]*(-E[I]+E[N+I]);
367: if ( TW == RW )
368: R += dp_hm(DF);
369: else if ( TW < RW )
370: break;
371: else
372: error("initial_part : cannot happen");
373: }
374: return dp_dtop(R,V);
375:
376: }
377:
1.1 noro 378: /* b-function of F ? */
379:
380: def bfct(F)
381: {
1.23 noro 382: /* XXX */
383: F = replace_vars_f(F);
384:
1.1 noro 385: V = vars(F);
386: N = length(V);
1.6 noro 387: D = newvect(N);
1.7 noro 388:
1.6 noro 389: for ( I = 0; I < N; I++ )
390: D[I] = [deg(F,V[I]),V[I]];
391: qsort(D,compare_first);
392: for ( V = [], I = 0; I < N; I++ )
393: V = cons(D[I][1],V);
1.1 noro 394: for ( I = N-1, DV = []; I >= 0; I-- )
395: DV = cons(strtov("d"+rtostr(V[I])),DV);
1.6 noro 396: V1 = cons(s,V); DV1 = cons(ds,DV);
1.7 noro 397:
398: G0 = indicial1(F,reverse(V));
399: G1 = dp_weyl_gr_main(G0,append(V1,DV1),0,1,0);
400: Minipoly = weyl_minipoly(G1,append(V1,DV1),0,s);
1.6 noro 401: return Minipoly;
402: }
403:
1.24 ! noro 404: /* called from bfct() only */
! 405:
! 406: def indicial1(F,V)
! 407: {
! 408: W = append([y1,t],V);
! 409: N = length(V);
! 410: B = [t-y1*F];
! 411: for ( I = N-1, DV = []; I >= 0; I-- )
! 412: DV = cons(strtov("d"+rtostr(V[I])),DV);
! 413: DW = append([dy1,dt],DV);
! 414: for ( I = 0; I < N; I++ ) {
! 415: B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
! 416: }
! 417: dp_nelim(1);
! 418:
! 419: /* homogenized (heuristics) */
! 420: G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
! 421: G1 = map(subst,G0,y1,1);
! 422: G2 = map(psi,G1,t,dt);
! 423: G3 = map(subst,G2,t,-s-1);
! 424: return G3;
! 425: }
! 426:
1.14 noro 427: /* b-function computation via generic_bfct() (experimental) */
428:
429: def bfct_via_gbfct(F)
430: {
431: V = vars(F);
432: N = length(V);
433: D = newvect(N);
434:
435: for ( I = 0; I < N; I++ )
436: D[I] = [deg(F,V[I]),V[I]];
437: qsort(D,compare_first);
438: for ( V = [], I = 0; I < N; I++ )
439: V = cons(D[I][1],V);
440: V = reverse(V);
441: for ( I = N-1, DV = []; I >= 0; I-- )
442: DV = cons(strtov("d"+rtostr(V[I])),DV);
443:
1.24 ! noro 444: B = [TMP_T-F];
1.14 noro 445: for ( I = 0; I < N; I++ ) {
1.24 ! noro 446: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
1.14 noro 447: }
1.24 ! noro 448: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
1.14 noro 449: W = newvect(N+1);
450: W[0] = 1;
1.21 noro 451: R = generic_bfct(B,V1,DV1,W);
1.14 noro 452:
453: return subst(R,s,-s-1);
454: }
455:
1.17 noro 456: /* use an order s.t. [t,x,y,z,...,dt,dx,dy,dz,...,h] */
457:
458: def bfct_via_gbfct_weight(F,V)
459: {
460: N = length(V);
461: D = newvect(N);
462: Wt = getopt(weight);
1.18 noro 463: if ( type(Wt) != 4 ) {
464: for ( I = 0, Wt = []; I < N; I++ )
465: Wt = cons(1,Wt);
466: }
467: Tdeg = w_tdeg(F,V,Wt);
468: WtV = newvect(2*(N+1)+1);
469: WtV[0] = Tdeg;
470: WtV[N+1] = 1;
471: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
472: for ( I = 1; I <= N; I++ ) {
473: WtV[I] = Wt[I-1];
474: WtV[N+1+I] = Tdeg-Wt[I-1]+1;
1.17 noro 475: }
1.18 noro 476: WtV[2*(N+1)] = 1;
477: dp_set_weight(WtV);
1.17 noro 478: for ( I = N-1, DV = []; I >= 0; I-- )
479: DV = cons(strtov("d"+rtostr(V[I])),DV);
480:
1.24 ! noro 481: B = [TMP_T-F];
1.17 noro 482: for ( I = 0; I < N; I++ ) {
1.24 ! noro 483: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
1.17 noro 484: }
1.24 ! noro 485: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
1.17 noro 486: W = newvect(N+1);
487: W[0] = 1;
1.18 noro 488: R = generic_bfct_1(B,V1,DV1,W);
489: dp_set_weight(0);
1.17 noro 490: return subst(R,s,-s-1);
491: }
492:
493: /* use an order s.t. [x,y,z,...,t,dx,dy,dz,...,dt,h] */
494:
495: def bfct_via_gbfct_weight_1(F,V)
496: {
497: N = length(V);
498: D = newvect(N);
499: Wt = getopt(weight);
1.18 noro 500: if ( type(Wt) != 4 ) {
501: for ( I = 0, Wt = []; I < N; I++ )
502: Wt = cons(1,Wt);
503: }
504: Tdeg = w_tdeg(F,V,Wt);
505: WtV = newvect(2*(N+1)+1);
506: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
507: for ( I = 0; I < N; I++ ) {
508: WtV[I] = Wt[I];
509: WtV[N+1+I] = Tdeg-Wt[I]+1;
1.17 noro 510: }
1.18 noro 511: WtV[N] = Tdeg;
512: WtV[2*N+1] = 1;
513: WtV[2*(N+1)] = 1;
514: dp_set_weight(WtV);
1.17 noro 515: for ( I = N-1, DV = []; I >= 0; I-- )
516: DV = cons(strtov("d"+rtostr(V[I])),DV);
517:
1.24 ! noro 518: B = [TMP_T-F];
1.17 noro 519: for ( I = 0; I < N; I++ ) {
1.24 ! noro 520: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
1.17 noro 521: }
1.24 ! noro 522: V1 = append(V,[TMP_T]); DV1 = append(DV,[TMP_DT]);
1.17 noro 523: W = newvect(N+1);
524: W[N] = 1;
1.21 noro 525: R = generic_bfct_1(B,V1,DV1,W);
1.19 noro 526: dp_set_weight(0);
527: return subst(R,s,-s-1);
528: }
529:
530: def bfct_via_gbfct_weight_2(F,V)
531: {
532: N = length(V);
533: D = newvect(N);
534: Wt = getopt(weight);
535: if ( type(Wt) != 4 ) {
536: for ( I = 0, Wt = []; I < N; I++ )
537: Wt = cons(1,Wt);
538: }
539: Tdeg = w_tdeg(F,V,Wt);
540:
541: /* a weight for the first GB computation */
542: /* [t,x1,...,xn,dt,dx1,...,dxn,h] */
543: WtV = newvect(2*(N+1)+1);
544: WtV[0] = Tdeg;
545: WtV[N+1] = 1;
546: WtV[2*(N+1)] = 1;
547: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
548: for ( I = 1; I <= N; I++ ) {
549: WtV[I] = Wt[I-1];
550: WtV[N+1+I] = Tdeg-Wt[I-1]+1;
551: }
552: dp_set_weight(WtV);
553:
554: /* a weight for the second GB computation */
555: /* [x1,...,xn,t,dx1,...,dxn,dt,h] */
556: WtV2 = newvect(2*(N+1)+1);
557: WtV2[N] = Tdeg;
558: WtV2[2*N+1] = 1;
559: WtV2[2*(N+1)] = 1;
560: for ( I = 0; I < N; I++ ) {
561: WtV2[I] = Wt[I];
562: WtV2[N+1+I] = Tdeg-Wt[I]+1;
563: }
564:
565: for ( I = N-1, DV = []; I >= 0; I-- )
566: DV = cons(strtov("d"+rtostr(V[I])),DV);
567:
1.24 ! noro 568: B = [TMP_T-F];
1.19 noro 569: for ( I = 0; I < N; I++ ) {
1.24 ! noro 570: B = cons(DV[I]+diff(F,V[I])*TMP_DT,B);
1.19 noro 571: }
1.24 ! noro 572: V1 = cons(TMP_T,V); DV1 = cons(TMP_DT,DV);
! 573: V2 = append(V,[TMP_T]); DV2 = append(DV,[TMP_DT]);
1.19 noro 574: W = newvect(N+1,[1]);
575: dp_weyl_set_weight(W);
576:
577: VDV = append(V1,DV1);
578: N1 = length(V1);
579: N2 = N1*2;
580:
581: /* create a non-term order MW in D<x,d> */
582: MW = newmat(N2+1,N2);
583: for ( J = 0; J < N1; J++ ) {
584: MW[0][J] = -W[J]; MW[0][N1+J] = W[J];
585: }
586: for ( J = 0; J < N2; J++ ) MW[1][J] = 1;
587: for ( I = 2; I <= N2; I++ ) MW[I][N2-I+1] = -1;
588:
589: /* homogenize F */
590: VDVH = append(VDV,[h]);
591: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,B,VDV)),VDVH);
592:
593: /* compute a groebner basis of FH w.r.t. MWH */
594: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
595:
596: /* dehomigenize GH */
597: G = map(subst,GH,h,1);
598:
599: /* G is a groebner basis w.r.t. a non term order MW */
600: /* take the initial part w.r.t. (-W,W) */
601: GIN = map(initial_part,G,VDV,MW,W);
602:
603: /* GIN is a groebner basis w.r.t. a term order M */
604: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
605:
606: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
607: for ( I = 0, T = 0; I < N1; I++ )
608: T += W[I]*V1[I]*DV1[I];
609:
610: /* change of ordering from VDV to VDV2 */
611: VDV2 = append(V2,DV2);
612: dp_set_weight(WtV2);
1.20 noro 613: for ( Pind = 0; ; Pind++ ) {
614: Prime = lprime(Pind);
615: GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-Prime,0);
616: if ( GIN2 ) break;
617: }
1.19 noro 618:
619: R = weyl_minipoly(GIN2,VDV2,0,T); /* M represents DRL order */
1.18 noro 620: dp_set_weight(0);
1.17 noro 621: return subst(R,s,-s-1);
622: }
623:
1.24 ! noro 624: /* minimal polynomial of s; modular computation */
! 625:
1.6 noro 626: def weyl_minipolym(G,V,O,M,V0)
627: {
628: N = length(V);
629: Len = length(G);
630: dp_ord(O);
631: setmod(M);
632: PS = newvect(Len);
633: PS0 = newvect(Len);
634:
635: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
636: PS0[I] = dp_ptod(car(T),V);
637: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
638: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
639:
640: for ( I = Len - 1, GI = []; I >= 0; I-- )
641: GI = cons(I,GI);
642:
643: U = dp_mod(dp_ptod(V0,V),M,[]);
1.17 noro 644: U = dp_weyl_nf_mod(GI,U,PS,1,M);
1.6 noro 645:
646: T = dp_mod(<<0>>,M,[]);
647: TT = dp_mod(dp_ptod(1,V),M,[]);
648: G = H = [[TT,T]];
649:
650: for ( I = 1; ; I++ ) {
1.14 noro 651: if ( dp_gr_print() )
652: print(".",2);
1.6 noro 653: T = dp_mod(<<I>>,M,[]);
654:
655: TT = dp_weyl_nf_mod(GI,dp_weyl_mul_mod(TT,U,M),PS,1,M);
656: H = cons([TT,T],H);
657: L = dp_lnf_mod([TT,T],G,M);
1.14 noro 658: if ( !L[0] ) {
659: if ( dp_gr_print() )
660: print("");
1.13 noro 661: return dp_dtop(L[1],[t]); /* XXX */
1.14 noro 662: } else
1.6 noro 663: G = insert(G,L);
664: }
665: }
666:
1.24 ! noro 667: /* minimal polynomial of s over Q */
! 668:
1.13 noro 669: def weyl_minipoly(G0,V0,O0,P)
1.6 noro 670: {
1.11 noro 671: HM = hmlist(G0,V0,O0);
1.13 noro 672:
673: N = length(V0);
674: Len = length(G0);
675: dp_ord(O0);
676: PS = newvect(Len);
677: for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
678: PS[I] = dp_ptod(car(T),V0);
679: for ( I = Len - 1, GI = []; I >= 0; I-- )
680: GI = cons(I,GI);
1.20 noro 681: PSM = newvect(Len);
1.13 noro 682: DP = dp_ptod(P,V0);
683:
1.20 noro 684: for ( Pind = 0; ; Pind++ ) {
685: Prime = lprime(Pind);
1.11 noro 686: if ( !valid_modulus(HM,Prime) )
687: continue;
1.20 noro 688: setmod(Prime);
689: for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
690: PSM[I] = dp_mod(dp_ptod(car(T),V0),Prime,[]);
1.13 noro 691:
692: NFP = weyl_nf(GI,DP,1,PS);
1.20 noro 693: NFPM = dp_mod(NFP[0],Prime,[])/ptomp(NFP[1],Prime);
694:
1.13 noro 695: NF = [[dp_ptod(1,V0),1]];
696: LCM = 1;
697:
1.20 noro 698: TM = dp_mod(<<0>>,Prime,[]);
699: TTM = dp_mod(dp_ptod(1,V0),Prime,[]);
700: GM = NFM = [[TTM,TM]];
701:
702: for ( D = 1; ; D++ ) {
1.14 noro 703: if ( dp_gr_print() )
704: print(".",2);
1.13 noro 705: NFPrev = car(NF);
706: NFJ = weyl_nf(GI,
707: dp_weyl_mul(NFP[0],NFPrev[0]),NFP[1]*NFPrev[1],PS);
708: NFJ = remove_cont(NFJ);
709: NF = cons(NFJ,NF);
710: LCM = ilcm(LCM,NFJ[1]);
1.20 noro 711:
712: /* modular computation */
713: TM = dp_mod(<<D>>,Prime,[]);
714: TTM = dp_mod(NFJ[0],Prime,[])/ptomp(NFJ[1],Prime);
715: NFM = cons([TTM,TM],NFM);
716: LM = dp_lnf_mod([TTM,TM],GM,Prime);
717: if ( !LM[0] )
718: break;
719: else
720: GM = insert(GM,LM);
1.13 noro 721: }
1.20 noro 722:
1.14 noro 723: if ( dp_gr_print() )
724: print("");
1.13 noro 725: U = NF[0][0]*idiv(LCM,NF[0][1]);
726: Coef = [];
727: for ( J = D-1; J >= 0; J-- ) {
728: Coef = cons(strtov("u"+rtostr(J)),Coef);
729: U += car(Coef)*NF[D-J][0]*idiv(LCM,NF[D-J][1]);
730: }
1.6 noro 731:
1.13 noro 732: for ( UU = U, Eq = []; UU; UU = dp_rest(UU) )
733: Eq = cons(dp_hc(UU),Eq);
734: M = etom([Eq,Coef]);
735: B = henleq(M,Prime);
736: if ( dp_gr_print() )
737: print("");
1.6 noro 738: if ( B ) {
1.13 noro 739: R = 0;
740: for ( I = 0; I < D; I++ )
741: R += B[0][I]*s^I;
742: R += B[1]*s^D;
1.6 noro 743: return R;
744: }
745: }
746: }
747:
748: def weyl_nf(B,G,M,PS)
749: {
750: for ( D = 0; G; ) {
751: for ( U = 0, L = B; L != []; L = cdr(L) ) {
752: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
753: GCD = igcd(dp_hc(G),dp_hc(R));
754: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
755: U = CG*G-dp_weyl_mul(CR*dp_subd(G,R),R);
756: if ( !U )
757: return [D,M];
758: D *= CG; M *= CG;
759: break;
760: }
761: }
762: if ( U )
763: G = U;
764: else {
765: D += dp_hm(G); G = dp_rest(G);
766: }
767: }
768: return [D,M];
769: }
770:
771: def weyl_nf_mod(B,G,PS,Mod)
772: {
773: for ( D = 0; G; ) {
774: for ( U = 0, L = B; L != []; L = cdr(L) ) {
775: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
776: CR = dp_hc(G)/dp_hc(R);
777: U = G-dp_weyl_mul_mod(CR*dp_mod(dp_subd(G,R),Mod,[]),R,Mod);
778: if ( !U )
779: return D;
1.1 noro 780: break;
1.6 noro 781: }
782: }
783: if ( U )
784: G = U;
785: else {
786: D += dp_hm(G); G = dp_rest(G);
1.1 noro 787: }
788: }
1.6 noro 789: return D;
1.1 noro 790: }
791:
792: def remove_zero(L)
793: {
794: for ( R = []; L != []; L = cdr(L) )
795: if ( car(L) )
796: R = cons(car(L),R);
797: return R;
798: }
799:
800: def z_subst(F,V)
801: {
802: for ( ; V != []; V = cdr(V) )
803: F = subst(F,car(V),0);
804: return F;
805: }
806:
807: def flatmf(L) {
808: for ( S = []; L != []; L = cdr(L) )
809: if ( type(F=car(car(L))) != NUM )
810: S = append(S,[F]);
811: return S;
812: }
813:
814: def intersection(A,B)
815: {
816: for ( L = []; A != []; A = cdr(A) )
817: if ( member(car(A),B) )
818: L = cons(car(A),L);
819: return L;
820: }
821:
822: def b_subst(F,V)
823: {
824: D = deg(F,V);
825: C = newvect(D+1);
826: for ( I = D; I >= 0; I-- )
827: C[I] = coef(F,I,V);
828: for ( I = 0, R = 0; I <= D; I++ )
829: if ( C[I] )
830: R += C[I]*v_factorial(V,I);
831: return R;
832: }
833:
834: def v_factorial(V,N)
835: {
836: for ( J = N-1, R = 1; J >= 0; J-- )
837: R *= V-J;
1.17 noro 838: return R;
839: }
840:
841: def w_tdeg(F,V,W)
842: {
843: dp_set_weight(newvect(length(W),W));
844: T = dp_ptod(F,V);
845: for ( R = 0; T; T = cdr(T) ) {
846: D = dp_td(T);
847: if ( D > R ) R = D;
1.23 noro 848: }
849: return R;
850: }
851:
852: def replace_vars_f(F)
853: {
854: return subst(F,s,TMP_S,t,TMP_T,y1,TMP_Y1,y2,TMP_Y2);
855: }
856:
857: def replace_vars_v(V)
858: {
859: V = replace_var(V,s,TMP_S);
860: V = replace_var(V,t,TMP_T);
861: V = replace_var(V,y1,TMP_Y1);
862: V = replace_var(V,y2,TMP_Y2);
863: return V;
864: }
865:
866: def replace_var(V,X,Y)
867: {
868: V = reverse(V);
869: for ( R = []; V != []; V = cdr(V) ) {
870: Z = car(V);
871: if ( Z == X )
872: R = cons(Y,R);
873: else
874: R = cons(Z,R);
1.17 noro 875: }
1.1 noro 876: return R;
877: }
878: end$
879:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to bfct CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / lib