# $Source: /u/maple/research/lib/algcurves/src/RCS/g_expand,v $
# $Notify: hoeij@sci.kun.nl $

# Purpose of this file: (was first called: ground_field)
# Procedures to allow code that works over different kinds of ground
# fields. So for example g_expand looks what kind of expression is given,
# and then decides what to do, expand, or evala(Expand()), or expanding
# with normalization of the lowest coefficient etcetera.

macro(
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	g_normal=`algcurves/g_normal`,
	g_expand=`algcurves/g_expand`,
	normal_tcoeff=`algcurves/normal_tcoeff`,
	g_evala=`algcurves/g_evala`,
	g_evala_rem=`algcurves/g_evala_rem`,
	g_zero_of=`algcurves/g_zero_of`,
	g_factors=`algcurves/g_factors`,
	rootof=`algcurves/rootof`,
	g_ext=`algcurves/g_ext`,
	g_ext_r=`algcurves/g_ext_r`,
	truncate=`algcurves/truncate`
):

	

g_conversion1:={}: # RootOf syntax -> rootof syntax
g_conversion2:={}: # rootof syntax -> RootOf syntax

# evala(Normal( ))
g_normal:=proc(aa)
	global g_conversion1,g_conversion2;
	local a;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	if has(aa,RootOf) then 
		RETURN(evala(Normal(aa),'independent'))
	fi;
	a:=subs(g_conversion2,aa);
	if has(a,RootOf) then
		subs(g_conversion1,evala(Normal(a),'independent'))
	else
		normal(a)
	fi
end:

# evala(Expand( ))
# Give 3'rd argument x for normalizing the lowest coefficient in x
# So then ldegree(result,x) gives the right answer
g_expand:=proc(a,ext)
	global g_conversion1,g_conversion2;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	if type(a,polynom) then
		g_evala(expand(a),ext)
	elif nargs=3 and type(a,polynom(anything,ext)) then
		normal_tcoeff(g_evala(expand(a),ext),args[3])
	elif nargs=3 then normal_tcoeff(procname(a,ext),args[3])
	else
		subs(g_conversion1,evala(Expand(subs(g_conversion2,a))
		 ,'independent'))
	fi
end:

# Normalize the lowest coefficient
normal_tcoeff:=proc(a,x)
	local c,i,r,d;
	# if not testeq(subs(g_conversion2,tcoeff(a,x))) then RETURN(a) fi;
	# Strangely enough this testeq line doesn't seem to speed things up in
	# the test examples
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	r:=a;
do
	d:=ldegree(r,x);
	c:=g_normal(coeff(r,x,d));
	if coeff(r,x,d)=c then RETURN(r) fi;
	r:=convert([c*x^d,seq(coeff(r,x,i)*x^i,i=d+1..degree(r,x))],`+`);
	if c<>0 then RETURN(r) fi
od
end:

# Input: A polynomial in the rootofs in ext
# Output: a evala normalization
g_evala:=proc(a,ext)
	local i,e;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	if nops(ext)=0 then RETURN(a)
	elif nops(ext)=1 then
		e:=ext[1];
		expand(convert([seq(coeff(a,e,i)*g_evala_rem(e^i)
		 ,i=0..degree(a,e))],`+`))
	else
	e:=g_evala(a,[op(2..nops(ext),ext)]);
	g_evala(expand(convert([seq(coeff(e,ext[1],i)*g_evala_rem(ext[1]^i)
	 ,i=0..degree(e,ext[1]))],`+`)),[op(2..nops(ext),ext)])
	fi
end:

g_evala_rem:=proc()
	global g_conversion1,g_conversion2;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	expand(subs(g_conversion1,evala(Expand(subs(g_conversion2,args))
	,'independent')))
end:

# Input : an irreducible polynomial kk in x, not necessarily monic
# Output: a root of k
# If an algebraic extension is needed it will be placed in ext.
g_zero_of:=proc(k,x,ext)
	global g_conversion1,g_conversion2;
	local a,v;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	if degree(k,x)=1 then
		ext:=NULL;
		RETURN(g_normal(-coeff(k,x,0)/coeff(k,x,1)))
	fi;
	a:=RootOf(subs(g_conversion2,k),x);
	if not member(_Z,indets(op(subs(g_conversion1,a)))) then
		a:=subs(g_conversion1,a);
		ext:=a;
		RETURN(a)
	fi;
	if g_conversion1={} then g_conversion1:=NULL fi;
	v:=nops(g_conversion2);
	g_conversion1:=a=rootof.v,g_conversion1;
	g_conversion2:={rootof.v=a,op(g_conversion2)};
	ext:=rootof.v
end:

# Factorization over the groundfield
# Both RootOf and rootof syntax
# ext: describes the groundfield
g_factors:=proc(f,x,ext)
	global g_conversion1,g_conversion2;
	local v,w,i;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	if not has(f,x) then RETURN([])
	elif ext<>[] and not has(ext,RootOf) then
		# rootof syntax input, so rootof syntax output
		RETURN(subs(g_conversion1,
		procname(op(subs(g_conversion2,[args])))))
	fi;
	userinfo(5,'algcurves',` Factorizing`,f);
	v:=traperror(factors(numer(f),{op(ext)}));
	if v=lasterror then
		# Now try evala(Factor()), it can handle more
		# cases than factors
		g_ext([args]);
		v:=subs(g_conversion2,factors(numer(subs(g_conversion1
		 ,evala(Factor(numer(f),{op(ext)}))))))
	fi;
	v:=v[2];
	userinfo(5,'algcurves',`Done factorization`);
	w:=NULL;
	for i in v do if has(i,x) then w:=w,i fi od;
	[w]
end:

# Gives the algebraic extensions appearing in aa.
g_ext:=proc(aa)
	global g_conversion1,g_conversion2;
	local v,i,result,ii,vv;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	v:=g_ext_r(subs(g_conversion2,aa));
	vv:=subs(g_conversion1,v);
	result:=NULL;
	for i from 0 to nops(g_conversion2) do if member(rootof.i,vv)
		then result:=rootof.i,result
	fi od;
	for i from nops(v) by -1 to 1 do if not member(subs(g_conversion1,v[i])
	 ,{seq(rootof.ii,ii=0..nops(g_conversion2))}) then
		if g_conversion1={} then g_conversion1:=NULL fi;
		vv:=nops(g_conversion2);
		g_conversion1:=v[i]=rootof.vv,g_conversion1;
		g_conversion2:={rootof.vv=v[i],op(g_conversion2)};
		result:=subs(g_conversion1,v[i]),result
	fi od;
	[result]
end:

# used by g_ext
g_ext_r:=proc(a)
local v,vv,i,tail;
options remember, 
    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	v:=indets(a,RootOf);
	if nops(v)=0 then RETURN([]) fi;
	vv:={};
	for i in v do vv:=vv union indets(op(1,i),RootOf) od;
	tail:=g_ext_r(vv);
	v:=[op(v minus vv)];
	[op(v),op(tail)]
end:

truncate:=proc(aa,n,y,ext)
local dummy,a;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	a:=collect(aa,y);
	a:=expand(convert([seq(y^dummy*coeff(a,y,dummy)
	 ,dummy=ldegree(a,y)..n-1)],`+`));
	g_evala(a,ext)
end:

#--------------------
# I moved the following procedures from the files puiseux,
#  integral_basis, integral_basis and singularities to here for
#  convenience with the diffop package.

macro(
	v_ext_m=`algcurves/v_ext_m`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	g_gcdex=`algcurves/gcdex`,
	degree_ext=`algcurves/degree_ext`
):

# Gives the zeros of the factors, their multiplicities and algebraic extensions
# The input of this procedure is the output of g_factors
v_ext_m:=proc(f,x)
	local ext,nulp,i,res;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	res:=NULL;
	for i in f do
		nulp:=g_zero_of(i[1],x,'ext');
		ext:=eval(ext);
		res:=res,[nulp,i[2],[ext],degree(i[1],x)]
	od;
	[res]
end:

# ext_to_coeffs does basically: coeffs(expression,RootOf( .. ))
# We need this procedure because coeffs doesn't always work this way.
ext_to_coeffs:=proc(a,ext)
	global g_conversion2;
	local i,aa;
	option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	aa:=(indets(a) minus indets(ext)) intersect
	{seq(rootof.i,i=0..nops(g_conversion2))};
	coeffs(a,[op(aa)])
end:

g_gcdex:=proc(a,b,c,x,ext)
	global g_conversion1,g_conversion2;
	local r,ss,tt,s,t;
	options remember, 
	    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	r:=subs(g_conversion2,[a,b,c]);
	if ext=[] then
		r:=gcdex(op(r),x,'ss','tt')
	else
		r:=subs(g_conversion1,evala(Gcdex(r[1],r[2],x,'ss','tt')));
		if g_normal(r-c)<>0 then ERROR(`found wrong gcd`) fi
	fi;
	s:=subs(g_conversion1,ss);
	t:=subs(g_conversion1,tt);
	[c,s,t]
end:

# Gives the algebraic degree of a over b
degree_ext:=proc(aa,bb)
global g_conversion2;
local a,b,v,i,all,d,var;
options remember, 
    `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`;
	a:=subs(g_conversion2,aa);
	b:=subs(g_conversion2,bb);
	v:=indets(a,RootOf) minus indets(b,RootOf);
	all:=[op(indets([a,b],RootOf))];
	all:={seq(all[i]=var[i],i=1..nops(all))};
	d:=1;
	for i in v do
		d:=d*degree(subs(all,op(1,i)),_Z)
	od;
	d
end:

#savelib('g_conversion1','g_conversion2','g_normal','g_expand',\
    'normal_tcoeff','g_evala','g_evala_rem','g_zero_of','g_factors',\
    'ext_to_coeffs','v_ext_m','degree_ext','g_gcdex',\
    'g_ext','g_ext_r','truncate','`algcurves/g_expand.m`'):