#	Commented out gegenbauer (currently line 48) in solver until I get it right
# 	Changed det to Jdet and findgauge to findgauge2
# 	and changed x's to _n's; tau's, t's, and T's to _tau's, t's, and T's (and protected _t, _T in my code)
# 	and made 'a' local and got rid of _Env...
_Env_LRE_x:=_n;
_Env_LRE_tau:=_tau;

macro(lclm = `LREtools/LCLM`, mult = `LREtools/mult`);


#########################################################################################
solver:=proc(LL)
	local gp,L, gexp,gexpq, r, l,gpm,i, A;	
	L := evala(Primpart(LL,_tau));
	gp:=gpmaxmin(L);
	gexp:=genexp(xtot(L),_tau,_t);
	gexpq:=genxq(gexp);
	r:=degree(gexpq[1], _T);
	l:=ldegree(gexpq[1], _T);
	gpm:={};
	if gp<>{} then 
	for i in gp do 
	if i[2]<>0 then 
	gpm:=gpm union {i}; 
	fi od fi;

	#Bessel functions
	if r=3 and l=2 and gpm={} then
		bessel(LL,gexpq[1]);

	#WhittakerM, Laguerre
	elif r=2 and l=0 and tcoeff(gexpq[1],_T)=1 then 
		A:=Laguerre(LL, gpm, gexpq[1]);
		if A=FAIL then
			WhitM(LL,gexpq[1],gpm) 
		else return A
		fi;


	#WhittakerW 
	elif r=1 and l=0 and tcoeff(gexpq[1],_t)=-3+2*sqrt(2) and nops(gpm)=1 or nops(gpm)=2 then
		WhitW(LL,gexpq[1],gpm);

	#Legendre, Gegenbauer
	elif r=0 and l=0 and nops(gpm)=1 then
		A:=Legendre(LL, gexpq[1]);
		if A=FAIL then
#			gegenbauer(LL, gpm, gexpq[1])
		else return A
		fi;	
	
	#Jacobi	
	#elif r=1 and l=0 then
	#	jacobi(LL, gpm, gexpq[1])

	else return FAIL
	fi;
end;
##########################################################################################
jacobi:=proc(LL, gp, g1)
	local Ljc, a, cdd, gp1, i, b, z1, sol, tsol;
	Ljc:=_tau^2-(1/2)*(2*_n+3+alpha+beta)*(alpha^2-beta^2+(2*_n+alpha+beta)*(2*_n+4+alpha+beta)*z)*_tau/((_n+2)*(_n+2+alpha+beta)*(2*_n+2+alpha+beta))-(_n+1+alpha)*(_n+alpha+beta)*(2*_n+4+alpha+beta)/((_n+2)*(_n+2+alpha+beta)*(2*_n+2+alpha+beta));
	
	if gp={[0,4]} then 
		cdd:=[[subs({alpha=0, beta=0}, Ljc),0,0,z]]
	elif gp={[0,3], [1/2, 1]} then
		cdd:=[[subs({alpha=3/2, beta=3/2}, Ljc),3/2,3/2,z]]
	elif gp={[0,1], [1/2, 1]} then
		cdd:=[[subs({alpha=1/2, beta=1/2}, Ljc),1/2, 1/2,z]]
	elif gp={[0,6]} then
		cdd:=[[subs({alpha=0, beta=2}, Ljc),0,2,z], [subs({alpha=1, beta=3}, Ljc),1,3,z]]
	elif gp={[0,4], [1/2, 2]} then
		cdd:=[[subs({alpha=0, beta=1}, Ljc),0,1,z], [subs({alpha=1, beta=0}, Ljc),1,0,z]]
	elif {gp[1,2], gp[2,2]}={2,3} then
		
		for i in gp do 
	if i[2]=3 then a:=-i[1] fi od;
		
		cdd:=[[subs({alpha=a, beta=1}, Ljc),a,1,z], [subs({alpha=a, beta=0}, Ljc),a,0,z]]
	elif {gp[1,2], gp[2,2]}={2,1} then
		
		for i in gp do 
		if i[2]=1 then a:=-i[1] fi od;
		
		cdd:=[[subs({alpha=a, beta=a}, Ljc),a,a,z]]
	elif nops(gp)=3 then 
		
		for i in gp do 
		if i[2]=1 then a:=-i[1]; gp1:= `minus`(gp, {i})[1,1]  fi od;
		b:=[ -gp1-a, -a-2*gp1];

		cdd:=[seq([subs({alpha=a, beta=i}, Ljc),a,i,z], `in`(i, b))]
	else return FAIL
	fi;

	
	z1:=[ solve(-2*z*sqrt(z^2+1)-2*z^2-1, z), solve(2*z*sqrt(z^2+1)-2*z^2-1, z)];
	cdd:=[seq( op(subs(z=i, cdd)), `in`(i, z1))];

	for i in cdd do
	if symg(i[1], LL)<>FAIL then 
		sol:=symg(i[1], LL);
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[1]*a(_n), a(_n), {}, output = basis)[1];
	return [plgop(sol[2],JacobiP(_n, i[2], i[3], i[4])*tsol)]
	fi od;
	FAIL
end;
		




##########################################################################################
gegenbauer:=proc(LL, gp, g1)
	local c , mm, i , sol, tsol, a;
	c:=[solve(-2*z*sqrt(z^2+1)-2*z^2-1=g1, z), solve(2*z*sqrt(z^2+1)-2*z^2-1=g1, z)];
	
	if gp[1,2]=2 then mm:= 1/2
	elif gp[1,2]=1 then mm:=-gp[1,1]/2
	elif gp[1,2]=0 then mm:=0
	fi;

	c:=[seq( [subs({ z=i, m=mm}, _tau^2-2*z*(m+_n+1)*_tau/(_n+2)-(2*m+_n)/(_n+2)), i], `in`(i,c))];

	for i in c do 
	if symg(i[1], LL)<>FAIL then
		sol:=symg(i[1], LL);
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[1]*a(_n), a(_n), {}, output = basis)[1];
	return [plgop(sol[2],GegenbauerC(_n, mm, i[2])*tsol)]
	fi od;
	FAIL
end;

##########################################################################################
Legendre:=proc(LL, g1)
	local c, z, i, sol, tsol, a;
	c:=[solve(2*z^2-2*z*sqrt(z^2-1)-1=g1, z), solve(2*z^2+2*z*sqrt(z^2-1)-1=g1, z) ];
	c:=[seq([subs(z = i, _tau^2-(2*_n+3)*z*_tau/(_n+2)+(_n+1)/(_n+2)), i], `in`(i, c))];

	for i in c do
	if symg(i[1], LL)<>FAIL then
		sol:=symg(i[1], LL);
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[1]*a(_n), a(_n), {}, output = basis)[1];
	return [plgop(sol[2],LegendreP(_n, i[2])*tsol),plgop(sol[2],LegendreQ(_n, i[2])*tsol)]
	fi od;
	FAIL
end;
##########################################################################################
Laguerre:=proc(LL, gp, g1)
	local c, sol, tsol, aa, a;
	c:=coeff(g1, _T, 1);
	c:=-c^2/4;
	
	if nops(gp)=2 then
		aa:=-(gp[1,1]+gp[2,1])
	elif nops(gp)=1 then
		aa:=0
	else return FAIL
	fi;
	
	sol:=symg( subs({alpha=aa, z=c}, _tau^2-(2*_n+3+alpha-z)*_tau/(_n+2)+(_n+1+alpha)/(_n+2)), LL);
	
	if  sol <> FAIL then
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[1]*a(_n), a(_n), {}, output = basis)[1];
		return [plgop(sol[2],LaguerreL(_n, aa, c)*tsol)]
	fi;
	FAIL
end;
##########################################################################################
WhitW:=proc(LL,gexpq,gpm)
	local d,z, li, i , j, cdd, sol, tsol, a;
	d:=evala(coeff(gexpq,_t,1)/coeff(gexpq,_t,0));
	z:=tcoeff(sqrt(2)*d,sqrt(2));
	#li:=[v,n]
	if nops(gpm)=1 then 
		li:=[[1/2-gpm[1,1],0], [1-gpm[1,1],1/2]];
	else
		li:=[[-(gpm[1,1]+gpm[2,1])/2, -(gpm[1,1]-gpm[2,1]-1)/2],[-(gpm[1,1]+gpm[2,1]+1)/2, -(gpm[1,1]-gpm[2,1])/2]];
	fi;
	cdd:=[];
	for i in li do
		cdd:=[op(cdd),[_tau^2+(z-2*i[1]-2*_n-2)*_tau-i[1]-_n-1/4-i[1]^2-2*i[1]*_n-_n^2+i[1]^2]];
	od;

	for j in cdd do
	if symg(j[1],LL)<>FAIL then 
		sol:=symg(j[1],LL);
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[1]*a(_n), a(_n), {}, output = basis)[1];
		return [plgop(sol[2],WhittakerW(j[2]+_n,j[3],z)*tsol), 1]
	fi od;
	FAIL
end;
##########################################################################################
WhitM:=proc(LL,gexpq,gpm)
	local z, li, i , j, cdd, sol, tsol, a;
	z:=-coeff(gexpq,_t,1)^2/4;
	#li:=[v,n]
	if nops(gpm)=1 then 
		li:=[[1/2-gpm[1,1],0], [1-gpm[1,1],1/2]];
	else
		li:=[[-(gpm[1,1]+gpm[2,1])/2, -(gpm[1,1]-gpm[2,1]-1)/2],[-(gpm[1,1]+gpm[2,1]+1)/2, -(gpm[1,1]-gpm[2,1])/2]];
	fi;
	cdd:=[];
	for i in li do
		cdd:=[op(cdd),[_tau^2*(2*i[2]+2*i[1]+3+2*_n)+(2*z-4*i[1]-4*_n-4)*_tau-2*i[2]+1+2*i[1]+2*_n,i[1],i[2]]];
	od;

	for j in cdd do
	if symg(j[1],LL)<>FAIL then 
		sol:=symg(j[1],LL);
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[1]*a(_n), a(_n), {}, output = basis)[1];
		return [plgop(sol[2],WhittakerM(j[2]+_n,j[3],z)*tsol), 1]
	fi od;
	FAIL
end;
##########################################################################################
#bessel function
bessel:=proc(LL,gexpq)
	local c,cdd,sol,i,j,vl,si,cx,d,z,tsol,a,v;
	c:=tcoeff(gexpq,_T);
	d:=lcoeff(gexpq, _T)/tcoeff(gexpq, _T);
	vl:=[-1/2*d-1/2,-1/2*(d+1)-1/2];
	if c > 0 then 
		z:=2*sqrt(c);
		cdd:=[seq([z*_tau^2-(2+2*v+2*_n)*_tau+z,v,z], v in vl), seq([z*_tau^2+(2+2*v+2*_n)*_tau+z,v,z], v in vl)];
	else
		z:=2*sqrt(-c);
		cdd:=[seq([z*_tau^2-(2+2*v+2*_n)*_tau-z,v,z], v in vl), seq([z*_tau^2+(2+2*v+2*_n)*_tau-z,v,z], v in vl)];
	fi;
	for j in cdd do
	if symg(j[1], LL)<> FAIL then 
		sol:=[j,symg(j[1],LL)];
		si:=lcoeff(sol[1,1],_tau)*tcoeff(sol[1,1],_tau);
		cx:=coeff(coeff(sol[1,1],_tau,1),_n);
		tsol:=`LREtools/hypergeomsols`(a(1+_n)-sol[2,1]*a(_n), a(_n), {}, output = basis)[1];
		if si > 0 and cx > 0 then
		sol:= [plgop(sol[2,2],BesselJ(sol[1,2]+_n,-sol[1,3])*tsol),plgop(sol[2,2],BesselY(sol[1,2]+_n,-sol[1,3])*tsol)];
		elif si > 0 and cx < 0 then
		sol:= [plgop(sol[2,2],BesselJ(sol[1,2]+_n,sol[1,3])*tsol),plgop(sol[2,2],BesselY(sol[1,2]+_n,sol[1,3])*tsol)];
		elif si < 0 and cx < 0 then
		sol:= [plgop(sol[2,2],BesselI(sol[1,2]+_n,-sol[1,3])*tsol),plgop(sol[2,2],BesselK(sol[1,2]+_n,sol[1,3])*tsol)];
		elif si < 0 and cx > 0 then
		sol:= [plgop(sol[2,2],BesselI(sol[1,2]+_n,sol[1,3])*tsol),plgop(sol[2,2],BesselK(sol[1,2]+_n,-sol[1,3])*tsol)];
		fi;
		sol:=[seq(collect(i,indets(i,function)),i in sol)];
		return sol; 
	fi od;
	FAIL
end;


##########################################################################################
#basic op
symprodx:=proc(L,M)
	local n,r,i,j;
	n:=degree(L,_tau);
	r:=M;
	coeff(L,_tau,n)*_tau^n+add(simplify(coeff(L, _tau, i)*(mul(subs(_n = _n+j, r), j = i .. n-1)))*_tau^i, i = 0 .. n-1)
end proc;

xtot:=proc(L)
	global _n, _t, _tau;
	collect(primpart(subs(_n=1/_t,L),_tau),_tau);
end;

with(LinearAlgebra); 
applygauge := proc (L, G) 
local rel, rec, a, b, ShouldBeZero, s, d, Y, A,trec,trel;
trec:=tautou(L);
trel:=v(_n)=tautou(G); 
rel := eval(trel, {_c[1] = 1, _c[2] = 1}); 
rec := trec; 
s := v(_n+2)+a(_n)*v(_n+1)+b(_n)*v(_n); 
ShouldBeZero := eval(s, {rel, subs(_n = _n+1, rel), subs(_n = _n+2, rel)}); 
solve(rec, u(_n+2)); 
solve(eval(rec, _n = _n+1), u(_n+3)); 
subs(u(_n+3) = %, u(_n+2) = `%%`, ShouldBeZero); 
collect(%, {u(_n), u(_n+1)}, factor); 
ShouldBeZero := {coeffs(%, {u(_n), u(_n+1)})}; 
d := solve(ShouldBeZero, {a(_n), b(_n)}); 
subs(d, s); 
vtotau(%);
end proc;

tautou:=proc(L)
local i;
	add(coeff(L,_tau,i)*u(_n+i),i=0..degree(L,_tau))
end;



tautov:=proc(L)
local i;
	add(coeff(L,_tau,i)*v(_n+i),i=0..degree(L,_tau))
end;



vtotau:=proc(L)
local i;
	add(coeff(L,v(_n+i))*_tau^i, i=0..2)
end;



Jdet:=proc(L)
	factor(coeff(L,_tau,0)/coeff(L,_tau,2))
end;



plgop:= proc (L, y)
local i; 
simplify(add(coeff(L, _tau, i)*subs(_n = _n+i, y), i = 0 .. degree(L, _tau))) 
end proc;



findgauge2:=proc(L,M)
	local G,eq0,eq1,eq2,eq3,eq4,c0_x1,eq5,eq6,c0_x0,eqc1,c,d,K,cfs,e,L1,L2;
	L1:=primpart(L,_tau);
	L2:=primpart(M,_tau);
	if L1=L2 then return 1 fi;
	cfs := indets([L1,L2],function);
	cfs := cfs union subs(_n=_n +1, cfs);
	cfs := cfs union subs(_n =_n +2, cfs);
	G := c1(_n )*_tau + c0(_n );
	K:=`LREtools/rightdivision`( `LREtools/mult`(L2, G), L1)[2];
	eq0, eq1 := coeff(%,_tau,0), coeff(%,_tau,1);
	if coeff(L1,_tau,1)=0 and coeff(L2,_tau,1)<>0 then
		c:=solve(eq1, c0(_n +1));
		eval(eq0, {c0(_n ) = subs(_n  = _n -1, %), c0(_n +2) = subs(_n  = _n +1, %)});
		collect(%, {c1(_n ),c1(_n +1),c1(_n +2),c1(_n +3),c1(_n +4)}, factor);
		d:=LREtools[ratpolysols](%, c1(_n ), {}, output = basis);
		if whattype(d)=exprseq then return FAIL fi;
		d:=d[1];
		eval(c, {c1(_n ) = d, c1(_n +2) = subs(_n =_n +2,d), c1(_n +1) = (_n =_n +1,d)});
		return collect(simplify(d*_tau+subs(_n =_n -1,%)),_tau)
	elif coeff(L1,_tau,1)<>0 and coeff(L2,_tau,1)=0 then
		c:=solve(eq1, c0(_n +2));
		eval(eq0, {c0(_n +2) = %, c0(_n ) = eval(%, _n  = _n -2)});
		d:=LREtools[ratpolysols](%, c1(_n ), {}, output = basis);
		if whattype(d)=exprseq then return FAIL fi;
		d:=d[1];
		eval(c, {c1(_n ) = %, c1(_n +2) = subs(_n =_n +2,%), c1(_n +1) = (_n =_n +1,%)});
		return collect(simplify(d*_tau+subs(_n =_n -2,%)),_tau)
	elif coeff(L1,_tau,1)=0 and coeff(L2,_tau,1)=0 then
		d:=LREtools[ratpolysols](eq1, c1(_n ), {}, output = basis);
		e:=LREtools[ratpolysols](eq0, c0(_n ), {}, output = basis);
		if whattype(d)=exprseq then d:=0
		else d:=d[1]
		fi;
		if whattype(e)=exprseq then e:=0
		else e:=simplify(e[1])
		fi;
		if d*_tau+e=0 then return FAIL fi;
		return d*_tau+e
	else
		eq2 := simplify(coeff(eq1, c0(_n +2) ) * eq0 - coeff(eq0, c0(_n +2) ) * eq1); # eq2 has no c0(_n +2) term.
		eq3 := simplify(subs(_n =_n +1, eq2));
		eq4 := simplify(coeff(eq1, c0(_n +2) ) * eq3 - coeff(eq3, c0(_n +2) )* eq1); # Has only c0(_n +1) term
		if has(indets(eq4),c0(_n +1)) then 		
			c0_x1 := solve( expand(eq4), c0(_n +1) );
			eq5 := simplify(coeff(eq2, c0(_n +1) ) * eq4 - coeff(eq4, c0(_n +1) ) * eq2);
			eq6 := simplify(coeff(eq0, c0(_n +2) ) * eq3 - coeff(eq3, c0(_n +2) ) * eq0);
			if not has(eq2, c0(_n )) then
			c0_x0 := solve(expand(eq6), c0(_n ))
			else
			c0_x0 := solve(expand(eq5), c0(_n ))
			fi;
			eqc1 := subs(_n =_n +1, c0_x0) - c0_x1;
		else eqc1:=eq4
		fi;
		d:=LREtools[ratpolysols](%, c1(_n ), {}, output = basis);
		if whattype(d)=exprseq then return FAIL
		else d:=d[1];
		#elif nops(d)>2 then
		#	error "input must be reducible"
		fi;
		c := subs({c1(_n )=d, c1(_n +1)=subs(_n =_n +1,d), c1(_n +2)=subs(_n =_n +2,d),c1(_n +3)=subs(_n =_n +3,d)}, c0_x0);
		c:=d*_tau + normal(c);
		collect(simplify(c),_tau)
	fi;
end;
##########################################################################################
#symg
#find symprod and gauge transformation that sends L1 to L2
symg:=proc(L1, L2)
	local A, LsA, cdd, i, d;
	d:=Jdet(L2)/Jdet(L1);
	A:=sheq(d);
	cdd:=[[symprodx(L1,A),A],[symprodx(L1,-A),-A]];
	for i in cdd do
	if findgauge2(i[1],L2)<> FAIL then
		return [i[2], findgauge2(i[1],L2)]
	fi od;
	FAIL
end;
	
vutau:=proc(L)
local i;
	subs({c[1]=1},rhs(L));
	add(coeff(%,u(_n+i))*_tau^i,i=0..1);
end;

#group shift eq factors together
pt:=proc(l)
	local d,e,lp,i ,c,s,re;
	if {op(l)}={} then return {} fi;
	s:={l[1]};
	c:=l[1,2];
	for i from 2 to nops(l) do
	if degree(l[1,1],_n)=degree(l[i,1],_n) then
		d:=degree(l[1,1],_n)-1;
		e:=Normalizer(coeff(l[1,1],_n,d)-coeff(l[i,1],_n,d));
		if type(e,integer) and e<>0 and irem(e,d+1)=0 and
		 Normalizer(l[1,1]-subs(_n=_n+e/(d+1),l[i,1]))=0 then 
		
			c:=c+l[i,2];
		 	s:=s union {l[i]};
		fi
	fi od;
	s:={op(l)} minus s;
	re:={[l[1,1],c]};
	re union procname([op(s)])
end;

#find shift eq of P
sheq:=proc(P)
	local li,c,a,i ;
	a:=1;
	li:=factors(P);
	c:=li[1];
	if c < 0 then c:=-c fi;
	li:=li[2];
	li:=pt(li);
	for i in li do
	if type(i[2], even) then 
		a:=a*i[1]^(i[2]/2);
	else return 1 
	fi od;
	sqrt(c)*a
end;
##########################################################################################
#gp

#returns [sing, gpmax-gpmin]
gpmaxmin:=proc(L)
	local F, ext, res, s, minmax,i;
	res:={};
	ext := indets([args], {'RootOf', 'radical', 'nonreal'});
	F := lcoeff(L, _tau)*tcoeff(L, _tau); 
	F := {seq(i[1], i = `LREtools/g_p/factors`(F, _n, ext))}; 	
	F := {seq(`LREtools/g_p/shorten_int`(i[1], _n), i = `LREtools/g_p/int_dif`(F, _n))}; 
	F:=PolynomialTools:-Sort([op(F)], _n);
	for s in F do
		minmax:=`LREtools/g_p`(primpart(L,_tau), RootOf(s,_n))[1];
		minmax := max(op(minmax))- min(op(minmax));
		res:= res union {[RootOf(s,_n),minmax]}
	od;
	res
end;
##########################################################################################
#gen exp
#quotient of two genexp of L of degree 2
genxq:=proc(gexp)
	local gexpq,tt, r,l, gg, a, b, d,i;
	gg:=gexp;
	r:=degree(rhs(gg[1,2]),_T);
	if nops(gg)=1 then 
		gg:= [allvalues(gg[1])];
		if nops(gg)=1 then return {gg[1,1]} fi;
	fi;

	gg:=[gg[1,1],gg[2,1]];
	a:= ldegree(gg[1], _T); b:=ldegree(gg[2], _T);

	if a=b then 
	gg:=[convert( series(gg[1]/gg[2], _T=0), polynom), convert( series(gg[2]/gg[1], _T=0), polynom)];
	gexpq:=[add( coeff( gg[1], _T, i )*_T^i , i=0..0+r ), add( coeff( gg[2], _T, i )*_T^i , i=0..r )];
	gexpq:={collect(gexpq[1], _T, evala), collect(gexpq[2], _T, evala)}

	elif a > b then
	gg:=convert( series(gg[1]/gg[2], _T=0), polynom);
	d:=ldegree(gg, _T);
	gexpq:=add( coeff( gg, _T, i )*_T^i , i=d..d+r );
	gexpq:={collect(gexpq, _T, evala)}

	else
	gg:=convert( series(gg[2]/gg[1], _T=0), polynom);
	d:=ldegree(gg, _T);
	gexpq:=add( coeff( gg, _T, i )*_T^i , i=d..d+r );
	gexpq:={collect(gexpq, _T, evala)}
	fi;
	
	return gexpq

end;


val:=proc(L, _tau, _t, r)
	local Ld, S, i,s,d;
	Ld:= collect(subs(_tau=delta+1, L),delta, simplify);
	d:=degree(Ld,delta);
	s:=min(seq(i+ldegree(coeff(Ld,delta,i),_t)/r ,i=0..d));
	S:={};
	for i from 0 to d do
		if i+ldegree(coeff(Ld,delta,i),_t)/r=s then
		S:= S union {i}
		fi;
	od;
	[s,S]
end;

Ind:=proc(L,_tau,_t, r)
	local Ld,v,M,P,i,k;
	Ld:= collect(subs(_tau=delta+1, L),delta, simplify);
	v:=val(L,_tau,_t, r);
	M:=v[2];
	P:=0;
	for i in M do
	P:=P+(-1)^i*mul(_n+k,k=0..i-1)*simplify(tcoeff(coeff(Ld,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 _tau ,_t;
	local d, P,j,i;
	d:=degree(L,_tau);
	P:=coeff(L,_tau,d)*_tau^d+add(simplify(coeff(L, _tau, i)*mul(subs(_t = Taut(_t,r,j,r*d+2), M), j = i .. d-1))*_tau^i, i = 0 .. d-1);
 	P:=primpart(P,_tau);
	P:=collect(P, _tau, simplify)
	end proc;

pwrcf:=proc(L,_tau)
local i;
	add(convert(series(coeff(L,_tau,i),_t=0),polynom)*_tau^i,i=0..degree(L,_tau))
end;
	

# L in C((_t))[_tau]
taupoly:=proc(L)
	local d,R,r,m,P,i,v,s,F, P1;
	d:=degree(L,_tau);
	if d=0 then return {} fi;
	for i from 0 to d do  
		v[i] := ldegree(coeff(L, _tau, i), _t) 
	end do;

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

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

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

	# Newton polynomial
	F:=[];
	P1:=add(tcoeff(coeff(L,_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 {[slope, NewtonPolynomial, c, ], ... }
	{seq([-s, factor(i[1]),i[2], denom(s)], `in`(i, P))} union procname(add(coeff(L, _tau, i)*_tau^i, i = 0 .. m))
end proc;

genexp:=proc(Lx)
	local Lc, Tp ,r, i, Gx, k, D, L ;
	L:=xtot(Lx);
	Tp:=taupoly(L);
	Gx:=[];
	for i in Tp do
	D:= gx(symprod(subs(_t=_t^i[4],L),1/(i[3]*_t^(i[1]*i[4])),i[4]), 1, 0, i[4] );
	Gx:=[op(Gx),  seq( [i[3]*_t^(i[1]*k[2])*k[1], k[2]] ,`in`(k, D)) ]
	od ;
	Gx:=[seq([subs(_t = _T, i[1]), _t = _T^i[2]], `in`(i, Gx))];
	
	return Gx	
end;

Lc:=proc(L)
	local Tp, cdd,i;
	cdd:=[];
	Tp:=taupoly(L);
	for i in Tp do 
	cdd:=[op(cdd),[symprod(subs(_t=_t^i[4],L),1/(i[3]*_t^(i[1]*i[4])), i[4]),i[4]]]
	od;
	return cdd
end;
	

gx:=proc(L, A, s, r)
	local Ld, Id, R, Dp, G, P, i;
	G:=[];
	R:=[];
	Ld:=symprod(L, 1/A, r);
	Id:=Ind(Ld,_tau,_t,r);
	if has(Id,_n) then 
		P:=factors(Id)[2];
		for i in P do 
		if has(i[1],_n) then 
		R:= [op(R), [A-RootOf(i[1],_n)*_t^r,r]]
		fi;
		od;
	return R
	fi;

	Dp:=deltapoly(Ld, s, r);
	procname(subs(_t=_t^Dp[1,4], L), subs(_t=_t^Dp[1,4], A)+Dp[1,3]*_t^(Dp[1,1]*Dp[1,4]),Dp[1,1]*Dp[1,4], r*Dp[1,4])


end;


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