Annotation of OpenXM_contrib2/asir2000/lib/bfct, Revision 1.22
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.22 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.21 2002/01/30 02:12:58 noro Exp $
1.10 noro 49: */
1.1 noro 50: /* requires 'primdec' */
1.22 ! noro 51:
! 52: extern LIBRARY_GR_LOADED$
! 53: extern LIBRARY_PRIMDEC_LOADED$
! 54:
! 55: if(!LIBRARY_GR_LOADED) load("gr"); else ; LIBRARY_GR_LOADED = 1$
! 56: if(!LIBRARY_PRIMDEC_LOADED) load("primdec"); else ; LIBRARY_PRIMDEC_LOADED = 1$
! 57:
! 58: /* toplevel */
! 59:
! 60: def bfunction(F)
! 61: {
! 62: V = vars(F);
! 63: N = length(V);
! 64: D = newvect(N);
! 65:
! 66: for ( I = 0; I < N; I++ )
! 67: D[I] = [deg(F,V[I]),V[I]];
! 68: qsort(D,compare_first);
! 69: for ( V = [], I = 0; I < N; I++ )
! 70: V = cons(D[I][1],V);
! 71: return bfct_via_gbfct_weight(F,V);
! 72: }
1.1 noro 73:
1.6 noro 74: /* annihilating ideal of F^s */
1.1 noro 75:
76: def ann(F)
77: {
78: V = vars(F);
79: N = length(V);
1.8 noro 80: D = newvect(N);
81:
82: for ( I = 0; I < N; I++ )
83: D[I] = [deg(F,V[I]),V[I]];
84: qsort(D,compare_first);
85: for ( V = [], I = N-1; I >= 0; I-- )
86: V = cons(D[I][1],V);
87:
1.1 noro 88: for ( I = N-1, DV = []; I >= 0; I-- )
89: DV = cons(strtov("d"+rtostr(V[I])),DV);
1.8 noro 90:
91: W = append([y1,y2,t],V);
1.1 noro 92: DW = append([dy1,dy2,dt],DV);
1.8 noro 93:
94: B = [1-y1*y2,t-y1*F];
1.1 noro 95: for ( I = 0; I < N; I++ ) {
96: B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
97: }
1.10 noro 98:
99: /* homogenized (heuristics) */
1.1 noro 100: dp_nelim(2);
1.10 noro 101: G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
1.1 noro 102: G1 = [];
103: for ( T = G0; T != []; T = cdr(T) ) {
104: E = car(T); VL = vars(E);
105: if ( !member(y1,VL) && !member(y2,VL) )
106: G1 = cons(E,G1);
107: }
1.12 noro 108: G2 = map(psi,G1,t,dt);
109: G3 = map(subst,G2,t,-1-s);
110: return G3;
1.1 noro 111: }
112:
1.10 noro 113: /*
114: * compute J_f|s=r, where r = the minimal integral root of global b_f(s)
115: * ann0(F) returns [MinRoot,Ideal]
116: */
117:
118: def ann0(F)
119: {
120: V = vars(F);
121: N = length(V);
122: D = newvect(N);
123:
124: for ( I = 0; I < N; I++ )
125: D[I] = [deg(F,V[I]),V[I]];
126: qsort(D,compare_first);
127: for ( V = [], I = 0; I < N; I++ )
128: V = cons(D[I][1],V);
129:
130: for ( I = N-1, DV = []; I >= 0; I-- )
131: DV = cons(strtov("d"+rtostr(V[I])),DV);
132:
133: /* XXX : heuristics */
134: W = append([y1,y2,t],reverse(V));
135: DW = append([dy1,dy2,dt],reverse(DV));
136: WDW = append(W,DW);
137:
138: B = [1-y1*y2,t-y1*F];
139: for ( I = 0; I < N; I++ ) {
140: B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
141: }
142:
143: /* homogenized (heuristics) */
144: dp_nelim(2);
145: G0 = dp_weyl_gr_main(B,WDW,1,0,6);
146: G1 = [];
147: for ( T = G0; T != []; T = cdr(T) ) {
148: E = car(T); VL = vars(E);
149: if ( !member(y1,VL) && !member(y2,VL) )
150: G1 = cons(E,G1);
151: }
1.12 noro 152: G2 = map(psi,G1,t,dt);
153: G3 = map(subst,G2,t,-1-s);
1.10 noro 154:
1.12 noro 155: /* G3 = J_f(s) */
1.10 noro 156:
157: V1 = cons(s,V); DV1 = cons(ds,DV); V1DV1 = append(V1,DV1);
1.12 noro 158: G4 = dp_weyl_gr_main(cons(F,G3),V1DV1,0,1,0);
159: Bf = weyl_minipoly(G4,V1DV1,0,s);
1.10 noro 160:
161: FList = cdr(fctr(Bf));
162: for ( T = FList, Min = 0; T != []; T = cdr(T) ) {
163: LF = car(car(T));
164: Root = -coef(LF,0)/coef(LF,1);
165: if ( dn(Root) == 1 && Root < Min )
166: Min = Root;
167: }
1.12 noro 168: return [Min,map(subst,G3,s,Min)];
1.10 noro 169: }
170:
1.7 noro 171: def indicial1(F,V)
1.6 noro 172: {
173: W = append([y1,t],V);
174: N = length(V);
175: B = [t-y1*F];
176: for ( I = N-1, DV = []; I >= 0; I-- )
177: DV = cons(strtov("d"+rtostr(V[I])),DV);
178: DW = append([dy1,dt],DV);
179: for ( I = 0; I < N; I++ ) {
180: B = cons(DV[I]+y1*diff(F,V[I])*dt,B);
181: }
182: dp_nelim(1);
1.10 noro 183:
184: /* homogenized (heuristics) */
1.7 noro 185: G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6);
1.6 noro 186: G1 = map(subst,G0,y1,1);
187: G2 = map(psi,G1,t,dt);
188: G3 = map(subst,G2,t,-s-1);
189: return G3;
190: }
191:
192: def psi(F,T,DT)
193: {
194: D = dp_ptod(F,[T,DT]);
195: Wmax = weight(D);
196: D1 = dp_rest(D);
197: for ( ; D1; D1 = dp_rest(D1) )
198: if ( weight(D1) > Wmax )
199: Wmax = weight(D1);
200: for ( D1 = D, Dmax = 0; D1; D1 = dp_rest(D1) )
201: if ( weight(D1) == Wmax )
202: Dmax += dp_hm(D1);
203: if ( Wmax >= 0 )
204: Dmax = dp_weyl_mul(<<Wmax,0>>,Dmax);
205: else
206: Dmax = dp_weyl_mul(<<0,-Wmax>>,Dmax);
207: Rmax = dp_dtop(Dmax,[T,DT]);
208: R = b_subst(subst(Rmax,DT,1),T);
209: return R;
210: }
211:
212: def weight(D)
213: {
214: V = dp_etov(D);
215: return V[1]-V[0];
216: }
217:
218: def compare_first(A,B)
219: {
220: A0 = car(A);
221: B0 = car(B);
222: if ( A0 > B0 )
223: return 1;
224: else if ( A0 < B0 )
225: return -1;
226: else
227: return 0;
228: }
229:
1.13 noro 230: /* generic b-function w.r.t. weight vector W */
231:
232: def generic_bfct(F,V,DV,W)
233: {
234: N = length(V);
235: N2 = N*2;
236:
1.16 noro 237: /* If W is a list, convert it to a vector */
238: if ( type(W) == 4 )
239: W = newvect(length(W),W);
1.15 noro 240: dp_weyl_set_weight(W);
241:
1.14 noro 242: /* create a term order M in D<x,d> (DRL) */
1.13 noro 243: M = newmat(N2,N2);
244: for ( J = 0; J < N2; J++ )
245: M[0][J] = 1;
246: for ( I = 1; I < N2; I++ )
247: M[I][N2-I] = -1;
248:
249: VDV = append(V,DV);
250:
251: /* create a non-term order MW in D<x,d> */
252: MW = newmat(N2+1,N2);
253: for ( J = 0; J < N; J++ )
254: MW[0][J] = -W[J];
255: for ( ; J < N2; J++ )
256: MW[0][J] = W[J-N];
257: for ( I = 1; I <= N2; I++ )
258: for ( J = 0; J < N2; J++ )
259: MW[I][J] = M[I-1][J];
260:
261: /* create a homogenized term order MWH in D<x,d,h> */
262: MWH = newmat(N2+2,N2+1);
263: for ( J = 0; J <= N2; J++ )
264: MWH[0][J] = 1;
265: for ( I = 1; I <= N2+1; I++ )
266: for ( J = 0; J < N2; J++ )
267: MWH[I][J] = MW[I-1][J];
268:
269: /* homogenize F */
270: VDVH = append(VDV,[h]);
271: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
272:
273: /* compute a groebner basis of FH w.r.t. MWH */
1.21 noro 274: dp_gr_flags(["Top",1,"NoRA",1]);
1.15 noro 275: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
1.21 noro 276: dp_gr_flags(["Top",0,"NoRA",0]);
1.13 noro 277:
278: /* dehomigenize GH */
279: G = map(subst,GH,h,1);
280:
281: /* G is a groebner basis w.r.t. a non term order MW */
282: /* take the initial part w.r.t. (-W,W) */
283: GIN = map(initial_part,G,VDV,MW,W);
284:
285: /* GIN is a groebner basis w.r.t. a term order M */
286: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
287:
288: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
289: for ( I = 0, T = 0; I < N; I++ )
290: T += W[I]*V[I]*DV[I];
1.14 noro 291: B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
1.13 noro 292: return B;
293: }
294:
1.18 noro 295: /* all term reduction + interreduce */
296: def generic_bfct_1(F,V,DV,W)
297: {
298: N = length(V);
299: N2 = N*2;
300:
301: /* If W is a list, convert it to a vector */
302: if ( type(W) == 4 )
303: W = newvect(length(W),W);
304: dp_weyl_set_weight(W);
305:
306: /* create a term order M in D<x,d> (DRL) */
307: M = newmat(N2,N2);
308: for ( J = 0; J < N2; J++ )
309: M[0][J] = 1;
310: for ( I = 1; I < N2; I++ )
311: M[I][N2-I] = -1;
312:
313: VDV = append(V,DV);
314:
315: /* create a non-term order MW in D<x,d> */
316: MW = newmat(N2+1,N2);
317: for ( J = 0; J < N; J++ )
318: MW[0][J] = -W[J];
319: for ( ; J < N2; J++ )
320: MW[0][J] = W[J-N];
321: for ( I = 1; I <= N2; I++ )
322: for ( J = 0; J < N2; J++ )
323: MW[I][J] = M[I-1][J];
324:
325: /* create a homogenized term order MWH in D<x,d,h> */
326: MWH = newmat(N2+2,N2+1);
327: for ( J = 0; J <= N2; J++ )
328: MWH[0][J] = 1;
329: for ( I = 1; I <= N2+1; I++ )
330: for ( J = 0; J < N2; J++ )
331: MWH[I][J] = MW[I-1][J];
332:
333: /* homogenize F */
334: VDVH = append(VDV,[h]);
335: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH);
336:
337: /* compute a groebner basis of FH w.r.t. MWH */
338: /* dp_gr_flags(["Top",1,"NoRA",1]); */
339: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
340: /* dp_gr_flags(["Top",0,"NoRA",0]); */
341:
342: /* dehomigenize GH */
343: G = map(subst,GH,h,1);
344:
345: /* G is a groebner basis w.r.t. a non term order MW */
346: /* take the initial part w.r.t. (-W,W) */
347: GIN = map(initial_part,G,VDV,MW,W);
348:
349: /* GIN is a groebner basis w.r.t. a term order M */
350: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
351:
352: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
353: for ( I = 0, T = 0; I < N; I++ )
354: T += W[I]*V[I]*DV[I];
355: B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */
356: return B;
357: }
358:
1.13 noro 359: def initial_part(F,V,MW,W)
360: {
361: N2 = length(V);
362: N = N2/2;
363: dp_ord(MW);
364: DF = dp_ptod(F,V);
365: R = dp_hm(DF);
366: DF = dp_rest(DF);
367:
368: E = dp_etov(R);
369: for ( I = 0, TW = 0; I < N; I++ )
370: TW += W[I]*(-E[I]+E[N+I]);
371: RW = TW;
372:
373: for ( ; DF; DF = dp_rest(DF) ) {
374: E = dp_etov(DF);
375: for ( I = 0, TW = 0; I < N; I++ )
376: TW += W[I]*(-E[I]+E[N+I]);
377: if ( TW == RW )
378: R += dp_hm(DF);
379: else if ( TW < RW )
380: break;
381: else
382: error("initial_part : cannot happen");
383: }
384: return dp_dtop(R,V);
385:
386: }
387:
1.1 noro 388: /* b-function of F ? */
389:
390: def bfct(F)
391: {
392: V = vars(F);
393: N = length(V);
1.6 noro 394: D = newvect(N);
1.7 noro 395:
1.6 noro 396: for ( I = 0; I < N; I++ )
397: D[I] = [deg(F,V[I]),V[I]];
398: qsort(D,compare_first);
399: for ( V = [], I = 0; I < N; I++ )
400: V = cons(D[I][1],V);
1.1 noro 401: for ( I = N-1, DV = []; I >= 0; I-- )
402: DV = cons(strtov("d"+rtostr(V[I])),DV);
1.6 noro 403: V1 = cons(s,V); DV1 = cons(ds,DV);
1.7 noro 404:
405: G0 = indicial1(F,reverse(V));
406: G1 = dp_weyl_gr_main(G0,append(V1,DV1),0,1,0);
407: Minipoly = weyl_minipoly(G1,append(V1,DV1),0,s);
1.6 noro 408: return Minipoly;
409: }
410:
1.14 noro 411: /* b-function computation via generic_bfct() (experimental) */
412:
413: def bfct_via_gbfct(F)
414: {
415: V = vars(F);
416: N = length(V);
417: D = newvect(N);
418:
419: for ( I = 0; I < N; I++ )
420: D[I] = [deg(F,V[I]),V[I]];
421: qsort(D,compare_first);
422: for ( V = [], I = 0; I < N; I++ )
423: V = cons(D[I][1],V);
424: V = reverse(V);
425: for ( I = N-1, DV = []; I >= 0; I-- )
426: DV = cons(strtov("d"+rtostr(V[I])),DV);
427:
428: B = [t-F];
429: for ( I = 0; I < N; I++ ) {
430: B = cons(DV[I]+diff(F,V[I])*dt,B);
431: }
432: V1 = cons(t,V); DV1 = cons(dt,DV);
433: W = newvect(N+1);
434: W[0] = 1;
1.21 noro 435: R = generic_bfct(B,V1,DV1,W);
1.14 noro 436:
437: return subst(R,s,-s-1);
438: }
439:
1.17 noro 440: /* use an order s.t. [t,x,y,z,...,dt,dx,dy,dz,...,h] */
441:
442: def bfct_via_gbfct_weight(F,V)
443: {
444: N = length(V);
445: D = newvect(N);
446: Wt = getopt(weight);
1.18 noro 447: if ( type(Wt) != 4 ) {
448: for ( I = 0, Wt = []; I < N; I++ )
449: Wt = cons(1,Wt);
450: }
451: Tdeg = w_tdeg(F,V,Wt);
452: WtV = newvect(2*(N+1)+1);
453: WtV[0] = Tdeg;
454: WtV[N+1] = 1;
455: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
456: for ( I = 1; I <= N; I++ ) {
457: WtV[I] = Wt[I-1];
458: WtV[N+1+I] = Tdeg-Wt[I-1]+1;
1.17 noro 459: }
1.18 noro 460: WtV[2*(N+1)] = 1;
461: dp_set_weight(WtV);
1.17 noro 462: for ( I = N-1, DV = []; I >= 0; I-- )
463: DV = cons(strtov("d"+rtostr(V[I])),DV);
464:
465: B = [t-F];
466: for ( I = 0; I < N; I++ ) {
467: B = cons(DV[I]+diff(F,V[I])*dt,B);
468: }
469: V1 = cons(t,V); DV1 = cons(dt,DV);
470: W = newvect(N+1);
471: W[0] = 1;
1.18 noro 472: R = generic_bfct_1(B,V1,DV1,W);
473: dp_set_weight(0);
1.17 noro 474: return subst(R,s,-s-1);
475: }
476:
477: /* use an order s.t. [x,y,z,...,t,dx,dy,dz,...,dt,h] */
478:
479: def bfct_via_gbfct_weight_1(F,V)
480: {
481: N = length(V);
482: D = newvect(N);
483: Wt = getopt(weight);
1.18 noro 484: if ( type(Wt) != 4 ) {
485: for ( I = 0, Wt = []; I < N; I++ )
486: Wt = cons(1,Wt);
487: }
488: Tdeg = w_tdeg(F,V,Wt);
489: WtV = newvect(2*(N+1)+1);
490: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
491: for ( I = 0; I < N; I++ ) {
492: WtV[I] = Wt[I];
493: WtV[N+1+I] = Tdeg-Wt[I]+1;
1.17 noro 494: }
1.18 noro 495: WtV[N] = Tdeg;
496: WtV[2*N+1] = 1;
497: WtV[2*(N+1)] = 1;
498: dp_set_weight(WtV);
1.17 noro 499: for ( I = N-1, DV = []; I >= 0; I-- )
500: DV = cons(strtov("d"+rtostr(V[I])),DV);
501:
502: B = [t-F];
503: for ( I = 0; I < N; I++ ) {
504: B = cons(DV[I]+diff(F,V[I])*dt,B);
505: }
506: V1 = append(V,[t]); DV1 = append(DV,[dt]);
507: W = newvect(N+1);
508: W[N] = 1;
1.21 noro 509: R = generic_bfct_1(B,V1,DV1,W);
1.19 noro 510: dp_set_weight(0);
511: return subst(R,s,-s-1);
512: }
513:
514: def bfct_via_gbfct_weight_2(F,V)
515: {
516: N = length(V);
517: D = newvect(N);
518: Wt = getopt(weight);
519: if ( type(Wt) != 4 ) {
520: for ( I = 0, Wt = []; I < N; I++ )
521: Wt = cons(1,Wt);
522: }
523: Tdeg = w_tdeg(F,V,Wt);
524:
525: /* a weight for the first GB computation */
526: /* [t,x1,...,xn,dt,dx1,...,dxn,h] */
527: WtV = newvect(2*(N+1)+1);
528: WtV[0] = Tdeg;
529: WtV[N+1] = 1;
530: WtV[2*(N+1)] = 1;
531: /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */
532: for ( I = 1; I <= N; I++ ) {
533: WtV[I] = Wt[I-1];
534: WtV[N+1+I] = Tdeg-Wt[I-1]+1;
535: }
536: dp_set_weight(WtV);
537:
538: /* a weight for the second GB computation */
539: /* [x1,...,xn,t,dx1,...,dxn,dt,h] */
540: WtV2 = newvect(2*(N+1)+1);
541: WtV2[N] = Tdeg;
542: WtV2[2*N+1] = 1;
543: WtV2[2*(N+1)] = 1;
544: for ( I = 0; I < N; I++ ) {
545: WtV2[I] = Wt[I];
546: WtV2[N+1+I] = Tdeg-Wt[I]+1;
547: }
548:
549: for ( I = N-1, DV = []; I >= 0; I-- )
550: DV = cons(strtov("d"+rtostr(V[I])),DV);
551:
552: B = [t-F];
553: for ( I = 0; I < N; I++ ) {
554: B = cons(DV[I]+diff(F,V[I])*dt,B);
555: }
556: V1 = cons(t,V); DV1 = cons(dt,DV);
557: V2 = append(V,[t]); DV2 = append(DV,[dt]);
558: W = newvect(N+1,[1]);
559: dp_weyl_set_weight(W);
560:
561: VDV = append(V1,DV1);
562: N1 = length(V1);
563: N2 = N1*2;
564:
565: /* create a non-term order MW in D<x,d> */
566: MW = newmat(N2+1,N2);
567: for ( J = 0; J < N1; J++ ) {
568: MW[0][J] = -W[J]; MW[0][N1+J] = W[J];
569: }
570: for ( J = 0; J < N2; J++ ) MW[1][J] = 1;
571: for ( I = 2; I <= N2; I++ ) MW[I][N2-I+1] = -1;
572:
573: /* homogenize F */
574: VDVH = append(VDV,[h]);
575: FH = map(dp_dtop,map(dp_homo,map(dp_ptod,B,VDV)),VDVH);
576:
577: /* compute a groebner basis of FH w.r.t. MWH */
578: GH = dp_weyl_gr_main(FH,VDVH,0,1,11);
579:
580: /* dehomigenize GH */
581: G = map(subst,GH,h,1);
582:
583: /* G is a groebner basis w.r.t. a non term order MW */
584: /* take the initial part w.r.t. (-W,W) */
585: GIN = map(initial_part,G,VDV,MW,W);
586:
587: /* GIN is a groebner basis w.r.t. a term order M */
588: /* As -W+W=0, gr_(-W,W)(D<x,d>) = D<x,d> */
589:
590: /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */
591: for ( I = 0, T = 0; I < N1; I++ )
592: T += W[I]*V1[I]*DV1[I];
593:
594: /* change of ordering from VDV to VDV2 */
595: VDV2 = append(V2,DV2);
596: dp_set_weight(WtV2);
1.20 noro 597: for ( Pind = 0; ; Pind++ ) {
598: Prime = lprime(Pind);
599: GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-Prime,0);
600: if ( GIN2 ) break;
601: }
1.19 noro 602:
603: R = weyl_minipoly(GIN2,VDV2,0,T); /* M represents DRL order */
1.18 noro 604: dp_set_weight(0);
1.17 noro 605: return subst(R,s,-s-1);
606: }
607:
1.6 noro 608: def weyl_minipolym(G,V,O,M,V0)
609: {
610: N = length(V);
611: Len = length(G);
612: dp_ord(O);
613: setmod(M);
614: PS = newvect(Len);
615: PS0 = newvect(Len);
616:
617: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
618: PS0[I] = dp_ptod(car(T),V);
619: for ( I = 0, T = G; T != []; T = cdr(T), I++ )
620: PS[I] = dp_mod(dp_ptod(car(T),V),M,[]);
621:
622: for ( I = Len - 1, GI = []; I >= 0; I-- )
623: GI = cons(I,GI);
624:
625: U = dp_mod(dp_ptod(V0,V),M,[]);
1.17 noro 626: U = dp_weyl_nf_mod(GI,U,PS,1,M);
1.6 noro 627:
628: T = dp_mod(<<0>>,M,[]);
629: TT = dp_mod(dp_ptod(1,V),M,[]);
630: G = H = [[TT,T]];
631:
632: for ( I = 1; ; I++ ) {
1.14 noro 633: if ( dp_gr_print() )
634: print(".",2);
1.6 noro 635: T = dp_mod(<<I>>,M,[]);
636:
637: TT = dp_weyl_nf_mod(GI,dp_weyl_mul_mod(TT,U,M),PS,1,M);
638: H = cons([TT,T],H);
639: L = dp_lnf_mod([TT,T],G,M);
1.14 noro 640: if ( !L[0] ) {
641: if ( dp_gr_print() )
642: print("");
1.13 noro 643: return dp_dtop(L[1],[t]); /* XXX */
1.14 noro 644: } else
1.6 noro 645: G = insert(G,L);
646: }
647: }
648:
1.13 noro 649: def weyl_minipoly(G0,V0,O0,P)
1.6 noro 650: {
1.11 noro 651: HM = hmlist(G0,V0,O0);
1.13 noro 652:
653: N = length(V0);
654: Len = length(G0);
655: dp_ord(O0);
656: PS = newvect(Len);
657: for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
658: PS[I] = dp_ptod(car(T),V0);
659: for ( I = Len - 1, GI = []; I >= 0; I-- )
660: GI = cons(I,GI);
1.20 noro 661: PSM = newvect(Len);
1.13 noro 662: DP = dp_ptod(P,V0);
663:
1.20 noro 664: for ( Pind = 0; ; Pind++ ) {
665: Prime = lprime(Pind);
1.11 noro 666: if ( !valid_modulus(HM,Prime) )
667: continue;
1.20 noro 668: setmod(Prime);
669: for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ )
670: PSM[I] = dp_mod(dp_ptod(car(T),V0),Prime,[]);
1.13 noro 671:
672: NFP = weyl_nf(GI,DP,1,PS);
1.20 noro 673: NFPM = dp_mod(NFP[0],Prime,[])/ptomp(NFP[1],Prime);
674:
1.13 noro 675: NF = [[dp_ptod(1,V0),1]];
676: LCM = 1;
677:
1.20 noro 678: TM = dp_mod(<<0>>,Prime,[]);
679: TTM = dp_mod(dp_ptod(1,V0),Prime,[]);
680: GM = NFM = [[TTM,TM]];
681:
682: for ( D = 1; ; D++ ) {
1.14 noro 683: if ( dp_gr_print() )
684: print(".",2);
1.13 noro 685: NFPrev = car(NF);
686: NFJ = weyl_nf(GI,
687: dp_weyl_mul(NFP[0],NFPrev[0]),NFP[1]*NFPrev[1],PS);
688: NFJ = remove_cont(NFJ);
689: NF = cons(NFJ,NF);
690: LCM = ilcm(LCM,NFJ[1]);
1.20 noro 691:
692: /* modular computation */
693: TM = dp_mod(<<D>>,Prime,[]);
694: TTM = dp_mod(NFJ[0],Prime,[])/ptomp(NFJ[1],Prime);
695: NFM = cons([TTM,TM],NFM);
696: LM = dp_lnf_mod([TTM,TM],GM,Prime);
697: if ( !LM[0] )
698: break;
699: else
700: GM = insert(GM,LM);
1.13 noro 701: }
1.20 noro 702:
1.14 noro 703: if ( dp_gr_print() )
704: print("");
1.13 noro 705: U = NF[0][0]*idiv(LCM,NF[0][1]);
706: Coef = [];
707: for ( J = D-1; J >= 0; J-- ) {
708: Coef = cons(strtov("u"+rtostr(J)),Coef);
709: U += car(Coef)*NF[D-J][0]*idiv(LCM,NF[D-J][1]);
710: }
1.6 noro 711:
1.13 noro 712: for ( UU = U, Eq = []; UU; UU = dp_rest(UU) )
713: Eq = cons(dp_hc(UU),Eq);
714: M = etom([Eq,Coef]);
715: B = henleq(M,Prime);
716: if ( dp_gr_print() )
717: print("");
1.6 noro 718: if ( B ) {
1.13 noro 719: R = 0;
720: for ( I = 0; I < D; I++ )
721: R += B[0][I]*s^I;
722: R += B[1]*s^D;
1.6 noro 723: return R;
724: }
725: }
726: }
727:
728: def weyl_nf(B,G,M,PS)
729: {
730: for ( D = 0; G; ) {
731: for ( U = 0, L = B; L != []; L = cdr(L) ) {
732: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
733: GCD = igcd(dp_hc(G),dp_hc(R));
734: CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD);
735: U = CG*G-dp_weyl_mul(CR*dp_subd(G,R),R);
736: if ( !U )
737: return [D,M];
738: D *= CG; M *= CG;
739: break;
740: }
741: }
742: if ( U )
743: G = U;
744: else {
745: D += dp_hm(G); G = dp_rest(G);
746: }
747: }
748: return [D,M];
749: }
750:
751: def weyl_nf_mod(B,G,PS,Mod)
752: {
753: for ( D = 0; G; ) {
754: for ( U = 0, L = B; L != []; L = cdr(L) ) {
755: if ( dp_redble(G,R=PS[car(L)]) > 0 ) {
756: CR = dp_hc(G)/dp_hc(R);
757: U = G-dp_weyl_mul_mod(CR*dp_mod(dp_subd(G,R),Mod,[]),R,Mod);
758: if ( !U )
759: return D;
1.1 noro 760: break;
1.6 noro 761: }
762: }
763: if ( U )
764: G = U;
765: else {
766: D += dp_hm(G); G = dp_rest(G);
1.1 noro 767: }
768: }
1.6 noro 769: return D;
1.1 noro 770: }
771:
772: def remove_zero(L)
773: {
774: for ( R = []; L != []; L = cdr(L) )
775: if ( car(L) )
776: R = cons(car(L),R);
777: return R;
778: }
779:
780: def z_subst(F,V)
781: {
782: for ( ; V != []; V = cdr(V) )
783: F = subst(F,car(V),0);
784: return F;
785: }
786:
787: def flatmf(L) {
788: for ( S = []; L != []; L = cdr(L) )
789: if ( type(F=car(car(L))) != NUM )
790: S = append(S,[F]);
791: return S;
792: }
793:
794: def member(A,L) {
795: for ( ; L != []; L = cdr(L) )
796: if ( A == car(L) )
797: return 1;
798: return 0;
799: }
800:
801: def intersection(A,B)
802: {
803: for ( L = []; A != []; A = cdr(A) )
804: if ( member(car(A),B) )
805: L = cons(car(A),L);
806: return L;
807: }
808:
809: def b_subst(F,V)
810: {
811: D = deg(F,V);
812: C = newvect(D+1);
813: for ( I = D; I >= 0; I-- )
814: C[I] = coef(F,I,V);
815: for ( I = 0, R = 0; I <= D; I++ )
816: if ( C[I] )
817: R += C[I]*v_factorial(V,I);
818: return R;
819: }
820:
821: def v_factorial(V,N)
822: {
823: for ( J = N-1, R = 1; J >= 0; J-- )
824: R *= V-J;
1.17 noro 825: return R;
826: }
827:
828: def w_tdeg(F,V,W)
829: {
830: dp_set_weight(newvect(length(W),W));
831: T = dp_ptod(F,V);
832: for ( R = 0; T; T = cdr(T) ) {
833: D = dp_td(T);
834: if ( D > R ) R = D;
835: }
1.1 noro 836: return R;
837: }
838: end$
839:
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