_Env_LRE_x:=x;
_Env_LRE_tau:=tau;

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


#########################################################################################
solver:=proc(LL)
	local gp,L, gexp,gexpq, r, l,gpm,i, Sol;	
	L := evala(Primpart(LL,_Env_LRE_tau));
	#L:=LL;
	gpm:=gpmaxmin(L);
	gexpq:=gen(L)[1];
	r:=degree(gexpq, T);
	l:=ldegree(gexpq, T);
	
	#Bessel functions
	if l=2 then
		Bessel(LL,gexpq);


	#Laguerre1, 2f0
	elif l=1 and r=2 then 
	Sol:=Twofzero(L, gpm, gexpq);
		if Sol=FAIL then 
		Sol:=Laguerre1(L, gexpq);

			if Sol=FAIL then return FAIL
			else return Sol fi;

		else return Sol fi;


	#Jacobi, Jacobi1, Gener, Legendre, Laguerre, Gegen, Kummer, 2f1, Hermite, WhittakerM
	elif l=0 then 
	Sol:=Legendre(L, gexpq);
		if Sol=FAIL then 
		Sol:=Laguerre(L, gpm, gexpq);
			if Sol=FAIL then 
			Sol:=Gegenbauer(L, gpm, gexpq);
				if Sol=FAIL then 
				Sol:=Jacobi1(L, gpm, gexpq);
					if Sol=FAIL then
					Sol:=Jacobi(L, gpm, gexpq);
						if Sol=FAIL then 
						Sol:=Gener(L, gpm, gexpq);
							if Sol=FAIL then 
							Sol:=Twofone(L, gpm, gexpq);
								if Sol=FAIL then 
								 Sol:=hermite(L, gexpq);
									if Sol=FAIL then 
									Sol:=Kummer(L, gpm, gexpq);
										if Sol=FAIL then 
										return WhitM(L, gpm, gexpq)

										else return Sol fi;
									else return Sol fi;		
								else return Sol fi;
							else return Sol fi;
						else return Sol fi;					
					else return Sol fi;
				else return Sol fi;
			else return Sol fi;
		else return Sol fi;


	
	#WhittakerW
	elif l=0 and r=1 then 
	return WhitW(L, gpm, gexpq)


	else return FAIL
	fi;


	
end;
##########################################################################################
Kummer:=proc(LL, gp, g1)
	local c, L, b, a, b1, a1, cdd, i, sol, tsol;
	#if degree(g1, T)<>1 then return FAIL fi;
	c:= coeff(g1, T, 1)^2/4;
		

	if nops(gp)=0 or nops(gp)>2 then return FAIL fi;

	#b=0 mod Z
	if nops(gp)=1 then 
	b:=0;
	a:=-gp[1,1];
	L:=(a+x+1)*tau^2+(-2*a-2*x-2+b-c)*tau+a+x+1-b;
	cdd:={[L,a,b]};
	
	#general case
	else
	b1:=gp[1,1]-gp[2,1];
	b1:=[b1, -b1];
	a1:=[ -gp[1,1]+b1[1], -gp[2,1]+b1[1], -gp[1,1]+b1[2], -gp[2,1]+b1[2] ];
	L:=(a+x+1)*tau^2+(-2*a-2*x-2+b-c)*tau+a+x+1-b;

	cdd:= {seq(op(subs(b = j, {seq(subs(a = i, [L,a,b]), i = a1)})), j = b1)};

	fi;

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

end;

##########################################################################################
Laguerre1:=proc(LL, g1)
	local z, c2, d, b1, cdd, Lgr1, i, sol, tsol;
	z:=coeff(g1, T, 1);
	if z=0 then return FAIL fi;
	b1:={coeff(g1, T, 2)/(2*z), (coeff(g1, T, 2)+1)/(2*z)};
	Lgr1:=(x+2)*tau^2+(x+z-b+1)*tau+z;
	cdd:={seq( [subs(b=i, Lgr1),i], `in`(i, b1) )};

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

end;
##########################################################################################
Twofone:=proc(LL, gp, g1)
	local z1, L, d, c, a, b1, c1, a1, cdd, i, j, k, l, sol, tsol;
	if degree(g1, T)<>1 then return FAIL fi;
	L:=(x+1+a)*(z-1)*tau^2-(z+z*a+z*x-2*a-2*x-2-z*b+c)*tau-x-1-a+c;
	z1:=[-coeff(g1, T, 0)+1, 1-1/coeff(g1, T, 0)];
	d:=coeff(g1, T, 1)/coeff(g1, T, 0);
	
	#if c in Z
	if nops(gp)=1 then 
	c:=0;
	a:=-gp[1,1];
	b1:=[ -1/2*d,-1/2*(d+1),1/2*d,1/2*(d+1)]; 
	cdd:=NULL;
	for i in z1 do for j in b1 do 
	cdd:= cdd, subs( z=i, b=j, [L, a, b, c, z])
	od od;
	cdd:=[cdd];

	
	#general case
	else
	c1:=gp[1,1]-gp[2,1];
	c1:=[c1, -c1];
	b1:=[seq(op(subs(c=i, {(1/2)*c+(1/2)*d, (1/2)*c+(1/2)*(d+1), (1/2)*c-(1/2)*d, (1/2)*c-(1/2)*(d+1)})), c=c1)];
	a1:=[ -gp[1,1]+c1[1], -gp[2,1]+c1[1], -gp[1,1]+c1[2], -gp[2,1]+c1[2] ];
	cdd:=NULL;
	for i in z1 do for j in a1 do for k in b1 do for l in c1 do 
	cdd:=cdd, subs( z=i, a=j, b=k, c=l, [L, a, b, c, z]) 
	od od od od;
	cdd:=[cdd];	

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

end;

##########################################################################################
Twofzero:=proc(LL, gp, g1)
	local c, b, a1, Lc, cdd, sol, i, tsol ;
	if ldegree(g1,T)<>1 and degree(g1, T)<>2 then return FAIL fi;
	c:=-1/coeff(g1, T, 1);
	b:=gp[1,1];
	a1:=[(coeff(g1, T, 2)*c -b)/(-2), (coeff(g1, T, 2)*c -b+1)/(-2)];
	Lc:=tau^2+(-c*b+c*x+c+c*a-1)*tau+c*(b-x-1);
	cdd:={seq( [subs( a=i, Lc), i] , `in`(i, a1) )};

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

##########################################################################################
Gener:=proc(LL, gp, g1)
	local Lge, a1, bd1, bd2, b1, a, b, d, d1, cdd, K, i, j, k, sol, tsol;
	Lge:=(x+2)*tau^2+(-a*b-d+(a+1)*(1+x))*tau+a*x-a*(b+d);

		
	if degree(g1, T)=1 and gp<>{[0,1]} then 

		# b neq d mod
		if nops(gp)=2 then 
			a1:=[tcoeff(g1, T), 1/tcoeff(g1, T)];
			
			#when a=1
			if a1=[1,1] then return FAIL fi;


			bd1:=coeff(g1, T, 1)/tcoeff(g1, T);
			bd2:=gp minus {[0,1]};
			bd2:=bd2[1,1];
			b1:=[(bd1+bd2)/2, (bd1+bd2+1)/2, (bd2-bd1)/1, (bd2-bd1+1)/2];
			d1:=[(bd2-bd1)/1, (bd2-bd1+1)/2, (bd1+bd2)/2, (bd1+bd2+1)/2];
			cdd := []; 
			for i in a1 do for j in b1 do for k in d1 do 
			cdd := [op(cdd), [subs(a = i, b = j, d = k, Lge),i,j,k] ]
			od od od;

		#b=-d mod Z
		elif nops(gp)=1 then 
			a1:=[tcoeff(g1, T), 1/tcoeff(g1, T)];

			#when a=1
			if a1=[1,1] then return FAIL fi;


			bd1:=coeff(g1, T, 1)/tcoeff(g1, T);
			b1:=[bd1/2, (bd1+1)/2, -bd1/2, -(bd1+1)/2];
			cdd := []; 
			for i in a1 do for j in b1 do
			cdd := [op(cdd), [subs(a = i, b = j, d = -j, Lge), i, j, -j] ]
			od od;

		else return FAIL

		fi;

	#when a=1
	elif g1=1 then return FAIL
	
	# b=d mod Z
	elif degree(g1, T)=0 and gp<>{[0,1]} then
		bd2:=gp minus {[0,1]};
		bd2:=bd2[1,1];
		b1:=[bd2/2, (bd2+1)/2];
		Lge:=subs(a=g1, Lge);
		cdd:=[seq( subs(b=i, d=i, [Lge, g1, b, d]), `in`(i, b1) )];
	
	else return FAIL

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

end;
	
	


##########################################################################################
hermite:=proc(LL, g1)
	local z1, c, i, sol, tsol ;
	if degree(g1, T)<>2 then return FAIL fi;
	z1:=[ coeff(g1, T, 1)/(I*sqrt(2)), -coeff(g1, T, 1)/(I*sqrt(2))];
	
	c:= seq([subs(z=i, 2+2*x-2*z*tau+tau^2),i], `in`(i, z1) );

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


##########################################################################################
#When z= +- 1
Jacobi1:=proc(LL, gp, g1)
	local Ljc,  cdd, gp1, i, k, sol, tsol, abz, s,t,u;
	Ljc:=(2*(x+1+a))*(x+1+b)*(2*x+4+a+b)+(-(2*x+3+a+b)*(a^2-b^2)-(2*x+a+b+2)*(2*x+3+a+b)*(2*x+4+a+b)*z)*tau+(2*(x+2))*(x+2+a+b)*(2*x+a+b+2)*tau^2;
	
		# only one of a or b is integer
	if nops(gp)=2 and [0,2] in gp then
		gp1:=gp minus {[0,2]};
		abz:=[ [-gp1[1,1], 0, 1] , [0, -gp1[1,1], -1] ];
		cdd:=[ seq( [subs(a=i[1], b=i[2], z=i[3], Ljc), i[1], i[2], i[3]], `in`(i, abz) ) ];	

	# both a,b are not integer but a=b mod integer
	elif nops(gp)=2 and [0,1]in gp then 
		gp1:=gp minus {[0,1]};
		abz:=[ [-gp1[1,1]/2, 1], [-(gp1[1,1]+1)/2+1, 1] ];	
		cdd:=[ seq( [subs(a=i[1], b=i[1], z=i[2], Ljc), i[1], i[1], i[2]], `in`(i, abz) ) ];

	# both a,b are not integer but a=-b mod integer
	elif nops(gp)=2 and gp[1,2]=1 and gp[2,2]=1 then
		abz:=[ [gp[1,1], gp[2,1], 1], [gp[2,1], gp[1,1], 1] ];
		cdd:=[ seq( [subs(a=i[1], b=i[2], z=i[3], Ljc), i[1], i[2],i[3]], `in`(i, abz) ) ];
		

	# both a,b are not integer
	elif nops(gp)=4 and [0,1] in gp then
	 	gp1:=gp minus {[0,1]};
		s:=gp1[1,1]; t:=gp1[2,1]; u:=gp1[3,1];
		
		if s=t+u or s=t+u-1 or s=t+u+1 then 
			abz:=[ [-u, -t, 1], [-t, -u, 1], [-u, -t, -1], [-t, -u, -1] ];
		elif t=s+u or t=s+u-1 or t=s+u+1 then
			abz:=[ [-s, -u, 1], [-u, -s, 1], [-s, -u, -1], [-u, -s, -1] ];
		elif u=s+t or u=s+t-1 or t=s+u+1 then
			abz:=[ [-s, -t, 1], [-t, -s, 1], [-s, -t, -1], [-t, -s, -1] ];
		fi;
		cdd:=[ seq( [subs(a=i[1], b=i[2], z=i[3], Ljc), i[1], i[2],i[3]], `in`(i, abz) ) ];
	

	elif gp={} then
		abz:=[ [1/2, 1], [1/2, -1] ];	
		cdd:=[ seq( [subs(a=i[1], b=i[1], z=i[2], Ljc), i[1], i[1], i[2]], `in`(i, abz) ) ];
	
	else return FAIL
	
	fi;

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

	FAIL
end;
	
##########################################################################################
Jacobi:=proc(LL, gp, g1)
	local Ljc,  cdd, gp1, i, z1, sol, tsol, abz, s,t,u;
	if degree(g1, T)<>0 then return FAIL fi;
	Ljc:=(2*(x+1+a))*(x+1+b)*(2*x+4+a+b)+(-(2*x+3+a+b)*(a^2-b^2)-(2*x+a+b+2)*(2*x+3+a+b)*(2*x+4+a+b)*z)*tau+(2*(x+2))*(x+2+a+b)*(2*x+a+b+2)*tau^2;
	
	z1:=[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)];


	# only one of a or b is integer
	if nops(gp)=2 and [0,2] in gp then
		abz:=[ [-gp[1,1], -gp[2,1], z1[1]], [-gp[2,1], -gp[1,1], z1[1]], [-gp[1,1], -gp[2,1], z1[2]], [-gp[2,1], -gp[1,1], z1[2]] ];
		cdd:=[ seq( [subs(a=i[1], b=i[2], z=i[3], Ljc), i[1], i[2], i[3]], `in`(i, abz) ) ];	

	# both a,b are not integer but a=b mod integer
	elif nops(gp)=2 and [0,1]in gp then 
		gp1:=gp minus {[0,1]};
		abz:=[ [-gp1[1,1]/2, z1[1]], [-gp1[1,1]/2, z1[2]], [-(gp1[1,1]+1)/2, z1[1]], [-(gp1[1,1]+1)/2+1, z1[2]] ];	
		cdd:=[ seq( [subs(a=i[1], b=i[1], z=i[2], Ljc), i[1], i[1], i[2]], `in`(i, abz) ) ];

	# both a,b are not integer but a=-b mod integer
	elif nops(gp)=2 and gp[1,2]=1 and gp[2,2]=1 then
		abz:=[ [gp[1,1], gp[2,1], z1[1]], [gp[2,1], gp[1,1], z1[1]], [gp[1,1], gp[2,1], z1[2]], [gp[2,1], gp[1,1], z1[2]] ];
		cdd:=[ seq( [subs(a=i[1], b=i[2], z=i[3], Ljc), i[1], i[2],i[3]], `in`(i, abz) ) ];
		

	# both a,b are not integer
	elif nops(gp)=4 and [0,1] in gp then
	 	gp1:=gp minus {[0,1]};
		s:=gp1[1,1]; t:=gp1[2,1]; u:=gp1[3,1];
		
		if s=t+u or s=t+u-1 or s=t+u+1 then 
			abz:=[ [-u, -t, z1[1]], [-t, -u, z1[1]], [-u, -t, z1[2]], [-t, -u, z1[2]] ];
		elif t=s+u or t=s+u-1 or t=s+u+1 then
			abz:=[ [-s, -u, z1[1]], [-u, -s, z1[1]], [-s, -u, z1[2]], [-u, -s, z1[2]] ];
		elif u=s+t or u=s+t-1 or t=s+u+1 then
			abz:=[ [-s, -t, z1[1]], [-t, -s, z1[1]], [-s, -t, z1[2]], [-t, -s, z1[2]] ];
		fi;
		cdd:=[ seq( [subs(a=i[1], b=i[2], z=i[3], Ljc), i[1], i[2],i[3]], `in`(i, abz) ) ];

	else return FAIL

	fi;
	

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

	FAIL
end;
		




##########################################################################################
Gegenbauer:=proc(LL, gp, g1)
	local c , mm, i , sol, tsol;
	if degree(g1, T)<>0 or nops(gp)<>1 then return FAIL fi;
	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)];
	
	mm:=[-gp[1,1]/2, -(gp[1,1]+1)/2];
	
	c := [seq([subs(z = i, tau^2-2*z*(m+x+1)*tau/(x+2)+(2*m+x)/(x+2)), i, m], `in`(i, c))]; 
	c := [seq(op(subs(m = i, c)), `in`(i, mm))];

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

##########################################################################################
Legendre:=proc(LL, g1)
	local c, z, i, sol, tsol;
	if degree(g1, T)<>0 then return FAIL fi;
	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*x+3)*z*tau/(x+2)+(x+1)/(x+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+x)-sol[1]*a(x), a(x), {}, output = basis)[1];
	return [plgop(sol[2],LegendreP(x, i[2])*tsol),plgop(sol[2],LegendreQ(x, i[2])*tsol)]
	fi od;
	FAIL
end;
##########################################################################################
Laguerre:=proc(LL, gp, g1)
	local c, sol, tsol, aa;
	c:=coeff(g1, T, 1);
	c:=-c^2/4;
	
	if gp={} then return FAIL fi;

	if nops(gp)=2 then
		aa:=-(gp[1,1]+gp[2,1])
	elif nops(gp)=1 then
		aa:=-gp[1,1];
	else return FAIL
	fi;

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

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

	

	for j in cdd do
	if symg(j[1],LL)<>FAIL then 
		sol:=symg(j[1],LL);
		tsol:=`LREtools/hypergeomsols`(a(1+x)-sol[1]*a(x), a(x), {}, output = basis)[1];
		return [plgop(sol[2],WhittakerM(j[2]+x,j[3],j[4])*tsol)]
	fi od;
	FAIL
end;
##########################################################################################
#bessel function
Bessel:=proc(LL,gexpq)
	local c,cdd,sol,i,j,vl,si,cx,d,z,tsol;
	if ldegree(gexpq, T)<>2 then return "Not Bessel" fi;
	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*x)*tau+z,v,z], v in vl), seq([z*tau^2+(2+2*v+2*x)*tau+z,v,z], v in vl)];
	else
		z:=2*sqrt(-c);
		cdd:=[seq([z*tau^2-(2+2*v+2*x)*tau-z,v,z], v in vl), seq([z*tau^2+(2+2*v+2*x)*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),x);
		tsol:=`LREtools/hypergeomsols`(a(1+x)-sol[2,1]*a(x), a(x), {}, output = basis)[1];
		if si > 0 and cx > 0 then
		sol:= [plgop(sol[2,2],BesselJ(sol[1,2]+x,-sol[1,3])*tsol),plgop(sol[2,2],BesselY(sol[1,2]+x,-sol[1,3])*tsol)];
		elif si > 0 and cx < 0 then
		sol:= [plgop(sol[2,2],BesselJ(sol[1,2]+x,sol[1,3])*tsol),plgop(sol[2,2],BesselY(sol[1,2]+x,sol[1,3])*tsol)];
		elif si < 0 and cx < 0 then
		sol:= [plgop(sol[2,2],BesselI(sol[1,2]+x,-sol[1,3])*tsol),plgop(sol[2,2],BesselK(sol[1,2]+x,sol[1,3])*tsol)];
		elif si < 0 and cx > 0 then
		sol:= [plgop(sol[2,2],BesselI(sol[1,2]+x,sol[1,3])*tsol),plgop(sol[2,2],BesselK(sol[1,2]+x,-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;
	n:=degree(L,tau);
	r:=M;
	coeff(L,tau,n)*tau^n+add(simplify(coeff(L, tau, i)*(mul(subs(x = x+j, r), j = i .. n-1)))*tau^i, i = 0 .. n-1)
end proc;

xtot:=proc(L)
	global x, t, tau;
	collect(primpart(subs(x=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)
	add(subs(x=n,coeff(L,tau,i))*u(n+i),i=0..degree(L,tau))
end;



tautov:=proc(L)
	add(subs(x=n,coeff(L,tau,i))*v(n+i),i=0..degree(L,tau))
end;



vtotau:=proc(L)
	subs(n=x,add(coeff(L,v(n+i))*tau^i, i=0..2))
end;



det:=proc(L)
	factor(coeff(L,tau,0)/coeff(L,tau,2))
end;



plgop:= proc (L, y) 
simplify(add(coeff(L, _Env_LRE_tau, i)*subs(_Env_LRE_x = x+i, y), i = 0 .. degree(L, _Env_LRE_tau))) 
end proc;

gen:=proc(L)
genxq(genexp(L))
end;

findgauge:=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 or L1=-L2 then return 1 fi;
	cfs := indets([L1,L2],function);
	cfs := cfs union subs(_Env_LRE_x=_Env_LRE_x +1, cfs);
	cfs := cfs union subs(_Env_LRE_x =_Env_LRE_x +2, cfs);
	G := c1(_Env_LRE_x )*_Env_LRE_tau + c0(_Env_LRE_x );
	K:=collect(`LREtools/rightdivision`( `LREtools/mult`(L2, G), L1)[2], tau);
	eq0, eq1 := coeff(%,_Env_LRE_tau,0), coeff(%,_Env_LRE_tau,1);
	if coeff(L1,_Env_LRE_tau,1)=0 and coeff(L2,_Env_LRE_tau,1)<>0 then
		c:=solve(eq1, c0(_Env_LRE_x +1));
		eval(eq0, {c0(_Env_LRE_x ) = subs(_Env_LRE_x  = _Env_LRE_x -1, %), c0(_Env_LRE_x +2) = subs(_Env_LRE_x  = _Env_LRE_x +1, %)});
		collect(%, {c1(_Env_LRE_x ),c1(_Env_LRE_x +1),c1(_Env_LRE_x +2),c1(_Env_LRE_x +3),c1(_Env_LRE_x +4)}, factor);
		d:=LREtools[ratpolysols](%, c1(_Env_LRE_x ), {}, output = basis);
		if whattype(d)=exprseq then return FAIL fi;
		d:=d[1];
		eval(c, {c1(_Env_LRE_x ) = d, c1(_Env_LRE_x +2) = subs(_Env_LRE_x =_Env_LRE_x +2,d), c1(_Env_LRE_x +1) = (_Env_LRE_x =_Env_LRE_x +1,d)});
		return collect(simplify(d*_Env_LRE_tau+subs(_Env_LRE_x =_Env_LRE_x -1,%)),_Env_LRE_tau)
	elif coeff(L1,_Env_LRE_tau,1)<>0 and coeff(L2,_Env_LRE_tau,1)=0 then
		c:=solve(eq1, c0(_Env_LRE_x +2));
		eval(eq0, {c0(_Env_LRE_x +2) = %, c0(_Env_LRE_x ) = eval(%, _Env_LRE_x  = _Env_LRE_x -2)});
		d:=LREtools[ratpolysols](%, c1(_Env_LRE_x ), {}, output = basis);
		if whattype(d)=exprseq then return FAIL fi;
		d:=d[1];
		eval(c, {c1(_Env_LRE_x ) = %, c1(_Env_LRE_x +2) = subs(_Env_LRE_x =_Env_LRE_x +2,%), c1(_Env_LRE_x +1) = (_Env_LRE_x =_Env_LRE_x +1,%)});
		return collect(simplify(d*_Env_LRE_tau+subs(_Env_LRE_x =_Env_LRE_x -2,%)),_Env_LRE_tau)
	elif coeff(L1,_Env_LRE_tau,1)=0 and coeff(L2,_Env_LRE_tau,1)=0 then
		d:=LREtools[ratpolysols](eq1, c1(_Env_LRE_x ), {}, output = basis);
		e:=LREtools[ratpolysols](eq0, c0(_Env_LRE_x ), {}, 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(_Env_LRE_x +2) ) * eq0 - coeff(eq0, c0(_Env_LRE_x +2) ) * eq1); # eq2 has no c0(_Env_LRE_x +2) term.
		eq3 := simplify(subs(_Env_LRE_x =_Env_LRE_x +1, eq2));
		eq4 := simplify(coeff(eq1, c0(_Env_LRE_x +2) ) * eq3 - coeff(eq3, c0(_Env_LRE_x +2) )* eq1); # Has no c0(_Env_LRE_x +1) term
		if has(indets(eq4),c0(_Env_LRE_x +1)) then 		
			c0_x1 := solve( expand(eq4), c0(_Env_LRE_x +1) );
			eq5 := simplify(coeff(eq2, c0(_Env_LRE_x +1) ) * eq4 - coeff(eq4, c0(_Env_LRE_x +1) ) * eq2);
			eq6 := simplify(coeff(eq0, c0(_Env_LRE_x +2) ) * eq3 - coeff(eq3, c0(_Env_LRE_x +2) ) * eq0);
			if not has(eq2, c0(_Env_LRE_x )) then
			c0_x0 := solve(expand(eq6), c0(_Env_LRE_x ))
			else
			c0_x0 := solve(expand(eq5), c0(_Env_LRE_x ))
			fi;
			eqc1 := subs(_Env_LRE_x =_Env_LRE_x +1, c0_x0) - c0_x1;
		else eqc1:=eq4
		fi;
		d:=LREtools[ratpolysols](%, c1(_Env_LRE_x ), {}, output = basis);
		if whattype(d)=exprseq then d:=0
		else d:=d[1];
		#elif nops(d)>2 then
		#	error "input must be reducible"
		fi;

		if d<>0 then 
		c := subs({c1(_Env_LRE_x )=d, c1(_Env_LRE_x +1)=subs(_Env_LRE_x =_Env_LRE_x +1,d), c1(_Env_LRE_x +2)=subs(_Env_LRE_x =_Env_LRE_x +2,d),c1(_Env_LRE_x +3)=subs(_Env_LRE_x =_Env_LRE_x +3,d)}, c0_x0);
		c:=d*tau + normal(c);
		return collect(simplify(c),tau)
		
		elif d=0 then
		c0_x0:=subs(c1(_Env_LRE_x)=0, c1(_Env_LRE_x+1)=0, c1(_Env_LRE_x+2)=0, eq2);
		c:=LREtools[ratpolysols](c0_x0, c0(_Env_LRE_x ), {}, output = basis);
			if whattype(c)=exprseq then return FAIL
			else return c[1]
			fi;
		fi;

	fi;
end;
##########################################################################################
#symg
#find symprod and gauge transformation that sends L1 to L2
symg:=proc(L1, L2)
	local A, LsA, cdd, i, d;
	d:=det(L2)/det(L1);
	A:=sheq(d);
	cdd:=[[symprodx(L1,A),A],[symprodx(L1,-A),-A]];
	for i in cdd do
	if findgauge(i[1],L2)<> FAIL then
		return [i[2], findgauge(i[1],L2)]
	fi od;
	FAIL
end;
	
vutau:=proc(L)
	subs({n=x,_c[1]=1},rhs(L));
	add(coeff(%,u(x+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],x)=degree(l[i,1],x) then
		d:=degree(l[1,1],x)-1;
		e:=Normalizer(coeff(l[1,1],x,d)-coeff(l[i,1],x,d));
		if type(e,integer) and e<>0 and irem(e,d+1)=0 and
		 Normalizer(l[1,1]-subs(x=x+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 evalf(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, gpm, i;
	res:={};
	ext := indets([args], {'RootOf', 'radical', 'nonreal'});
	F := lcoeff(L, _Env_LRE_tau)*tcoeff(L, _Env_LRE_tau); 
	F := {seq(i[1], i = `LREtools/g_p/factors`(F, _Env_LRE_x, ext))}; 	
	F := {seq(`LREtools/g_p/shorten_int`(i[1], _Env_LRE_x), i = `LREtools/g_p/int_dif`(F, _Env_LRE_x))}; 
	F:=PolynomialTools:-Sort([op(F)], _Env_LRE_x);
	for s in F do
		minmax:=`LREtools/g_p`(primpart(L,_Env_LRE_tau), RootOf(s,_Env_LRE_x))[1];
		minmax := max(op(minmax))- min(op(minmax));
		res:= res union {[RootOf(s,_Env_LRE_x),minmax]}
	od;
	
	gpm:={};
	if res<>{} then 
	for i in res do 
	if i[2]<>0 then 
	gpm:=gpm union {i}; 
	fi od fi;
	
	gpm

end;
##########################################################################################
#gen exp
#quotient of two genexp of L of degree 2
genxq:=proc(gexp)
	local gexpq,tt, r,l, gg, a, b, d;
	gg:=gexp;
	r:=degree(rhs(gg[1,2]),T);
	if nops(gg)=1 then 
		gg:= [allvalues(gg[1])];
		if nops(gg)=1 then return {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 ,n;
	Ld:= collect(subs(tau=delta+1, L),delta, simplify);
	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,tau,t, r)
	local Ld,v,M,P,i;
	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 n, P;
	n:=degree(L,tau);
	P:=coeff(L,tau,n)*tau^n+add(simplify(coeff(L, tau, i)*mul(subs(t = Taut(t,r,j,r*n+2), M), j = i .. n-1))*tau^i, i = 0 .. n-1);
 	P:=primpart(P,tau);
	P:=collect(P, tau, simplify)
	end proc;

pwrcf:=proc(L,tau)
	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 n,R,r,m,P,i,v,s,F, P1;
	n:=degree(L,tau);
	if n=0 then return {} fi;
	for i from 0 to n do  
		v[i] := ldegree(coeff(L, tau, i), t) 
	end do;

	# 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,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))];
	P:=allvalues( P ); 

	# 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, n,  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);
	for i in Dp do 
	procname(subs(t=t^i[4], L), subs(t=t^i[4], A)+RootOf(i[2],t)*t^(i[1]*i[4]),i[1]*i[4], r*i[4])
	od;

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;