# Compute a "generating" set, or compute all subfields.
#
# This implementation is meant to be an illustration of the idea and
# is not optimized for speed.
#
# Mark van Hoeij,   August 2004.


# B = basefield
# f is irreducible in B[x]
#
# The output is of the form
#
#   [E1,E2,..,Ek], vars, g
#
# such that gi := subs(Ei, g) is a generic element
# of some subfield Li of K := B[x]/(f).
# What is meant by that is that:
#   [coeffs(gi, indets(gi) intersect vars)]
# is a B-vectorspace-basis of that subfield Li.
#
# These subfields L1,L2,..,Lk of K := B[x]/(f) have the
# following property: Every proper subfield of K containing B
# is an intersection of one or more Li's.
GeneratingSet := proc(f,x,B)
	local n,R,V,i,c,ansatz,vars,res,g,E;
	n := degree(f,x);
	ansatz := add(c[i]*x^i,i=0..n-1);
	vars := {seq(c[i],i=0..n-1)};
	if n<4 then
		# There can be no non-trivial subfields.
		return [],vars,ansatz
	fi;
	R := RootOf(f,x);
	# The factors of f in K[x] except x-R:
	# We also throw away any factor of degree >= n/2.
	V := FactorOverK(f,x,R,B,n);
	res := NULL;
	for g in V do
		# Get the equations. Perhaps this could be done faster?
		E := subs(x=R, ansatz) - evala(Rem(ansatz,g,x));
		E := evala({coeffs(evala(Expand(E)),{x,R})});
		# Solve the equations.
		E := SolveTools:-Linear(E,vars);
		# If subfield=B then next:
		if degree(subs(E,ansatz),x)=0 then next fi;
		res := res, E
	od;
	[res], vars, ansatz
end:


# Factor f/(x-R) over the field K.
# This is a preliminary implementation, there are tricks to compute
# this much faster.
FactorOverK := proc(f,x,R,B,n)
	local v,i;
	v := evala(Factors(numer(normal(f)), [R, op(B)]))[2];
	# Throw away polys of degree >= n/2 and polys that don't
	# contain x.
	[op( {seq(`if`(has(i,x) and degree(i[1],x) < n/2,
	 collect(i[1]/lcoeff(i[1],x),x,evala), NULL),i=v)} minus {x-R} )]
end:


# Notation the same as in GeneratingSet.
# The output is a list [L1,L2,..Ls] consisting of
# all subfields of K that contain B.
# Each Li is given in the form [b.i1,b.i2,...] where the b.ij's
# form a B-vectorspace-basis of Li.
AllSubfields := proc(f,x,B)
	local L,vars,ansatz,n,k,res,curf,p,num;
	L, vars, ansatz := GeneratingSet(args);
	n := degree(f,x);
	k := nops(L);
	res := ansatz;
	curf[0] := ansatz;
	p := 1; num[p] := -1;
	# The equations of the current field:
	curf[0] := {};
	# Find all intersections with a back-tracking search:
	while p>0 do
		if p>k or num[p]=1 then
			p:=p-1; next
		fi;
		num[p]:=num[p]+1;
		if num[p]=0 then
			curf[p] := curf[p-1];
		fi;
		if num[p]=0 and Contains(L[p], curf[p]) then
			num[p] := 1
		elif num[p]=1 then
			curf[p] := IntersectFields(L[p], curf[p], vars);
			# if curf[p] equals B or is contained in some L[i]
			# with num[i]=0 then
			if CheckCondition(num, p, L, curf[p]) then
				p := p-1; next
			fi;
			res := res, collect(subs(curf[p], ansatz),vars)
		fi;
		p := p+1;  num[p] := -1
	od;
	[seq({coeffs(p, vars intersect indets(p))}, p=[res]), {1}]
end:


# Output is true if the field L contains the field F,
# where both L and F are given in the form of equations.
Contains := proc(L, F)
	local i;
	for i in subs(F,L) do if evala(Normal(rhs(i)-lhs(i)))<>0 then
		return false
	fi od;
	true
end:


# Intersect two fields given by equations.
# Output is again given by equations.
IntersectFields := proc(L1,L2,vars)
	SolveTools:-Linear(L1 union L2, vars)
end:


CheckCondition := proc(num, p, L, F)
	local i;
	if add(`if`(rhs(i)=0,0,1),i=F)=1 then
		# now F = basefield
		return true
	fi;
	for i to p do if num[i]=0 and Contains(L[i], F) then
		# F is contained in L[i]
		return true
	fi od;
	false
end: