# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/expsols,v $
# $Notify: hoeij@sci.kun.nl $

# Procedures for users
#    ratbeke, expsols
# Procedures for use only in the rest of the diffop package:
#    none
# Procedures for use only in this file:
#    try_combinations

# macro's for g_ext:
macro(
	infosubs=`DEtools/diffop/infosubs`,
	modulus=`DEtools/diffop/modulus`,
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	xDF=`DEtools/diffop/xDF`
):

macro(
	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	rootof=`algcurves/rootof`,
	degree_ext=`algcurves/degree_ext`,
	g_evala=`algcurves/g_evala`,
	g_evala_rem=`algcurves/g_evala_rem`,
	g_zero_of=`algcurves/g_zero_of`,
	v_ext_m=`algcurves/v_ext_m`,
	ext_to_coeffs=`algcurves/e_to_coeff`,
	g_ext_r=`algcurves/g_ext_r`,
	g_gcdex=`algcurves/gcdex`,
	truncate=`algcurves/truncate`,
	g_factors=`algcurves/g_factors`
);
macro(
	new_laurent2=`DEtools/diffop/new_laurent2`,
	differentiate=`DEtools/diffop/differentiate`,
	g_ext=`DEtools/diffop/g_ext`,
	g_factors=`DEtools/diffop/g_factors`,
	g_normal=`DEtools/diffop/g_normal`,
	set_laurents=`DEtools/diffop/set_laurents`,
	description_laurent=`DEtools/diffop/description_laurent`,
	modm=`DEtools/diffop/modm`,
	g_solve=`DEtools/diffop/g_solve`,
	g_quotient=`DEtools/diffop/g_quotient`,
	g_rem=`DEtools/diffop/g_rem`,
	g_quo=`DEtools/diffop/g_quo`,
	truncate_x=`DEtools/diffop/truncate_x`,
	g_expand=`DEtools/diffop/g_expand`
):

macro(
	gen_exp=readlib(`DEtools/diffop/gen_exp`),
	substitute=readlib(`DEtools/diffop/substitute`),
	solve=readlib(`solve/linear`),
	LCLM=readlib(`DEtools/diffop/LCLM`),
	rightdivision=readlib(`DEtools/diffop/rightdivision`),
	readlibs=readlib(`DEtools/diffop/readlibs`),
	l_p=`DEtools/diffop/l_p`,
	nt=`DEtools/diffop/nt`,
	DFactor=`DEtools/diffop/DFactor`,

	ratbeke=`DEtools/diffop/ratbeke`,
	try_combinations=`DEtools/diffop/try_combinations`,
	expsols=`DEtools/diffop/expsols`
):

# Compute all exponential solutions (i.e. first order right hand factors) that
# do not involve algebraic extensions over the base field, by Beke's method (cf. for
# example section 3.4 in my PhD thesis).
# RootOf syntax.

# If the last argument has a DF in it, then ratbeke will assume that the
# singularities of that operator are the only non-apparant singularities of f.
# In some applications this can speed up the computation.
ratbeke:=proc(f,ext)
	global x,DF;
	local singularities,p,m,i,N,opt;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if nargs<2 or not(type(ext,list)) then
		# for the options remember to work:
		RETURN(procname(f,g_ext([args]),args[2..nargs]))
	fi;
	if has([args],rational) then
		# Compute only integer generalized exponents, i.e. only
		# rational solutions, no exponential solutions, will be
		# computed.
		opt:=integer
	elif has([args],radical) then
		opt:=rational
	else
		opt:=NULL
	fi;
	N:=`if`(has(args[nargs],DF),args[nargs],f);
	singularities:={seq(i[1],i=v_ext_m(g_factors(denom(
	 g_normal(N/lcoeff(N,DF))),x,ext),x)),infinity};
	singularities:=subs(g_conversion2,singularities);
for p in singularities do
	# Compute generalized exponents, through away the ones that
	# are not minimal, that have ramification index <> 1 and those
	# that involve an algebraic extension:
	m[p]:=[seq(`if`(degree_ext(i,[p,args])=1,
		# args has to be sent to gen_exp, cause if you don't you get a bug
		# in test example
		# de41=98*x^2+931392+(x^3+28908*x)*DF+296*x^2*DF^2+x^3*DF^3
	 i[1],NULL),i=gen_exp(f,p,x,'minimal','ramification1',opt,args))];
	# Is the option type=operator specified?
	N:=subs({seq(`if`(type(i,equation),i,NULL),i=args)},type);
	# If yes, select the proper exponential part:
	if degree(N,DF)=1 then
	 m[p]:=[seq(`if`(type(evala(Normal(subs(g_conversion2,
		nt(coeff(l_p(N,p),xDF,0),1)-i
	  ))),integer),i,NULL),i=m[p])]
	fi;
	if m[p]=[] then
		# there are no exponential solution that don't involve
		# algebraic extensions
		userinfo(2,'diffop',`0 possible combinations to check`);
		RETURN([])
	fi
od;
	N:=convert([seq(nops(m[p]),p=singularities)],`*`);
	userinfo(2,'diffop',N,`possible combinations to check`);
	if has([args],`not too hard`) and N>100 then
		RETURN(`giving up`)
	fi;

	# Optimization that can be done: suppose the dimension of those
	# solutions that we found so far that have exponential part e at p
	# is \mu_{e,p}(f), then we can skip e at p from the list, thereby reducing
	# the number of combinations that remains to be checked.
	# Perhaps I'll implement this (not very important) optimization later.

	# build up the input for try_combinations:
N:=try_combinations(
  [seq(
	[seq(map(evala,map(Trace,
	 # Take the Trace of these two things:
	 [subs(x=x-p,i/x),coeff(i,x,0)],
	 # over the following field extension:
	 op(subs(g_conversion2,[g_ext([p,args]),g_ext([args])]))
	)),i=m[p])]
  ,p=singularities minus {infinity})]
,m[infinity],f,DF);
	userinfo(2,'diffop',nops(N),
		`types of rational exponential solutions found`);
	N
end:

# recursive procedure for checking all combinations, doing the corresponding
# substitutions, and trying to find exponential solutions (which are polynomial
# solutions after the substitutions).
# V is for the point infinity
# v other points
# R1 is of order 1, I collect the substitutions in there.
# RootOf syntax.
try_combinations:=proc(v,V,f,R1)
	global x,DF;
	local a,m,ansatz,i,s,res,F,S,j,C;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	res:=NULL;
if nops(v)=0 then
	# Now look at V, the gen. exp. at infinty:

	# skip the ones that don't have an integer as
	# the constant term of the generalized exponent at infinity
	for i in V do m:=-g_normal(coeff(i,x,0)); if type(m,integer) and m>=0 then
		S:=collect(subs(x=1/x,i+m)/x,x);
		if S<>0 then
			# after this substitution the irregular singular
			# part of this generalized exponent is gone, what remains
			# is the exponent -m
			F:=substitute(DF-S,f)
		else
			F:=f
		fi;
		# Make an ansatz for the polynomial solutions, find
		# linear equations, solve, and return the answer.
		ansatz:=convert([seq(C[j]*x^j,j=0..m)],`+`);
		a:=ansatz;
		s:=a*coeff(F,DF,0);
		for j from 1 to m do
			a:=diff(a,x);
			s:=s+a*coeff(F,DF,j)
		od;
		m:=indets(ansatz) minus {x};
		# the following are the polynomial solutions:
		a:=subs(solve({coeffs(collect(numer(normal(s)),x),x)},m)
		 ,ansatz);
		# Now make a basis of solutions and return the answer:
		m:=m intersect indets(a);
		s:=solve(m,m);
		# basis of the polynomial solutions:
		s:=[seq(subs(j=1,s,a),j=m)];
		if s<>[] then
			res:=res,[s,substitute(DF+S,R1)]
		fi
	fi od;
else
	for i in v[1] do
		# i[1]=the subs that must be done
		# i[2]=the effect of this subs on the gen. exp. at infinity
		res:=res,op(procname(v[2..nops(v)],[seq(j+i[2],j=V)],
		 substitute(DF+i[1],f),substitute(DF-i[1],R1)))
	od
fi;
	[res]
end:

# Input: a differential operator f and the non-homogeneous part h
#   (h used to have to be of the form exp(int(something,x))), but now
#   is allowed to be anything).
# Output: the exponential solutions, cf. the help page for expsols.
# The option rational will cause that for the homogeneous equation only
# the rational solutions are computed.
expsols:=proc(f,h)
	global x,DF;
	local i,j,L,N,R,V,W,ext,sol;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if nargs=1 or args[2]=rational then
		RETURN(procname(f,0,args[2..nargs]))
	fi;
	# trivial cases first:
	if degree(f,DF)=0 then
		if h<>0 then
			RETURN([[],h/f])
		else
			RETURN([])
		fi
	elif degree(f,DF)=1 then
		sol:=Int(-coeff(f,DF,0)/coeff(f,DF,1),x);
                L := `if`(h=0, [exp(sol)],
			[[exp(sol)],exp(sol)*Int(h*exp(-sol)/coeff(f,DF,1),x)]);
		RETURN(L)
	fi;

	# load the necessary code:
	readlibs(ratbeke);

	# First: rational exponential solutions:
	V:={op(ratbeke(f,args[3..nargs]))};
	# Dimension of the rational exponential solutions:
	N:=convert(map(i -> nops(i[1]),[op(V)]),`+`);
	R:=1;
	# solutions of R:
	W:=[];
if N<degree(f,DF)-1 and not has([args],rational) then
	# see if there are more exponential solutions:
	ext:=g_ext([args]);
	if N=0 then
		readlibs(DFactor);
		N:=DFactor(f,`one step`);
		if nops(N)>1 then
			R:=LCLM(N[2],`and conjugates`,ext);
			W:=procname(N[2],0,ext)
		fi
	else
		# remove the right hand factors to obtain the
		# following right hand factor of f:
		R:=LCLM(seq(seq(i[2]-diff(j,x)/j,j=i[1]),i=V))
	fi;
	if R<>1 then
		# Now compute the left hand factor:
		L:=rightdivision(f,R);
		if L[2]<>0 then ERROR() fi;
		L:=L[1];
		# apply recursion:
		L:=procname(L,0,ext);
		for i in L do if degree_ext(i,ext)>1 then
		  # We skip the ones with degree_ext(i,ext)=1 because
		  # those have already been caught by ratbeke.
		  # Now check if we don't already have this type:
		    j:=type=DF-g_normal(diff(i,x)/i);
		    if R=1 or ratbeke(R,j,R)=[] then
			V:=V union {op(ratbeke(f,
			 type=DF-g_normal(diff(i,x)/i),f))}
		    fi
		fi od
	fi
fi;
	V:=[op(W),seq(seq(exp(Int(DF-i[2],x))*j,j=i[1]),i=V)];
	if h<>0 then
		L:=substitute(DF+g_normal(diff(h,x)/h),f);
# The following line is not so efficient, one should
# compute the exponents from the gen.exp of f (since these
# were already known).
		L:=readlib(`DEtools/ratsols`)(
		  [seq(coeff(L,DF,i),i=0..degree(L,DF))],1,x);
		V:=[V,seq(g_normal(i*h),i=L[2..nops(L)])]
	fi;
	V
end:

#savelib('ratbeke','try_combinations','`DEtools/diffop/ratbeke.m`'):
#savelib('expsols','`DEtools/diffop/expsols.m`'):