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

# Procedures for users
#    none
# Procedures for use only in the rest of the diffop package:
#    eval_laurent, new_laurent, nmterms_laurent, nt,
#    nterms_expression, nthterm_laurent, subsvalueslaurents, valuation_laurent
# Procedures for use only in this file:
#    new_laurent2, series_val

# The purpose of this file is computation with infinite series.

# Global variables: (these are macro'ed)
# x
# LL
# accuracy_laurent
# description_laurent
# set_laurents
# value_laurent
#
# LL.i are indeterminates which stand for Laurent series. The value of a
# Laurent series is given by:
# value_laurent[LL.i]. This value is determined modulo x^accuracy_laurent[LL.i]
#
# description_laurent[LL.i]: consists of: [procedure_name, arguments]
# such that procedure_name(arguments,accuracy) produces the value of LL.i
# modulo x^accuracy.
#
#
# Procedures:
# differentiate           : gives the derivative of an expression
# eval_laurent            : converts an expression to a Laurent series.
# new_laurent             : used for introducing new Laurent series.
# new_laurent2            : used by eval_laurent
# nmterms_laurent         : computes a laurent series mod x^n
# nt                      : computes an expression mod x^n
# nterms_expression       : used by eval_laurent
# nthterm_laurent         : gives the coefficient of x^n in a Laurent series
# subsvalueslaurents      : name says it all
# series_val              : used by nterms_expression
# valuation_laurent       : the valuation (minimum n with coefficient x^n <>0)

# 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`
):

# Adapted to Maple 5.4 syntax for []
new_laurent:=proc()
	global set_laurents,accuracy_laurent,value_laurent,description_laurent;
	local i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	i:=nops(set_laurents)+1;
	set_laurents:=set_laurents union {LL.i};
	if nargs=3 then
		accuracy_laurent[LL.i]:=args[1];
		value_laurent[LL.i]:=args[2];
		description_laurent[LL.i]:=
		 [op(args[3][1..min(1,nops(args[3]))]),LL.i
		 ,`if`(nops(args[3])>1,op(args[3][2..nops(args[3])]),NULL)]
	else
		accuracy_laurent[LL.i]:=NULL;
		value_laurent[LL.i]:=NULL;
		description_laurent[LL.i]:=NULL
	fi;
	LL.i
end:

subsvalueslaurents:=proc(a)
	local i,res;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	res:=a;
	for i in indets(a) intersect set_laurents do
		res:=subs(i=value_laurent[i],res)
	od;
	res
end:

# Adapted to Maple 5.4 syntax for []
# Computes a laurent series modulo x^n
# Input: laurent series, upperbound n, (optional) lower bound m
nmterms_laurent:=proc()
	global accuracy_laurent,value_laurent;
	local g,procedure,dummy;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if accuracy_laurent[args[1]]<args[2] then
		g:=description_laurent[args[1]];
		procedure:=g[1];
		g:=procedure(op(g[2..nops(g)]),args[2]);
		if g=NULL then
			# procedure has handled everything, checking accuracy:
			if accuracy_laurent[args[1]]<args[2] then
				RETURN(nmterms_laurent(args))
			fi
		else
			accuracy_laurent[args[1]]:=args[2];
			value_laurent[args[1]]:=g
		fi
	fi;
	g:=value_laurent[args[1]];
	convert([seq(coeff(g,x,dummy)*x^dummy,dummy=max(args[3..nargs],
	 ldegree(g,x))..args[2]-1)],`+`)
end:

# Gives the n'th coefficient
nthterm_laurent:=proc(l,n)
	global x;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if not member(l,set_laurents) then RETURN(coeff(l,x,n)) fi;
	if accuracy_laurent[l]<=n then
		coeff(nmterms_laurent(l,n+1,n),x,n)
	else
		coeff(value_laurent[l],x,n)
	fi
end:

# Ramification of an element in k((x)). Needed because the
# procedure ramification_of allows only degree(f,xDF)>0.
ramification_laur:=proc(L,r,ext)
	global x;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if r=x then RETURN(L) fi;
	new_laurent(max0_val(L)*degree(r,x),0,[lift_ramification_laur,L,r,ext])
end:

lift_ramification_laur:=proc(laur,L,r,ext,ac)
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	modm(g_expand(subs(x=r,nmterms_laurent(L,ceil(ac/degree(r,x))
	,ceil(accuracy_laurent[laur]/degree(r,x)))),ext))+
	value_laurent[laur]
end:

# Converts an expression to a Laurent series.
# Arguments: f, optional ext,
# optional var (then f is a polynomial in var)
eval_laurent:=proc()
	global x;
	local v,f,ext,i;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	f:=args[1];
	if member(f,set_laurents) then RETURN(f) fi;
	if nargs=1 then ext:=g_ext(f) else ext:=args[2] fi;
	if nargs=3 then
		# Consider the expression as a polynomial over Laurent series
		v:=args[3];
		if lcoeff(f,v)=1 then
			# Don't apply eval_laurent on the constant 1.
			RETURN(expand(convert([seq(eval_laurent(coeff(f,v,i)
			 ,ext)*v^i,i=0..degree(f,v)-1),v^degree(f,v)],`+`)))
		else
			RETURN(expand(convert([seq(eval_laurent(coeff(f,v,i)
			 ,ext)*v^i,i=0..degree(f,v))],`+`)))
		fi
	fi;
# The following is to save something on the number of laurent series used.
for i in indets(f) intersect set_laurents do
if   description_laurent[i][1]=truncate_x and
type(description_laurent[i][3],rational) then
	f:=subs(i=description_laurent[i][3],f)
fi od;
if f<>args[1] then
	RETURN(procname(f,args[2..nargs]))
fi;
	if has(f,log(x)) then for i in indets(f,anything) do
	 if i=log(x) then
		RETURN(expand(convert([seq(i^v*eval_laurent(coeff(f,i,v),ext)
		,v=ldegree(f,i)..degree(f,i))],`+`)))
	fi od fi;
	v:=indets(f) intersect set_laurents;
	if v={} then if type(f,ratpoly(rational,[x,op(ext)])) then
		f:=[g_expand(numer(f),ext,x),g_expand(denom(f),ext,x)];
		if type(f[2],polynom(rational,x))
		and ldegree(f[2],x)=degree(f[2],x) then
			f:=expand(f[1]/f[2]);
			RETURN(new_laurent(0,truncate_x(f,f,0),[truncate_x,f]))
		fi;
		RETURN(new_laurent2([seq(new_laurent(0,truncate_x(i,i,0),
		 [truncate_x,i]),i=f)],`/`,ext))
	elif type(f,procedure) then
		RETURN(new_laurent2(f,procedure,ext))
	else
		RETURN(new_laurent2(evala(Normal(subs(g_conversion2,f))),'maple',ext))
	fi fi;
	if type(f,polynom(rational,v)) and degree(f,v)<=2 and type(f,`+`)
	and coeffs(f,v)=1 then
		RETURN(new_laurent2([f,v],`deg 2 polynom`,ext))
	fi;
	v:=[op(v)];
	if not type(f,polynom(anything,v)) then
		v:=g_normal(f);
		RETURN(new_laurent2(map(eval_laurent,[numer(v),denom(v)])
		,`/`,ext))
	fi;
	v:=v[1];
	if f=x^degree(f,x)*v then
		RETURN(new_laurent2([degree(f,x),v],`*x^n`,ext))
	fi;
	f:=collect(f,v);
	if f=v^degree(f,v) then
		f:=eval_laurent(f/v);
		RETURN(new_laurent2([f,v],`*`,ext))
	fi;
	new_laurent2([eval_laurent(coeff(f,v,0)),new_laurent2(
	 [eval_laurent(expand((f-coeff(f,v,0))/v)),v],`*`,ext)],`+`,ext)
end:


# For use with eval_laurent.
new_laurent2:=proc()
	local l;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if member(`deg 2 polynom`,[args]) then
		l:=lowerbound_val([args]);
		l:=min(l,2*l);
		new_laurent(l,0,[nterms_expression,args])
	elif nargs>1 and args[2]=`/` then
		new_laurent(max0_val(args[1][1])-valuation_laurent(args[1][2],infinity)
		,0,[nterms_expression,args])
	elif nargs>1 and args[2]=`*` then
		new_laurent(max0_val(args[1][1])+max0_val(args[1][2])
		,0,[nterms_expression,args])
	else
		new_laurent(-10,nterms_expression(0,args,-10),[nterms_expression,args])
	fi
end:

# For use with eval_laurent
# This procedure computes an expression mod x^n.
nterms_expression:=proc(lau,aa,what,ext,n)
	global x;
	local a,s,d,v1,v2;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if what='maple' then
	s:=series(aa,x=0,max(n+series_val(denom(aa)),
	 series_val(denom(aa))+2));
	if order(s)<>infinity then
	d:=0;
	while order(s)<>infinity and order(s)<n and s<>aa do
		d:=d+1;
		s:=series(aa,x=0,max(n+series_val(denom(aa))
		 ,series_val(denom(aa))+2)+d)
	od
	fi;
	s:=evala(Expand(convert(s,polynom)));
	normal_coeffs(subs(g_conversion1,convert(
	 [seq(coeff(s,x,d)*x^d,d=ldegree(s,x)..n-1)],`+`)),x)
elif what=`+` then
	if lau<>0 and value_laurent[lau]<>0 then
		convert(map(nmterms_laurent,aa,n,accuracy_laurent[lau]),`+`)
		 +value_laurent[lau]
	else
		normal_coeffs(convert(map(nmterms_laurent,aa,n),`+`),x)
	fi
elif what=`*` then
	v2:=max0_val(aa[2]);
	a:=nmterms_laurent(aa[1],n-v2);
	modm(g_expand(value_laurent[lau]+convert([seq(coeff(a,x,d)*x^d
	 *nmterms_laurent(aa[2],n-d,accuracy_laurent[lau]-d)
	 ,d=ldegree(a,x)..degree(a,x))],`+`),ext))
elif what=`*x^n` then
	expand(x^aa[1]*nmterms_laurent(aa[2],n-aa[1]))
elif what=`/` then
	v1:=max0_val(aa[1]);
	v2:=valuation_laurent(aa[2],infinity);
	if lau=0 or value_laurent[lau]=0 then
		g_quotient(nmterms_laurent(aa[1],v2+n),nmterms_laurent(
		 aa[2],2*v2+n-v1),n,ext)
	else
		a:=nmterms_laurent(aa[2],2*v2+n-v1);
		value_laurent[lau]+g_quotient(nmterms_laurent(aa[1],n+v2
		 ,accuracy_laurent[lau]+v2)-
		convert([seq(convert([seq(x^s*coeff(value_laurent[lau],x,s)
		 ,s=accuracy_laurent[lau]-d+v2..n-d-1+v2)],`+`)*x^d
		 *coeff(a,x,d),d=ldegree(a,x)..degree(a,x))],`+`),a,n,ext)
	fi
elif what=differentiate then
	diff(nmterms_laurent(aa,n+1),x)
elif what=`deg 2 polynom` then
	# Now aa[1] is a sum
	normal_coeffs(
	convert(map(nterms_expression,[op(aa[1])],accuracy_laurent[lau],
	`deg 2 1 term`,ext,n),`+`),x)
elif what=`deg 2 1 term` then
	if type(lau,`*`) then
		v2:=max0_val(op(2,lau));
		a:=nmterms_laurent(op(1,lau),n-v2);
		modm(g_expand(convert([seq(coeff(a,x,d)*x^d
		 *nmterms_laurent(op(2,lau),n-d,aa-d)
		 ,d=ldegree(a,x)..degree(a,x))],`+`),ext))
	else
		nmterms_laurent(lau,n,aa)
	fi
elif what=procedure then
	convert([value_laurent[lau],seq(aa(d)*x^d
	,d=accuracy_laurent[lau]..n)],`+`)
else bug()
fi
end:

# valuation of a. It looks at most as far as bound
valuation_laurent:=proc(a,bound)
	global x,value_laurent;
	local i;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	i:=value_laurent[a];
	if i<>0 then
		if g_normal(tcoeff(i,x))=0 then
			value_laurent[a]:=collect(i,x,g_normal);
			RETURN(procname(args))
		fi;
		i:=ldegree(i,x);
		if accuracy_laurent[a] < min(i+1,bound) then
			# Increase accuracy_laurent[a]:
			nmterms_laurent(a,min(i+1,bound)+1);
			RETURN(procname(args))
		fi;
		RETURN(min(bound,i))
	fi;
	i:=accuracy_laurent[a]-1;
	while (bound=infinity or i<bound) and nmterms_laurent(a,i+1)=0 do
		i:=i+1
	od;
	if i>=bound then bound else procname(args) fi
end:

# This procedure computes the valuation of an expression.
series_val:=proc(a)
	global x;
	local d,n;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if a=0 then ERROR(`division by zero`) fi;
	d:=1;
	n:=0;
	while n=0 do
		d:=d+1;
		n:=evala(Expand(convert(series(a,x=0,d),polynom)))
	od;
	ldegree(n,x)
end:

# Returns a lower bound for the valuation of the Laurent series in f.
lowerbound_val:=proc(f)
	local i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	min(seq(
	max0_val(i),
	 i=indets(f) intersect set_laurents))
end:

# A lower bound (not much too low) for the valuation
max0_val:=proc(a)
	global x;
	local v;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	v:=value_laurent[a];
	if v=0 then
		v:=accuracy_laurent[a];
		if v<0 then valuation_laurent(a,0) else v fi
	else
		ldegree(v,x)
	fi
end:

# nt (n terms), compute f mod x^n
# Only for expressions linear in Laurent series, like the standard form
# for local differential operators.
nt:=proc(f,n)
	global x;
	local a,i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if has(f,log(x)) then RETURN(subs(a=log(x),nt(subs(log(x)=a,f),n))) fi;
	if type(f,list) then RETURN(map(nt,f,n)) fi;
	map(nmterms_laurent,indets(f) intersect set_laurents,n);
	truncate(subsvalueslaurents(f),n,x,[])
end:

#savelib('new_laurent','subsvalueslaurents','nmterms_laurent', \
  'nthterm_laurent','ramification_laur','lift_ramification_laur','eval_laurent', \
  'new_laurent2','nterms_expression','valuation_laurent', \
  'series_val','lowerbound_val','max0_val','nt','`DEtools/diffop/eval_laurent.m`'):