=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v retrieving revision 1.1 retrieving revision 1.19 diff -u -p -r1.1 -r1.19 --- OpenXM_contrib2/asir2000/lib/bfct 2000/06/05 04:59:34 1.1 +++ OpenXM_contrib2/asir2000/lib/bfct 2002/01/29 02:03:41 1.19 @@ -1,90 +1,729 @@ -/* $OpenXM$ */ +/* + * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED + * All rights reserved. + * + * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited, + * non-exclusive and royalty-free license to use, copy, modify and + * redistribute, solely for non-commercial and non-profit purposes, the + * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and + * conditions of this Agreement. For the avoidance of doubt, you acquire + * only a limited right to use the SOFTWARE hereunder, and FLL or any + * third party developer retains all rights, including but not limited to + * copyrights, in and to the SOFTWARE. + * + * (1) FLL does not grant you a license in any way for commercial + * purposes. You may use the SOFTWARE only for non-commercial and + * non-profit purposes only, such as academic, research and internal + * business use. + * (2) The SOFTWARE is protected by the Copyright Law of Japan and + * international copyright treaties. If you make copies of the SOFTWARE, + * with or without modification, as permitted hereunder, you shall affix + * to all such copies of the SOFTWARE the above copyright notice. + * (3) An explicit reference to this SOFTWARE and its copyright owner + * shall be made on your publication or presentation in any form of the + * results obtained by use of the SOFTWARE. + * (4) In the event that you modify the SOFTWARE, you shall notify FLL by + * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification + * for such modification or the source code of the modified part of the + * SOFTWARE. + * + * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL + * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND + * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS + * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES' + * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY + * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY. + * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT, + * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY + * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL + * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES + * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES + * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY + * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF + * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART + * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY + * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE, + * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE. + * + * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.18 2002/01/28 02:42:27 noro Exp $ + */ /* requires 'primdec' */ -/* annihilating ideal of F^s ? */ +/* annihilating ideal of F^s */ def ann(F) { V = vars(F); - W = append([y1,y2,t],V); N = length(V); - B = [1-y1*y2,t-y1*F]; + D = newvect(N); + + for ( I = 0; I < N; I++ ) + D[I] = [deg(F,V[I]),V[I]]; + qsort(D,compare_first); + for ( V = [], I = N-1; I >= 0; I-- ) + V = cons(D[I][1],V); + for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); + + W = append([y1,y2,t],V); DW = append([dy1,dy2,dt],DV); + + B = [1-y1*y2,t-y1*F]; for ( I = 0; I < N; I++ ) { B = cons(DV[I]+y1*diff(F,V[I])*dt,B); } - ctrl("do_weyl",1); + + /* homogenized (heuristics) */ dp_nelim(2); - G0 = dp_gr_main(B,append(W,DW),0,0,6); + G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6); G1 = []; for ( T = G0; T != []; T = cdr(T) ) { E = car(T); VL = vars(E); if ( !member(y1,VL) && !member(y2,VL) ) G1 = cons(E,G1); } - G2 = map(subst,G1,dt,1); - G3 = map(b_subst,G2,t); - G4 = map(subst,G3,t,-1-s); - ctrl("do_weyl",0); - return G4; + G2 = map(psi,G1,t,dt); + G3 = map(subst,G2,t,-1-s); + return G3; } -/* b-function of F ? */ +/* + * compute J_f|s=r, where r = the minimal integral root of global b_f(s) + * ann0(F) returns [MinRoot,Ideal] + */ -def bfct(F) +def ann0(F) { - G4 = ann(F); - - ctrl("do_weyl",1); V = vars(F); N = length(V); + D = newvect(N); + + for ( I = 0; I < N; I++ ) + D[I] = [deg(F,V[I]),V[I]]; + qsort(D,compare_first); + for ( V = [], I = 0; I < N; I++ ) + V = cons(D[I][1],V); + for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); - N1 = 2*(N+1); + /* XXX : heuristics */ + W = append([y1,y2,t],reverse(V)); + DW = append([dy1,dy2,dt],reverse(DV)); + WDW = append(W,DW); - M = newmat(N1+1,N1); - for ( J = N+1; J < N1; J++ ) + B = [1-y1*y2,t-y1*F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+y1*diff(F,V[I])*dt,B); + } + + /* homogenized (heuristics) */ + dp_nelim(2); + G0 = dp_weyl_gr_main(B,WDW,1,0,6); + G1 = []; + for ( T = G0; T != []; T = cdr(T) ) { + E = car(T); VL = vars(E); + if ( !member(y1,VL) && !member(y2,VL) ) + G1 = cons(E,G1); + } + G2 = map(psi,G1,t,dt); + G3 = map(subst,G2,t,-1-s); + + /* G3 = J_f(s) */ + + V1 = cons(s,V); DV1 = cons(ds,DV); V1DV1 = append(V1,DV1); + G4 = dp_weyl_gr_main(cons(F,G3),V1DV1,0,1,0); + Bf = weyl_minipoly(G4,V1DV1,0,s); + + FList = cdr(fctr(Bf)); + for ( T = FList, Min = 0; T != []; T = cdr(T) ) { + LF = car(car(T)); + Root = -coef(LF,0)/coef(LF,1); + if ( dn(Root) == 1 && Root < Min ) + Min = Root; + } + return [Min,map(subst,G3,s,Min)]; +} + +def indicial1(F,V) +{ + W = append([y1,t],V); + N = length(V); + B = [t-y1*F]; + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + DW = append([dy1,dt],DV); + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+y1*diff(F,V[I])*dt,B); + } + dp_nelim(1); + + /* homogenized (heuristics) */ + G0 = dp_weyl_gr_main(B,append(W,DW),1,0,6); + G1 = map(subst,G0,y1,1); + G2 = map(psi,G1,t,dt); + G3 = map(subst,G2,t,-s-1); + return G3; +} + +def psi(F,T,DT) +{ + D = dp_ptod(F,[T,DT]); + Wmax = weight(D); + D1 = dp_rest(D); + for ( ; D1; D1 = dp_rest(D1) ) + if ( weight(D1) > Wmax ) + Wmax = weight(D1); + for ( D1 = D, Dmax = 0; D1; D1 = dp_rest(D1) ) + if ( weight(D1) == Wmax ) + Dmax += dp_hm(D1); + if ( Wmax >= 0 ) + Dmax = dp_weyl_mul(<>,Dmax); + else + Dmax = dp_weyl_mul(<<0,-Wmax>>,Dmax); + Rmax = dp_dtop(Dmax,[T,DT]); + R = b_subst(subst(Rmax,DT,1),T); + return R; +} + +def weight(D) +{ + V = dp_etov(D); + return V[1]-V[0]; +} + +def compare_first(A,B) +{ + A0 = car(A); + B0 = car(B); + if ( A0 > B0 ) + return 1; + else if ( A0 < B0 ) + return -1; + else + return 0; +} + +/* generic b-function w.r.t. weight vector W */ + +def generic_bfct(F,V,DV,W) +{ + N = length(V); + N2 = N*2; + + /* If W is a list, convert it to a vector */ + if ( type(W) == 4 ) + W = newvect(length(W),W); + dp_weyl_set_weight(W); + + /* create a term order M in D (DRL) */ + M = newmat(N2,N2); + for ( J = 0; J < N2; J++ ) M[0][J] = 1; - for ( J = 0; J < N+1; J++ ) - M[1][J] = 1; -#if 0 - for ( I = 0; I < N+1; I++ ) - M[I+2][N-I] = -1; + for ( I = 1; I < N2; I++ ) + M[I][N2-I] = -1; + + VDV = append(V,DV); + + /* create a non-term order MW in D */ + MW = newmat(N2+1,N2); + for ( J = 0; J < N; J++ ) + MW[0][J] = -W[J]; + for ( ; J < N2; J++ ) + MW[0][J] = W[J-N]; + for ( I = 1; I <= N2; I++ ) + for ( J = 0; J < N2; J++ ) + MW[I][J] = M[I-1][J]; + + /* create a homogenized term order MWH in D */ + MWH = newmat(N2+2,N2+1); + for ( J = 0; J <= N2; J++ ) + MWH[0][J] = 1; + for ( I = 1; I <= N2+1; I++ ) + for ( J = 0; J < N2; J++ ) + MWH[I][J] = MW[I-1][J]; + + /* homogenize F */ + VDVH = append(VDV,[h]); + FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH); + + /* compute a groebner basis of FH w.r.t. MWH */ + dp_gr_flags(["Top",1,"NoRA",1]); + GH = dp_weyl_gr_main(FH,VDVH,0,1,11); + dp_gr_flags(["Top",0,"NoRA",0]); + + /* dehomigenize GH */ + G = map(subst,GH,h,1); + + /* G is a groebner basis w.r.t. a non term order MW */ + /* take the initial part w.r.t. (-W,W) */ + GIN = map(initial_part,G,VDV,MW,W); + + /* GIN is a groebner basis w.r.t. a term order M */ + /* As -W+W=0, gr_(-W,W)(D) = D */ + + /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */ + for ( I = 0, T = 0; I < N; I++ ) + T += W[I]*V[I]*DV[I]; + B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */ + return B; +} + +/* all term reduction + interreduce */ +def generic_bfct_1(F,V,DV,W) +{ + N = length(V); + N2 = N*2; + + /* If W is a list, convert it to a vector */ + if ( type(W) == 4 ) + W = newvect(length(W),W); + dp_weyl_set_weight(W); + + /* create a term order M in D (DRL) */ + M = newmat(N2,N2); + for ( J = 0; J < N2; J++ ) + M[0][J] = 1; + for ( I = 1; I < N2; I++ ) + M[I][N2-I] = -1; + + VDV = append(V,DV); + + /* create a non-term order MW in D */ + MW = newmat(N2+1,N2); + for ( J = 0; J < N; J++ ) + MW[0][J] = -W[J]; + for ( ; J < N2; J++ ) + MW[0][J] = W[J-N]; + for ( I = 1; I <= N2; I++ ) + for ( J = 0; J < N2; J++ ) + MW[I][J] = M[I-1][J]; + + /* create a homogenized term order MWH in D */ + MWH = newmat(N2+2,N2+1); + for ( J = 0; J <= N2; J++ ) + MWH[0][J] = 1; + for ( I = 1; I <= N2+1; I++ ) + for ( J = 0; J < N2; J++ ) + MWH[I][J] = MW[I-1][J]; + + /* homogenize F */ + VDVH = append(VDV,[h]); + FH = map(dp_dtop,map(dp_homo,map(dp_ptod,F,VDV)),VDVH); + + /* compute a groebner basis of FH w.r.t. MWH */ +/* dp_gr_flags(["Top",1,"NoRA",1]); */ + GH = dp_weyl_gr_main(FH,VDVH,0,1,11); +/* dp_gr_flags(["Top",0,"NoRA",0]); */ + + /* dehomigenize GH */ + G = map(subst,GH,h,1); + + /* G is a groebner basis w.r.t. a non term order MW */ + /* take the initial part w.r.t. (-W,W) */ + GIN = map(initial_part,G,VDV,MW,W); + + /* GIN is a groebner basis w.r.t. a term order M */ + /* As -W+W=0, gr_(-W,W)(D) = D */ + + /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */ + for ( I = 0, T = 0; I < N; I++ ) + T += W[I]*V[I]*DV[I]; + B = weyl_minipoly(GIN,VDV,0,T); /* M represents DRL order */ + return B; +} + +def initial_part(F,V,MW,W) +{ + N2 = length(V); + N = N2/2; + dp_ord(MW); + DF = dp_ptod(F,V); + R = dp_hm(DF); + DF = dp_rest(DF); + + E = dp_etov(R); + for ( I = 0, TW = 0; I < N; I++ ) + TW += W[I]*(-E[I]+E[N+I]); + RW = TW; + + for ( ; DF; DF = dp_rest(DF) ) { + E = dp_etov(DF); + for ( I = 0, TW = 0; I < N; I++ ) + TW += W[I]*(-E[I]+E[N+I]); + if ( TW == RW ) + R += dp_hm(DF); + else if ( TW < RW ) + break; + else + error("initial_part : cannot happen"); + } + return dp_dtop(R,V); + +} + +/* b-function of F ? */ + +def bfct(F) +{ + V = vars(F); + N = length(V); + D = newvect(N); + for ( I = 0; I < N; I++ ) - M[I+2+N+1][N1-1-I] = -1; -#elif 1 - for ( I = 0; I < N1-1; I++ ) - M[I+2][N1-I-1] = 1; -#else - for ( I = 0; I < N1-1; I++ ) - M[I+2][I] = 1; -#endif + D[I] = [deg(F,V[I]),V[I]]; + qsort(D,compare_first); + for ( V = [], I = 0; I < N; I++ ) + V = cons(D[I][1],V); + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); V1 = cons(s,V); DV1 = cons(ds,DV); - G5 = dp_gr_main(cons(F,G4),append(V1,DV1),0,0,M); - for ( T = G5, G6 = []; T != []; T = cdr(T) ) { - E = car(T); - if ( intersection(vars(E),DV1) == [] ) - G6 = cons(E,G6); + + G0 = indicial1(F,reverse(V)); + G1 = dp_weyl_gr_main(G0,append(V1,DV1),0,1,0); + Minipoly = weyl_minipoly(G1,append(V1,DV1),0,s); + return Minipoly; +} + +/* b-function computation via generic_bfct() (experimental) */ + +def bfct_via_gbfct(F) +{ + V = vars(F); + N = length(V); + D = newvect(N); + + for ( I = 0; I < N; I++ ) + D[I] = [deg(F,V[I]),V[I]]; + qsort(D,compare_first); + for ( V = [], I = 0; I < N; I++ ) + V = cons(D[I][1],V); + V = reverse(V); + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + + B = [t-F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+diff(F,V[I])*dt,B); } - ctrl("do_weyl",0); - G6_0 = remove_zero(map(z_subst,G6,V)); - F0 = flatmf(cdr(fctr(dp_gr_main(G6_0,[s],0,0,0)[0]))); - for ( T = F0, BF = []; T != []; T = cdr(T) ) { - FI = car(T); - for ( J = 1; ; J++ ) { - S = map(srem,map(z_subst,idealquo(G6,[FI^J],V1,0),V),FI); - for ( ; S != [] && !car(S); S = cdr(S) ); - if ( S != [] ) + V1 = cons(t,V); DV1 = cons(dt,DV); + W = newvect(N+1); + W[0] = 1; + R = generic_bfct_1(B,V1,DV1,W); + + return subst(R,s,-s-1); +} + +/* use an order s.t. [t,x,y,z,...,dt,dx,dy,dz,...,h] */ + +def bfct_via_gbfct_weight(F,V) +{ + N = length(V); + D = newvect(N); + Wt = getopt(weight); + if ( type(Wt) != 4 ) { + for ( I = 0, Wt = []; I < N; I++ ) + Wt = cons(1,Wt); + } + Tdeg = w_tdeg(F,V,Wt); + WtV = newvect(2*(N+1)+1); + WtV[0] = Tdeg; + WtV[N+1] = 1; + /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */ + for ( I = 1; I <= N; I++ ) { + WtV[I] = Wt[I-1]; + WtV[N+1+I] = Tdeg-Wt[I-1]+1; + } + WtV[2*(N+1)] = 1; + dp_set_weight(WtV); + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + + B = [t-F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+diff(F,V[I])*dt,B); + } + V1 = cons(t,V); DV1 = cons(dt,DV); + W = newvect(N+1); + W[0] = 1; + R = generic_bfct_1(B,V1,DV1,W); + dp_set_weight(0); + return subst(R,s,-s-1); +} + +/* use an order s.t. [x,y,z,...,t,dx,dy,dz,...,dt,h] */ + +def bfct_via_gbfct_weight_1(F,V) +{ + N = length(V); + D = newvect(N); + Wt = getopt(weight); + if ( type(Wt) != 4 ) { + for ( I = 0, Wt = []; I < N; I++ ) + Wt = cons(1,Wt); + } + Tdeg = w_tdeg(F,V,Wt); + WtV = newvect(2*(N+1)+1); + /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */ + for ( I = 0; I < N; I++ ) { + WtV[I] = Wt[I]; + WtV[N+1+I] = Tdeg-Wt[I]+1; + } + WtV[N] = Tdeg; + WtV[2*N+1] = 1; + WtV[2*(N+1)] = 1; + dp_set_weight(WtV); + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + + B = [t-F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+diff(F,V[I])*dt,B); + } + V1 = append(V,[t]); DV1 = append(DV,[dt]); + W = newvect(N+1); + W[N] = 1; + R = generic_bfct(B,V1,DV1,W); + dp_set_weight(0); + return subst(R,s,-s-1); +} + +def bfct_via_gbfct_weight_2(F,V) +{ + N = length(V); + D = newvect(N); + Wt = getopt(weight); + if ( type(Wt) != 4 ) { + for ( I = 0, Wt = []; I < N; I++ ) + Wt = cons(1,Wt); + } + Tdeg = w_tdeg(F,V,Wt); + + /* a weight for the first GB computation */ + /* [t,x1,...,xn,dt,dx1,...,dxn,h] */ + WtV = newvect(2*(N+1)+1); + WtV[0] = Tdeg; + WtV[N+1] = 1; + WtV[2*(N+1)] = 1; + /* wdeg(V[I])=Wt[I], wdeg(DV[I])=Tdeg-Wt[I]+1 */ + for ( I = 1; I <= N; I++ ) { + WtV[I] = Wt[I-1]; + WtV[N+1+I] = Tdeg-Wt[I-1]+1; + } + dp_set_weight(WtV); + + /* a weight for the second GB computation */ + /* [x1,...,xn,t,dx1,...,dxn,dt,h] */ + WtV2 = newvect(2*(N+1)+1); + WtV2[N] = Tdeg; + WtV2[2*N+1] = 1; + WtV2[2*(N+1)] = 1; + for ( I = 0; I < N; I++ ) { + WtV2[I] = Wt[I]; + WtV2[N+1+I] = Tdeg-Wt[I]+1; + } + + for ( I = N-1, DV = []; I >= 0; I-- ) + DV = cons(strtov("d"+rtostr(V[I])),DV); + + B = [t-F]; + for ( I = 0; I < N; I++ ) { + B = cons(DV[I]+diff(F,V[I])*dt,B); + } + V1 = cons(t,V); DV1 = cons(dt,DV); + V2 = append(V,[t]); DV2 = append(DV,[dt]); + W = newvect(N+1,[1]); + dp_weyl_set_weight(W); + + VDV = append(V1,DV1); + N1 = length(V1); + N2 = N1*2; + + /* create a non-term order MW in D */ + MW = newmat(N2+1,N2); + for ( J = 0; J < N1; J++ ) { + MW[0][J] = -W[J]; MW[0][N1+J] = W[J]; + } + for ( J = 0; J < N2; J++ ) MW[1][J] = 1; + for ( I = 2; I <= N2; I++ ) MW[I][N2-I+1] = -1; + + /* homogenize F */ + VDVH = append(VDV,[h]); + FH = map(dp_dtop,map(dp_homo,map(dp_ptod,B,VDV)),VDVH); + + /* compute a groebner basis of FH w.r.t. MWH */ + GH = dp_weyl_gr_main(FH,VDVH,0,1,11); + + /* dehomigenize GH */ + G = map(subst,GH,h,1); + + /* G is a groebner basis w.r.t. a non term order MW */ + /* take the initial part w.r.t. (-W,W) */ + GIN = map(initial_part,G,VDV,MW,W); + + /* GIN is a groebner basis w.r.t. a term order M */ + /* As -W+W=0, gr_(-W,W)(D) = D */ + + /* find b(W1*x1*d1+...+WN*xN*dN) in Id(GIN) */ + for ( I = 0, T = 0; I < N1; I++ ) + T += W[I]*V1[I]*DV1[I]; + + /* change of ordering from VDV to VDV2 */ + VDV2 = append(V2,DV2); + dp_set_weight(WtV2); + GIN2 = dp_weyl_gr_main(GIN,VDV2,0,-1,0); + + R = weyl_minipoly(GIN2,VDV2,0,T); /* M represents DRL order */ + dp_set_weight(0); + return subst(R,s,-s-1); +} + +def weyl_minipolym(G,V,O,M,V0) +{ + N = length(V); + Len = length(G); + dp_ord(O); + setmod(M); + PS = newvect(Len); + PS0 = newvect(Len); + + for ( I = 0, T = G; T != []; T = cdr(T), I++ ) + PS0[I] = dp_ptod(car(T),V); + for ( I = 0, T = G; T != []; T = cdr(T), I++ ) + PS[I] = dp_mod(dp_ptod(car(T),V),M,[]); + + for ( I = Len - 1, GI = []; I >= 0; I-- ) + GI = cons(I,GI); + + U = dp_mod(dp_ptod(V0,V),M,[]); + U = dp_weyl_nf_mod(GI,U,PS,1,M); + + T = dp_mod(<<0>>,M,[]); + TT = dp_mod(dp_ptod(1,V),M,[]); + G = H = [[TT,T]]; + + for ( I = 1; ; I++ ) { + if ( dp_gr_print() ) + print(".",2); + T = dp_mod(<>,M,[]); + + TT = dp_weyl_nf_mod(GI,dp_weyl_mul_mod(TT,U,M),PS,1,M); + H = cons([TT,T],H); + L = dp_lnf_mod([TT,T],G,M); + if ( !L[0] ) { + if ( dp_gr_print() ) + print(""); + return dp_dtop(L[1],[t]); /* XXX */ + } else + G = insert(G,L); + } +} + +def weyl_minipoly(G0,V0,O0,P) +{ + HM = hmlist(G0,V0,O0); + + N = length(V0); + Len = length(G0); + dp_ord(O0); + PS = newvect(Len); + for ( I = 0, T = G0, HL = []; T != []; T = cdr(T), I++ ) + PS[I] = dp_ptod(car(T),V0); + for ( I = Len - 1, GI = []; I >= 0; I-- ) + GI = cons(I,GI); + DP = dp_ptod(P,V0); + + for ( I = 0; ; I++ ) { + Prime = lprime(I); + if ( !valid_modulus(HM,Prime) ) + continue; + MP = weyl_minipolym(G0,V0,O0,Prime,P); + D = deg(MP,var(MP)); + + NFP = weyl_nf(GI,DP,1,PS); + NF = [[dp_ptod(1,V0),1]]; + LCM = 1; + + for ( J = 1; J <= D; J++ ) { + if ( dp_gr_print() ) + print(".",2); + NFPrev = car(NF); + NFJ = weyl_nf(GI, + dp_weyl_mul(NFP[0],NFPrev[0]),NFP[1]*NFPrev[1],PS); + NFJ = remove_cont(NFJ); + NF = cons(NFJ,NF); + LCM = ilcm(LCM,NFJ[1]); + } + if ( dp_gr_print() ) + print(""); + U = NF[0][0]*idiv(LCM,NF[0][1]); + Coef = []; + for ( J = D-1; J >= 0; J-- ) { + Coef = cons(strtov("u"+rtostr(J)),Coef); + U += car(Coef)*NF[D-J][0]*idiv(LCM,NF[D-J][1]); + } + + for ( UU = U, Eq = []; UU; UU = dp_rest(UU) ) + Eq = cons(dp_hc(UU),Eq); + M = etom([Eq,Coef]); + B = henleq(M,Prime); + if ( dp_gr_print() ) + print(""); + if ( B ) { + R = 0; + for ( I = 0; I < D; I++ ) + R += B[0][I]*s^I; + R += B[1]*s^D; + return R; + } + } +} + +def weyl_nf(B,G,M,PS) +{ + for ( D = 0; G; ) { + for ( U = 0, L = B; L != []; L = cdr(L) ) { + if ( dp_redble(G,R=PS[car(L)]) > 0 ) { + GCD = igcd(dp_hc(G),dp_hc(R)); + CG = idiv(dp_hc(R),GCD); CR = idiv(dp_hc(G),GCD); + U = CG*G-dp_weyl_mul(CR*dp_subd(G,R),R); + if ( !U ) + return [D,M]; + D *= CG; M *= CG; break; + } } - BF = cons([FI,J],BF); + if ( U ) + G = U; + else { + D += dp_hm(G); G = dp_rest(G); + } } - return BF; + return [D,M]; } +def weyl_nf_mod(B,G,PS,Mod) +{ + for ( D = 0; G; ) { + for ( U = 0, L = B; L != []; L = cdr(L) ) { + if ( dp_redble(G,R=PS[car(L)]) > 0 ) { + CR = dp_hc(G)/dp_hc(R); + U = G-dp_weyl_mul_mod(CR*dp_mod(dp_subd(G,R),Mod,[]),R,Mod); + if ( !U ) + return D; + break; + } + } + if ( U ) + G = U; + else { + D += dp_hm(G); G = dp_rest(G); + } + } + return D; +} + def remove_zero(L) { for ( R = []; L != []; L = cdr(L) ) @@ -138,6 +777,17 @@ def v_factorial(V,N) { for ( J = N-1, R = 1; J >= 0; J-- ) R *= V-J; + return R; +} + +def w_tdeg(F,V,W) +{ + dp_set_weight(newvect(length(W),W)); + T = dp_ptod(F,V); + for ( R = 0; T; T = cdr(T) ) { + D = dp_td(T); + if ( D > R ) R = D; + } return R; } end$