=================================================================== RCS file: /home/cvs/OpenXM_contrib2/asir2000/lib/bfct,v retrieving revision 1.5 retrieving revision 1.6 diff -u -p -r1.5 -r1.6 --- OpenXM_contrib2/asir2000/lib/bfct 2000/12/11 02:00:42 1.5 +++ OpenXM_contrib2/asir2000/lib/bfct 2000/12/13 05:37:31 1.6 @@ -45,11 +45,11 @@ * 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.4 2000/12/08 08:26:09 noro Exp $ + * $OpenXM: OpenXM_contrib2/asir2000/lib/bfct,v 1.5 2000/12/11 02:00:42 noro Exp $ */ /* requires 'primdec' */ -/* annihilating ideal of F^s ? */ +/* annihilating ideal of F^s */ def ann(F) { @@ -77,68 +77,290 @@ def ann(F) return G4; } +def indicial1(F) +{ + V = vars(F); + 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); + G0 = dp_weyl_gr_main(B,append(W,DW),0,0,6); + G1 = map(subst,G0,y1,1); + Mat = newmat(2,2,[[-1,1],[0,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; +} + /* b-function of F ? */ def bfct(F) { - G4 = ann(F); - + G4 = indicial1(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); for ( I = N-1, DV = []; I >= 0; I-- ) DV = cons(strtov("d"+rtostr(V[I])),DV); + V1 = cons(s,V); DV1 = cons(ds,DV); + G0 = dp_weyl_gr_main(G4,append(V1,DV1),0,1,0); + Minipoly = weyl_minipoly(G0,append(V1,DV1),0,s); + return Minipoly; +} - N1 = 2*(N+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); - M = newmat(N1+1,N1); - for ( J = N+1; J < N1; 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 = 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 - V1 = cons(s,V); DV1 = cons(ds,DV); - dp_nelim(0); -/* G4 = dp_weyl_gr_main(G4,append(V1,DV1),0,0,10); */ - for ( PrimeIndex = 0; ; PrimeIndex++ ) { - Prime = lprime(PrimeIndex); - dp_nelim(0); /* XXX */ - Success = dp_weyl_gr_main(cons(F,G4),append(V1,DV1),0,Prime,10); - if ( !Success ) - continue; - dp_nelim(N+1); - G5 = dp_weyl_gr_main(cons(F,G4),append(V1,DV1),0,-Prime,10); - if ( G5 ) - break; + 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,[]); + + T = dp_mod(<<0>>,M,[]); + TT = dp_mod(dp_ptod(1,V),M,[]); + G = H = [[TT,T]]; + + for ( I = 1; ; I++ ) { + 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] ) + return dp_dtop(L[1],[V0]); + else + G = insert(G,L); } - for ( T = G5, G6 = []; T != []; T = cdr(T) ) { - E = car(T); - if ( intersection(vars(E),DV1) == [] ) - G6 = cons(E,G6); +} + +def weyl_minipoly(G0,V0,O0,V) +{ + for ( I = 0; ; I++ ) { + Prime = lprime(I); + MP = weyl_minipolym(G0,V0,O0,Prime,V); + for ( D = deg(MP,V), TL = [], J = 0; J <= D; J++ ) + TL = cons(V^J,TL); + dp_ord(O0); + NF = weyl_gennf(G0,TL,V0,O0)[0]; + + LHS = weyl_nf_tab(-car(TL),NF,V0); + B = weyl_hen_ttob(cdr(TL),NF,LHS,V0,Prime); + if ( B ) { + R = ptozp(B[1]*car(TL)+B[0]); + return R; + } } - 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 != [] ) +} + +def weyl_gennf(G,TL,V,O) +{ + N = length(V); Len = length(G); dp_ord(O); PS = newvect(Len); + for ( I = 0, T = G, HL = []; T != []; T = cdr(T), I++ ) { + PS[I] = dp_ptod(car(T),V); HL = cons(dp_ht(PS[I]),HL); + } + for ( I = 0, DTL = []; TL != []; TL = cdr(TL) ) + DTL = cons(dp_ptod(car(TL),V),DTL); + for ( I = Len - 1, GI = []; I >= 0; I-- ) + GI = cons(I,GI); + T = car(DTL); DTL = cdr(DTL); + H = [weyl_nf(GI,T,T,PS)]; + + T0 = time()[0]; + while ( DTL != [] ) { + T = car(DTL); DTL = cdr(DTL); + if ( dp_gr_print() ) + print(".",2); + if ( L = search_redble(T,H) ) { + DD = dp_subd(T,L[1]); + NF = weyl_nf(GI,dp_weyl_mul(L[0],dp_subd(T,L[1])),dp_hc(L[1])*T,PS); + } else + NF = weyl_nf(GI,T,T,PS); + NF = remove_cont(NF); + H = cons(NF,H); + } + print(""); + TNF = time()[0]-T0; + if ( dp_gr_print() ) + print("gennf(TAB="+rtostr(TTAB)+" NF="+rtostr(TNF)+")"); + return [H,PS,GI]; +} + +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 weyl_hen_ttob(T,NF,LHS,V,MOD) +{ + if ( length(T) == 1 ) + return car(T); + T0 = time()[0]; M = etom(weyl_leq_nf(T,NF,LHS,V)); TE = time()[0] - T0; + T0 = time()[0]; U = henleq(M,MOD); TH = time()[0] - T0; + if ( dp_gr_print() ) { + print("(etom="+rtostr(TE)+" hen="+rtostr(TH)+")"); + } + return U ? vtop(T,U,LHS) : 0; +} + +def weyl_leq_nf(TL,NF,LHS,V) +{ + TLen = length(NF); + T = newvect(TLen); M = newvect(TLen); + for ( I = 0; I < TLen; I++ ) { + T[I] = dp_ht(NF[I][1]); + M[I] = dp_hc(NF[I][1]); + } + Len = length(TL); INDEX = newvect(Len); COEF = newvect(Len); + for ( L = TL, J = 0; L != []; L = cdr(L), J++ ) { + D = dp_ptod(car(L),V); + for ( I = 0; I < TLen; I++ ) + if ( D == T[I] ) + break; + INDEX[J] = I; COEF[J] = strtov("u"+rtostr(J)); + } + if ( !LHS ) { + COEF[0] = 1; NM = 0; DN = 1; + } else { + NM = LHS[0]; DN = LHS[1]; + } + for ( J = 0, S = -NM; J < Len; J++ ) { + DNJ = M[INDEX[J]]; + GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD; + S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J]; + DN *= CS; + } + for ( D = S, E = []; D; D = dp_rest(D) ) + E = cons(dp_hc(D),E); + BOUND = LHS ? 0 : 1; + for ( I = Len - 1, W = []; I >= BOUND; I-- ) + W = cons(COEF[I],W); + return [E,W]; +} + +def weyl_nf_tab(A,NF,V) +{ + TLen = length(NF); + T = newvect(TLen); M = newvect(TLen); + for ( I = 0; I < TLen; I++ ) { + T[I] = dp_ht(NF[I][1]); + M[I] = dp_hc(NF[I][1]); + } + A = dp_ptod(A,V); + for ( Z = A, Len = 0; Z; Z = dp_rest(Z), Len++ ); + INDEX = newvect(Len); COEF = newvect(Len); + for ( Z = A, J = 0; Z; Z = dp_rest(Z), J++ ) { + D = dp_ht(Z); + for ( I = 0; I < TLen; I++ ) + if ( D == T[I] ) + break; + INDEX[J] = I; COEF[J] = dp_hc(Z); + } + for ( J = 0, S = 0, DN = 1; J < Len; J++ ) { + DNJ = M[INDEX[J]]; + GCD = igcd(DN,DNJ); CS = DNJ/GCD; CJ = DN/GCD; + S = CS*S + CJ*NF[INDEX[J]][0]*COEF[J]; + DN *= CS; + } + return [S,DN]; } def remove_zero(L)