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

# 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's for g_ext:
macro(
	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`
);
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's for eval_laurent:
macro(	accuracy_laurent=`DEtools/diffop/accuracy_laurent`,
	value_laurent=`DEtools/diffop/value_laurent`,
	LL=_LL
):
macro(
	lowerbound_val=`DEtools/diffop/lowerbound_val`,
	max0_val=`DEtools/diffop/max0_val`,
	new_laurent=`DEtools/diffop/new_laurent`,
	nmterms_laurent=`DEtools/diffop/nmterms_laurent`,
	nterms_expression=`DEtools/diffop/nterms_expression`,
	nthterm_laurent=`DEtools/diffop/nthterm_laurent`,
	ramification_laur=`DEtools/diffop/ramification_laur`,
	lift_ramification_laur=`DEtools/diffop/lift_ramification_laur`,
	eval_laurent=`DEtools/diffop/eval_laurent`,
	subsvalueslaurents=`DEtools/diffop/subsvalueslaurents`,
	series_val=`DEtools/diffop/pseries_val`,
	valuation_laurent=`DEtools/diffop/valuation_laurent`,
	nt=`DEtools/diffop/nt`
):

# macro's for formal_sol:
macro(
	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`
): # :)

# Should already be loaded by readlibs
# macro(
# g_ext=readlib(`DEtools/diffop/g_ext`),
# l_p=readlib(`DEtools/diffop/l_p`),
# factor_op=readlib(`DEtools/diffop/factor_op`)
# ):
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 bug() 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','`DEtools/diffop/formal_sol.m`'):