# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/DFactorLCLM,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

# Algorithm: DFactorLCLM.
#
# Input: operator L.
#
# Output:
#	operators [L1,L2,..Ln] such that
#	1) L = LCLM(L1,L2,..Ln)
#	2) order(L)=sum(orders(L.i)).
#	3) n is maximal with the above properties. This does not imply that
#	   the L.i are irreducible, however, it does imply that the L.i
#	   can not be split further in this way.
#
# Notes:
#	n equals the largest number of eigenvalues of an endomorphism. A basis
#	of the endomorphisms is computed with DEtools[eigenring].
#
#	If L is what is called "completely reducible" (this means:
#	irreducible or an LCLM of irreducibles) then the L.i will
#	be irreducible. However, this procedure will not determine if
#	this is the case or not, use DFactor to see if L.i can be
#	factored further.
#
# Side-effects:
#	DFactorLCLM(L.i) gets the value [L.i].
#
#	If L.i has no nontrivial endomorphisms then eigenring(L.i)
#	gets the value [1]. If L.i does have nontrivial endomorphisms
#	then L.i is reducible, and DFactor(L.i,`one step`) will be
#	assigned a value.
#
# Purpose:
#	If you have a bases {y.i1, y.i2, ...} of the solution spaces
#	of the L.i, then union of these bases is a basis for L.
#	Condition 1) above guarantees that y.ij span the solutions
#	of L, and condition 2) guarantees linear independence.
#	Whenever n = nops(output) is greater than 1, this simplifies
#	solving L.
#					Mark van Hoeij, July 1999.

macro(
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	eigenring=`DEtools/diffop/eigenring`,
	endomorphism_charpoly=`DEtools/diffop/endomorphism_charpoly`,
	infosubs=`DEtools/diffop/infosubs`,
	GCRD=`DEtools/diffop/GCRD`,
	mult=`DEtools/diffop/mult`,
	rightdivision=`DEtools/diffop/rightdivision`,
	adjoint=`DEtools/diffop/adjoint`,
	substitute=`DEtools/diffop/substitute`,
	rational_solutions=`DEtools/diffop/rational_solutions`,
	LCLM=`DEtools/diffop/LCLM`,
	ratbeke=`DEtools/diffop/ratbeke`,
	g_factors=`DEtools/diffop/g_factors`,
	DFactorLCLM=`DEtools/diffop/DFactorLCLM`,
	er_works=`DEtools/diffop/er_works`
):

er_works:={}:

# Note: use `DEtools/DFactorLCLM`. It does the necessary conversions and
# then calls this procedure.
#
# ext gives the basefield. Can be omitted, in which case ext=indets(L,RootOf).
DFactorLCLM:=proc(L,ext)
	global x,DF,er_works;
	local i,n,v,nv,j,cp,R,res;
 options `Copyright (c) 1999 Waterloo Maple Inc. All rights reserved. \
 Author: M. van Hoeij`;

	i:=collect(L,DF,Normalizer);
	n:=degree(i,DF);
if n<1 then
	error "order should be > 0"
elif i<>L or lcoeff(i,DF)<>1 or nargs=1 then
	procname(collect(i/lcoeff(i,DF),DF,Normalizer),indets([args],RootOf))
elif n=1 then
	[L]
elif nargs=2 then
	# First ratbeke, then for each type try to put in on the
	# left. If that succeeds, then divide out and apply recursion
	# with an option not to do this a second time.
	v:=ratbeke(L);
	if add(nops(i[1]),i=v)=n then
		# L is an LCLM of first order operators that do
		# not involve algebraic extensions.
		[seq(seq(i[2]-diff(j,x)/j,j=i[1]),i=v)];
	else
		res:=NULL;
		R:=L;
		for i in v do
			cp:=adjoint(R);
			nv:=rational_solutions(substitute(i[2],cp));
			if nops(nv)<>nops(i[1]) then
				next
			fi;
			cp:=rightdivision(cp,LCLM(1,seq(
			 DF-coeff(i[2],DF,0)-diff(j,x)/j,j=nv)));
			if cp[2]<>0 then
				error "should not happen"
			fi;
			nv:=seq(i[2]-diff(j,x)/j,j=i[1]);
			cp:=adjoint(cp[1]);
			if GCRD(cp,LCLM(1,nv))<>1 then
				next
			fi;
			R:=cp;
			res:=res,nv
		od;
		R:=collect(R/lcoeff(R,DF),DF,Normalizer);
		if degree(R,DF)<2 then
			error "should not happen"
		fi;
		[res,op(procname(R,ext,`skip ratbeke`))]
	fi
elif _Env_DFactorLCLM=`no eigenring` then
	# skip computation, let dsolve use DFactor instead which is
	# usually faster.
	[L]
else
	v:=eigenring(L);
	nv:=nops(v);
	userinfo(2,'diffop',`Dimension endomorphisms of`,infosubs(L),`is:`,nv);
  if nv=1 then
	[L]
  elif nv=n^2 then
	error "should not happen"
  else
	v:=sort([op({op(v)} minus {1})],(S,T) -> evalb(length(S)<length(T)));
	v:=[op(v),seq(add(rand(-j..j)()*i,i=v),j=[2,9,99,10^4,10^9])];
	# We added a few "generic" ones to the list v. The assumption is that
	# one of those will have a maximal number of eigenvalues, which is highly
	# likely but not certain.
	i,j := 0,0;
	while i<nops(v) and j=0 do
		i:=i+1;
		cp:=endomorphism_charpoly(L,v[i]);
		j:=Normalizer((x+coeff(cp,x,n-1)/n)^n-cp)
	od;
    if j=0 then
	# The following is to tell DFactor to use the eigenring because
	# it works, and it has already been computed.
	er_works:=er_works union {L};
	# All endomorphisms have only 1 eigenvalue (with high probability),
	# so the result is:
	[L]
    else
	v:=v[i];
	res:=NULL;
	cp:=g_factors(cp,x,[op(ext)]);
	for i in cp do
		R,j := 1,0;
		while degree(R,DF)<i[2] do
			R:=GCRD(mult(R,v-RootOf(i[1],x)),L);
			j:=j+1;
			if j>i[2] then
				error "should not happen"
			elif j>1 then
				er_works:=er_works union {R}
			fi
		od;
		res:=res,R
	od;
	res:=[res];
	i:=add(degree(i[1],x),i=cp);
	userinfo(2,'diffop',`Used endomorphism with`,i,`eigenvalues.`);
	if i>nv then
		error "should not happen"
	elif i=nv then
		for j in res do
			eigenring(j):=[1]
		od
	elif i>1 then
		res:=[seq(op(procname(j,indets([j,ext],RootOf))),j=res)]
	fi;
	res
    fi
  fi
fi
end:

#savelib('DFactorLCLM','er_works','`DEtools/diffop/DFactorLCLM.m`'):