# Solver for differential equations of order 2 with rational function coeffs
#
# Finds all solutions that can be expressed in terms of  0F1 and 1F1 functions 
#such as Bessel{I,J,K,Y}(nu, sqrt (f)) or Airy{A,B}(f)
# or Whittaker{M,W}(mu, nu, f) or Kummer{M,U}(mu, nu, f) with mu, nu
# constants and f an arbitrary rational function.
#
# Version 1.0 --- September 23, 2011
#
# Copyright: Quan Yuan, Mark van Hoeij and Ruben Debeerst.


macro(
	EXP_INT = `DEtools/kovacicsols/e_int`,
	EXP_SOLS = `DEtools/Expsols`,

	NO_SOL = "No 0F1 or 1F1 solutions",
	KUMMER = `dsolveBessel/KummerScore`
):

BesselSolver := proc()
	local ans;
	# use frontend to replace non-algebraic constants by names
	ans := traperror(frontend(BesselSolver_doit,[args],
		[{`+`,`*`,RootOf,radical,nonreal},
		indets([args[2..-1]],function)]));
	if ans=lasterror then
		if ans=NO_SOL then
			ans := 0
		else
			error ans
		fi
	fi;
	if _EnvListOutput = true and ans=0 then ans := [] fi;
	ans
end:




# Input:
#    i) a differential operator L and (optional) the domain y
#   ii) a differential equation L and the dependent variable y
# output: 0F1 or 1F1 type solution, i.e Bessel, Whittaker, or Kummer type solution if such solution exists
#         otherwise return either 0 or NO_SOL.
BesselSolver_doit := proc(L,y)
	#option trace;
	local Dx,x,i, ext,res, Sirr, Sreg, mayB, caseB, caseW, T, sol,Sa;
	i := NULL;
	if nargs>1 then
		if type(y,list(name)) and nops(y)=2 then
			_Envdiffopdomain:=y
		elif type(y,function) then
			x := op(y);
			_Envdiffopdomain:=[Dx,x];
			return procname(DEtools['de2diffop'](L,y))
		else
			error "wrong number or type of arguments"
		fi
	fi;
	if not assigned(_Envdiffopdomain) then
		error "domain is not assigned"
	fi;
	Dx,x := op(_Envdiffopdomain);
	if degree(L, Dx)<>2 then
		error "input should be of order 2"
	fi;
	ext := indets([args], {radical, nonreal});  # All algebraics not in RootOf notation
	if ext<>{} then
		i := [args];
		ext := radfield(ext);  # p[1] translates to RootOf notation, p[2] translates back.
		return   op(eval( [procname(op(eval(i, ext[1])))],  ext[2]))

	fi;
	#Convert everything into RootOf notation to avoid 
	#that sometimes Maple command doesn't work properly in radical case. 
	i:={Dx,x};
	ext := indets(L, {name,RootOf}) minus i;
	if ext = {} then
		Normalizer := normal;
	else
		Normalizer := evala;
	fi;
	Sirr, Sreg, Sa,mayB, caseB, caseW := singInfo(L,ext,T); 
	#collect information about exponent difference and possible cases from singularities.  
	userinfo(1,dsolveBessel,"Exponents difference", Sirr,Sreg);
	if mayB then		
		sol := findBessel(L,Sirr,Sreg,Sa,caseB, ext,x,T);
		if sol <> 0 then 
			return sol;			
		fi;
		
	fi;
	if caseW>0 then
		findWhittaker(L, Sirr, Sreg, Sa,ext, T, x, caseB, caseW)
	else
		error NO_SOL
	fi
		

end:



findWhittaker := proc(L, oldSirr, Sreg,Sa,ext,T,x, oldcaseB, caseW)
	#option trace;
	local vf,v,c,f,Sirr,mult,i,k,s,mulist,ans,mu, vflist, j,Srest,caseB,Sir;
	Sir, Srest,caseB:= updateSin(oldSirr,Sreg,Sa,ext,T,x, oldcaseB);
	if caseW = 3 then
		for i in Sir do
			if type(i[3],list) then
				v := coeff(i[3,1],T,0)/ldegree(i[3,1],T)/2;
				v := evala(v^2*i[3,2])
			else
				v := coeff(i[3],T,0)/ldegree(i[3],T)/2;
				v := x^2 - poly_d(v, x, ext)
			fi;
			if assigned(c) and evala(c-v)<>0 then error NO_SOL fi;
			c := v
		od;
		Sirr := SqrtConst(Sir, ext, T, x, c)[1,1]
	else
		c := 1;
		Sirr := Sir;
	fi;
	vflist := NULL;
	for i in [findBesselvf(caseB, Sirr,Srest,ext,T)] do
		if nops(i)=3 then
			# The logarithmic case
			for j in i[1..2], [abs(i[1]-1/2), i[2]] do
				if j[1]=i[3] then
					vflist := j, vflist
				else
					vflist := vflist, j
				fi
			od
		else
			vflist := vflist, i
		fi;
	od;
	for vf in [vflist] do
		v, f := op(vf);
		userinfo(2,dsolveBessel,"Trying v,c,f = ", v,c,f);
		ans := [1,-1];
		for i from 1 to nops(Sirr) while c=1 do
			mult := -ldegree(Sirr[i,3],T);
			userinfo(3,dsolveBessel,"[p,mult]",[Sirr[i,1],mult]);
			k := coeff(Sirr[i,3],T,0);
			k := [seq(k+j,j=0..mult-1)];
			mulist := k/(2*mult);
			if i=1 then
				ans := mulist
			else
				ans := [seq(`if`(member(true,
					  [seq(seq(type(2*evala(j-k*s),integer)
					   ,k=mulist),s={1,-1})]
					),j,NULL),j=ans)]
			fi;
			userinfo(4,dsolveBessel,"mu's so far:",ans);
			if ans = [] then break fi;
		od;
		if ans<>[] then f := convert(f,sqrfree,x) fi;

		userinfo(1,dsolveBessel,"remaining mu's:",ans);
		for mu in ans do
			userinfo(1,dsolveBessel,"testing mu",mu);
			s := sign(length([-mu,-2*f]) -  length([mu,2*f]));
			if s=0 then s := `if`(abs_value(mu)=mu, 1, -1) fi;
			i := whittakerequiv(L, s*mu, v, s*2*f, c);
			if i<>0 then
				return SimplifyAnswer(i)
			fi
		od;
	od;
	0


end:

# Input: list of singularities Sirr and Srest, extensions ext, name T
# Output: list of possible pairs [v,f]
findBesselvf := proc(caseB, Sirr, Srest, ext, T)
	local f, inf;

	f := besselsubst(Sirr,T,ext);
	userinfo(1,dsolveBessel,"Possibilities for f ",f);
	inf := member(infinity, map(i->i[1],Sirr));

	if caseB >1 then
		# find bessel function with v \nin \Q
		findBesselvfirrat(Srest,f,inf,ext)
	elif caseB=0 then
		# find bessel function with v=0
		findBesselvfln(Srest,f,ext)
	elif caseB=1 then
		# find bessel function with v\in\Q
		findBesselvfrat(Srest,f,ext)
	else
		# find bessel function with v\in\Q where all singularity
		# differences are integers
		findBesselvfint(f,inf,ext,Srest)
	fi;
end:

# Input: set of singularities SC with exp differences not in K,
#   set of  possible rational functions f, a boolean inf=(\inf\in\Sirr), 
#   extensions ext
# Output: list of possible pairs [v,f] with v\nin K
#   such that exp-gauge-transformations of B(v,f(x)) can be a solution
#   of the differential operator belonging to SC
findBesselvfirrat := proc(SC,f,inf,ext)
	local Dx,x,z,i,g,gg,h,hh,d,r,a,b,l,v,sol,vf,ordeRoot,hasx;
	Dx,x := op(_Envdiffopdomain);
	
	vf := NULL;
	#gen-exp differences are sqrt(.)+rational
	#make minpolys for sqrt(.)
	hasx := 0;
	for i from 1 to nops(SC) do
		d[i] := SC[i,3];
		if has(d[i],x) then hasx := hasx+1 fi;
	od;
	if hasx = 0 then
		return findBesselvfK(args)
	elif hasx < nops(SC) then
		error NO_SOL
	fi;
	userinfo(3,dsolveBessel,"minpolys for gen-exp differences",
	 	[seq(d[i],i=1..nops(SC))] );
	
	# all differences are a rational multiple of d[1]
	# compute these rational coefficients
	a := coeff(d[1],x,0);
	r[1] := 1;
	for i from 2 to nops(SC) do
		b := coeff(d[i],x,0);
		r[i] := sqrt(Normalizer(b/a));
		if not type(r[i],rational) then return NULL fi;
	od;
	userinfo(4,dsolveBessel,"ratios for gen-exp differences",
	 	[seq(r[i],i=1..nops(SC))] );
	l := ilcm( seq(denom(r[i]),i=1..nops(SC)) );
	h := mul(SC[i,4]^(r[i]*l),i=1..nops(SC));
	userinfo(2,dsolveBessel,"numerator must be a power of",h);
	
	# test all possibilities of f (try +/- each polar part)
for i from 1 to 2^(nops(f)-1) do
	g := Normalizer(possibility(f,i));
	userinfo(3,dsolveBessel,"Testing possibility", g);
	# get c from first singularity in SC 
	g := changeconstant(g,x,SC[1,1],ext);
	if g=NULL then next; fi;
	# test if other singularities in SC are zeros
	if not testzeros(g,x,map(j->j[1],SC)) then next; fi;

	gg:=numer(g);
	if degree(h,x)=0 then
		a:=1;
	else
		a:=degree(gg,x)/degree(h,x);
	fi;
	# now h^a=g
	if not type(a,integer) then
		next
	fi;
	gg:=collect(gg/lcoeff(gg,x),x,Normalizer);
	hh:=collect(h/lcoeff(h,x),x,Normalizer);
	if Normalizer(gg-hh^a)=0 then
		if SC[1,1]=infinity then
			ordeRoot := degree(denom(g),x)-degree(gg,x)
		else
			# h^a contains i'th factor a*r[i]*l times, and r[1]=1
			ordeRoot := a*l
		fi;
		v := sqrt(factor(-coeff(d[1],x,0)), 'symbolic');
		v := v / (2*ordeRoot);
		vf := vf, [v,g];
	fi;
od;
	vf
end:

# Input: set of singularities SC with non-rational exp differences in K,
#   set of  possible rational functions f, a boolean inf=(\inf\in\Sirr), 
#   extensions ext
# Output: list of possible pairs [v,f]
#   such that exp-gauge-transformations of B(v,f(x)) can be a solution
#   of the differential operator belonging to SC
findBesselvfK := proc(SC, f, inf, ext)
	local Dx,x,z,i,g,gg,h,hh,r0,r1,a,l,v,vv,k,sol,vf,ordeRoot,flag;
	Dx,x := op(_Envdiffopdomain);
	
	vf := NULL;	
	# all differences are a rational multiple of d[1]
	# compute these rational coefficients
	r0[1] := 0;
	r1[1] := 1;
	for i from 2 to nops(SC) do
		r0[i],r1[i] := compare(SC[i,3],SC[1,3]);
	od;
	userinfo(4,dsolveBessel,"ratios for exp differences",
	 	[seq(r1[i],i=1..nops(SC))] );
	l := ilcm( seq(denom(r1[i]),i=1..nops(SC)) );
	h := mul(SC[i,4]^(r1[i]*l),i=1..nops(SC));
	userinfo(2,dsolveBessel,"numerator must be a power of",h);
	
	# test all possibilities of f
	for i from 1 to 2^(nops(f)-1) do
		g := Normalizer(possibility(f,i));
		userinfo(3,dsolveBessel,"Testing possibility", g);
		# get c from first singularity in SC 
		g := changeconstant(g,x,SC[1,1],ext);
		if g=NULL then next; fi;
		# test if other singularities in SC are zeros
		if not testzeros(g,x,map(j->j[1],SC)) then next; fi;

		gg:=numer(g);
		if degree(h,x)=0 then
			a:=1;
		else
			a:=degree(gg,x)/degree(h,x);
		fi;
		if not type(a,integer) then next fi;
		gg:=collect(gg/lcoeff(gg,x),x,Normalizer);
		hh:=collect(h/lcoeff(h,x),x,Normalizer);
		if Normalizer(gg-hh^a) <> 0 then next fi;
		# now h^a=gg
		
		if SC[1,1]=infinity then
			ordeRoot := degree(denom(g),x)-degree(gg,x)
		else
			# h^a contains the i'th factor a*l times
			ordeRoot := a*l
		fi;
		v:=SC[1,3]/(2*ordeRoot);
		vv := seq(v+k/(2*ordeRoot),k=0..(2*ordeRoot-1));
		userinfo(2,dsolveBessel,"set for v", vv);
		
		# test v with other zeros
		for v in vv do
			flag := true;
			for k from 2 to nops(SC) do
				ordeRoot := `if`(SC[k,1]=infinity,
					degree(denom(g),x) - degree(gg,x),
					r1[k]*a*l);
				if not type(Normalizer(v*2*ordeRoot-SC[k,3])
				,integer) then
					flag := false; break
				fi;
			od;
			if flag then vf := vf, [v,g] fi;
		od;
	od;
	vf
end:

# Output: a new variable
newvar := proc()
	local t;
	t
end:

# Input: A,B in k
# Output: r0,r1 such that A=r0+r1*B
compare := proc(A,B)
	local iszero,E,sol;
	iszero := A- (r0+r1*B);
	iszero := evala(Expand(numer(Normalizer(iszero))));
	iszero := convert(iszero,RootOf);
	do
		E := indets(iszero,RootOf);
		if E = {} then break fi;
		iszero := subs(E[1]=newvar(),iszero);
	od;
	sol := SolveTools:-Linear({coeffs(expand(iszero),
		indets(iszero,{name,symbol}) minus {r0,r1})},{r0,r1});
	if not type(sol,set) or not type(map(rhs,sol),set(rational)) then
		error NO_SOL
	fi;
	op(subs(sol, [r0,r1]))
end:


# Input: set of singularities Sln with exp differences in Z and log in 
#   formal solution, set of  possible rational functions f, a boolean 
#   inf=(\inf\in\Sirr), extensions ext
# Output: list of possible pairs [v,f] with v=0
#   such that exp-gauge-transformations of B(v,f(x)) can be a solution
#   of the differential operator belonging to SC
findBesselvfln := proc(Sln,f,ext)
	local Dx,x,i,j,g,sol,v,vf;
	Dx,x := op(_Envdiffopdomain);
	vf := NULL;
	# test all possibilities of f
	for i from 1 to 2^(nops(f)-1) do
		g := Normalizer(possibility(f,i));
		userinfo(3,dsolveBessel,"Testing possibility", g);
		# compute constant
		g := changeconstant(g,x,Sln[1,1],ext);
		if g=NULL then next; fi;
		# test for other zeros
		if not testzeros(g,x,map(j->j[1],Sln)) then next; fi;
		v := 2*max(degree(numer(g),x),degree(denom(g),x));
		v := add(j[3]*j[5], j=Sln)/v;
		vf := vf, [ceil(v),g,v]
	od;
	vf
end:

# Input: set of singularities SQ with exp differences in Q,
#   set of  possible rational functions f, extensions ext
# Output: list of possible pairs [v,f] with v in Q
#   such that exp-gauge-transformations of B(v,f(x)) can be a solution
#   of the differential operator belonging to SC
findBesselvfrat := proc(SQ,f,ext)
	local Dx,x,i,j,k,g,h,fctrs,mult,v,vv,c,sol, newh,vf;
	Dx,x := op(_Envdiffopdomain);
	vf := NULL;
	for i from 1 to 2^(nops(f)-1) do
		g := possibility(f,i);
		userinfo(3,dsolveBessel,"Testing possibility", g);
		# get constant from first singularity
		g := changeconstant(g,x,SQ[1,1],ext);
		if g=NULL then next; fi;
		if not testzeros(g,x,map(j->j[1],SQ)) then next; fi;
		
		h := numer(g);
		for j from 1 to nops(SQ) do
			# look for multiplicities of the zeros
			mult[j] := 0;
			if SQ[j,1] = infinity then
				mult[j] := degree(denom(g))-degree(numer(g));
			else
				mult[j] := 0;
				while evala(Rem(h, SQ[j,4], x, 'newh')) = 0 do
					mult[j] := mult[j] + 1;
					h := newh
				od;
				# mult[j] \neq 0 since we already know
				# that this is a zero
			fi;
			# possibilities for v
			v := [seq(op([SQ[j,3]-i,SQ[j,3]+i]),i=0..mult[j])]
				/ (2*mult[j]);
			if j=1 then
				vv := [seq(`if`(MatchInt([0,1/2,op(v[1..i-1])],
					v[i]), NULL, abs(v[i])), i=1..nops(v))]
			else
				vv := [seq(`if`(MatchInt(v,i),i,NULL),i=vv)]
			fi;
			if nops(vv) = 0 then break; fi;
		od;
		userinfo(3,dsolveBessel,"multiplicities of factors",
			[seq(mult[j],j=1..nops(SQ))]);
		userinfo(2,dsolveBessel,
			"possibilities for v", vv);
		# if we get here, then we have a function g and a parameter v
		# where for each singularity [x,t,d] in SQ there is a
		# corresponding factor t in g, and it's multiplicity m is
		# such that d=2*v*m. The rest-polynomial h has now to be a
		# p-power such that 2*v*p \in \Z
		for v in vv do
			# rest must be a p-power
			c:=lcoeff(h,x);
			newh:=ispower(h/c,x,denom(2*v));
			if Normalizer(h/c-newh^denom(2*v)) = 0 then
				vf := vf, [v,g];
			fi;
		od;		
	od;
	vf
end:


# Input: set of  possible rational functions f, a boolean inf=(\inf\in\Sirr), 
#   extensions ext
# Output: list of possible pairs [v,f] with v in Q
#   such that exp-gauge-transformations of B(v,f(x)) can be a solution
#   of a differential operator with no half-regular singularities
findBesselvfint := proc(f,inf,ext, Sr)
	local Dx,x,ff,n,p,i,g,lc,a,k,v,sol,C,c,cc,fctrs,vf,s;
	Dx,x := op(_Envdiffopdomain);
	vf:= NULL;
	for i from 1 to 2^(nops(f)-1) do
		ff := Normalizer(possibility(f,i)+C);
		userinfo(3,dsolveBessel,"Testing possibility", ff);
		if inf then
			n := degree(numer(ff));
		else
			n := degree(denom(ff));
		fi;
		for p from 2 to n do if irem(n,p)=0 then
			userinfo(3,dsolveBessel,
			  "Test whether numerator is a p-power", p);
			# if inf then test p-power and get c
			# if not then get c from p power or c is zero
			g := numer(ff);
			lc := lcoeff(g,x);
			a := ispower(collect(g/lc,x,Normalizer),x,p);
			userinfo(5,dsolveBessel, "Is a power of", a);
			s := numer(Normalizer(evala(Content(g-lc*a^p,x))));
			s := evala(Factors(s,ext));
			cc := {seq(`if`(degree(k[1],C)=1, 
			-coeff(k[1],C,0)/coeff(k[1],C,1),NULL),k=s[2])};

			if not inf then
				# test c=0
				fctrs := sqrfree(subs(C=0,g),x);
				if add(modp(k[2],p), k=fctrs[2])=0 then
					# the multiplicity of each factor of g
					# is a multiple of p
					cc:=cc union {0};
				fi;
			fi;

			userinfo(2,dsolveBessel, "possibilities for c", cc);
			for c in cc do
				g := subs(C=c,ff);
				for k from 1 to p-1 do
					v := k/(2*p);
					if denom(2*v)<p then next fi;
					if Sr<>[] and UpdateV(v,g,x,Sr,'v') then
						vf := [v,g], vf
					else
						vf := vf, [v,g]
					fi;
				od;
			od;
		fi od;
	od;
	vf
end:

# Update v in such a way to make sure that no gauge transformation
# is used when none is needed. This procedure is not needed for
# correctness or completeness, but only to make the output look nicer.
UpdateV := proc(v,g,x,Sr,newv)
	local s,m,mg,h,N,n1,n2;
	m := {op(map(s -> s[3], Sr))};
	m := min(op(m));
	for s in Sr while s[3]>m do od;
	if s[1]=infinity then
		mg := degree(denom(g),x)-degree(numer(g),x)
	else
		mg := 0;  h := numer(g);
		while evala(Rem(h, s[4], x, 'h'))=0 do mg := mg+1 od
	fi;
	if mg=0 then return false fi;
	n1 := solve( (v+N)*2*mg = m, N);
	if type(n1,integer) then
		newv := v+n1;
		return true
	fi;
	n2 := solve( (N-v)*2*mg = m, N);
	if type(n2,integer) then
		newv := n2-v;
		return true
	fi;
	if n1>1 then
		if abs(round(n1)-n1) <= abs(round(n2)-n2) then
			newv := v+round(n1)
		else
			newv := round(n2)-v
		fi
	fi;
	false
end:

# Input: set of singularities Sirr of a differential operator L with 
#   non-constant generalized exponent
# Output: the set {f1,...fn} such that if L has a Bessel solution, then it is 
#   of the form r1*E(x)*B(v,f(x))+r2*(E(x)*B(v,f(x)))', where 
#   f=(-1)^k1*f1+...+(-1)^kn*fn+C, r1,r2\in\C(x), ki\in\Z and E(x) is a
#   solution of a first order differential equation
besselsubst := proc(Sirr,T,ext)
	local i,d,f,ff,j;
	f := NULL;
	for i from 1 to nops(Sirr) do
		d := 1/2 * Sirr[i,3];
		ff := add( coeff(d,T,j)/j*Sirr[i,2]^j, j=ldegree(d,T) .. -1);
		f := f, evala(Trace(ff,ext,ext))
	od;
	[f]
end:

# Input: a list f=[f1,...fn] and an integer 1<=k<=2^n
# Output: the sum +-f1+-...+-fn, where the k-th possibility
#   of signs is choosen
possibility := proc(f,k)
	local b,n,a,i;
	n := nops(f);
	if k>2^n then error "there are just 2^n possibilities" fi;
	b := k-1;
	for i from n by -1 to 1 do
		b := iquo(b,2,a[i])
	od;
	add((-1)^a[i]*f[i],i=1..n)
end:

# Input: a function f(x) and a point p
# Output: f(x)+c such that f(p)+c=0,
#    if p=\infty then f(x) is returned
changeconstant := proc(f,x,p,ext)
	local c,pp,ff,C;
	if p=infinity then
		return f;
	fi;
	if type(p, RootOf) then
		pp:=evala(Norm(x-p,ext,ext));
		ff:=numer(Normalizer(f+C));
		c:=SolveTools:-Linear({coeffs(evala(Rem(ff,pp,x)),x)},{C});
		if c=NULL then
			NULL
		else
			Normalizer(f+subs(c[1],C))
		fi;
	else
		Normalizer(f-subs(x=p,f))
	fi;
end:

# Input: a function f(x) and a set of points pp
# Output: true, if all points in pp are a zero of the numerator
#   and false otherwise
#   for \infty the condition changes to
#   degree(denom(f))>degree(numer(f))
testzeros := proc(f,x,pp)
	local p,g;
	g:=numer(f);
	for p in pp do
		if p=infinity then
			if not degree(g,x)<degree(denom(f),x) then
				return false;
			fi;
		else
			if not evala(subs(x=p,g))=0 then
				return false;
			fi;
		fi;
	od;
	true
end:






whittakerequiv := proc(L,mu,v,f,c)
	local Dx,x,eq,t,L2,Y,g,wc, s,w;
	global KUMMER;

	Dx,x := op(_Envdiffopdomain);
	if c=1 then
		# See if the answer is shorter in terms of Kummer functions:
		g := seq(seq([1/2-s*mu+w, 1+2*w,s*f] ,w=[v,-v]), s=[1,-1]);
		g := sort([g],length)[1];
		if length(g) < length([mu,v,f]) + 5 * KUMMER then
			return Kummerequiv(L,op(g),x)
		fi
	fi;

	# Whittaker equation after subs:  (x,mu) -> sqrt(c) * (x,mu)
	eq :=diff(Y(x),x,x)-(4*evala(v^2)+c*x^2-4*c*mu*x-1)/4/x^2*Y(x);
	eq := PDEtools['dchange']({x=subs(x=t,f)},eq,[t],
		params=(indets([args],{name,symbol}) minus {x}) );
	L2 := DEtools['de2diffop'](eq,Y(t)):
	L2 := s_equiv(L2,L);
	if L2 <> 0 then
		wc := sqrt(c, 'symbolic');
		t:=1/sqrt(1/c,'symbolic');if length(t)<length(wc) then wc:=t fi;
		L2, [WhittakerM(wc*mu, v, wc*f), WhittakerW(wc*mu, v, wc*f)]
	else
		0
	fi;
end:



# Input: The list of irregular singularities Sirr.
# Let K denote the field given by ext, and let T be the name used in 
# gen_exp, and x the independent variable.
#
# Each i in Sirr is a list with
#  i[1] = point p in C union {infinity}
#  i[4] = minpoly of i[1] over K, whenever i[1] is not infinity.
#  i[5] = degree of this minpoly
#  i[2] = local parameter (not used in this procedure)
#  i[3] = Either the gen.exp-difference D in K(p)[1/T],
#         or a list [E,D] such that gen.exp-difference is E * sqrt(D)
#         and where E in K(p)[1/T] and D in K(p).
#
# Output: a set of pairs [Updated(Sirr), c].
# Here an updated Sirr is one in which i[3] is never a list [E,D], but
# always an element D in K(p)[1/T], and where the gen.exp-difference i[3]
# (which was E*sqrt(D) or D depending on whether i[3] was a list or not)
# has been divided by sqrt(c) for all i in Sirr.
# This means that if we had an i[3] = [E,D] in the input, then for 
# this D it must have been true that D/c is a square in K(p). And if we
# had an i[3] = D in the input, then for this i, we must have had sqrt(c)
# in K(p) where p denotes i[1].
#
# In most cases the c can be read off from one of the i[3]'s. In the
# remaining case we'll have to find candidate c's by computing Subfields.
#
# The c's coming out of this procedure are used as follows: When we compute
# the f for which the solution is in terms of BesselI(v, f) then this f
# should still be multiplied by sqrt(c). If this c is negative
# then replace c by -c and replace BesselI by BesselJ.
#
SqrtConst := proc(Sirr, ext, T, x)
	local Sodd, i, nl, ol, S, disc, ans, c, flag, res, extp, pol, wc, i3;

	# The [Updated(Sirr), c]'s found so far:
	ans := NULL;

if nargs=5 then
	# c is already given in the optional input args[5]
	S := {args[5]}
else	
	# Singularities for which K(p) has odd degree over K:
	Sodd := [seq(`if`(irem(i[5],2)=1, i, NULL), i=Sirr)];

	# The number of i's for which i[3] is a list:
	nl := add(`if`(type(i[3],list),1,0),i=Sirr);
	ol := add(`if`(type(i[3],list),1,0),i=Sodd);

	# c can only be 1 if no i[3] is a list (i.e. if nl=0).
	if nl=0 then
		ans := [Sirr, 1]
	fi;

	# If any odd sing does not have a list, then no cases other than
	# perhaps c = 1 need to be considered:
	if ol < nops(Sodd) then
		return {ans}
	fi;
	# At this point i[3] must be a list for all i in Sodd.

	# Find all singularities where c is easy to read from:
	S := NULL;
	for i in Sirr do
		if not type(i[3],list) then
			# Such i is not in Sodd, so the degree is even. If
			# the degree is 2, then c is easy to read.
			if i[5]=2 then
				S := S, Normalizer(discrim(i[4],x))
			fi;
			next
		fi;
		disc := i[3][2];
		if indets(disc,RootOf) minus ext <> {} then
			if member(i, Sodd) then
				disc := evala(Norm(disc, ext, ext))
			else
				next
			fi
		fi;
		S := S, Normalizer(disc)
	od;
	
	#if S = NULL then
		# We found no points where c was easy to read, this
		# means we'll have to resort to Subfields to get the
		# potential values of c.
		#S := seq(`if`(type(i[3],list),
		# sqrfree(evala(Norm(x^2-i[3][2],ext,ext)),x)[2,1,1], i[4])
		#,i=Sirr);
		S := seq(`if`(type(i[3],list),
		subs(_Z=x,op(1, lhs(
		evala(Primfield({`if`(degree(i[4],x)>0, RootOf(i[4],x), NULL), RootOf(sqrfree(evala(Norm(x^2-i[3][2],ext,ext)),x)[2,1,1],x) }))[1,1]
		)))
		, i[4])
		,i=Sirr);
		
		S := evala(Subfields({S}, 2, ext, x));
		S := {seq(evala(Norm(x-i,ext,ext)), i=S)};
		userinfo(2,dsolveBessel,"subfields",S);
	#else
	#	# Found at least some point where c could be read, so
	#	# just pick the first one and try to simplify it a bit:
	#	S := sort([S],length)[1];
	#	S := {PolynomialTools:-Shorten(x^2 - S,x)}
	#fi;
	# The set of all possible c's other than c=1.
	S := {seq(Normalizer(discrim(i,x)/4), i=S)};
	userinfo(2,dsolveBessel,"possible c's other than c=1 in f=sqrt(c)*g",S);
fi;
	for c in S do
		# adjust Sirr if possible
		flag := true;
		res := NULL;
		for i in Sirr do
			extp := ext union indets(i[1],RootOf);
			pol := `if`(type(i[3],list), i[3,2], 1)*x^2 - c;
			pol := evala(Primpart(pol,x));
			pol := evala(Factors(pol, extp))[2];
			if nops(pol)=1 then
				flag := false;  break
			fi;
			pol := collect(pol[1,1],x);
			wc := coeff(pol,x,0)/coeff(pol,x,1);
			i3 := `if`(type(i[3],list), i[3,1], i[3]) / wc;
			res := res, subsop(3 = collect(i3,T,evala), i)
		od;
		if flag then
			ans := ans, [[res], c]
		fi
	od;
	{ans}
end:


updateSin := proc (oldSirr,Sreg,Sa,ext,T,x,oldCaseB)
	#option trace;
	local Sirr, Srest,Sr, s,extp, R,caseB,d,c,R1;
	Sirr :=NULL;
	Srest := NULL;
	Sr := Sreg;
	caseB:=oldCaseB;
	if Sr = [] then 
		Sr := Sa;
		caseB:= -1	
	 fi;
	for s in oldSirr do
		extp := indets(s[1],RootOf) union ext;
		R := indets(s[3],RootOf) minus extp;
		if nops(R) >0 then			
			if nops(R)<>1 then
				error "should not happen"
			fi;
			R := R[1];
			d := evala(discrim( evala(Norm(x-R,extp,extp)), x));
			c := evala(Factors(x^2-d,{R} union extp))[2][1][1];
			c := coeff(c,x,0);
			c := s[3]/c;
			s[3] := [collect(c,T,evala),d];	;
		fi;
		Sirr := Sirr,s;			
	od;
	for s in Sr do
		extp := indets(s[1],RootOf) union ext;
		R := indets(s[3],RootOf) minus extp;
		R1 := indets(s[3],RootOf) minus ext;
		if nops(R) >0 then
			s[3] := x^2-s[3];
		elif nops(R1)>0 then
			s[3] := poly_d(s[3], x, ext)
		else
			s[3] := abs_value(s[3])
		fi;
		Srest := Srest,s;

	od;
	[Sirr],[Srest], caseB
end:

Kummerequiv := proc(L,mu,nu,f,x)
	local eq,t,L2,Y;

	# Kummer equation
	eq := x * diff(Y(x),x,x) + (nu-x)*diff(Y(x),x) - mu*Y(x);
	eq := PDEtools['dchange']({x=subs(x=t,f)},eq,[t],
		params=(indets([args],{name,symbol}) minus {x}) );
	L2 := DEtools['de2diffop'](eq,Y(t)):
	L2 := s_equiv(L2,L);
	if L2 <> 0 then
		L2, [KummerM(mu, nu, f), KummerU(mu, nu, f)]
	else
		0
	fi
end:


# Minpoly of c should be an integer shift of x^2-const, if so, return x^2-const
poly_d := proc(c, x, ext)
	local t,d;
	d := sqrfree(evala(Norm(x-c,ext,ext)),x)[2,1,1];
	d := collect(d/lcoeff(d,x),x,Normalizer);
	t := coeff(d,x,1)/2;
	if degree(d,x)<>2 or not type(t,integer) then
		error NO_SOL
	fi;
	collect(subs(x=x-t,d),x,Normalizer)
end:

# Return -a if a looks negative, otherwise return a.
abs_value := proc(a)
	local i;
	i := traperror(sign(a));
	if i=lasterror then
		i := traperror(sign(Re(evalf[10](a))))
	fi;
	if i=-1 then -a else a fi
end:


# Input: differential operator L, extensions ext, name T
# Output:
#1) a list of singularities S (seperate by regular singularities and irregular singularities) such that
#   for each singularity s of L at p, S contains [p,t,d,pol,deg], where
#   - t is local parameter, ie x-p or 1/x
#   - d is exponent difference in terms of T, if D is exponent differnce
#	-- d = D for unramified case
#       -- d = [D, aT^2=b]
#   - pol is polynomial over K with zero p
#   - deg = degree(pol)
# 2) mayB, a boolean value indicate whether Bessel case is possible.
# 3) caseB, possible Bessel case list:
#	-0: easy case;
#	-1: lograthemic case;
#	-2: rational case;
#	-3: irrational case;
# 4) caseW,
# 	-0: no Whittaker solutions
# 	-1: when mu is rational
#	-2: when mu is irrational but in K
# 	-3: when mu is not in K	
  
singInfo := proc(L,ext,T)
	#option trace;
	local Dx,x,p,S,Sirr,Sreg,mayB,mayK,caseB,caseW, muD, A, q, g, t, d,c, R, extp, ae,Sa;
	global KUMMER; # Used to see if we should use Kummer functions.

	mayB := true; # true means L may have Bessel solutions
	Sirr := NULL;
	Sreg := NULL;
	Sa := NULL;
	caseB, caseW, muD, KUMMER := {}, {}, 0, 0;
	Dx,x := op(_Envdiffopdomain);

	A := evala(Primpart(L,Dx)); A := [seq(coeff(A,Dx,i), i=0..2)];
	p := A[3];
	userinfo(3,dsolveBessel,"singularities are zeros of",p);
	

	# Each finite non-apparant singularity is a root of this p:
	#
	# p := evala(Gcd(p, denom(Normalizer(
	#	DEtools[LCLM](L,Dx-rand(10..100)()) )) ));
	#
	# The following line computes such p much faster:
	#
	p := seq(`if`(A[2]=0 or q[2]>1, q[1], Normalizer(q[1] /
	   gcd( gcd(q[1], A[2] + diff(A[3],x)),
		diff(A[1],x)*A[2]-A[1]*diff(A[2],x)-A[1]^2)
	  )), q=sqrfree(A[3],x)[2]);
	p := mul(q, q=[p]);
	S := [infinity, seq(i[1],i=evala(Factors(p,ext))[2])];
	

	# Every non-apparant singularity is in S
	for q in S do
		R := {};
		if q=infinity then 
			p := infinity
		else
			if not has(q,x) then next fi;
			p := RootOf(q,x)
		fi;
		userinfo(3,dsolveBessel,"processing singularity",p);
		extp := indets(p,RootOf) union ext;
		g := DEtools[gen_exp](L,[Dx,x],T,x=p);
		#
		# The above line works OK, however, kovacicsols has
		# a gen_exp implementation optimized for order 2.
		# So we'll use that instead:
		#
		# g := [`DEtools/kovacicsols/gen_exp`(A,T,x,p,extp)];
		userinfo(3,dsolveBessel,"point and gen_exp",[p,g]);
		t := `if`(p=infinity, 1/x, x-p); # t = local parameter at p
		if nops(g)=1 and nops(g[1]) = 2 then
			# gen_exp involves an algebraic extension
			if lhs(g[1,-1])<>T then				 
				d := add(coeff(g[1,1],T,i) * 2 * mods(i,2) * T^i, i=ldegree(g[1,1],T)..0);
				d := [d, g[1,2]]
			else
				d := evala(Trace(g[1, 1], ext, ext), x)-2*g[1, 1];
				if has(d,t) then
					ae := true;
					R := indets(d,RootOf) minus extp;
					if nops(R)>1 then
						error "should not happen";
					fi;
				else
					ae := false;
				fi;
								
			fi;
		else
			d := g[1,1] - `if`(nops(g)=2, g[2,1], g[1,2]);
			ae := false;
		fi;
		d := collect(d, T, evala);		
		p := [p,t,d,op(`if`(q=infinity, [1,1], [q,degree(q,x)]))];
		userinfo(2,dsolveBessel,"gen_exp-array", p);
		if has(d,T) then
			if type(d,list) then
				caseW := {0,1};
				#squareroot case, no Whittaker solutions.			
				
			elif ae  then #has algebraic extenion but not ramified
				t := coeff(d,T,0);
				if t=0 then
					# No Whittaker solutions:
					caseW := {0,1}
				else
					mayB := false;
					# 0: no Whittaker solutions
					# 1: when mu is rational
					# 2: when mu is irrational but in K
					# 3: when mu is not in K
					caseW := caseW union {3}
				fi

			else		
				t := coeff(d,T,0);				
				mayB := mayB and type(t,integer);
				if type(t,rational) then
					muD := igcd(muD, denom(t)*ldegree(d,T));
					# denom(2*mu) divides muD, so if muD=1
					# then 2*mu\in Z, which is Bessel case
					caseW := caseW union {1,`if`(muD=1,0,1)}
				elif indets(t,RootOf) minus ext = {} then
					caseW := caseW union {2}
				else
					caseW := caseW union {3}
				fi;
				# if generalized exponents look like Kummer
				# after change of variable, then increase the
				# Kummer score; if it looks like whittaker
				# then decrease the Kummer score:
				c := {seq(min(0,ldegree(t[1],T)),t=g)};
				KUMMER := KUMMER + p[5]*
				 `if`(nops(c)=1, op(c), abs(c[1]-c[2]))
			fi;
			if not mayB then
				userinfo(1,dsolveBessel,
				"There are no Bessel solutions")
			fi;
			if nops(caseW)>1 then
				userinfo(1,dsolveBessel,
				"There are no Whittaker solutions");
			fi;
			Sirr := Sirr, p
		elif R <> {} then
			caseB := caseB union {3};
			Sreg := Sreg, p;
		elif not type(d,rational) then
			caseB := caseB union {2};
			Sreg := Sreg, p;
		elif not type(d,integer) then
			caseB := caseB union {1};
			Sreg := Sreg, p;
		elif d=0 or ( DEtools['formal_sol'](L,`has logarithm?`,x=p[1]) ) then
			# d is an integer and has a logarithm iff
			# there is a gauge operation such that d=0
			caseB :=caseB union {0};
			Sreg := Sreg,p;
		elif abs(d)>2 then
			# Apparant singularities with large d
			Sa := Sa,p				
		fi;			
		if nops(caseB)>1 or (nops(caseW)>1 and not mayB) then
			error NO_SOL
		fi;
	od;
	if Sirr = NULL then
		error NO_SOL
	elif caseB={} then
		caseB:={1};
		#Sreg is empty set, it is rational case. 
	fi;
	[Sirr], [Sreg], [Sa], mayB, op(caseB), `if`(nops(caseW)>1, 0, op(caseW))
end:

#input: the operataor L and information of exponents difference
#output: the possible list of pairs (f,nu), f represents the 
#	change of variable and nu is Bessel parameter. 

findBessel := proc(L,oldSirr,oldSreg,Sa, caseB, oldext,x,T)
	#option trace;
	local Sirr, Sreg, easy, dA, B, vflist, vf,c,ai,needext,extp,s,temp, uf, ext;
	Sirr, Sreg,B,dA, easy := singSeries(oldSirr,oldSreg, ext,x,T);
	#find the new information we need to find the possible pairs
	ext := oldext;
	needext := {op(map(s -> s[5], Sirr))};
	needext := min(op(needext));
	if needext>1 and (not easy) and caseB=1 then
		_Env_original_field := ext;
		
		for s in Sirr while s[5]>needext do od;
		extp := indets(s[1],RootOf) union ext;
		ext := extp;
		temp := singInfo(L,extp,T); # return Sirr, Sreg, Sa,mayB, caseB, caseW
		temp := singSeries(temp[1], temp[2],extp,x, T); #return Sirr, Sreg,B,dA, easy 
		Sirr := temp[1]; Sreg:=temp[2];
	fi;
	vflist := NULL;
	
	if easy then
		vflist:= sqrtEasy(Sirr, Sreg,caseB,B,dA,ext,x);
		if (indets(vflist, {name}) minus {x}<>{}) then vflist :=NULL; fi;
	fi;
	if (vflist=NULL) then
		if caseB=0 then
			vflist := sqrtLog(Sirr, Sreg,B,dA,ext,x);
		elif caseB=1 then
			vflist := sqrtRat(L, Sirr, Sreg, Sa,B,dA,ext,x,T);
		elif caseB>1 then
			vflist := sqrtIrrat(Sirr, Sreg,B,dA,ext,x);
		fi;
	fi;
	for vf in vflist do
		userinfo(2,dsolveBessel,"testing",vf);
		
		ai := testAiry(vf[2],vf[1],ext);
		if ai[1] then 
			B := Airyequiv(L,ai[2]);
			if B<>0 then
				 return SimplifyAnswer(B) 
			fi;
		fi;
		c := sign_value(lcoeff(numer(vf[2]),x));
		B:=BesselSqrtequiv(L,oldext,c,op(vf));
		if B<>0 then
			return SimplifyAnswer(B)
		fi			
	od;
	return 0;
	
end:

# Input: differential operator L, a rational function f and a parameter v
# Output: M, y
#   y is the general bessel solution B(v,f)
#   M is a differential operator such that M(y) is a solution of L
#   if such a solution is not found by the equiv-program, then 0 is 
#   returned
BesselSqrtequiv := proc(L,ext, cpos,v,f)
	#option trace;
	local Dx,x,eq,t,L2,Y,g,modifiedbessel;
	Dx,x := op(_Envdiffopdomain);
	modifiedbessel := evalb(cpos = 1);
	# The BesselI equation after change of variables x --> sqrt(x) 
	eq := (-(1/4)*x-(1/4)*v^2)*Y(x)+x*(diff(Y(x), x))+x^2*(diff(diff(Y(x), x), x));
	eq := PDEtools['dchange']({x=subs(x=t,f)},eq,[t],
		params=(indets([args],{name,symbol}) minus {x}));
	L2 := DEtools['de2diffop'](eq,Y(t)):
	L2 := s_equiv(L2,L);
	#if L2 <> 0 then
	t:=0;
	if type(v,rational) and nargs=5 and L2[2]<>1 and denom(v)>2 then
		g := length(L2[2]);
		t := procname(L, ext,cpos,abs(v+1-2*frac(v)),  f, g);
	fi;
	if L2=0 then
		if t<>0 then
			return t;
		else
			return 0;
		fi;
	else
	
		if length(t[2])<g and t<>0 then
				userinfo(2,dsolveBessel,"found smaller answer");
				return t;
		fi;
			
	
		g := compute_sqrt(collect(cpos*f,x,evala),x,ext);
		if modifiedbessel then
			L2, [BesselI(v,g),BesselK(v,g)]
		else
			L2, [BesselJ(v,g),BesselY(v,g)]
		fi;
	#else
	#	
	fi;
end:

sign_value := proc(a)
	local i;
	i := traperror(sign(a));
	if i=lasterror then
		i := traperror(sign(Re(evalf[10](a))))
	fi;
	i	
end:




compute_sqrt := proc(f,x,ext)
	local v,p1,p2,i,w,f1;
	#option trace;
	v := evala(Normal(f,expanded));
	v := evala(Sqrfree(v,x));
	p1 := 1;
	p2 := 1;
	for i in v[2] do
		if i[2] mod 2 =0 then
			p1 := p1*evala(Factor(i[1],ext))^(i[2]/2);
		else			
			p2 := p2*evala(Factor(i[1],ext))^i[2];
		fi;

	od;
	p2 := p2^(1/2);
	if has(p1,x) then
		sqrt(v[1])*p1*p2;
	else
		if type(v[1],`rational`) then return sqrt(v[1])*p1 * p2 fi;
		try
			w:=evala(Factors(x^2-v[1],{op(ext)}));
		catch :
			w:=evala(Factors(numer(normal(x^2-v[1])),{op(ext)}))
		end try;
		w:=select(has,w[2],x);
		if nops(w)=1 then
			sqrt(v[1]) * p1*p2
		else
			-coeff(w[1][1],x,0)/coeff(w[1][1],x,1) *p1* p2
		fi
	fi;
end:

testAiry:=proc(f,v,ext)
	#option trace;
	local fas,fa,T, g;
	if (not type(v,rational)) or denom(v)<> 3  then 
		return [false,0]; 
	fi;
	fas := evala(Factors(numer(f), ext))[2];
	g:=1;
	for fa in fas do
		if irem(fa[2],3) <> 0 then 
			return [false,0];
		fi;
		g:=g*fa[1]^(fa[2]/3);
	od;
	fas := evala(Factors(denom(f), ext))[2];
	for fa in fas do
		if irem(fa[2],3) <> 0 then 
			return [false,0];
		fi;
		g:=g/(fa[1]^(fa[2]/3));
	od;
	[true, g];	
end:

Airyequiv := proc(L,f)
	local eq,t,L2,Y;

	# Airy equation
	eq := diff(Y(x),x,x) -x*Y(x);
	eq := PDEtools['dchange']({x=subs(x=t,f)},eq,[t],
		params=(indets([args],{name,symbol}) minus {x}) );
	L2 := DEtools['de2diffop'](eq,Y(t)):
	L2 := s_equiv(L2,L);
	if L2 <> 0 then
		L2, [AiryAi(f), AiryBi(f)]
	else
		0
	fi



end:


# Apply operator exp(Int(d,x)) * B1 on the functions in B2.
SimplifyAnswer := proc(d, B1, B2)
	local c, c2, C, Dx, x, i, de, Y, f, sol, p;
	global _C1, _C2;
	userinfo(1,dsolveBessel,"Operator and functions",exp(Int(d,x))*B1, B2);
	Dx, x := op(_Envdiffopdomain);
	de := DEtools['diffop2de'](B1,Y(x));
	f := eval(de, Y(x) = B2[1]);
	i := select(has, indets(f,function), op(0,B2[1]));

	# f := collect(primpart(f, i, 'c'), i, factor);
	#
	# Bug, when f := A( 3^(1/2)*(2*x+4) ); i := {f};
	# then primpart(f, i, 'c'); returns 1 instead of f.
	#
	# Fix:
	p := proc(f,v) local w,i,j;
		w := {seq([v[i]=w[i],w[i]=v[i],w[i]], i=1..nops(v))};
		w := seq(map(i -> i[j], w), j=1..3);
		i := primpart(eval(f,w[1]), w[3], 'j');
		eval(collect(i,w[3],factor),w[2]), j
	end:
	f, c := p(f, i);

	C := EXP_INT(Normalizer(d + diff(c,x)/c),x);
	sol := C * f;
	f := eval(de, Y(x) = B2[2]);
	i := select(has, indets(f,function), op(0,B2[2]));

	# f := collect(primpart(f, i, 'c2'), i, factor);
	# Same bug as above, fix:
	f, c2 := p(f, i);

	sol := [sol, Normalizer(c2/c) * C * f];
	if _EnvListOutput = true then
		sol
	else
		_C1 * sol[1] + _C2 * sol[2]
	fi
end:


singSeries:=proc(oldSirr, OldSreg,ext,x, T)
	#option trace;
	local Sirr, Sreg, dA, B, neq, i,p, d, j,l, acc, temp,easy;
	Sirr := NULL;
	Sreg, easy, dA, B :=  OldSreg, false, 0,1;
	#here we  also check if it is easy case. 
	neq := add(i[5], i in Sreg); # use neq to track the possible equations we can get.
	for p in oldSirr do
		if type(p[3], list) then
			d:=p[3][1];
			acc := -ldegree(d,T);
			if p[1] <> infinity then 
				B := B*p[4]^acc;
			else 
				dA := acc; 
			fi;
			acc := ceil((1+acc)/2);
			d := collect(d,T,evala);
			if type(d, `+`) then l := [op(d)] else l := [d] fi;
			d := add(i/ldegree(i, T), i = l); 
			d := collect(d^2, T, evala);
			temp := evala(rhs(p[3,2])/coeff(lhs(p[3,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); 
			if j=FAIL then error NO_SOL fi;
			d := 0; 
			for i in l do if degree(i, T) < j then d := d+i fi od; 
			d := subs(T = temp, d);
			# addjust coefficient and take squre to get series after squre.
			neq := neq + acc*p[5];	
			Sirr := Sirr, [p[1],p[2],[d,acc], p[4],p[5]]
		else
			d:=p[3];
			acc:=-ldegree(d,T);
			if p[1] <> infinity then 
				B := B*p[4]^(2*acc); 
			else 
				dA := 2*acc; 
			fi;
			d := d-coeff(d,T,0);
			d := collect(d,T,evala);
			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 = p[2], d));
			neq := neq +acc*p[5]; 
			Sirr := Sirr, [p[1],p[2],[d,acc], p[4],p[5]]
		fi;
	od;
	dA := degree(B, x)+dA;
	if neq >dA then easy:=true fi;
	# if has(Sreg, infinity) then dA := dA-1 fi; #adjust degree for infinity to reduce some computation.
	[Sirr],Sreg,B,dA, easy;
end:

sqrtEasy:=proc(Sirr, Sreg, caseB,B,dA, ext, x)
	#option trace;
	local a,f,fn,p,i,j,t,linearr,temp,res,d,d1,nonpole,und;
	
	userinfo(3,myname,`Now it is the easy case`);					  
	f := add(a[i]*x^i, i = 0 .. dA)/B;
	und := indets(f) minus {x} ;
	if has(Sreg, infinity) then 
		f:=subs(a[dA]=0,f);
	fi;
	fn := numer(f); 
	for p in Sreg do 
		fn := fn/p[4]; 
	od; 
	linearr := evala(Rem(numer(fn), denom(fn), x)); 
	linearr := coeffs(linearr, x); 
	temp := SolveTools:-Linear({linearr},und); 
	f := subs(op(temp), f);
	#each p in Sreg should be a zero. 
	if has(Sreg, infinity) then 
		f:=subs(a[dA]=0,f);
	fi;
	d := 0; 
	j := 0; # use j indicate that if infinity is irregularsigularity
	for p in Sirr do 
		
		if p[1] = infinity then 
			d := d+ p[3][1];
			j := p[3][2]; 
		else
			d := d+evala(Trace(p[3][1], ext, ext), x); 
		fi; 
	od; 
	nonpole:=1;
	for p in Sirr do
		nonpole:= nonpole*p[4]^p[3][2];
	od;
	nonpole := Normalizer((f-d)*nonpole); 
	temp := evala(Rem(numer(nonpole), denom(nonpole), x));
	linearr := coeffs(temp, x);
	#each irregular singularity represents pole part of f. 
	if j<>0 then #the equations at infinity
		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 := SolveTools:-Linear({linearr},und);
	temp := op(temp); 
	if temp <> NULL then 
		f := subs(temp, f); 
		f := op(simplef({f}));
	else 
		f := NULL;
	fi;
	if caseB =0 then
		return findnuLog(f,Sreg,x);
	elif caseB =1 then
		return findnueasyrat(f,Sreg,B,x);
	else 
		return findnueasyIrrat(f,Sreg,x);
	fi;
end:

findnueasyIrrat:= proc(f,Sreg,x)
	#option trace;
	local mult,h,res,i,s, smin, multmin;
	mult := 0;
	h := numer(f);
	res := {};
	multmin := infinity;
	for s in Sreg do
		if s[1] = infinity then
			mult := degree(denom(f),x)-degree(h,x);
		else
			while evala(Rem(h, s[4], x, 'h')) = 0 do mult := mult+1 od
		fi;
		if mult<multmin then
			smin:=s;
			multmin := mult;
		fi;
	od;
	for i from 0 to floor(multmin/2) do  
		res := res union {[(smin[3]+i)/multmin, f]}
	od;
	return res;
end:

findnueasyrat:=proc(f,Sreg,B,x)
	#option trace;
	local j,h,v, vv, vf, mult, newh, c;
	vf:={};
	h := numer(f);
	for j from 1 to nops(Sreg) do
			# look for multiplicities of the zeros
			mult[j] := 0;
			if Sreg[j,1] = infinity then
				mult[j] := degree(denom(f),x)-degree(numer(f),x);
			else
				mult[j] := 0;
				while evala(Rem(h, Sreg[j,4], x, 'newh')) = 0 do
					mult[j] := mult[j] + 1;
					h := newh
				od;
				# mult[j] \neq 0 since we already know
				# that this is a zero
			fi;
			# possibilities for nu
			v := [seq(op([Sreg[j,3]-i,Sreg[j,3]+i]),i=0..mult[j])]
				/ (mult[j]);
			if j=1 then
				vv := [seq(`if`(MatchInt([0,1/2,op(v[1..i-1])],
					v[i]), NULL, abs(v[i])), i=1..nops(v))]
			else
				vv := [seq(`if`(MatchInt(v,i),i,NULL),i=vv)]
			fi;
	od;
	# if we get here, then we have a function g and a parameter v
	# where for each singularity [x,t,d] in SQ there is a
	# corresponding factor t in g, and it's multiplicity m is
	# such that d=v*m. The rest-polynomial h has now to be a
	# p-power such that v*p \in \Z
	for v in vv do
		# rest must be a p-power
		c:=lcoeff(h,x);
		newh:=ispower(h/c,x,denom(v));
		if Normalizer(h/c-newh^denom(v)) = 0 then
			vf := vf union {[v,f]};
		fi;
	od;	
	vf

end:

MatchInt := proc(L, i)
	local j;
	for j in L do
		if type(i-j,integer) or type(i+j,integer) then
			return true
		fi
	od;
	false
end:



sqrtLog := proc(Sirr, Sreg,B,dA,ext,x)
	#option trace;
	local res,knl, kn,f, j, nu,g, d, c, p, nonpole, temp, linearr, d1;
	userinfo(3,myname,`Now it is the logarithmic case`);
	res:=NULL;
	knl := searchKnLog(Sreg,dA);
	f := {};
	for kn in knl do
		d := 0; 
		j := 0; # use j indicate that if infinity is irregularsigularity
		for p in Sirr do 
		
			if p[1] = infinity then 
				d := d+ p[3][1];
				j := p[3][2]; 
			else
				d := d+evala(Trace(p[3][1], ext, ext), x); 
			fi; 
		od; 
		nonpole:=1;
		for p in Sirr do
			nonpole:= nonpole*p[4]^p[3][2];
		od;
		nonpole := Normalizer((c*kn/B-d)*nonpole); 
		
		temp := evala(Rem(numer(nonpole), denom(nonpole), x));
		temp := collect(temp, x, Normalizer);
		linearr := coeffs(temp, x);
		#each irregular singularity represents pole part of f. 
		if j<>0 then #the equations at infinity
			temp := evala(Quo(numer(nonpole), denom(nonpole),x));
			temp := collect(temp, x, Normalizer);
			d1:=degree(temp, x);			 
			temp := seq(coeff(temp, x, i), i = d1-j+1 .. d1); 
			linearr := linearr, temp;
		fi; 
		temp := SolveTools:-Linear({linearr},{c});
		if temp <> NULL then f := f union {subs(temp, c*kn/B)} fi; 
	od;
	if f<>{} then
		f := simplef(f);
	fi;
	for g in f do
		res:= res union findnuLog(g,Sreg,x);
	od;
end:

findnuLog := proc(g,Sreg,x)
	local nu,j;
	nu := max(degree(numer(g),x),degree(denom(g),x));
	nu := ceil(add(j[3]*j[5], j=Sreg)/nu);
	{[abs(nu),g]}
end:


searchKnLog :=proc(Sreg,d)
	#option trace;
	local res,temp,s,deg,t,k;
	res:=NULL;
	if nops(Sreg)=1 then 
		if d mod Sreg[1][5] =0 then 
			if Sreg[1][1]= infinity then 
				res:=res,1;
			else
				res:= res, Sreg[1][4]^(d/Sreg[1][5]);
			fi; 
		fi;
	else
		deg := Sreg[1][5];
		for k from 1 to iquo(d,deg) do
			temp := searchKnLog([seq(Sreg[i],i=2..nops(Sreg))],d-k*deg);
			if temp<>NULL then
				if Sreg[1][1]<>infinity then
					for t in temp do
						res:= res,Sreg[1][4]^k*t;
					od;
				else
					for t in tmep do
						res:=res,t;
					od;
				fi;
			fi;
		od;
			
	fi;
	[res];
end:

# Input: A monic polynomial f\in K[x] and a positive integer p\in\N
# Output: g\in K[x], s.t.
#   if there exists a solution to y^p=f, then y=g is such a solution
ispower := proc(f,x,p)
	local d,n,A,Ap,a,i,b;
	if p=1 or degree(f,x)=0 then return f fi;
	d:=degree(f,x);
	n:=d/p;
	if not type(n,integer) then return NULL fi;
	A:=add(a[i]*x^i,i=0..n-1)+x^n;
	Ap:=collect(A^p,x);
	for i from 1 to n do
	  # Solve a linear equation in one unknown:
	  b[n-i]:=solve(subs({seq(a[n-j]=b[n-j],j=1..i-1)},coeff(Ap,x,d-i))
			-coeff(f,x,d-i),a[n-i] );
	od;
	Ap:=subs({seq(a[n-j]=b[n-j],j=1..n)},Ap);
	subs({seq(a[n-j]=b[n-j],j=1..i-1)},A)
end:



sqrtIrrat := proc(Sirr, Sreg,B,dA,ext,x)
	#option trace;
	local ira, A, c, sum, i,s,d, mult, acc, f,vf,j, p, nonpole, temp, linearr, d1,IND, smin,multmin, minpo, r, a,b;
	userinfo(3,myname,`Now it is the irrational case`);
	f := NULL;
	vf := {};
	ira := indets(Sreg[1][3],{RootOf, name})[1];
if (ira in ext) then
	IND := [op(indets([args],{RootOf,name}))];	
	A :=c;
	sum := add(sign(i[3],IND)*i[3]*i[5], i=Sreg);
	sum := collect(sum, ira, evala);
	sum := lcoeff(sum, ira);
	multmin := infinity;	
	for s in Sreg do
		mult := collect(sign(s[3],IND)*s[3], ira, evala);
		mult := lcoeff(mult, ira);
		mult := mult/sum*dA;
		if type(mult, `integer`) then
			A := A*s[4]^mult;
			if mult<multmin then
				multmin := mult;
				smin := s;
			fi;
		else
			return NULL; #degree doesn't mathc
		fi;
	od;
else
	for i from 1 to nops(Sreg) do
		minpo[i] := poly_d(Sreg[i][3],x, ext);
		if degree(minpo[i],x)<>2 then return NULL fi;
	od;
	a:=coeff(minpo[1],x,0);
	r[1]:=1;
	multmin:=1;
	smin:=Sreg[1];
	for i from 2 to nops(Sreg) do
		b:=coeff(minpo[i],x,0);
		r[i] := sqrt(Normalizer(b/a));
		if not type(r[i],rational) then return NULL fi;
		if r[i]<multmin then 
			multmin :=r[i]; 
			smin:= Sreg[i] 
		fi;
	od;
	sum:=add(r[i],i=1..nops(Sreg));
	sum:=dA/sum;
	multmin := multmin *sum;
	if type(sum, integer) then
		A:=c*mul(Sreg[i][4]^(r[i]*sum),i=1..nops(Sreg));
	else
		return NULL;
	fi;
fi;
	d := 0; 
	j := 0; # use j indicate that if infinity is irregularsigularity
	for p in Sirr do 
		if p[1] = infinity then 
			d := d+ p[3][1];
			j := p[3][2]; 
		else
			d := d+evala(Trace(p[3][1], ext, ext), x); 
		fi; 
	od; 
	nonpole:=1;
	for p in Sirr do
		nonpole:= nonpole*p[4]^p[3][2];
	od;
	#nonpole := mul(p[4]^(p[3][2]),p=Sirr);
	nonpole := Normalizer((A/B-d)*nonpole);
	temp := evala(Rem(numer(nonpole), denom(nonpole), x));
	linearr := coeffs(temp, x);
	#each irregular singularity represents pole part of f. 
	if j<>0 then #the equations at infinity
		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 := SolveTools:-Linear({linearr},{c});
	if temp <> NULL then 
		f := subs(temp, A/B);
		f := op(simplef({f}));
		for i from 0 to floor(multmin/2) do
			vf:= vf union {[(smin[3]+i)/multmin,f]};
		od
	fi; 
	
	
	vf

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,s,l;
	l:=L[1];
	try
	   	v:=evala(Factors(X^d -l,ext));
	catch :
		v:=evala(Factors(numer(normal(X^d - l)),{op(ext)}))
	end try;
	v:=select(has,v[2],X);
	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 := L1[3]*L2[3];
	[evala(L1[1]*L2[1]),L1[2]+L2[2],
	  convert(series(f,tp=0,acc),polynom),
	  acc]
end:

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

sqrtRat := proc(L,Sirr, Sreg,Sa,B,dA,ext,x,T)
	#option trace;
	local  s, l, dnu, m, dv, d1,d2, res, i, knl, j, k,f,ff;
	res:={};
	if Sreg=[] then 
		dnu := {};
		m := dA;
		for i from 3 to floor(m/2) do
			if irem(m,i)=0 then dnu:= dnu union {i};fi;			
		od;
		dnu := dnu union {m};
	else
		dnu := {};
		dv := denom(Sreg[1][3]);
		for s in Sreg do
			dv:=lcm(dv,denom(s[3]));
		od;
		for i from 1 to dA do
			d1:=0;
			d2:=0;
			for s in Sreg do
				d1 := d1+i*dv/denom(s[3])*s[5];
				d2 := gcd(d2, i*dv/denom(s[3]));
			od;
			if d1<=dA and irem(dA,d2)=0 then
				 dnu := dnu union {i*dv}; 
												
			fi;
		od;
	fi;
	userinfo(3,myname,`The possible list of the denominator of nu is`,dnu);
	for i in dnu do
	 	knl :=searchA1k(Sreg,dA,i);
		if knl <> [] then
			ff := findSqrtf(L,Sirr,B,i,knl,ext,x,T);
			ff := simplef(ff);
			if assigned(_Env_original_field) then
				f := {};
				for k in ff do
					if indets(k, RootOf) minus _Env_original_field ={} then
						f:=f union {k};
					fi;
				od;
			else
				f:= ff;
			fi;
			
		fi;
		for j in f do
			res := res union findnuRat(j,i,Sreg,Sa,B,x);
			#for k from 1 to floor(i/2) do
			#	if j<>0  and gcd(k,i)=1 then res := res union {[k/i,j]} fi;
			#od;
		od;
	od;

	res;
end:

findnuRat :=proc(f,d,Sreg,Sa,B,x)
	#option trace;
	local mult, vf, h, s,k,n1,n2,v,p,newh;
	mult:=0;
	h := numer(f);
	vf :={};
	if nops(Sreg) <> 0 then 
		s:=Sreg[1]; 
		if s[1]=infinity then
			mult := degree(denom(f),x)-degree(numer(f),x);
		else
			
			while evala(Rem(h, s[4], x, 'newh')) = 0 do
					mult := mult + 1;
					h := newh
			od;
		fi;
	else
		for p in Sa do
			if p[1]=infinity then
				mult := degree(denom(f),x)-degree(numer(f),x);
			else
			
				while evala(Rem(h, p[4], x, 'h'))=0 do mult := mult+1 od;
			fi;
			if mult<>0 then
				s := p;
				break;
			fi;
		od;
	fi; 

	for k from 1 to floor(d/2) do
		if gcd(k,d)=1 then
			v := k/d;
			if mult=0 then
				vf:=vf union {[v,f]};
			else
				n1 := solve( (v+N)*mult = s[3], N);
				if type(n1,integer) then
					vf := vf union {[abs(v+n1),f]};
					break;
				fi;
				n2 := solve( (N-v)*mult = s[3], N);
				if type(n2,integer) then					
					vf := vf union {[abs(n2-v),f]};
					break;
				fi;
				if n1>1 then
					if abs(round(n1)-n1) <= abs(round(n2)-n2) then
						vf := vf union {[abs(v+round(n1)),f]};
						break;
					else
						vf := vf union {[abs(round(n2)-v), f]};
						break;
					fi;
				fi;
				vf:=vf union {[v,f]};
			fi;

		fi;
	od;
	vf
	
end:

simplef := proc(flist)
	#option trace;
	local f,res;
	res := {};
	for f in flist do
		f := evala(Normal(f, expanded));
		f := evala(Sqrfree(f,x));
		f := f[1]*mul(i[1]^i[2],i=f[2]);
		res := res union {f};
	od;
	res;
end:



searchA1k:=proc(Sreg,dA,d)
	#option trace;
	local de,s,dr,k,res,temp;
	res :=NULL;
	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(dA,d) do
		dr :=dA-k*d;
		temp :=searchA1(Sreg,dr,d);
		for s in temp do res := res, [s,k] od;
	od;
	{res};	
end:

searchA1 :=proc(Sreg,de,d)
	#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[3]);
	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 := searchA1([seq(Sreg[i],i=2..nops(Sreg))],de-k*i,d);
		if temp<>NULL then 
			for t in temp do
				res := res,s[4]^(i*k)*t;
			od;	
		fi;		
	od;	
	[res];
end:

findSqrtf := proc(L,Sirr,B,d,knl,ext,x,T)
	#option trace;
	local res,kn,l,acc,ep, Bs, temp, ser,lp, uk, und,i, j, lr,f, extp, ukn,uki,flist,nlist;
	userinfo(3,myname,`Now it is the rational case and the denominator of nu is`, d);
	res := {};		
	if knl <> {} then 
		for kn in knl do
			for i from 1 to nops(Sirr) do				
				acc := Sirr[i][3][2]; 
				ep := poly_to_powerseries(kn[1],x,Sirr[i][1],T,acc);
				ep := pow_reciprocal(ep,T);
				Bs := poly_to_powerseries(B,x,Sirr[i][1],T,acc);
				temp := pow_time(ep, Bs, T);
				ep := poly_to_powerseries(Sirr[i][3][1],x,Sirr[i][1],T,acc);
				ser[i] := pow_time(temp, ep,T);
				lp := ser[i][1];
				und := indets(lp,RootOf) minus ext;
				if und={} then
					if not assigned(l) then
						l:=lp;
					else
						if length(lp)<length(l) then 
							l := lp;
						fi;
					fi;
				fi;
			od;
			if not assigned(l) then
				return {};
			fi;
			userinfo(3,myname,`Let g=CA_1A_2^d`);
			userinfo(3,myname,`Fix C as`,l);
			userinfo(3,myname,`Fix A_1 as`,kn[1]);	
			uk := add(a[i]*x^i,i=0..kn[2]);
			und := indets(uk) minus {x};
			for i from 1 to nops(Sirr) do
				temp := ser[i];
				if temp[2] <> 0 and Sirr[i][1]<>infinity then error "we made an error" fi;
				temp[1] := temp[1]/l;
				extp := indets(Sirr[i][1],RootOf) union ext;
				ep := dthroot_powseries(temp,T,d,extp);
				ukn := NULL;								 
				for uki in {uk} do 
					temp := poly_to_powerseries( uki,x,Sirr[i][1],T,Sirr[i][3][2]);
					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(Sirr[i][1],RootOf) and not member(Sirr[i][1],ext) then
							lr := map(evala@Expand, lr);
							lr := map(coeffs, lr, Sirr[i][1])
						fi;
						lr := SolveTools:-Linear(lr,und);
						if lr <> NULL then
							f:= subs(op(lr),uki);
							ukn := ukn,f;
						fi;
					od;
				od;
				uk := ukn;
				if uk = NULL then break fi;
			od;
			userinfo(3,myname,`Compute A_2 equal`,uk);
			if uk <> NULL then res := res union {seq(evala(l*i^d*kn[1]/B),i=[uk])} fi;
		od;
	fi;
	
	res;
end:



# Input:
#    i) Two differential operators M1 and M2 and optional the environment d
#   ii) Two differential equations M1 and M2 and the indeterminate function d
# Output: A differential operator M3 such that M3(y) is a solution
#   of M2 for each solution y of M1
#
# Answer is split into 2 parts.
s_equiv := proc(M1,M2)
	local Dx,x,T1,L1,T2,L2,C0,D0,S,C1,D1,d,v,v1,i,w;
	if not assigned(_Envdiffopdomain) then
		error "domain is not assigned"
	fi;
	Dx,x := op(_Envdiffopdomain);
	if not type([M1,M2],list(polynom(ratpoly(anything,x),Dx))) then
		error "wrong arguments or dependent variable not specified"
	elif degree(M1,Dx)<>2 or degree(M2,Dx)<>2 then
		error "case not handled yet"
	fi;
	T1:=Normalizer(coeff(M1,Dx,1)/coeff(M1,Dx,2)/2);
	L1:=DEtools['symmetric_product'](M1,Dx-T1);
	T2:=Normalizer(coeff(M2,Dx,1)/coeff(M2,Dx,2)/2);
	L2:=DEtools['symmetric_product'](M2,Dx-T2);
	C0:=Normalizer(coeff(L1,Dx,0)/coeff(L1,Dx,2));
	D0:=Normalizer(coeff(L2,Dx,0)/coeff(L2,Dx,2));
	if not Same_p_curvature(C0,D0,x) then
		userinfo(2,dsolveBessel,"case discarded by p-curvature test");
		return 0
	fi;
	D0:=Normalizer(D0-C0);
	if D0=0 then
		userinfo(1,dsolveBessel,"No gauge transformation needed");
		d := 0;
		S := 1
	else
		userinfo(3,dsolveBessel,"Computing gauge transformation");
		C1:=Normalizer(diff(C0,x));
		D1:=Normalizer(diff(D0,x));
		d:=Normalizer(D1/D0);
		# v = coeff_list of symmetric_product(L1, L2)
		v:=map(Normalizer,[2*diff(C1,x)+diff(D1,x)-2*d*C1-D1*d+D0^2,
		6*C1+D1-4*d*C0, 2*(D0+2*C0), -d, 1]);

		v, d := RemoveAppSingOrd4(v,x);

		# This _Env passes information on apparant singularities
		# to the expsols code, which then runs faster:
		_Env_DF_sq := table([x, denom(C0)*denom(D0)]);

		# The options tell Expsols to only look at exponential
		# solutions with generalized exponents in 1/2 * Z  (i.e. only
		# exponential solutions whose square is a rational function)
		v:=EXP_SOLS(v,0,x,`use Int`,`no algext`,radical,denom,2);
		
		if nops(v)=0 then
			return 0
		elif nops(v)>1 then
			userinfo(1,'dsolve',"input is reducible");
			return 0
		fi;
		v1 := `if`(type(v[1],`*`),[op(v[1])],v);
		v := 1;
		for i in v1 do
			if not has(i,x) then next fi;
			if type(i,function) and op(0,i)=exp
			and type(op(i),function) and op(0,op(i))=Int
			and op(2,op(i))=x then
				d := d + op(1,op(i))
			else
				v := v * i
			fi
		od;
		w[0] := Normalizer(v);
		for i to 3 do
			w[i] := Normalizer(diff(w[i-1],x) + d*w[i-1])
		od;
		S:=collect(w[0]*Dx +
		  (w[3]+(4*C0+5*D0)*w[1]+w[0]*(2*C1+D1))/2/D0-2*w[1],
		  Dx,Normalizer);
		S := DEtools['symmetric_product'](S, Dx+T1);
		S := evala(Primpart(S, Dx, 'i'));
		d := d + diff(i,x)/i;
		S := frontend(collect,[S,Dx,factor])
	fi;
	d := Normalizer(d + T1 - T2);
	d, S # This represents the operator exp(Int(d,x)) * S
end:

# Try to remove some apparant singularities we may have introduced.
RemoveAppSingOrd4 := proc(v, x)
	local d,i,Dx,L;
	d := mul(`if`(i[2]=1,i[1],1),i=sqrfree(denom(v[-2]),x)[2]);
	if not has(d,x) then return v, 0 fi;
	i := denom(Normalizer(v[-2]-2*diff(d,x)/d));
	d := denom(Normalizer(i/d));
	if not has(d,x) then return v, 0 fi;
	i := 2*v[-3] - 3*(v[-2]*diff(d,x)+2*diff(d,x,x))/d + (3*diff(d,x)/d)^2;
	d := denom(Normalizer( denom(Normalizer(i))/d ));
	if not has(d,x) then return v, 0 fi;
	d := 1/2 * diff(d,x)/d;
	L := DEtools['symmetric_product'](add(v[i+1]*Dx^i,i=0..4), Dx-d, [Dx,x]);
	[seq(Normalizer(coeff(L,Dx,i)),i=0..4)], -d
end:

# Operators Dx^2 + a, Dx^2 + b are compared mod 3
# If the reduction mod 3 fails, then compare mod 5 instead.
Same_p_curvature := proc(a,b, x)
	local v,i,s,s3,T;
	if a=b then return true fi;
	v := traperror(`mod/ReduceField`([a,b],{x},3));
	if v = lasterror then
		return Same_5_curvature(args)
	fi;
	for i to 2 do
		s := Normal(v[i]) mod 3;
		s3 := subs({seq(j=j^3, j=indets(s, {name,symbol,RootOf}))}, s);
		s := Normal(diff(s^2,x,x)) mod 3;
		T[i] := s3 + s
	od;
	evalb(Normal(T[1] - T[2]) mod 3 = 0)
end:

# Operators Dx^2 + a, Dx^2 + b are compared mod 5
Same_5_curvature := proc(a,b, x)
	local v,i,s,s5,j,T;
	v := traperror(`mod/ReduceField`([a,b],{x},5));
	if v = lasterror then
		return true
	fi;
	for i to 2 do
		s := Normal(v[i]) mod 5;
		s5 := subs({seq(j=j^5, j=indets(s, {name,symbol,RootOf}))}, s);
		s := Normal( s * (2*s^2 + diff(s,x,x)) ) mod 5;
		for j to 4 do
			s := Normal( diff(s,x) ) mod 5
		od;
		T[i] := s + s5
	od;
	evalb(Normal(T[1] - T[2]) mod 5 = 0)
end: