# Help page of this implementation can be found at https://sites.google.com/site/yongjaecha/code/hom. 
##########################################################################################################
Hom := proc(L1, L2)
	local ind, G, r, i;
	if coeff(L1, _Env_LRE_tau,0)=0 or coeff(L2, _Env_LRE_tau, 0)=0 then 
		return "input should be a normal operator"
	fi;
	ind := indets({L1, L2}) minus {_Env_LRE_x, _Env_LRE_tau};
	G := adjustR(RatPolySym(adjointOperator(L1), L2, ind), L1);

	# Check if G is correct:
	r := numer(`LREtools/rightdivision`(`LREtools/mult`(L2, G), L1)[2]);
	if r<>0 then
		lprint("Discarding wrong solutions");  # Explanation, see [**] below.
		G:=subs(CheckZero(r, ind), G)
	fi;
	ind := indets(G) minus {op(ind), _Env_LRE_x, _Env_LRE_tau};
	[seq(collect(coeff(G, i), _Env_LRE_tau, normal), i=ind )]
end:


CheckZero:=proc(r, ind)
	local tosolve, ind1, i, j;
	tosolve := seq(coeff(r, _Env_LRE_tau, i), i=ldegree(r, _Env_LRE_tau)..degree(r, _Env_LRE_tau));
	tosolve := seq( seq( coeff( i,_Env_LRE_x, j), j=ldegree(i,_Env_LRE_x)..degree(i,_Env_LRE_x) ), `in`(i, {tosolve}) );
	ind1 := indets(r) minus {_Env_LRE_tau,_Env_LRE_x} minus ind;
	solve( {tosolve}, ind1)
end:

##########################################################################################################
# change of basis from v=d[n-1] to d[i]'s
adjustR := proc(mat,L)
	local action, V, b, i, n, bd, v;
	n := degree(L, _Env_LRE_tau);
	action := subsDual(L, n, bd);
	V[0] := bd[n-1];
	for i to n-1 do 
		V[i] := normal(subs(action, subs(_Env_LRE_x = _Env_LRE_x+1, V[i-1]))) 
	od;
	V := SolveTools:-Linear({seq(V[i]-v[i], i = 0 .. n-1)}, {seq(bd[i], i = 0 .. n-1)});
	b := subs(V, add(bd[i]*_Env_LRE_tau^i, i = 0 .. n-1));
	collect(subs({seq( v[i-1]=mat[i], i=1..n )}, b), _Env_LRE_tau, normal)
end:
adjointOperator := proc(L) local n, i;
	n := degree(L, _Env_LRE_tau); 
	add(normal(subs(_Env_LRE_x = _Env_LRE_x+i-1, coeff(L, _Env_LRE_tau, n-i)/lcoeff(L,_Env_LRE_tau))) * _Env_LRE_tau^i, i=0..n)
end:
subsDual := proc(L, n, bd) 
	local j;
	{seq(bd[j] = -coeff(L, _Env_LRE_tau, j+1)*bd[0]/coeff(L, _Env_LRE_tau, 0)+`if`(j = n-1, 0, bd[j+1]), j = 0 .. n-1)} 
end:

######################################################################################################
# Find ratsol of symmetric product of L1 and L2;
RatPolySym:=proc(L1, L2, ind)
	local Ll, Lh, sw;
	if degree(L1, _Env_LRE_tau) > degree(L2, _Env_LRE_tau) then
		Ll, Lh := L2, L1;
		# Switched the two operators so that the first operator has the lowest order.
		sw := true
	else
		Ll, Lh := L1, L2;
		sw := false
	fi;	
	FindCoeff(RatCdd(Ll, Lh), Ll, Lh, ind, sw);
end:


# To make FindCoeff faster, need to limit number of eqns in tosolve  [**]
FindCoeff:=proc(ll, Ll, Lh, ind, sw)
	local W, dh, dl, cff, i, j, relbackLl, tosolve, solved, sol;
	dh := degree(Lh, _Env_LRE_tau);
	dl := degree(Ll, _Env_LRE_tau);
	cff := indets(ll) minus {op(ind), _Env_LRE_x};
	sol := ll;
	for i from 0 to dl-1 do
		W[0, -i] := subs(_Env_LRE_x = _Env_LRE_x-i, ll[i+1])
	od;

	# rel back with Ll to W[-1,_]
	relbackLl := collect(subs(_Env_LRE_x = _Env_LRE_x-1, -Ll/tcoeff(Ll, _Env_LRE_tau) + 1 ), _Env_LRE_tau, normal);
	
	for j from 0 to dh while nops(cff) > 1 do
	    W[0,1+j] := normal(subs(_Env_LRE_x=_Env_LRE_x+1, add(coeff(relbackLl, _Env_LRE_tau, dl-i)*subs(_Env_LRE_x=_Env_LRE_x+dl-1-i, W[0,j-dl+1+i]), i=0..dl-1)));
	    if dl+j >= dh then 
		tosolve := numer(add( subs(_Env_LRE_x=_Env_LRE_x+j-dh+1,coeff( Lh, _Env_LRE_tau, i))*W[0,1+j-dh+i], i=0..dh));
		tosolve := seq( coeff( tosolve, _Env_LRE_x, i), i=0..degree(tosolve,_Env_LRE_x));
		solved := SolveTools:-Linear( {tosolve}, cff );

		for i from -dl+1 to 1+j do
			W[0, i] := subs(solved, W[0,i])
		od;
		sol := subs(solved, sol);
		cff := {seq( op(indets(W[0,i])), i=-dl+1..1+j)} minus {op(ind), _Env_LRE_x}

	    fi
	od;
	if sw then
		collect([seq(W[0,i], i=0..dh-1)], cff)
	else
		collect(sol, cff)
	fi
end:


# Cdd of Rational solution of L1 tensor L2
RatCdd := proc(Ll, Lh)
	local gxp2, Gp1, Gp2, Denomi, DenomDeg, Numbdd, i, dh, j, c, k;
	dh := degree(Lh, _Env_LRE_tau);
	gxp2 := genexp(Lh);
	Gp1 := plist(Ll);
	Gp2 := plist(Lh);
	Denomi := ListOfDenomBound(Gp1, Gp2, dh);
	DenomDeg := [seq(degree( numer( Denomi[i]), _Env_LRE_x)- degree( denom( Denomi[i]), _Env_LRE_x), i=1..dh )];
	Numbdd := [seq(max(seq(CompareGenExp(i, subs(_Env_LRE_x=_Env_LRE_x+j, Ll), T), i=gxp2)), j=0..dh-1)];
	[seq(Denomi[j]*add(c[j-1,k]*_Env_LRE_x^k, k=0..Numbdd[j]-DenomDeg[j]), j=1..dh)]
end:


# Denom bound for tau^i(u)*v, u \in V(L1), v \in V(L2), i.e W[0,0], .., W[i, 0]
ListOfDenomBound:=proc(l1, l2, n)
	local i;
	[seq(DenomBound(Shiftl1(l1, i), l2 ), i=0..n-1)];
end:
Shiftl1:=proc(l2,k)
	local i,j;
	{seq([seq( i[j], j=1..7), i[-2]-k, i[-1]-k], i=l2)}
end:


# compare a genexp of L1 and the set of genexp of L2 and 
# return when there product is in the form of 1+c*T^d
CompareGenExp:=proc(g, L, T)
	local rg, i, cdd, Ls, Res,eq;
	rg:=degree(rhs(g[2]), T);
	Res:=NULL;	
	Ls:=subs(t=t^(rg), xtot(L));
	eq:=Ind(symprod(Ls, subs(T=t, g[1]), rg), t, rg);
	cdd:=solve(eq, _n);
	for i in cdd do if type(i*rg, integer) then
		Res := Res, -floor(i)
	fi od;
   Res;
end;

#denominator bound for ratsol of L1 tensor L2
#l1=plist(L1), l2=plist(L2),
DenomBound:=proc( l1, l2)
	local  GGp,i,j, Gp, ll2, denomi, maxlist, dem, leftbound, rightbound, st; 

	#group into same sing.pt. 
	GGp:=NULL;
	ll2:=l2;
	for i in l1 do
	Gp[i[1]]:= i;
	 for j in ll2 do
		if j[1]=i[1] then
		Gp[i[1]]:=Gp[i[1]], j;
		ll2:= ll2 minus {j};		
		fi;
	 od; 
	GGp:=GGp, [Gp[i[1]]];
	
	od; 
	for i in ll2 do
	Gp[i[1]]:= i;
	GGp:=GGp, [Gp[i[1]]];
	od;

	GGp:=[GGp];

	#generate denominator bound	
	denomi:=1;
	for i in GGp do
	if nops(i)=1 then
		maxlist:= [seq(max(i[1,7, j], i[1,6, j]), j = 1 .. nops(i[1,6]))];
		st:=i[1,-2]-i[1,2]
	else
		leftbound:=MatchBound(i[1,-4], i[1,-2]-i[1,2], i[2,-4], i[2, -2]-i[2,2],
		i[1,-4,1], i[1,-4,1]+min(i[1,3]), i[2,-4,1], i[2,-4,1]+min(i[2,3]) );

		rightbound:=MatchBound(i[1,-3], i[1,-2]-i[1,2], i[2,-3], i[2, -2]-i[2,2],
		i[1,-3,-1]-max(i[1,3]), i[1,-3,-1], i[2,-3,-1]-max(i[2,3]), i[2,-3,-1] );

		maxlist:= [seq(max( leftbound[1,j], rightbound[1,j] ), j = 1 .. nops(leftbound[1]))];
		st:=leftbound[2];
	fi;
	
	dem:=mul(evala(Norm(_Env_LRE_x-(i[1,1]+st+j-1)))^maxlist[j], j = 1 .. nops(maxlist));
	denomi:=denomi*dem; 
	od;

denomi
end:

# match two bounds
# l1, l2:= bounds, pl1, pl2:= left starting point
# returns matched bounds of two bounds and starting point, st, of the bounds
MatchBound:=proc( l1, pl1, l2, pl2, first_l1, last_l1, first_l2, last_l2)
	local mn, mx, r1, r2, ll1, ll2, i;
	mn:=min(pl1, pl2);
	r1:=pl1+nops(l1);
	r2:=pl2+nops(l2);
	mx:=max(r1, r2);
	ll1:=[ first_l1$(pl1-mn), op(l1), last_l1$(mx-r1)];
	ll2:=[ first_l2$(pl2-mn), op(l2), last_l2$(mx-r2)];
	[seq( ll1[i]+ll2[i], i=1..nops(ll1)  )], mn
end:



######################################################################################################
xtot:=proc(L)
	global t;
	collect(primpart(subs(_Env_LRE_x=1/t,L),_Env_LRE_tau),_Env_LRE_tau);
end;

sing := proc(L) 
	local d, i, M, ext; 
	d := degree(L, _Env_LRE_tau); 
	if ldegree(L, _Env_LRE_tau) <> 0 then 
	error "not yet implemented" end if; 
	for i from 5 while degree(M, _Env_LRE_tau) <> d+1 do 
	M := evala(Primpart(`LREtools/LCLM`(L, _Env_LRE_tau-i), _Env_LRE_tau)) end do;

	ext := indets([args], {'RootOf', 'radical', 'nonreal'});  
	M := {gcd(lcoeff(L, _Env_LRE_tau), eval(lcoeff(M, _Env_LRE_tau), _Env_LRE_x = _Env_LRE_x - 1)), gcd(tcoeff(L, _Env_LRE_tau), tcoeff(M, _Env_LRE_tau))};
	M := {seq( op(`LREtools/g_p/factors`(i, _Env_LRE_x, ext)) , i=M)};
	M := {seq(i[1], i = M)}; 	
	M := {seq(`LREtools/g_p/shorten_int`(i[1], _Env_LRE_x), i = `LREtools/g_p/int_dif`(M, _Env_LRE_x))}; 
	PolynomialTools:-Sort([op(M)], _Env_LRE_x)

end proc;


#[ sing.pt, degree of L, gp]
plist:=proc(L)
	local singg,res,s,minmax, Lp, n;
	Lp:=primpart(L, _Env_LRE_tau);
	n:=degree(Lp, _Env_LRE_tau);
	singg:=sing(Lp);
	res := {};
	for s in singg do
		minmax:=`LREtools/g_p`(Lp, RootOf(s,_Env_LRE_x));
		res := res union { [RootOf(s,_Env_LRE_x), n, minmax] }
	od;
	res
end;



##########################################################################################################
#gen exp

val:=proc(Ld, t, r, delta)
	local  S, i,s ,n;
	n:=degree(Ld,delta);
	s:=min(seq(i+ldegree(coeff(Ld,delta,i),t)/r ,i=0..n));
	S:={};
	for i from 0 to n do
		if i+ldegree(coeff(Ld,delta,i),t)/r=s then
		S:= S union {i}
		fi;
	od;
	[s,S]
end;

Ind:=proc(L, t, r)
	local Ldd,v,M,P,i,delta,k;
	Ldd:= collect(subs(_Env_LRE_tau = delta+1, L), delta, evala);
	v:=val(Ldd,t, r, delta);
	M:=v[2];
	P:=0;
	for i in M do
	P:=P+(-1)^i*mul(_n+k,k=0..i-1)*simplify(tcoeff(coeff(Ldd,delta,i),t))
	od;
	collect(P,_n)
end;



Taut := proc(t, r, j, N)
	convert(series(t/(1+j*t^r)^(1/r),t=0,N), polynom);
end:

#sym prod of L and tau-M in terms of t
symprod:=proc(L,M,r)
	global t;
	local n, P, i,j;
	n:=degree(L,_Env_LRE_tau);
	P:=coeff(L,_Env_LRE_tau,n)*_Env_LRE_tau^n+add(simplify(coeff(L, _Env_LRE_tau, i)*mul(subs(t = Taut(t,r,j,r*n+2), M), j = i .. n-1))*_Env_LRE_tau^i, i = 0 .. n-1);
 	P:=primpart(P,_Env_LRE_tau);
	P:=collect(P, _Env_LRE_tau, simplify)
	end proc;


# L in C((t))[tau]
taupoly:=proc(L)
	local n,r,m,P,i,v,s,F, P1;
	n:=degree(L,_Env_LRE_tau);
	if n=0 or L=0 then return {} fi;
	for i from 0 to n do  
		v[i] := ldegree(coeff(L, _Env_LRE_tau, i), t) 
	od;

	# s = the maximal slope
	s := max(seq((v[n]-v[i])/(n-i), i = 0 .. n-1));

	# r = {points on the maximal slope}
	r := {n, seq(`if`( (v[n]-v[i]) / (n-i) = s, i, NULL), i=0..n-1)};

	# m = the first point on the maximal slope
	m:=min(op(r));

	# Newton polynomial
	F:=[];
	P1:=add(tcoeff(coeff(L,_Env_LRE_tau,i), t)*t^(i-m), i=r);
	P:=factors(P1)[2];
	for i in P do
	if has(i[1],t) then F:=[op(F),i[1]] fi od;
	P:=[ seq([P1, RootOf(i,t)], `in`(i, F)) ];

	# We return the set {[slope, NewtonPolynomial, c, ramification]}
	{seq([-s, factor(i[1]),i[2], denom(s)], `in`(i, P))} union procname(add(coeff(L, _Env_LRE_tau, i)*_Env_LRE_tau^i, i = 0 .. m))
end:

genexp:=proc(Lx)
	local Tp, i, Gx, k, D, L ;
	L:=xtot(Lx):
	Tp:=taupoly(L);
	Tp:=[seq([i[3]*t^(i[1]*i[4]), i[4]], `in`(i, Tp))];
	Gx:=[];
	
	for i in Tp do

	D:= gx(symprod(subs(t=t^i[2],L), 1/i[1], i[2]), [[1, 0, 0]],i[2] );
	Gx:=[op(Gx),  seq([subs(t=t^(k[2]/i[2]),i[1])*k[1], k[2]] ,`in`(k, D)  ) ];

	od;
	
	[op({seq([subs(t = T, i[1]), t = T^i[2]], `in`(i, Gx))})]
end:

# A :=list of [so far found, degree of so far found, c::check value]
gx:=proc(L, A,r)
	
	local Ld, R, Dp, P, j, k, i, rr;  
	R:=[];
	for i in A do 

		# skip the ones which are already done
		if i[3]=1 then 
		R:=[op(R), i];

		elif i[3]=0 then
		# this is to avoid doing symp twice.
		if i[1]=1 then
		Ld:=L;

		else
		Ld:=symprod(L,1/i[1],r);
		fi;
		Dp:=deltapoly(Ld, i[2], r);
		for k in Dp do  
			# when we reach the final part.
			#if k[2]=r then
			if degree(k[1], t)=k[2]*r or k[1]=0 then      
			P:=factors(Ind(Ld,t,r))[2];     
				for j in P do     
					if has(j[1],_n) then     
					R:= [op(R), [i[1]-RootOf(j[1],_n)*t^r, r , 1]]    
					fi;
				od;          
	
			# when need to continue. 
			else R:=[op( R ), [subs(t=t^k[2],i[1])+k[1],degree(k[1],t),0] ];
				if k[2] <>1 then				
				rr:=r*k[2];
				Ld:=subs(t=t^k[2], L);
				
				else
				rr:=r;
				Ld:=L;
				fi;	
	
			fi; 
		od; 

		fi;

	od;
	if mul(i[3], `in`(i, R))=1 then return R fi;

	procname(Ld, R, rr)
end:


#find delta poly with slope between r1 and r2
deltapoly:=proc(L, r1, r2)
	local Ld, R, Tp, i;
	R:={};
	Ld:=subs(_Env_LRE_tau=_Env_LRE_tau+1,L):
	Tp:=taupoly(Ld);
	for i in Tp do
	if i[1]<=r2 and r1<i[1] then 
	R:= R union {i};
	fi; od;

	if R={} then R:=[[0, r2]]; 
	#return [the term found by deltapoly with ramification modified, new ramification from this deltapoly]
	else R:=[seq([i[3]*t^(i[1]*i[4]), i[4]], `in`(i, R))];
	fi;

	return R
end;



xtot:=proc(L)
	collect(primpart(subs(_Env_LRE_x=1/t,L),_Env_LRE_tau),_Env_LRE_tau);
end;