macro(
	# tau=_EnvLREdomain[1],
	# x=_EnvLREdomain[2],

	adjoint=`LREtools/adjoint`,
	mult=`LREtools/mult`,
	rightdivision=`LREtools/rightdivision`,
	GCRD=`LREtools/GCRD`,
	LRE2op=`LREtools/LRE2op`,
	op2LRE=`LREtools/op2LRE`
):

adjoint:=proc(L)
	local i;
	global x,tau;
	collect(add(mult(tau^i,subs(x=-x,coeff(L,tau,i)))
	 ,i=ldegree(L,tau)..degree(L,tau)),[tau,x])
end:
#savelib('adjoint'):

LRE2op:=proc(eq::{algebraic,algebraic=algebraic},fn::function(name))
	local R,i;
	global x,tau;
	R:=op(4,LREtools[REcreate](args));
	if not (R['linear'] and R['polycoeffs']) then
	 error `Expecting linear recurrence relation with rational functions coeffs`
	fi;
	add(R['coeffs'][R[1+i-'min_shift']]*tau^i,i=R['min_shift']..R['max_shift'])
end:
#savelib('LRE2op'):

op2LRE:=proc(L,fn::function(name))
	local i,u,n;
	global x,tau;
	if nargs=1 then
		procname(L,u(n))
	else
		LREtools['REcreate'](subs(x=op(fn),
		 add(op(0,fn)(op(fn)+i)*coeff(L,tau,i),i=ldegree(L,tau)..degree(L,tau))
		),fn,{})
	fi
end:
#savelib('op2LRE'):

mult:=proc(L1,L2)
	local i;
	global x,tau;
	if nargs=0 then
		1
	elif nargs=1 then
		L1
	elif nargs=2 then
		sort(collect(add(coeff(L1,tau,i)*
		 subs(x=x+i,L2)*tau^i,i=0..degree(L1,tau))
		,tau,Normalizer),tau)
	else
		procname(L1,procname(args[2..-1]))
	fi
end:
#savelib('mult'):

# result: [a,b] such that a*right+b=f
rightdivision:=proc(f,right)
	local a,b,t;
	global x,tau;
	a:=0;
	b:=f;
	if lcoeff(right,tau)<>1 then
		a:=procname(f,collect(right/lcoeff(right,tau),tau,Normalizer));
		return [mult(a[1],1/lcoeff(right,tau)),a[2]]
	elif right=1 then
		return [f,0]
	fi;
	userinfo(5,'LREtools',`right division`,f,right);
	while degree(b,tau)>=degree(right,tau) do
		a:=a+lcoeff(b,tau)*tau^(degree(b,tau)-degree(right,tau));
		b:=collect(b-lcoeff(b,tau)*mult(seq(tau
		 ,t=1..degree(b,tau)-degree(right,tau)),right),tau,Normalizer)
	od;
	userinfo(5,'LREtools',`done right division`);
	[a,b]
end:
#savelib('rightdivision'):

GCRD:=proc(f,g)
	global tau;
	local a;
if f=0 or g=0 then
	f+g
else
	a:=collect(g/lcoeff(g,tau),tau,Normalizer);
	procname(a,rightdivision(f,a)[2])
fi
end:
#savelib('GCRD'):


macro(
	LCLM=`LREtools/LCLM`,
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	degree_ext=`algcurves/degree_ext`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	g_ext=`DEtools/diffop/g_ext`,
	g_normal=`DEtools/diffop/g_normal`,
	g_expand=`algcurves/g_expand`,
	solve=`solve/linear`
):

# Code is the same as in diffop (I should make only 1 procedure
# to handle both cases).
#
# Input: sequence of operators
# or:
# sequence of operators,`and conjugates`,ext
# ext: optional, if not present then ext=[]
# Output: the LCLM of the operators and their conjugates over Q(ext)
# RootOf syntax
LCLM:=proc(L)
	local ext,j,i,ind,sol,ansatz,a;
	global x,tau;
	ext:=g_ext([args]);
if has([args],`and conjugates`) then
	if not type(args[nargs],list) then
		# ground field is Q
		RETURN(procname(args,[]))
	fi;
	if nargs>3 then
	 return procname(seq(procname(args[j],args[nargs-1..nargs]),j=1..nargs-2))
	fi;
	j:=subs(g_conversion1,args[3]);
	ansatz:=convert([seq(a[i]*tau^i,i=0..degree(L,tau)
	 *degree_ext(ext,j))],`+`);
	sol:={seq(coeffs(i,tau),i=[ext_to_coeffs(g_expand(subs(g_conversion1,
	 numer(rightdivision(ansatz,L)[2])),ext),j)])}
else
	if type(args[nargs],list) then
		# ext should not be given
		return procname(args[1..nargs-1])
	elif nargs=1 then
		return collect(L/lcoeff(L,tau),tau,g_normal)
	fi;
	ansatz:=convert([seq(a[i]*tau^i,i=0..convert([seq(degree(j,tau),j=args)]
	,`+`))],`+`);
	sol:={seq(coeffs(g_expand(numer(rightdivision(
	 ansatz,subs(g_conversion1,i))[2]),ext),tau),i=args)}
fi;
	ind:=indets(ansatz) minus {tau};
	sol:=solve(subs(g_conversion2,sol),ind);
	i:=degree(ansatz,tau);
	while nops(indets(subs(sol,ansatz)) intersect ind)>1 do
		sol:=solve({a[i]=0,op(sol)},ind);
		i:=i-1
	od;
	sol:=solve({a[i]=1,op(sol)},ind);
	collect(subs(sol,ansatz),tau,g_normal)
end:
#savelib('LCLM'):