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

# Procedures for users
#    none
# Procedures for use only in the rest of the diffop package:
#    pade2, lift_pade2 (until now only used in DFactor)
# Procedures for use only in this file:
#    smaller_pade2, smallest_pade2, valuation_pade2, wipe_pade2

#####################################
#    Pade Hermite approximations    #-----------------------------------------
#####################################		diffop_pade2

# The following is a variant on Harm Derksen's (hderksen@sci.kun.nl) algorithm
# for computing Pade expansions. You can find his algorithm in the share
# library under the name pade2. G. Labahn and him found indepent of each other
# a good method for computing these expansions.

# I made the following changes:
# 1) This pade2 accepts my syntax for elements of k((x))  (laur1, laur2 etc.)
# 2) Harm pointed out to me that pade2 can also work with vectors of power
#  series. I replaced the variable z in his pade2 by z[1], z[2], etc. to 
#  achieve this.
# 3) point is always x=0
# 4) ext is the coefficients field
# Note that it is useful to have the original pade2 in order to understand
# how the method works.

# 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`,
	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'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(	lift_pade2=`DEtools/diffop/lift_pade2`,
	pade2=`DEtools/diffop/pade2`,
	smaller_pade2=`DEtools/diffop/smaller_pade2`,
	smallest_pade2=`DEtools/diffop/smallest_pade2`,
	valuation_pade2=`DEtools/diffop/valuation_pade2`,
	wipe_pade2=`DEtools/diffop/wipe_pade2`
):

# functionlist is a list of lists of elements in k((x))
# pade2 looks for a k[x]-linear relation between these lists (interpret these
# lists as vectors). If all these lists contain only 1 element, then this is
# basically Harm Derksen's pade2 algorithm.
pade2:=proc(functionlist,accuracy,ext)
	local i,j,n,y,z,yvars,zvars,fl,appr,degrees;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	n:=min(seq(seq(valuation_laurent(j,0),j=i),i=functionlist));
	if n<0 then
		RETURN(pade2([seq([seq(eval_laurent(x^(-n)*j)
		 ,j=i)],i=functionlist)],accuracy,ext))
	fi;
	n:=nops(functionlist);
	degrees:=[seq(0,i=1..n)]; # Not used yet, for later generalizations
	zvars:=[seq(z[i],i=1..nops(functionlist[1]))];
	yvars:=[seq(y[i],i=1..n)];
	fl:=0,[seq(convert([seq(functionlist[i][j]*z[j],j=1..nops(zvars))]
	 ,`+`),i=1..n)],n,degrees,yvars,zvars,ext,0;
	if subs(infinity=0,accuracy)<>accuracy then RETURN([fl]) fi;
	appr:=lift_pade2(fl,accuracy);
	smallest_pade2(op(appr))
end:

# This procedure lifts the list appr
# Output: new value of appr
# appr, acc, n and degrees is the same as in Harm Derksen's pade2
# functionlist: same as Harm's but multiplied by z
# y is a list of variables [y[1],y[2],..] are the same as in Harm's algorithm
# zvars is a list of variables [z[1],z[2],..] which play the role of Harm's
# variable z.
# old_acc: 0 on the first call, and the previous acc otherwise
# appr: 0 on the first call, and the previous result of lift_pade2 otherwise
lift_pade2:=proc(appr,functionlist,n,degrees,y,zvars,ext,old_acc,extra_acc)
local app,i,k,j,z,acc;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	userinfo(9,'diffop','lifting',extra_acc,'steps');
	acc:=old_acc+extra_acc;
	if old_acc=0 then
		# Since I can't change appr, I'll call it app here.
		app:=[seq(expand(x^degrees[i]*y[i]+nt(
		 functionlist[i],acc)),i=1..n)]
	else
		app:=map(g_expand,subs({seq(y[i]=y[i]+x^(-degrees[i])
		 *(nt(functionlist[i],acc)-nt(functionlist[i],old_acc))
		 ,i=1..n)},appr),ext)
	fi;
	app:=modm(app);
	for i from old_acc to acc-1 do for z in zvars do
		k:=0;
		for j to n do if valuation_pade2(app[j],z,x,acc)=i
		 and (k=0 or smaller_pade2(app[j],app[k],x,y)) then
			k:=j
		fi od;
		if k>0 then
			app:=[seq(wipe_pade2(app[j],app[k],z,x,i,ext),j=1..k-1),
			 expand(x*app[k]),seq(wipe_pade2(app[j],app[k],
			 z,x,i,ext),j=k+1..n)]
		fi
	od od;
	userinfo(9,'diffop',`done lifting`);
	[app,functionlist,n,degrees,y,zvars,ext,acc]
end:

smallest_pade2:=proc(appr,dummy,n,degrees,y,zvars,ext,dummy2)
local smallest,i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	smallest:=1;
	appr;
	for i from 2 to n do if smaller_pade2(appr[i],appr[smallest],x,y) then
		smallest:=i;
	fi od;
	[seq(g_expand(coeff(appr[smallest],y[i],1)/x^degrees[i],ext),i=1..n)]
end:

valuation_pade2:=proc(f,z,y,acc)
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if coeff(f,z,1)=0 then acc else ldegree(coeff(f,z,1),y) fi
end:

wipe_pade2:=proc(f,g,z,y,i,ext)
	local ff,gg;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	ff:=coeff(coeff(f,z,1),y,i);
	gg:=g_normal(ff/coeff(coeff(g,z,1),y,i));
	modm(g_expand(f-gg*g,ext))
end:

smaller_pade2:=proc(f,g,x,vars)
local ff,gg,i,df,dg,n;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
n:=nops(vars);
ff:=[seq(degree(coeff(f,vars[i],1),x),i=1..n)];df:=max(op(ff));
gg:=[seq(degree(coeff(g,vars[i],1),x),i=1..n)];dg:=max(op(gg));
if df<dg then RETURN(true)
elif df>dg then RETURN(false)
fi;
for i to n do
if ff[i]=df and coeff(f,vars[i],1)<>0 then RETURN(false) fi;
if gg[i]=dg and coeff(g,vars[i],1)<>0 then RETURN(true);fi;
od;
end:

#savelib('pade2','lift_pade2','smaller_pade2','smallest_pade2',\
  'valuation_pade2','wipe_pade2','`DEtools/diffop/pade2.m`'):