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

# Procedures for users
#    formal_sol, but I think I want to write a user interface for that, so then
#    this procedure would fall in the next category:
# Procedures for use only in the rest of the diffop package:
#    none
# Procedures for use only in this file:
#    solve_semiregular, integrate_logs, lift_integral, sol_1order_eq,
#    lift_sol1

###############################################
#  Formal solutions of differential equations #--------------------------------
###############################################		diffop_formalsol

# This part of diffop computes solutions of f=xDF^n+laur*xDF^(n-1)+...
# in algclos( k((x)) )[log(x),exp( .. )]
# All solutions are given up to conjugation over k((x))
# The result is a list which contains a list of solutions for every exponential
# part e. Such a list for 1 exp part looks like:
# [ [a1,a2,..] , e] where a.i in k((x))[e,log(x)]
# The corresponding solutions are: a.i*ExpInt(e) where ExpInt(e)=exp(int(e/x,x))
# If any ramifications (then x is not really x but a root of a*t^n=x) have been
# used then the list is longer and looks like:
# [ [a1,a2,..] , e,ext,ram,`alg over k((x))` ] where ram=a*x^n
# The same form is used when algebraic extensions were needed.

macro(
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	xDF=`DEtools/diffop/xDF`,
	g_ext=`DEtools/diffop/g_ext`,
	set_laurents=`DEtools/diffop/set_laurents`,
	modm=`DEtools/diffop/modm`,
	g_expand=`DEtools/diffop/g_expand`,
	accuracy_laurent=`DEtools/diffop/accuracy_laurent`,
	value_laurent=`DEtools/diffop/value_laurent`,
	new_laurent=`DEtools/diffop/new_laurent`,
	nmterms_laurent=`DEtools/diffop/nmterms_laurent`,
	eval_laurent=`DEtools/diffop/eval_laurent`,

	formal_sol=`DEtools/diffop/formal_sol`,
	solve_semiregular=`DEtools/diffop/solve_semiregular`,
	integrate_logs=`DEtools/diffop/integrate_logs`,
	lift_integral=`DEtools/diffop/lift_integral`,
	sol_1order_eq=`DEtools/diffop/sol_1order_eq`,
	lift_sol1=`DEtools/diffop/lift_sol1`
	# test_result=`DEtools/diffop/test_result`
): # :)

macro(
	l_p=`DEtools/diffop/l_p`,
	factor_op=`DEtools/diffop/factor_op`
):

# Adapted to Maple 5.4 syntax for []
# Input: an operator f of order >= 1
# Output: the formal solutions
formal_sol:=proc(f,ext)
	global x,DF;
	local v,i,res,c,lp;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if nargs=1 then return procname(f,g_ext(f)) fi;
	if has(f,DF) then return procname(l_p(f,0),ext) fi;
	v:=factor_op(f,'semireg',ext);
	res:=NULL;
	for i in v do
		c:=solve_semiregular(factor_op(i[1],`split over k((x))`,i[2]),i[2]);
		lp:=[c,i[4]];
		if not (i[2]=ext and i[3]=x) then
			lp:=[op(lp),op(i[2..3]),`alg over k((x))`]
		fi;
		# test_result(formal_sol,f,ext,op(lp)[1..2],i);
		res:=res,lp
	od;
	[res]
end:

# Input: a list of factors [xDF+r1,...xDF+rn] where
# the 0'th coefficients of r.i are in Z.
# Output: the solutions
solve_semiregular:=proc(v,ext)
	global x,xDF;
	local s,i;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if v=[] then return [] fi;
	s:=procname([seq(xDF+eval_laurent(v[i]-v[nops(v)]
	 ,ext),i=1..nops(v)-1)],ext);
	s:=[op(map(integrate_logs,s,ext)),1];
	[seq(eval_laurent(expand(i*sol_1order_eq(coeff(v[nops(v)],xDF,0)
	 ,ext)),ext,log(x)),i=s)]
end:

# used by solve_semiregular
# determine g such that x*d/dx (g)=ff where ff in k((x))[log(x)]
integrate_logs:=proc(ff,ext)
	global x;
	local i,d,l,f;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	d:=degree(ff,log(x));
if d=0 then
	if has(ff,set_laurents) then
		f:=eval_laurent(ff,ext);
		# l becomes the integral of f/x:
		nmterms_laurent(f,1,0)*log(x)+new_laurent(1,expand(int(
		 nmterms_laurent(f,0)/x,x)),[lift_integral,f,ext])
	else
		g_expand(int(ff/x,x),ext)
	fi
else
	# Bugfix (oct 94)
	l:=procname(coeff(ff,log(x),d),ext);
	l:=coeff(l,log(x),0)+log(x)*coeff(l,log(x),1)/(d+1);
	expand(l*log(x)^d+procname(convert([seq(coeff(ff,log(x),i)
	*log(x)^i,i=0..d-1)],`+`)-d*log(x)^(d-1)*coeff(l,log(x),0),ext))
fi
end:

# Used by integrate_logs
# lifts LL which is int(aa/x , x)
lift_integral:=proc(LL,aa,ext,acc)
	global x;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	value_laurent[LL]+expand(int(
	nmterms_laurent(aa,acc,accuracy_laurent[LL])/x,x))
end:

# used by solve_semiregular
# Input: a Laurent series a in Z + x k[[x]]
# Output: a solution L of a*L + x*d/dx L = 0
sol_1order_eq:=proc(a,ext)
	global x;
	local n;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	n:=-nmterms_laurent(a,1,0); # n := coefficient of x^0 in a
	if not type(n,integer) then error "should not happen" fi;
	new_laurent(n+1,x^n,[lift_sol1,a,n,ext])
end:

# used by sol_1order_eq
lift_sol1:=proc(LL,aa,n,ext,acc)
	global x;
	local L,a,t,i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	L:=value_laurent[LL];
	a:=nmterms_laurent(aa,acc-n);
	for t from accuracy_laurent[LL] to acc-1 do
		# divide the coefficient of x^t in a*L by n-t:
		L:=L+modm(g_expand(convert([seq(coeff(L,x,i)*coeff(a,x,t-i)
		 ,i=ldegree(L,x)..degree(L,x))],`+`)/(n-t)*x^t,ext))
	od;
	L
end:

#savelib('formal_sol','solve_semiregular','integrate_logs','lift_integral',\
  'sol_1order_eq','lift_sol1'):