# Author:   Mark van Hoeij,  July 2004.
#
# As an example, type this:
#
#    read Hom:
#    l1,l2 := Dx^2-Dx-1/x^2, Dx^2+x:
#    m1,m2 := seq(rightdivision(LCLM(i,Dx-x),Dx-x)[1],i=[l1,l2]):
#    L1 := LCLM(l1, m2);  L2 := LCLM(l2, m1);
#    v := Hom(L1,L2);
#    seq(rightdivision(mult(L2,i),L1)[2],i=v); # To check correctness.
#
# Make sure you first put this file, with the filename Hom, in a folder
# where Maple's read command can find it.



# Load some code:
with(DEtools):


_Envdiffopdomain := [Dx,x]:
# This tells Maple we use the notation Dx to represent the operator d/dx.


macro(
	APPARANT1 = [infinity,1],
	APPARANT2 = [infinity,2],
	degree_ext = `algcurves/degree_ext`,
	subsDual = `Hom/subsDual`,
	CompareGenExp = `Hom/CompareGenExp`,
	CharProperties = `Hom/CharProperties`,
	GenExpMinPol = `Hom/GenExpMinPol`,
	IntShift = `Hom/IntShift`,
	AppCheck = `Hom/AppCheck`,
	adjustR = `Hom/adjustR`,
	adjustS = `Hom/adjustS`,
	newvar = `Hom/newvar`,
	get_eqns = `Hom/get_eqns`,
	RatFuncBounds = `Hom/RatFuncBounds`
):

# Let D = C(x)[Dx] and let M be the D-module D/DL with
# basis b[0] .. b[n-1] where b[i] = image of Dx^i in M.
# Now denote the corresponding basis of the dual module
# by bd[0] .. bd[n-1].
# This procedure describes the action of Dx on this basis.
subsDual := proc(L, bd, Dx)
	local j,n;
	if lcoeff(L,Dx)<>1 then
	  return procname( collect(L/lcoeff(L,Dx),Dx,Normalizer), bd, Dx)
	fi;
	n := degree(L,Dx);
	{seq(bd[j] = coeff(L,Dx,j)*bd[n-1] - `if`(j=0,0,bd[j-1]),j=0..n-1)}
end:

# Given one gen.exp of L1, and the list of gen.exp of L2,
# we search for the gen.exp of L2 that has integer/r sum
# with the given gen.exp of L1 where r=ramification.
# The name of the local parameter is the third entry T
# and the base field is given in the last entry ext.
#
# If we find any, then we return the "integer/r sum",
# as well as the valuation (in {0,-1,-2,..}/r) of the gen.exp.
CompareGenExp := proc(g, v, T, ext)
	local r,c,S,i,lambda,j;
	c  := CharProperties(g, T, ext);
	r := c[1]; # ramification index.
	S := NULL;
	for i in v do if c = CharProperties(i, T, ext) then
		S := S, i
	fi od;
	if S=NULL then return S fi;
	c := GenExpMinPol(g, T, ext, lambda);
	c := collect((-1)^degree(c,lambda)*subs(lambda = -lambda,c),lambda);
	for i in [S] do
		j := IntShift(GenExpMinPol(i, T, ext, lambda), c, lambda, r);
		if j<>NULL then
			return [j,min(0,ldegree(g[1],T))/r]
		fi
	od;
	NULL
end:

# Return a list of some properties of a generalized exponent g.
CharProperties := proc(g, T, ext)
	local i;
	[degree(lhs(g[-1]),T), degree_ext(g, ext), degree_ext(g[-1],ext),
	seq([	`if`(coeff(g[1],T,i)=0,0,i),
		degree_ext(coeff(g[1],T,i), [ext,g[-1]])
	    ],i=min(0,ldegree(g[1],T))..0)]
end:

GenExpMinPol := proc(g, T, ext, lambda)
	local R;
	R := RootOf(numer(Normalizer(lhs(g[-1])-rhs(g[-1]))),T);
	evala(Norm(lambda - subs(T=R, g[1]), ext,ext))
end:

# Check if two polynomials are an integer/r shift of each other.
IntShift := proc(f,g,x,r)
	local n,s;
	n := degree(f,x);
	s := Normalizer( (coeff(g,x,n-1)-coeff(f,x,n-1))/n );
	if not type(r*s,integer) or Normalizer(subs(x=x+s,f)-g)<>0 then
		NULL
	else
		s
	fi
end:

# If output=true then p is guaranteed to be an apparant singularity of L.
# If output=false then we don't know if p is apparant or not.
AppCheck := proc(L,p,Dx,x)
	local R,LL;
	R := Dx - rand(10..100)();
	LL := DEtools['LCLM'](L,R);
	#
	# This LL is always monic, so the singularities of LL can be
	# found by looking at the denominator of LL. If the singularity p
	# of L is no longer a singularity of LL, then p is an apparant
	# singularity of L. In this case p will not contribute 
	# denominators in our matrix mat[i,j] in Hom.
	# This is a significant improvement in denominator bound
	# because (the minpoly of) an apparant singularity often has
	# high degree so if we can discard such p from the denominator
	# bound then that reduces the number of unknowns quite a bit.
	#
	has(denom(Normalizer(denom(Normalizer(LL))/p)),x)
end:

# Find a basis of operators R in C(x)[Dx] that map solutions of L1
# to solutions of L2.
Hom := proc(L1, L2, domain)
	local ext,T,i,j,sing,L1a,p,v1,v2,b,n1,n2,vars,pr,d,mat,solved,R,E,k
	 ,Dx,x;
	if nargs>2 and type(domain,list(name)) and nops(domain)=2 then
		_Envdiffopdomain := domain;
		return procname(L1,L2)
	elif not type(_Envdiffopdomain,list(name)) then
		error "differential algebra not specified"
	fi;
	Dx, x := op(_Envdiffopdomain);
	# lprint("entry time:", time());
	ext := indets([L1,L2],{'RootOf','radical','nonreal'});
	Normalizer:=`if`(ext={},normal,evala);
	# Make input monic.
	if lcoeff(L1,Dx)<>1 or lcoeff(L2,Dx)<>1 then
	  return procname(seq(collect(i/lcoeff(i,Dx),Dx,Normalizer),i=[L1,L2]))
	fi;
	n1 := degree(L1,Dx);
	n2 := degree(L2,Dx);
	sing := PolynomialTools:-Sort(map(i -> i[1],
	 `DEtools/kovacicsols/factors`(denom(L1)*denom(L2), x, ext)), x);
	L1a := DEtools['adjoint'](L1);
	pr := 1;
	for i in [infinity,op(sing)] do
	# lprint(SING,i,time());
		if i=infinity then p:=i else p:=RootOf(i,x) fi;
		v1 := DEtools['gen_exp'](L1a,T,x=p,'groundfield'=ext);
		v2 := DEtools['gen_exp'](L2 ,T,x=p,'groundfield'=ext);
		b[i] := [seq(CompareGenExp(j,v2,T,ext),j=v1)];
		if b[i]=[] then return 0 fi;
		if p<>infinity then
		  # Special treatment of apparant singularities.
		  if nops({op(v2[1][1..-2])})=n2 and v2[1][1]=0
		  and AppCheck(L2,i,Dx,x) then
			b[i] := [op(b[i]),APPARANT2]
		  fi;
		  if (member(APPARANT2,b[i]) and L1=L2)
		  or (nops({op(v1[1][1..-2])})=n1 and type(v1[1][1],integer)
		  and AppCheck(L1,i,Dx,x)) then
			b[i] := [op(b[i]),APPARANT1];
			if v1[1][1]<>0 then
				pr := pr*i^v1[1][1]
			fi;
		  fi
		fi
	od;
	if pr<>1 then
		# lprint(time(),AP_PROD,1/pr);
		b[infinity] := [seq([i[1]+degree(numer(pr),x)-degree(denom(pr)
		 ,x),i[2]],i=b[infinity])];
		L1a := collect(mult(L1a,pr)/pr,Dx,Normalizer)
	fi;
	# Now we have enough info for the bounds and we have to form
	# the system. We form entries of the matrix mat[i,j]. This
	# matrix corresponds to an element of L1a tensor L2.
	# Now this is not exactly the usual "Hom system" but it is
	# isomorphic to it (at the end, I'll have to do a conversion
	# using the procedure adjustR). The reason for using this
	# L1a tensor L2 system rather than the usual Hom system is
	# that this makes it easier to figure out the bounds. The
	# generalized exponents for L1a tensor L2 are simply the
	# sums of the generalized exponents of L1a and those of L2.
	# The procedure CompareGenExp figures out which of those
	# sums are integers (divided by the ramification index) and
	# returns information to let us find the smallest sum.
	#
	# The entries of mat[i,j] satisfy relations (which are
	# either used to compute the next entry, or are used to
	# get equations on the unknown coefficients). Also, one has
	# a denominator bound on each entry (if some entry fails this
	# denominator bound, then we compute linear equations from
	# that too).
	vars := {};
	solved := {};
	for i from 0 to n1 do for j from 0 to n2 do
		R,d := RatFuncBounds(sing,b,i,j,x);
		# lprint(i,j,d,time());
		if j<n2 and i=0 then
			p := add(newvar()*x^i,i=0..d);
			vars := vars union (indets(p) minus {x});
			# lprint("new vars:",d+1, "total:", nops(vars));
		else
			if j<n2 then
				p := diff(mat[i-1,j],x)-mat[i-1,j+1]
			else
				p := -add(coeff(L2,Dx,k)*mat[i,k],k=0..n2-1);
			fi;
			p := Normalizer(subs(solved,p)/R);
			E,p := get_eqns(p, x, vars, d);
			solved := adjustS(solved,E)
		fi;
		mat[i,j] := R*p;
		if i=n1 then # Get additional equations
			p := add(coeff(L1a,Dx,k)*mat[k,j],k=0..n1);
			p := numer(Normalizer(subs(solved,p)));
			E,p := get_eqns(p, x, vars, -1); # Equate p to 0.
			solved := adjustS(solved,E)
		fi;
	od od;
	# lprint(DONE_SOLVING, time());
	mat := subs(solved,
	 [seq(mat[i,0],i=0..n1-1)]);
	vars := indets(mat) intersect vars;
	if pr=1 then T := L1 else
		T := mult(L1,1/pr)
	fi;
	# adjust the syntax to the syntax one would have gotten with
	# the usual Hom system. Then sort the output so that the
	# smallest ones appear first.
	PolynomialTools:-Sort( adjustR(mat,n1,T,vars, pr, Dx,x) ,Dx)
end:

# Convert from my Hom syntax to the usual one.
adjustR := proc(mat, n1, L1, vars, pr, Dx,x)
	local b,bd,V,v,i,action,j;
	action := subsDual(L1,bd,Dx);
	V[0] := bd[n1-1];
	for i to n1-1 do
	  V[i] := Normalizer(diff(V[i-1],x)+subs(action, V[i-1]))
	od;
	b := seq(bd[i],i=0..n1-1);
	V := SolveTools:-Linear({seq(V[i]-v[i],i=0..n1-1)}, {b});
	b := subs(V,add(bd[i]*Dx^i,i=0..n1-1));
	V := solve(vars);
	[seq(collect( mult(subs({seq(v[i-1]=Normalizer(subs(j=1,V,mat[i])),
	 i=1..n1)},b),pr), Dx, Normalizer),j=vars)]
end:

# Update the set of solved equations.
adjustS := proc(sol,E)
	if E={} then return sol fi;
	Normalizer(subs(E,sol)) union E
end:

newvar := proc()
	local v;
	v;
end:

# Input p is a rational function in x, which we want to turn
# into a polynomial of degree <= d.
# If p is not a polynomial and has a denominator, then
# this procedure computes linear conditions on the coefficients
# in p equivalent to denom(p) cancelling. If p is a polynomial
# of degree >d, then again produce linear equations. Then solve
# the equations, return the solutions of those equations as well
# as the updated p, which now can only be a polynomial of degree <= d.
get_eqns := proc(p,x,vars,d)
	local e,a,i;
	e := NULL;
	if has(denom(p),x) then
		e := coeffs(collect(rem(numer(p),denom(p),x),x),x);
		# lprint(HAS_DENOM, "number of eqns:", nops([e]));
	fi;
	a := collect(numer(p),x);
	e := {e,seq(coeff(a,x,i),i=degree(denom(p),x)+d+1 .. degree(a,x))};
	# i:=nops(e minus {0});
	# if i>0 then lprint("total number of eqns:", i); fi;
	e := SolveTools:-Linear(e,indets(e) intersect vars);
	e, collect(Normalizer(subs(e,p)),x)
end:

# Rational function bounds for mat[i,j]. If i=j=0 then we add
# the generalized exponents and take the minimum that is an
# integer divided by the ramification (this is already done in
# procedure CompareGenExp) and then round above to an integer
# with the procedure ceil. If i or j is not 0, then we have
# differentiated i+j times, which lowers the valuation bound
# except when p is apparant (in which case we do not lower
# the bound).
RatFuncBounds := proc(sing, b, i, j, x)
	local p,d,v,pr,S,s;
	for p in [infinity,op(sing)] do
		s := `if`(member(APPARANT1,b[p]),0,i)
		   + `if`(member(APPARANT2,b[p]),0,j);
		d := min(seq(ceil( v[1] + s*(v[2]-1) ), v=b[p]));
		if p=infinity then
			pr := 1;
			S := d+2*(i+j)
		else
			pr := pr*p^d;
			S := S + d*degree(p,x)
		fi
	od;
	pr, -S
end: