bessellist:= proc(L)
	#ption trace;
	local Dx,x,Sreg, Sirr, dA,B,A,ext,res,i,neq,und,dinf,mu,m,fac,mun,f,t,j,s,ind,d1,d2,d,dv,knl,l,k;
	Dx, x := op(_Envdiffopdomain); 
	ext := indets(L, {RootOf, nonreal, radical,name}) minus {x,Dx};
	Sirr,Sreg,B,dA,dinf := handleSingularity(L,Dx,x,ext);
	neq:=0;
	for i in Sirr do 
		if i[1]= infinity then 
			neq:=neq+i[4];
		else
			neq:=neq+i[4]*degree(evala(Norm(i[2],ext,ext),x))
		fi;
	od;
	res := {};
	for i in Sreg do
		neq:=neq+degree(evala(Norm(x-i[1],ext,ext),x));
	od;
	if neq>dA then
		f := easycase(Sirr,Sreg,B,dA,ext,Dx,x);
		if f<>NULL then
			if Sreg = [] then 
				mu := seq((dinf+i)/degree(B,x),i=0..degree(B,x)-1);
				mu := {seq((numer(i) mod denom(i))/denom(i),i=mu)};
			else
				t := 0;
				for s in Sreg do
					m := 0;
					fac := evala(Norm(x-s[1],ext,ext),x);
					while evala(Rem(numer(f), fac^(m+1), x)) = 0 do
						m := m+1;
					od;
					if t=0 then 
						mu := seq((s[2]+i)/m,i=0..m-1);
						mu := seq((numer(i) mod denom(i))/denom(i),i=[mu]);
						mu := {mu};
					else
						mun := seq((s[2]+i)/m,i=0..m-1);
						mun:= seq((numer(i) mod denom(i))/denom(i),i=[mun]);
						mun := {mun};
						mu := mu intersect mun;
					fi;
					t := t+1;
				od;
				if dinf<>"f" then 
					m := degree(B,x)-degree(numer(f),x);
					mun := seq((dinf+i)/m,i=0..m-1);
					mun:= seq((numer(i) mod denom(i))/denom(i),i=mun);
					mun := {mun};
					mu := mu intersect mun;
				fi;
			
				 
			fi;
			
			for i in mu do
				res := res union {[compute_sqrt(f,x,ext),i]};
			od; 
		fi;
		
	else
		if Sreg=[] and dinf="f" then 
			mu := {};
			m := dA;
			for i from 3 to floor(m/2) do
				if irem(m,i)=0 then mu:= mu union {i};fi;			
			od;
			mu := mu union {m};
			res := {};
			for i in mu do
				knl :=findukn(Sreg,dA,m,i,dinf,ext);
				f,und:= rationalcase(Sirr,B,i,knl,ext,Dx,x);
				while f=0 and und<>{} do 
					ext := ext union und;
					Sirr,Sreg,B,dA,dinf := handleSingularity(L,Dx,x,ext);
					knl :=findukn(Sreg,dA,degree(B,x),i,dinf,ext);
					f,und := rationalcase(Sirr,B,i,knl,ext,Dx,x);
				od;
				for j in f do
					for k from 1 to floor(i/2) do
						if j<>0  and gcd(k,i)=1 then res := res union {[j,k/i]} fi;
					od;
				od;
			od;
		else
			if dinf<>"f" then 
				t := [infinity,dinf]; 
				d :=degree(B,x); 
			else 
				t := Sreg[1]; 
				d:=dA;
			fi;
			if type(t[2],`integer`) then
				if dinf<>"f" then Sreg :=[seq(l,l=Sreg),[infinity,dinf]];fi;
				knl := searchknlog(Sreg,d,ext);
				f:= logcase(Sirr,B,knl,ext,Dx,x);
				for j in f do
					res:= res union {[j,0]};
				od;
			elif type(t[2],`rational`) then
				mu := {};
				dv := denom(t[2]);
				
				for s in Sreg do
					dv:=lcm(dv,denom(s[2]));
					
				od;
				for i from 1 to max(1,iquo(d,dv)) do
					d1:=0;
					d2:=0;
					if dinf<>"f" then 
						d1 := i*dv/denom(dinf);
						d2 := dv/denom(dinf);
					fi;
					for s in Sreg do
						d1 := d1+i*dv/denom(s[2])*degree(evala(Norm(x-s[1],ext,ext),x));
						d2 := gcd(d2, i*dv/denom(s[2]));
					od;
					if d1<=d and irem(d,d2)=0 then
						 mu := mu union {i*dv}; 
												
					fi;
				od;					
				res := {};
				for i in mu do
					if dinf<>"f" then Sreg :=[seq(l,l=Sreg),[infinity,dinf]];fi;
					knl :=findukn(Sreg,dA,degree(B,x),i,dinf,ext);
					f,und:= rationalcase(Sirr,B,i,knl,ext,Dx,x);
					while f=0 and und<>{} do 
						ext:= ext union und;
						Sirr,Sreg,B,dA,dinf := handleSingularity(L,Dx,x,ext);
						if dinf<>"f" then Sreg :=[seq(l,l=Sreg),[infinity,dinf]];fi;
						knl :=findukn(Sreg,dA,degree(B,x),i,dinf,ext union und);
						f,und := rationalcase(Sirr,B,i,knl,ext,Dx,x);
					od;
					for j in f do
						for k from 1 to floor(i/2) do
							if j<>0  and gcd(k,i)=1 then res := res union {[j,k/i]} fi;
						od;
					od;
				od;
			else
				ind := indets(t[2],{RootOf, nonreal, radical,name});
				ind := ind[1];
				m := 0;
				for s in Sreg do m:= m+degree(evala(Norm(x-s[1],ext,ext),x))*abs(coeff(s[2],ind,1)/coeff(t[2],ind,1)) od;
				if dinf="f" then
					m := dA/m ;
				else
					m:= m + 1;
					m := degree(B,x)/m;
					
				fi;
				mu := {seq((t[2]+i)/m,i=0..m-1)};
				
				
				res := {};
				for i in mu do
					f := irrationalcase(Sirr,Sreg,B,dA,i,ext,Dx,x);
					if f<>NULL then res := res union {[f,i]};fi;
				od;

			fi;
		fi;
		
	fi;
	res;
end:

handleSingularity := proc(L,Dx,x,ext)
	#option trace;
	local S,Sirr,Sreg,B,dA,temp,T,d,l,i,j,acc,p,g,t,dinf;
	dinf := "f";
	S := NULL;
	temp := evala(Factors(denom(L),ext))[2]; 
	for i in temp do if has(i[1],x) then S :=S,i[1] fi;od;
	temp := evala(Factors(coeff(L, Dx, 2), ext))[2];
 	for i in temp do if has(i[1],x) then S :=S,i[1] fi;od;
	S := [infinity, seq(RootOf(i,x), i = S)] ;
	dA:=0;
	B:=1;
	Sirr:=NULL;
	Sreg:=NULL;
	for p in S do 
		g := DEtools[gen_exp](L, [Dx, x], T, x = p); 
		t := `if`(p = infinity, 1/x, x-p);		 
		if `and`(nops(g) = 1, nops(g[1]) = 2) then 
			
			if lhs(g[1, -1]) <> T then 
				if p <> infinity then 
					B := B*evala(Norm((x-p)^(-ldegree(g[1][1], T)), ext, ext), x); 
				else 
					dA := -ldegree(g[1][1], T); 
				fi;
				
				d := 2*(g[1, 1]-coeff(g[1, 1], T, 0));
				acc := (1-ldegree(d,T))/2;
				if type(d, `+`) then l := [op(d)] else l := [d] fi;
				d := add(i/ldegree(i, T), i = l); 
				d := expand(d^2); 
				temp := evala(rhs(g[1,2])/coeff(lhs(g[1,2]),T,2)); 
				if type(d, `+`) then l := [op(d)] else l := [d] fi; 
				d := add(i*T^(-(1/2)*ldegree(i, T)), i = l); 
				if type(d, `+`) then l := [op(d)] else l := [d] fi; 
				j := (1/2)*ldegree(d, T); 
				d := 0; 
				for i in l do if degree(i, T) < j then d := d+i fi od; 
				d := subs(T = temp, d); 
				# d := evala(Trace(d, ext, ext), x); 
				Sirr := Sirr, [p,t, d,acc]; 
			else 
				if p <> infinity then 
					B := B*evala(Norm((x-p)^(-2*ldegree(g[1][1], T)), ext, ext), x);
				else 
					dA := -2*ldegree(g[1][1], T); 
				fi;
				d := evala(Trace(g[1, 1], ext, ext), x)-2*g[1, 1];
				acc := -ldegree(d,T); 			
				d := expand(d-coeff(d,T,0)); 
				if type(d, `+`) then l := [op(d)] else l := [d] fi;
				d := add(i/(2*ldegree(i, T)), i = l); 
				d := expand(d^2); 
				d := evala(subs(T = t, d)); 
				# d := evala(Trace(d, ext, ext), x);
				# d := poly_to_powerseries(d,x,p,T,acc);
				Sirr := Sirr, [p,t, d,acc]; 
			fi; 
		else 
			d := g[1, 1]-`if`(nops(g) = 2, g[2, 1], g[1, 2]); 
			acc:=-ldegree(d,T);			 
			if has(d, T) then 
				if p <> infinity then 
					B := B*evala(Norm((x-p)^(-2*ldegree(g[1][1], T)), ext, ext), x); 
				else 
					dA := -2*ldegree(g[1][1], T); 
				fi;
				d := expand(d-coeff(d,T,0)); 
				if type(d, `+`) then l := [op(d)] else l := [d] fi; 
				d := add(i/(2*ldegree(i, T)), i = l); 
				d := d^2;
				d := evala(subs(T = t, d)); 
				# d := evala(Trace(d, ext, ext),x); 
				# d := poly_to_powerseries(d,x,p,T,acc);
				Sirr := Sirr, [p,t,d,acc] 
			else 
				if not type(d, integer) or d = 0 or DEtools['formal_sol'](L, `has logarithm?`, x = p) then 
					if p = infinity then 
						dA := -1; 
						dinf := d;
					else 
						Sreg := Sreg, [p,d] ;
					fi; 
				fi; 
			fi; 
		fi; 
	od; 
	dA := degree(B, x)+dA;
	[Sirr],[Sreg],B,dA,dinf
end:

rationalcase := proc(Sirr,B,d,knl,ext,Dx,x)
	#option trace;
	local res,T,s,uk,i,kn,j,lr,acc,ep,temp,ukn,Bs,extp,und,f,l,pot;
	res := {};		
	if knl <> {} then 
		for kn in knl do
			l,pot := computel(Sirr,B,kn[1],T,ext);
			if l=0 then return l,pot; fi;
				
			uk := add(a[i]*x^i,i=0..kn[2]);
			und := indets(uk) minus {x};
			und := [op(und)];
			for s in Sirr do 
				acc := s[4];
				ep := poly_to_powerseries(kn[1],x,s[1],T,acc);
				ep := pow_reciprocal(ep,T);
				Bs := poly_to_powerseries(B,x,s[1],T,acc);
				temp := pow_time(ep, Bs, T);
				ep := poly_to_powerseries(s[3],x,s[1],T,acc);
				temp := pow_time(temp, ep,T);
				if temp[2] <> 0 and s[1]<>infinity then error "we made an error" fi;
				temp[1] := temp[1]/l;
				extp := indets(s[1],RootOf) union ext;
				ep := dthroot_powseries(temp,T,d,extp);
				ukn := NULL;								 
				for i in {uk} do 
					temp := poly_to_powerseries( i,x,s[1],T,acc);
					for j in ep do
						lr := temp[3]-j[3];
						lr := coeffs(lr,T);
						lr := {lr,temp[1]-j[1]};
						lr := map(numer,lr);
						if type(s[1],RootOf) and not member(s[1],ext) then
							lr := map(evala@Expand, lr);
							lr := map(coeffs, lr, s[1])
						fi;
						lr := solve(lr,{op(und)});
						if lr <> [] and lr <> NULL then
							f:= subs(op(lr),i);
							ukn := ukn,f;
						fi;
					od;
				od;
				uk := ukn;
				if uk = NULL then break fi;
			od;
			if uk <> NULL then res := res union {seq(compute_sqrt(l*i^d*kn[1]/B,x,ext),i=[uk])} fi;
		od;
	fi;
	res,{};
end:



computel := proc(Sirr,B,kn,T,ext)
	#option trace;
	local l,s,ep,und,Bs,temp,acc;
	l := 0;
	und := {};
	for s in Sirr do
		acc := s[4]; 
		ep := poly_to_powerseries(kn,x,s[1],T,acc);
		ep := pow_reciprocal(ep,T);
		Bs := poly_to_powerseries(B,x,s[1],T,acc);
		temp := pow_time(ep, Bs, T);
		ep := poly_to_powerseries(s[3],x,s[1],T,acc);
		temp := pow_time(temp, ep,T);
		und := indets(temp[1],RootOf) minus ext;
		if und={} then 
			l:= temp[1];
			break;
		fi;
	od;
	l,und;
end:

irrationalcase := proc (Sirr,Sreg,B,dA,mu,ext,Dx,x)
	#option trace;
	local res,i,m,s,mul,A,c,acc,ep,Bs,temp,lr;	
	m := nops(Sreg);
	A :=c;
	if m <>0 then 
		for s in Sreg do
			for i from 1 to dA-m+1 do
				if type((i*mu-s[2]),`integer`) or type ((i*mu+s[2]),`integer`) then 
					A := A*evala(Norm(x-s[1], ext, ext),x)^i;
					break;
				fi;
			od;
		od;
	fi;
	if degree(A,x)>dA then error "degree doesn't march" fi;
	for s in Sirr do 
		acc := s[4];
		Bs := poly_to_powerseries(B,x,s[1],T,acc);
		ep := poly_to_powerseries(A,x,s[1],T,acc);
		temp:=poly_to_powerseries(s[3],x,s[1],T,acc);
		temp := pow_time(Bs,temp,T);
		lr := temp[3]-ep[3];
		lr := coeffs(lr,T);
		lr := {lr,temp[1]-ep[1]};
		if type(s[1],RootOf) and not member(s[1],ext) then
				lr := map(evala@Expand, lr);
				lr := map(coeffs, lr, s[1])
		fi;
		lr := solve(lr,c);
		if lr = NULL then
			res := NULL;
			break;
		else
			A:=subs(lr,A);
			res := compute_sqrt(A/B,x,ext);
		fi;
	od;
	res;	
end:

easycase := proc(Sirr,Sreg,B,dA,ext,Dx,x)
	#option trace;
	local f,f1,p,i,j,t,linearr,temp,res,d,d1,nonpole,und;					  
	f := add(a[i]*x^i, i = 0 .. dA)/B;
 	und := indets(f) minus {x} ;
	und := [op(und)];
	f1 := numer(f); 
	for p in Sreg do 
		t := evala(Norm(x-p[1], ext, ext),x); 
		f1 := f1/t; 
	od; 
	linearr := evala(Rem(numer(f1), denom(f1), x)); 
	linearr := coeffs(linearr, x); 
	temp := solve({linearr},und); 
	f := subs(op(temp), f); 
	d := 0; 
	j := 0;
	for p in Sirr do 
		
		if p[1] = infinity then 
			d := d+ p[3];
			j := p[4]; 
		else
			d := d+evala(Trace(p[3], ext, ext), x); 
		fi; 
	od; 
	nonpole := f-d; 
	for p in Sirr do
		if p[1]<> infinity then 
			t := evala(Norm(x-p[1], ext, ext), x); 
			nonpole := normal(nonpole*t^p[4]);
		fi; 
	od; 
	temp := evala(Rem(numer(nonpole), denom(nonpole), x));
	linearr := coeffs(temp, x);
	if j<>0 then
		temp := evala(Quo(numer(nonpole), denom(nonpole),x));
		d1:=degree(temp, x); 
		temp := seq(coeff(temp, x, i), i = d1-j+1 .. d1); 
		linearr := linearr, temp;
	fi; 
	temp := solve({linearr},und);
	temp := op(temp); 
	if temp <> NULL then 
		f := subs(temp, f); 
		res := f; 
	else 
		res := NULL;
	fi;
	res;
end:			
poly_to_powerseries := proc(P, x, p, T, acc)
	# option trace;
	local d, l, i, S;
	if P=0 then error "only treat non-zero powseries" fi;
	if p = infinity then
		d, l:= degree(P,x), lcoeff(P,x);
		[l, -d,
		 add(evala(coeff(P,x,d-i)/l)*T^i, i=0..acc-1), acc]
	else
		S := collect(eval(P, x = T + p), T, evala);
		d, l := ldegree(S,T), tcoeff(S,T);
		# Note: the evala in the line S := .. is
		# needed to ensure that d is correct.
		[l, d,
		 add(evala(coeff(S,T,d+i)/l)*T^i, i=0..acc-1), acc]
	fi
end:

dthroot_powseries := proc(L, tp, d, ext)
	#option trace;
	local X,v,res,i;
	v := evala(Factors(X^d - L[1], ext))[2];
	res := NULL;
	for i in v do
		if degree(i[1],X) = 1 then
			res := res, -coeff(i[1],X,0)/coeff(i[1],X,1)
		fi
	od;
	v := series(L[3]^(1/d), tp=0, L[4]);
	v := collect(convert(v,polynom),tp,evala);
	res :={seq([i, L[2]/d, v, L[4]],  i = [res])};
	res;
end:

pow_time := proc(L1,L2,tp)
	local acc,f,ind;
	acc := min(L1[4],L2[4]);
	f := expand(L1[3]*L2[3]); 
	ind := [op(indets(f))]; 
	if nops(ind) = 2 then f := subs(ind[1] = tp, f); f := subs(ind[2] = tp, f) fi;	
	[L1[1]*L2[1],L1[2]+L2[2],
	  convert(series(f,tp=0,acc),polynom),
	  acc]
end:

pow_reciprocal := proc(L1,tp)
	[1/L1[1],-L1[2],
	 convert(series(1/L1[3],tp=0,L1[4]),polynom),
	L1[4]]
end:

findukn :=proc(Sreg,dA,dB,d,dinf,ext)
	#option trace;
	local de,s,dr,k,res,temp;
	res :=NULL;
	if dinf<>"f" then 		
		de:=dB;
	else		
		de:=dA;
	fi;
	if nops(Sreg)=0 then 
		if irem(dA,d)<>0 then 
			return {};
		else
			return {[1,dA/d]};
		fi;
	fi;
	for k from 0 to iquo(de,d) do
		dr :=de-k*d;
		temp :=searchkn(Sreg,dr,d,ext);
		for s in temp do res := res, [s,k] od;
	od;
	{res};	
end:

searchkn :=proc(Sreg,de,d,ext)
	#option trace;
	local s,res,i,temp,k,t;	
	res :=NULL;
	if nops(Sreg)=0 then 
		if de=0 then 
			return 1; 
		else 
			return NULL;
		fi;
	fi;
	s :=Sreg[1];
	i :=d/denom(s[2]);
	if (de-i)<0 then return NULL; fi;
	for k from 1 to iquo(de-nops(Sreg)+1,i) do
		if i*k>=d then break; fi;
		temp := searchkn([seq(Sreg[i],i=2..nops(Sreg))],de-k*i,d,ext);
		if temp<>NULL then 
			if s[1]<> infinity then 
				for t in temp do
					res := res,evala(Norm(x-s[1], ext, ext),x)^(i*k)*t;
				od;
			else
				for t in temp do
					res := res, t;
					
				od;
			fi;
		fi;		
	od;	
	[res];
end:

searchknlog :=proc(Sreg,d,ext)
	#option trace;
	local res,temp,s,deg,t,k;
	res:=NULL;
	if nops(Sreg)=1 then 
		deg:= degree(evala(Norm(x-Sreg[1][1], ext, ext),x),x);
		if d mod deg =0 then 
			if Sreg[1][1]= infinity then 
				res:=res,1;
			else
				res:= res, evala(Norm(x-Sreg[1][1], ext, ext),x)^(d/deg);
			fi; 
		fi;
	else
		deg := degree(evala(Norm(x-Sreg[1][1], ext, ext),x),x);
		for k from 1 to iquo(d,deg) do
			temp := searchknlog([seq(Sreg[i],i=2..nops(Sreg))],d-k*deg,ext);
			if temp<>NULL then
				if Sreg[1][1]<>infinity then
					for t in temp do
						res:= res,evala(Norm(x-Sreg[1][1], ext, ext),x)^k*t;
					od;
				else
					for t in tmep do
						res:=res,t;
					od;
				fi;
			fi;
		od;
			
	fi;
	[res];
	
end:
	
logcase := proc(Sirr,B,knl,ext,Dx,x)
	#option trace;
	local s,c,kn,res,T,temp,ep,acc;
	res := NULL;	
	for kn in knl do
		c:=NULL;
		for s in Sirr do
			acc := s[4];
			ep := poly_to_powerseries(kn,x,s[1],T,acc);
			ep := pow_reciprocal(ep,T);
			temp := poly_to_powerseries(B,x,s[1],T,acc);
			temp := pow_time(ep, temp, T);
			ep := poly_to_powerseries(s[3],x,s[1],T,acc);
			temp := pow_time(temp, ep,T);
			if temp[3]<>1 then 
				c:=NULL;
				break;				
			fi;
			if c=NULL then 
				c:=temp[1];
			elif c<>temp[1] then
				c:=NULL;
				break;
			fi;			
		od;
		if c<>NULL then res:=res,compute_sqrt(c*kn,x,ext);fi;
	od;
	[res];
end:

compute_sqrt := proc(f,x,ext)
	local v,p1,p2,i,w,f1;

	f1:=expand(numer(f))/expand(denom(f));
	v := evala(Sqrfree(f1,x));
	p1 := mul(i[1]^(i[2] mod 2), i=v[2]);
	p1 := factor(evala(Expand( numer(p1) ))/evala(Expand(denom(p1))));
	p2 := factor(collect(mul(i[1]^( floor(i[2]/2) ), i=v[2]), x, evala));
	# Now f = v[1] * p1 * p2^2   with v[1] constant,
	#   p1 square-free monic poly, p2 rat funct.
	if has(p1,x) then
		sqrt(v[1]*p1) * p2
	else
		#w := evala(Factors(x^2-v[1], ext))[2];
		#if nops(w)=1 then
		#	sqrt(v[1]) * p2
		#else
		#	-coeff(w[1][1],x,0)/coeff(w[1][1],x,1) * p2
		#fi
		sqrt(v[1]*p1) * p2
	fi
end: