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$