# Desingularization (i.e. removal of apparant singularities) for
# differential operators. Note: the singularity at infinity is ignored,
# only finite apparant singularities are removed.
#
# See the paper by:   Abramov, Barkatou, van Hoeij.


# Written for Maple 8.
#
# Implemented by:   Mark van Hoeij, April 26, 2003.
#
desing := proc(LL)
	local Dx,x,lc,L,ext,La,p,w,fm,v,res,s,y,mi,a,sd,t,LE;
	if not assigned(_Envdiffopdomain) then
		_Envdiffopdomain := args[2]
	fi;
	if not type(_Envdiffopdomain,list) then
		error "_Envdiffopdomain should be a list"
	fi;
	Dx, x := op(_Envdiffopdomain);
	lc := lcoeff(LL,Dx);
	L := collect(LL/lc, Dx, evala);
	ext := indets(L, {'nonreal', 'radical', 'algext'});
	La := evala(L);
	p := `DEtools/kovacicsols/factors`(denom(La), x, ext);
	w := NULL;
	for fm in p do
		v := is_desing_f(L, op(fm));
		if [v][1] then
			w := w, [op(fm),v]
		fi
	od;
	if nops([w])=1 and nops(w)=5 then
	    res := w[5]
	else
	    w := {w};
	    res := NULL;
	    while w<>{} do
		w := [seq(exp_and_acc(L,i,Dx,x),i=w)];
		s := mul(i[1]^i[5],i=w);
		mi := mul(i[1],i=w);
		a := mul(i[1]^i[6],i=w);
		s := eval(subs(y(x)=s, DEtools['diffop2de'](numer(La),y(x))));
		s := evala(Normal( s/denom(La) ));
		s := evala(Rem( numer(s), a, x ))/denom(s);
		s := evala(mi * diff(s,x)/s);
		sd := 'sd';
		t := 't';
		if evala(Gcdex(denom(s), a, x, sd, t))<>1 then error fi;
		LE := Dx-evala(Rem( eval(sd)*numer(s), a, x))/mi;
		res := LE, res;
		L := mult(LE, L);
		La := evala(L);
		w := {seq(`if`(i[4][-1]=nops(i[4])-1,NULL,i[1..4]),i=w)}
	    od
	fi;
	res := [res];
	op(res[1..-2]), mult(res[-1], 1/lc)
end:

exp_and_acc := proc(L,v,Dx,x)
	local i;
	for i while v[4][i]=i-1 do od;
	[op( subsop(4 = [op(v[4][1..i-1]), i-1, op(v[4][i..-1])], v[1..4])),
	 i-1, # lowest missing exponent
	 v[4][-1] - i + 2 # accuracy we need to work with
	]
end:


# f is irred. poly
# m = multiplicity of f as a factor of the denominator of L.
is_desing_f := proc(L, f, m)
	local Dx,x,p,q,v,n,T,LL,so;

	Dx, x := op(_Envdiffopdomain);
	if lcoeff(L,Dx)<>1 then
		error "input should be monic"
	fi;
	p := denom(evala(L));

	if m=0 then
		error "f should be in the denominator"
	fi;
	p := RootOf(f, x);
	n := degree(L,Dx);
	v := DEtools['gen_exp'](L, T, x=p);
	if not (nops(v)=1 and type(v[1][1..-2],list(nonnegint))
	and nops({op(v[1][1..-2])}) = n) then
		return false
	fi;

	v := v[1][1..-2];
	# if v[n-m] <> n-m-1 or v[n-m+1] = n-m then
	#	return false, v
	# fi;
	# Double check that this is OK.

	if v[n]=n and v[n-1]=n-2 then
		# This is a special case that occurs very
		# frequently. So I cover this case specially,
		# all other cases are treated with a general
		# method below.

		LL := Dx + diff(f,x)/f + add(c[i]*x^i,i=0..degree(f,x)-1);
		so := evala(Rem( numer(evala(DEtools['mult'](LL,L)*f)) ,f,x));
		so := SolveTools:-Linear( map(coeffs,map(collect,{coeffs(so,x)}
		 ,Dx),Dx),{seq(c[i],i=0..degree(f,x)-1)});
		if so=NULL then
			false, v
		else
			true, v, Dx + evala(subs(so,coeff(LL,Dx,0)))
		fi
	else
		# The general case
		evalb(not DEtools['formal_sol'](L, `has logarithm?`, x=p)), v
	fi
end: