# Remove at least all apparant singularities
#  (and sometimes but not always more).
#
# Note that this is a small modification of the algorithm given in our
# paper. Namely, in the code below we compute a linear combination of L
# and all L2's that occur during the computation, whereas in the paper we
# compute a linear combination of L and just the last L2. As a
# consequence, the code below removes as least as many singularities as
# the program in the paper, but might remove more. However, it still does
# not always remove all singularities that can be removed (that can be
# done with the linear algebra based desingularization program which you
# can find at this URL:
# http://www.math.fsu.edu/~hoeij/papers/desing/DifferenceCase/desing
#
tdesing := proc(L, domain)
	local tau,x,a0,ad,L2,L3,i,s,t,n;
	tau,x := op(domain);
	a0 := tcoeff(L,tau);
	ad := lcoeff(L,tau);
	L3 := L;
	n := LREtools['dispersion'](ad,a0,x,'maximal');
	if not type(n,integer) then n := 0 fi;
	L2 := L/a0;
	for i to n do
		L2 := collect(L2, tau, Normalizer);
		# Clear the tau^i coefficient of L2:
		L2 := L2 - coeff(L2,tau,i)*
			subs(x=x+i,L/a0)*tau^i;
		L2 := primpart(L2, tau);
		gcdex(tcoeff(L3,tau),coeff(L2,tau,0),x,'s','t');
		L3 := collect(s*L3+t*L2, tau, expand)
	od;
	sort(collect(L3, tau, factor), domain)
end: