File: [local] / OpenXM_contrib2 / asir2000 / lib / de (download)
Revision 1.5, Tue Aug 9 05:34:10 2005 UTC (19 years, 1 month ago) by noro
Branch: MAIN
CVS Tags: R_1_3_1-2, RELEASE_1_3_1_13b, RELEASE_1_2_3_12, KNOPPIX_2006, HEAD, DEB_REL_1_2_3-9 Changes since 1.4: +4 -4
lines
Fixed a bug.
|
module de;
localf split_lexgb;
localf sort_lex_dec,sort_lex_inc;
localf inverse_or_split, linear_dim;
localf dp_monic_mod,monic_gb;
localf membership_test;
localf dp_chrem,intdptoratdp,intdpltoratdpl;
localf comp_by_ht,dp_gr_mod,gr_chrem;
localf construct_sqfrbasis;
/*
* G : a 0-dim lex gb, reduced
* if there exists a non-monic g = g(x')x^k+...
* then g(x') is a zero divisor and Id(G) splits
*/
def split_lexgb(G,V)
{
G = map(ptozp,G);
G = sort_lex_dec(G,V);
TG = G;
for ( ; TG != []; TG = cdr(TG) ) {
F = car(TG);
for ( TV = V; !deg(F,car(TV)); TV = cdr(TV) );
VF = car(TV);
DF = deg(F,VF);
CF = coef(F,DF,VF);
if ( type(CF) == 2 ) {
L = inverse_or_split(V,G,CF);
R = map(split_lexgb,L,V);
return append(R[0],R[1]);
}
}
return [G];
}
/*
* V : a variable list
* Id : a 0-dim radical ideal represented by a lex basis
* F : a poly
* if F is a unit of Q[V]/Id, then 1/F is returned
* else F is a zero divisor and Id = (Id:F)\cap(Id+<F>)
* [GB(Id:F),GB(Id+<F>)] is returned.
*/
def inverse_or_split(V,Id,F)
{
Id = map(ptozp,Id);
N = length(V);
dp_ord(2);
set_field(Id,V,2);
DF = dptodalg(dp_ptod(F,V));
Ret = inv_or_split_dalg(DF);
if ( type(Ret) == 1 ) {
/* Ret = 1/F */
Dp = dalgtodp(Ret);
return dp_dtop(Dp[0],V)/Dp[1];
} else {
/* Ret = GB(Id:F) */
/* compute GB(Id+<f>) */
Gquo = append(map(ptozp,map(dp_dtop,Ret,V)),Id);
/* inter-reduction */
Gquo = nd_gr_postproc(Gquo,V,0,2,0);
DTotal = linear_dim(Id,V,2);
Dquo = linear_dim(Gquo,V,2);
Drem = DTotal-Dquo;
B = cons(F,Id);
Grem = gr_chrem(B,V,2,Drem);
return [map(ptozp,Gquo),map(ptozp,Grem)];
}
}
/* add F(X,V) to Id(B) */
/* returns a list of splitted ideals */
/* B should be a triangular basis */
def construct_sqfrbasis(F,X,B,V)
{
if ( type(F) == 1 )
return [];
B = sort_lex_dec(B,V);
V1 = cons(X,V);
F = nd_nf(F,reverse(B),cons(X,V),2,0);
D = deg(F,X);
H = coef(F,D,X);
if ( type(H) == 2 ) {
Ret = inverse_or_split(V,B,H);
if ( type(Ret) == 4 ) {
/* H != 0 on Id_nz, H = 0 on Id_z */
B0=construct_sqfrbasis(F,X,Ret[0],V);
B1=construct_sqfrbasis(F,X,Ret[1],V);
return append(B0,B1);
} else
F = nd_nf(F*Ret,reverse(B),cons(X,V),2,0);
}
B1 = cons(F,B);
/* F is monic */
M = minipoly(B1,V1,2,X,zzz);
S = sqfr(M); S = cdr(S);
if ( length(S) == 1 && car(S)[1] == 1 )
return [cons(F,B)];
else {
R = [];
for ( T = S; T != []; T = cdr(T) ) {
G = nd_gr_trace(cons(subst(car(T)[0],zzz,X),B1),V1,1,1,2);
R1 = split_lexgb(G,V1);
R = append(R1,R);
}
return R;
}
}
def sort_lex_dec(B,V)
{
dp_ord(2);
B = map(dp_ptod,B,V);
B = vtol(qsort(ltov(B),comp_by_ht));
B = map(dp_dtop,B,V);
return reverse(B);
}
def sort_lex_inc(B,V)
{
dp_ord(2);
B = map(dp_ptod,B,V);
B = vtol(qsort(ltov(B),comp_by_ht));
B = map(dp_dtop,B,V);
return B;
}
def linear_dim(G,V,Ord)
{
dp_ord(Ord);
MB = dp_mbase(map(dp_ptod,G,V));
return length(MB);
}
def membership_test(B,G,V,O)
{
B = map(ptozp,B);
G = map(ptozp,G);
for ( T = B; T != []; T = cdr(T) )
if ( nd_nf(car(T),G,V,O,0) ) return 0;
return 1;
}
def gr_chrem(B,V,O,Dim)
{
B = map(ptozp,B);
G = []; HS = []; Mod = 1;
for ( I = 0; ; I++ ) {
P = lprime(I);
GM = nd_gr(B,V,P,O);
if ( linear_dim(GM,V,O) != Dim ) continue;
L = monic_gb(GM,V,O,P); GM = L[0]; HSM = L[1];
G = dp_chrem(G,HS,Mod,GM,HSM,P);
Mod *= P;
if ( G != [] )
HS = HSM;
M = idiv(isqrt(2*Mod),2);
R1 = intdpltoratdpl(G,Mod,M);
if ( R1 ) {
if ( Found && R == R1
&& (GB=nd_gr_postproc(map(dp_dtop,R,V),V,0,O,1))
&& membership_test(B,GB,V,O) )
break;
else {
R = R1; Found = 1;
}
}
}
return GB;
}
def comp_by_ht(A,B)
{
HA = dp_ht(A); HB = dp_ht(B);
if ( HA > HB )
return 1;
else if ( HA < HB )
return -1;
else
return 0;
}
def intdpltoratdpl(G,Mod,M)
{
for ( R = []; G != []; G = cdr(G) ) {
T = intdptoratdp(car(G),Mod,M);
if ( !T )
return 0;
else
R = cons(T,R);
}
R = reverse(R);
return vtol(qsort(newvect(length(R),R),comp_by_ht));
}
def intdptoratdp(F,Mod,M)
{
for ( T = F, N = 0; T; T = dp_rest(T), N++ );
C = newvect(N);
for ( I = 0, T = F; I < N; T = dp_rest(T), I++ ) {
L = inttorat(dp_hc(T),Mod,M);
if ( !L )
return 0;
else
C[I] = (L[0]/L[1])*dp_ht(T);
}
for ( R = 0, I = N-1; I >= 0; I-- )
R += C[I];
return R;
}
def dp_chrem(G,HS,Mod,GM,HSM,P)
{
if ( G == [] )
return GM;
if ( HS != HSM )
return [];
R = [];
M1 = inv(Mod,P);
ModM1 = Mod*M1;
for ( ; G != []; G = cdr(G), GM = cdr(GM) ) {
E = car(G); EM = car(GM);
E1 = E+(EM-E)*ModM1;
R = cons(E1,R);
}
return reverse(R);
}
def monic_gb(G,V,O,P)
{
dp_ord(O); setmod(P);
D = map(dp_ptod,G,V);
D = map(dp_monic_mod,D,P);
D = vtol(qsort(newvect(length(D),D),comp_by_ht));
return [D,map(dp_ht,D)];
}
def dp_monic_mod(F,P)
{
FP = dp_mod(F,P,[]);
return dp_rat(FP/dp_hc(FP));
}
endmodule;
end$