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

# Procedures for users
#    eigenring, endomorphism_charpoly, gen_exp
# Procedures for use only in the rest of the diffop package:
#    none
# Procedures for use only in this file:
#    bounds_frlf, local_frlf0, lift_frlf0, liftstep_frlf0, eigenring_lineq


#####################################
#  The mixed differential equation  #----------------------------------------
#####################################		diffop_mixed_eqn

# Purpose of this file: computing the endomorphisms of the solution space,
# c.f. my ISSAC'96 paper.

# macro's for g_ext:
macro(
	infosubs=`DEtools/diffop/infosubs`,
	modulus=`DEtools/diffop/modulus`,
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	xDF=`DEtools/diffop/xDF`
):

macro(
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	rootof=`algcurves/rootof`,
	degree_ext=`algcurves/degree_ext`,
	g_evala=`algcurves/g_evala`,
	g_evala_rem=`algcurves/g_evala_rem`,
	g_zero_of=`algcurves/g_zero_of`,
	v_ext_m=`algcurves/v_ext_m`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	g_ext_r=`algcurves/g_ext_r`,
	g_gcdex=`algcurves/gcdex`,
	truncate=`algcurves/truncate`,
	g_factors=`algcurves/g_factors`
);
macro(
	new_laurent2=`DEtools/diffop/new_laurent2`,
	differentiate=`DEtools/diffop/differentiate`,
	g_ext=`DEtools/diffop/g_ext`,
	g_factors=`DEtools/diffop/g_factors`,
	g_normal=`DEtools/diffop/g_normal`,
	set_laurents=`DEtools/diffop/set_laurents`,
	description_laurent=`DEtools/diffop/description_laurent`,
	modm=`DEtools/diffop/modm`,
	g_solve=`DEtools/diffop/g_solve`,
	g_quotient=`DEtools/diffop/g_quotient`,
	g_rem=`DEtools/diffop/g_rem`,
	g_quo=`DEtools/diffop/g_quo`,
	truncate_x=`DEtools/diffop/truncate_x`,
	g_expand=`DEtools/diffop/g_expand`
):

# macro's for eval_laurent:
macro(	accuracy_laurent=`DEtools/diffop/accuracy_laurent`,
	value_laurent=`DEtools/diffop/value_laurent`,
	LL=_LL
):
macro(
	lowerbound_val=`DEtools/diffop/lowerbound_val`,
	max0_val=`DEtools/diffop/max0_val`,
	new_laurent=`DEtools/diffop/new_laurent`,
	nmterms_laurent=`DEtools/diffop/nmterms_laurent`,
	nterms_expression=`DEtools/diffop/nterms_expression`,
	nthterm_laurent=`DEtools/diffop/nthterm_laurent`,
	ramification_laur=`DEtools/diffop/ramification_laur`,
	lift_ramification_laur=`DEtools/diffop/lift_ramification_laur`,
	eval_laurent=`DEtools/diffop/eval_laurent`,
	subsvalueslaurents=`DEtools/diffop/subsvalueslaurents`,
	series_val=`DEtools/diffop/pseries_val`,
	valuation_laurent=`DEtools/diffop/valuation_laurent`,
	nt=`DEtools/diffop/nt`
):

# macro's for factor_op and substitute
macro(
	Newtonpolygon=`DEtools/diffop/Newtonpolygon`,
	coefs_operator=`DEtools/diffop/coefs_operator`,
	factor_newton=`DEtools/diffop/factor_newton`,
	factor_newton2=`DEtools/diffop/factor_newton2`,
	factor_op=`DEtools/diffop/factor_op`,
	factor_riccati=`DEtools/diffop/factor_riccati`,
	lift_newton=`DEtools/diffop/lift_newton`,
	lift_ramification=`DEtools/diffop/lift_ramification`,
	lift_substitute=`DEtools/diffop/lift_substitute`,
	nm_mult=`DEtools/diffop/nm_mult`,
	nm_block=`DEtools/diffop/nm_block`,
	op_with_slope=`DEtools/diffop/op_with_slope`,
	ram_laur=`DEtools/diffop/ram_laur`,
	ramification_of=`DEtools/diffop/ramification_of`,
	rem_lift=`DEtools/diffop/rem_lift`,
	skipped_factors=`DEtools/diffop/skipped_factors`,
	substitute=`DEtools/diffop/substitute`
):

macro(
	mult=readlib(`DEtools/diffop/mult`),
	l_p=readlib(`DEtools/diffop/l_p`),
	rightdivision=readlib(`DEtools/diffop/rightdivision`)
):

macro(
	gen_exp=`DEtools/diffop/gen_exp`,

	bounds_frlf0=`DEtools/diffop/bounds_frlf0`,
	eigenring=`DEtools/diffop/eigenring`,
	endomorphism_charpoly=`DEtools/diffop/endomorphism_charpoly`,
	local_frlf0=`DEtools/diffop/local_frlf0`,
	lift_frlf0=`DEtools/diffop/lift_frlf0`,
	liftstep_frlf0=`DEtools/diffop/liftstep_frlf0`,
	eigenring_lineq=`DEtools/diffop/eigen_lineq`
):

# Compute linear equations for coefficients of r, from the condition that
# the r(V(f)) \subset V(f) mod x^(n+N)
eigenring_lineq:=proc(f,r)
	local a,i,S,s,Si,n,N,var;
	global DF,x;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if g_normal(subs(x=0,denom(f)))=0 then
	RETURN(`not a regular point`)
fi;
	n:=degree(f,DF);
	N:=2; # We'll obtain N*n equations
	# Solution of f:
	S:=convert([seq(a[i]*x^i,i=0..2*n-1+N)],`+`);
	# Apply f:
	s:=0;
	Si:=S;
	for i from 0 to n do
		s:=s+convert(taylor(coeff(f,DF,i),x=0,2*n+N),polynom)*Si;
		Si:=diff(Si,x)
	od;
	s:=collect(s,x);
	S:=subs(readlib(`solve/linear`)(evala({seq(coeff(s,x,i),i=0..n-1+N)})
	 ,{seq(a[i],i=n..2*n-1+N)}),S);
	# Apply r on the solution S:
	s:=0;
	Si:=S;
	for i from 0 to n do
		s:=s+convert(taylor(coeff(r,DF,i),x=0,n+N),polynom)*Si;
		Si:=diff(Si,x)
	od;
	s:=collect(s,x);
	var:={seq(a[i]=0,i=0..n-1)};
	s:=collect(s-convert([seq(subs(var,subs(a[i]=1,S))*coeff(s,x,i)
	 ,i=0..n-1)],`+`),x);
	var:=indets(var);
	S:=[seq(coeffs(expand(coeff(s,x,i)),var),i=n..n-1+N)];
	{seq(g_normal(numer(i)),i=S)}
end:

# Compute the ring of endomorphisms E_D(f), in other words the solutions r
# of the mixed equation f*r+l*f=0.
eigenring:=proc(F)
	local a,B,n,f0,L,xDFpow,i,values,zeros,r,ac,j,ne,f;
	global DF,x,xDF;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;

	f:=collect(F,DF,g_normal);
if degree(f,DF)=0 then
	ERROR(`degree input should be > 0`)
elif degree(f,DF)=1 then
	RETURN([1])
fi;
	f:=collect(f/lcoeff(f,DF),DF,g_normal);
	userinfo(2,'diffop',`Computing endomorphism ring of`,infosubs(f));
	if subs(x=0,denom(f))=0 then RETURN(map(collect,subs(x=x+1,
		procname(subs(x=x-1,f),args[2..nargs]
		 ,`DEtools/ subs`=f)),DF,g_normal))
	fi;
	# Now we can assume that x=0 is regular
	L:=select(has,[args],`DEtools/ subs`);
	if L=[] then
		B:=bounds_frlf0(f)
	else
		# it is more efficient to compute the bounds using the
		# original operator because it's quite likely that the
		# local arithmetic has already been done for that operator.
		B:=bounds_frlf0(subs(L[1],`DEtools/ subs`))
	fi;
	userinfo(2,'diffop',`bounds are:`,infosubs(B));
	if L<>[] then
		B:=map(evala@Normal,subs(x=x-nops(L),B),expanded)
	fi;
	n:=degree(f,DF);
	# nc:=convert([seq(i[2]+1,i=B)],`+`); # Number of coeffs to be determined
	f0:=l_p(f,0);
	L:=local_frlf0(f0);
	# The following suffices for expressing the candidate global solutions
	# in n^2 variables, using the local solutions up to this accuracy.
	ac:=max(seq(B[i][2]+2-i,i=1..n));
	# If we take ac we have the following number of equations
	ne:=convert([seq(ac-B[i][2]-2+i,i=1..n)],`+`);
	# We need at least n^2 equations:
	ac:=ac+max(0,ceil((n^2-ne)/n));
	L:=subs(g_conversion2,lift_frlf0(f0,L,ac,`only r`));
	# Now transform to d/dx syntax:
	xDFpow:=[1];
	for i from 1 to n-1 do
		xDFpow:=[op(xDFpow),mult(x*DF,xDFpow[nops(xDFpow)])]
	od;
	values:=NULL;
	zeros:=[seq(0,i=1..n)];
	for i from 1 to n^2 do
		r:=collect(L[i],xDF);
		# convert to DF syntax as follows
		r:=expand(convert([seq(coeff(r,xDF,j-1)
		 *xDFpow[j],j=1..n)],`+`));
		# The accuracy of coeff(r,DF,j) is ac+j
		# Now divide by B[j+1][1]
		r:=convert([seq(convert(series(
		 coeff(r,DF,j)/B[j+1][1],x=0,ac+j),polynom)*DF^j,j=0..n-1)],`+`);
		# The high part should be cancel, the low part contributes to
		# the value of the solution
		zeros:=[seq(zeros[j]+a[i]*quo(coeff(r,DF,j-1)*x^(-min(0,B[j][2]+1)),x^max(0,B[j][2]+1),x)
		,j=1..n)];
		values:=values,r
	od;
	values:=[values];
	userinfo(3,'diffop',`solving linear equations in`,n^2,`variables...`);
	j:={seq(a[i],i=1..n^2)};
	r:=readlib(`solve/linear`)({seq(coeffs(collect(i,x),x),i=zeros)},j);
	userinfo(3,'diffop',`done solving`);
	r:=collect(convert([seq(subs(r,a[i])*values[i],i=1..n^2)],`+`),DF);
	r:=convert([seq(g_normal(coeff(r,DF,i)*B[i+1][1])*DF^i,i=0..n-1)],`+`);
	i:=0;
	while zeros<>{0} do
		i:=i+1;
		zeros:=eigenring_lineq(op(subs(x=x-i,[f,r])));
		if not member(zeros,{`not a regular point`,{0}}) then
			userinfo(3,'diffop',`solving extra linear equations`);
			r:=subs(`solve/linear`(zeros,j),r)
		fi
	od;
	j:=indets(r) intersect j;
if member('verify'=true,[args]) then
lprint(`Found `,nops(j),`linearly independent solutions. Checking correctness...`);
	zeros:=[coeffs(collect(rightdivision(numer(mult(f,r)),f)[2],DF),DF)];
	zeros:={seq(coeffs(collect(numer(i),x),x),i=zeros)};
	r:=subs(`solve/linear`(zeros,j),r);
	j:=indets(r) intersect j;
lprint(`Verified correctness. Found `,nops(j),`independent solutions`)
fi;
	[seq(collect(subs(i=1,solve(j),r),DF,g_normal),i=j)]
end:

# Bounds for the mixed equation f*r+l*f=0
bounds_frlf0:=proc(f)
	local v,ext,n,singularities,p,i,E,r,vp,s,m,B,D,N;
	global g_conversion1,g_conversion2,x,xDF,DF;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if nargs=1 then
		v:=g_ext(f);
		RETURN(subs(g_conversion2,
		procname(subs(g_conversion1,f),v)))
	fi;
	ext:=args[2];
	n:=degree(f,DF);
	if n<2 then ERROR(`wrong arguments`) fi;
	singularities:=[seq(i[1],i=v_ext_m(g_factors(denom(
	 g_normal(f/lcoeff(f,DF))),x,ext),x)),infinity];
	userinfo(3,'diffop',`Computing the bounds...`);
	for i from 0 to n-1 do D[i]:=1; N[i]:=0 od;
for p in singularities do
	B[p,[]]:=0;
	for i from 0 to n-1 do B[p,i]:=infinity od;
	v:=factor_op(l_p(f,p),'semireg');
	for E in v do
		if not type(E,list) then bug() fi;
		r:=degree(E[3],x); # ramification index
		vp:=min(0,ldegree(E[4],x)/r); # the number v' in lemma 5
		s:=E[1]; # semi-regular part
		s:=[solve(nt(s,1),xDF)]; # list of exponents of s
		for i in s do if not type(i,integer) then bug() fi od;
		# Smallest difference in the canonical list that is an
		# integer divided by the ramification index:
		m:=(min(op(s))-max(op(s)))/r;
		for i from 0 to n-1 do # lemma 5, order(r) <= n-1
			B[p,i]:=min(B[p,i],ceil(m+(n-1-i)*vp))
		od;
		if p<>infinity and type(E[4],integer) and
		 min(seq(E[4]+i,i=s))=0 then
			# Can improve the bound a bit for this case:
			i:=0;
			while member(i-E[4],s) do i:=i+1 od;
			B[p,[]]:=i
		fi
	od;
	for i from 0 to n-1 do
	if p<>infinity then
		if g_ext([f,p])<>ext then
		 r:=subs(_Z=x,subs(g_conversion1,op(subs(g_conversion2,p))))
		else
		 r:=x-p
		fi;
		B[p,i]:=B[p,i]+i; # From the xDF -> DF conversion
		# If 0,1,..,B[p,[]]-1 are exponents and there are no negative
		# integer exponents then the bound can be tightened a bit,
		# using the following argument.
# Suppose r=r_(n-1)*DF^(n-1)+..+r_0*DF^0
# Suppose that at the point x=p (where p is finite) we have lower bounds
# B[i],i=0..n-1 for the valuation of r_i. Because of the xDF -> DF conversion
# we always have B[0]<B[1]<B[2]<...
# Now suppose that 0 is the smallest integer exponent. Take a solution y with
# degl=0 that has exponent 0. Then apply r, and you see that the valuation
# of r(y) is B[0] unless the coefficient of r_0 has valuation > B[0].
# However, the valuation of r(y) must be >=0 because 0 was the smallest integer
# exponent and r(y) must be a (possibly 0) solution of f. If now
# B[0]<0, then we can replace B[0] by B[0]+1. If then there is also an
# exponent 1, apply same argument, and can increase B[1] (but then also B[0]),
# if B[1]<0.
# Use the property: if y in C((x))[log(x)] with valuation v>=0 and degl=0,
# then the n'th derivative of y has valuation >=v-n+1 (for the +1 degl=0 is
# needed).
		if B[p,[]]=n then
			B[p,i]:=0
		elif i<B[p,[]] then
			B[p,i]:=max(B[p,i],min(0,B[p,B[p,[]]]+B[p,[]]))
		fi;
		D[i]:=D[i]/r^B[p,i];
		N[i]:=N[i]-degree(r,x)*B[p,i]
	else
		# the +i is from the xDF -> DF conversion (gives -i) followed
		# by the DF (at infinity) -> DF conversion (gives + 2*i)
		N[i]:=N[i]-B[p,i]+i
	fi
	od;
od;
	userinfo(3,'diffop',`Done computing the bounds`);
	[seq([1/D[i],N[i]],i=0..n-1)]
end:

# Compute a basis of the local solutions of the mixed equation f*r + l*f = 0
# at the point x=0, assuming that x=0 is regular.
# Input: f a local differential operator
local_frlf0:=proc(f)
	local res,f0,n,a,g,j,r,l;
	global x,xDF;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	res:=NULL;
	f0:=nt(f,1);
	n:=degree(f,xDF);
for a from 1-n to n-1 do
	g:=normal(subs(xDF=xDF+a,f0)/f0);
	g:=[expand(numer(g)),expand(denom(g))];
	for j from 0 to n-1-abs(a) do
		r:=x^a*xDF^j*g[2];
		l:=-x^a*xDF^j*g[1];
		res:=res,[expand(l),expand(r),a+1] # l,r, and the accuracy
	od
od;
	[res]
end:

# v is a list of solutions l,r,accuracy
# The output is this list up to accuracy a (or higher if the accuracy was
# already more than a)
lift_frlf0:=proc(f,v,a)
	local f0,n,res,ext,i,j;
	global xDF;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	f0:=nt(f,1);
	n:=degree(f,xDF);
	res:=NULL;
	ext:=g_ext(f);
	for i in v do
		j:=i;
		while j[3]<a do
			j:=liftstep_frlf0(f,f0,n,ext,op(j))
		od;
		if args[nargs]=`only r` then j:=j[2] fi;
		res:=res,j
	od;
	[res]
end:

# f0 is the Newton polynomial of f, so f0=xDF*(xDF-1)*..
# n is the order of f
# a is the accuracy of l and r
# This will be increased by 1
liftstep_frlf0:=proc(f,f0,n,ext,l,r,a)
	local F,c,d,l_extra,r_extra,frlf,i;
	global g_conversion1,g_conversion2,x,xDF;
	# Compute a sufficient number of terms of f:
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	F:=nt(f,a+n);
	# Compute f*r + l*f:
	frlf:=subs(x=1,nm_mult(F,r,a,a,ext)+nm_mult(l,F,a,a,ext));
	# Uncomment for debug mode:
	# if expand(nm_mult(F,r,a-1,a-1,ext)+nm_mult(l,F,a-1,a-1,ext))<>0
	# then ERROR()
	# fi;
	userinfo(10,'diffop',`Computing the`,a,`th term of a local endomorphism`);
	if a<n then
		# Add (x^a * the following) to l and r:
		l_extra:=convert([seq(c[i]*xDF^i,i=0..n-1)],`+`);
		# The basis element for r have the form x^j (j=0,1,..) times
		# something of degree abs(a), so modulo that basis we
		# go up to degree abs(a)-1.
		r_extra:=convert([seq(d[i]*xDF^i,i=0..abs(a)-1)],`+`);
		frlf:=collect(frlf+subs(xDF=xDF+a,f0)*r_extra +f0*l_extra,xDF);
		F:=modm(subs(readlib(`solve/linear`)(
		  {coeffs(subs(g_conversion2,frlf),xDF)}
		  ,indets([l_extra,r_extra]))
		 ,g_conversion1,[l_extra,r_extra]));
		if (indets(F) intersect indets([l_extra,r_extra]))
		   minus {xDF} <> {} then ERROR('failed') fi;
		[l+x^a*F[1],r+x^a*F[2],a+1]
	else
		gcdex(subs(xDF=xDF+a,f0),f0,-subs(g_conversion2,frlf),xDF
		 ,'r_extra','l_extra');
		subs(g_conversion1,[l+x^a*modm(l_extra),r+x^a*modm(r_extra),a+1])
	fi
end:

macro(
	solve=readlib(`solve/linear`)
):

# f: monic differential operator
# r: element of E_D(f,f), so r is a k linear map r: V(f) -> V(f)
# RootOf syntax, parameters allowed.
# Output: characteristic polynomial of the k-linear action of r on V(f)
endomorphism_charpoly:=proc(f,r)
	local a,i,S,s,Si,j,n;
	global DF,x;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if evala(Normal(subs(x=0,denom(f)*lcoeff(f,DF))))=0 then
	# x=0 is not a regular point
	RETURN(procname(seq(collect(subs(x=x+1,i),DF,evala),i=[f,r])))
fi;
	n:=degree(f,DF);
	# Solution of f:
	S:=convert([seq(a[i]*x^i,i=0..2*n-1)],`+`);
	# Apply f:
	s:=0;
	Si:=S;
	for i from 0 to n do
		s:=s+convert(taylor(coeff(f,DF,i),x=0,2*n),polynom)*Si;
		Si:=diff(Si,x)
	od;
	s:=collect(s,x);
	S:=subs(solve(evala({seq(coeff(s,x,i),i=0..n-1)})
	 ,{seq(a[i],i=n..2*n-1)}),S);
	# Apply r on the solution S:
	s:=0;
	Si:=S;
	for i from 0 to n do
		s:=s+convert(taylor(coeff(r,DF,i),x=0,n),polynom)*Si;
		Si:=diff(Si,x)
	od;
	s:=collect(s,x);
	S:={seq(a[i]=0,i=0..n-1)};
	collect(linalg[charpoly](linalg[matrix](n,n,[seq(seq(subs(S,
	 subs(a[j]=1,coeff(s,x,i))),i=0..n-1),j=0..n-1)]),x),x,evala)
end:

macro(
	solve=solve
):

# Input: operator L and point p, and (optional) a sequence of options
# If the third entry is a variable then the
# Option minimal -> give only the minimal generalized exponents.
gen_exp:=proc(f,p,T)
	global x,xDF,DF;
	local v,E,r,s,res,i;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if has([args],{integer,rational}) then
	# compute only the integer or rational exponents
	v:=Newtonpolygon(l_p(args),`include the Newton polynomials`);
	if v[1][3]<>0 then RETURN([]) fi;
	v:=factors(collect(v[1][4]/lcoeff(v[1][4],x),x,g_normal))[2];
	v:=[seq(`if`(type(x-i[1],rational),(x-i[1])$i[2],NULL),i=v)];
	if has([args],integer) then
		v:=[seq(`if`(type(i,integer),i,NULL),i=v)]
	fi;
	if has([args],'minimal') then
		# give only the minimal ones
		[[min(op(v)),T=x]]
	else
		[sort(v),T=x]
	fi
else
	# not just: l_p(f,p) because args might contain additional
	# algebraic numbers:
	v:=factor_op(l_p(args),'semireg',g_ext([args]));
	res:=NULL;
	for E in v do
		# if not type(E,list) then bug() fi;
		r:=degree(E[3],x); # ramification index
		s:=E[1]; # semi-regular part
		# list of exponents of s (these are integers)
		s:=sort(map(i -> `if`(has(i[1],xDF),solve(i[1],xDF)$i[2],NULL),
		 factors(nt(s,1))[2]));
		if has([args],'minimal') then
			# give only the minimal ones
			s:=[min(op(s))]
		fi;
		# now add E[4], then divide by r
		s:=[seq((i+E[4])/r,i=s)];
		if has([args],'ramification1') and r>1 then
			s:=NULL
		else
			s:=[op(subs(x=args[3],s)),subs(x=args[3],E[3])=x]
		fi;
		res:=res,s
	od;
	subs(g_conversion2,[res])
fi
end:

#savelib('gen_exp','`DEtools/diffop/gen_exp.m`'):
#savelib('endomorphism_charpoly','`DEtools/diffop/endomorphism_charpoly.m`'):
#savelib('eigenring_lineq','eigenring','bounds_frlf0','local_frlf0',\
  'lift_frlf0','liftstep_frlf0','`DEtools/diffop/eigenring.m`'):