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

# factor_op.m:
# Procedures for users
#    none. If I want to export some of the procedures in here I think
#    it would be better to place them in another file, so that this file
#    will remain "internal tools only".
# Procedures for use only in the rest of the diffop package:
#    nm_mult, op_with_slope, nm_block, factor_op, Newtonpolygon,
#    coefs_operator, ramification_of
# Procedures for use only in this file:
#    factor_newton, factor_newton2, lift_newton, factor_riccati, lift_ramification
# Things I intend to delete:
#    skipped_factors

# make_rightfactor.m
# Procedures for users
#    none
# Procedures for use only in the rest of the diffop package:
#    faster_riccati_split, make_rightfactor
# Procedures for use only in this file:
#    lift_rightfactor, lift_rsplit, lift_rsplit2, indeterminates_op,
#    indeterminate, nm_block2

# substitute.m and subs_local.m
# Procedures for users
#    substitute
# Procedures for use only in the rest of the diffop package:
#    subs_local
# Procedures for use only in this file:
#    lift_substitute

##############################
# Computation in k((x))[xDF] #------------------------------------------------
##############################			diffop_local

# 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'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 factor_op and substitute
macro(
	Newtonpolygon=`DEtools/diffop/Newtonpolygon`,
	coefs_operator=`DEtools/diffop/coefs_operator`,
	factor_newton=`DEtools/diffop/factor_newton`,
	factor_newton2=`DEtools/diffop/factor_newton2`,
	factor_op=`DEtools/diffop/factor_op`,
	factor_riccati=`DEtools/diffop/factor_riccati`,
	lift_newton=`DEtools/diffop/lift_newton`,
	lift_ramification=`DEtools/diffop/lift_ramification`,
	lift_substitute=`DEtools/diffop/lift_substitute`,
	nm_mult=`DEtools/diffop/nm_mult`,
	nm_block=`DEtools/diffop/nm_block`,
	op_with_slope=`DEtools/diffop/op_with_slope`,
	ram_laur=`DEtools/diffop/ram_laur`,
	ramification_of=`DEtools/diffop/ramification_of`,
	rem_lift=`DEtools/diffop/rem_lift`,
	skipped_factors=`DEtools/diffop/skipped_factors`,
	substitute=`DEtools/diffop/substitute`
):

# macro's for make_rightfactor
macro(
	faster_riccati_split=`DEtools/diffop/faster_riccati_split`,
	indeterminate=`DEtools/diffop/indeterminate`,
	indeterminates_op=`DEtools/diffop/indeterminates_op`,
	lift_rightfactor=`DEtools/diffop/lift_rightfactor`,
	lift_rsplit=`DEtools/diffop/lift_rsplit`,
	lift_rsplit2=`DEtools/diffop/lift_rsplit2`,
	make_rightfactor=`DEtools/diffop/make_rightfactor`,
	nm_block2=`DEtools/diffop/nm_block2`
):
macro(
	rightdivision=readlib(`DEtools/diffop/rightdivision`),
	substitute=readlib(`DEtools/diffop/subs_local`),
	make_rightfactor=readlib(`DEtools/diffop/make_rightfactor`)
):

# This part of diffop deals with
# the local factorization. The syntax for local operators is:
# f=xDF^n + LL.? * xDF^(n-1) + ..
# So f is monic, and has elements of set_laurents as coefficients.
# Furthermore these operators must use my own syntax for algebraic numbers,
# and must use xDF (differentiation followed by a multiplication by x)
# as a syntax for operators, instead of DF.

# Gives terms of a multiplication in k[x,xDF] from x^low to x^high
# So the accuracy is high+1
nm_mult:=proc(l,r,low,high,ext)
	global x,xDF;
	local j,i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	modm(g_expand(convert([seq(coeff(r,x,j)*x^j*subs(xDF=xDF+j
	,convert([seq(x^i*coeff(l,x,i),i=max(ldegree(l,x),low-j)..min(degree(
	l,x),high-j))],`+`)),j=max(ldegree(r,x),low-degree(l,x))..min(high
	-ldegree(l,x),degree(r,x)))],`+`),ext))
end:

# Adapted to Maple 5.4 syntax for []
# Returns an operator having a given slope.
op_with_slope:=proc(order,slope,shift,d)
global description_laurent,xDF;
local i,result;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	result:=xDF^args[1];
	for i from 0 to args[1]-1 do
		result:=result+new_laurent(ceil((i-args[1])*
		 slope+shift),0,[])*xDF^i # Note: shift <= 0.
	od;
	if d=[] then RETURN(result)
	elif member(d[1],{lift_newton,lift_rsplit}) then
		result:=[result,op_with_slope(d[2],slope,0,[])]
	fi;
	for i in indets(result) minus {xDF} do
		description_laurent[i]:=
		 [d[1],i,args[1..2],op(d[2..nops(d)]),result]
	od;
	result
end:

# Gives all terms from x^n to x^m of f if slope=0.
nm_block:=proc(f,n,m,slope,round_off)
	global xDF;
	local i,res;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if round_off(n)<=0 and 0<round_off(m)
		then res:=xDF^degree(f,xDF)
	else
		res:=0
	fi;
	for i from degree(f,xDF)-1 by -1 to 0 do
		res:=res+nmterms_laurent(coeff(f,xDF,i),
		 round_off(m+slope*(i-degree(f,xDF))),
		 round_off(n+slope*(i-degree(f,xDF))))*xDF^i
	od;
	modm(expand(res))
end:

#########################################
#       Factorization in k((x))[xDF]    #
#########################################

# Input: a monic operator f in xDF, in Laurent form
# a string what:
# `split over k((x))` gives a factorization over k((x))
# `all right factors` gives all safe right factors over k((x))
# `all alg factors` gives all safe right factors over alg. closure of k((x))
# `alg factor` gives one safe right factor over alg. closure k((x))
# ext: gives the constant field
factor_op:=proc()
	global xDF;
	local dummy;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	userinfo(6,'diffop','factorizing',infosubs(args));
	if nargs=1 then
		# shouldn't use factor_op with 1 argument
		RETURN(factor_op(eval_laurent(args[1],[],xDF)
		 ,`split over k((x))`,[]))
	elif nargs=2 then RETURN(factor_op(args,g_ext([args])))
	fi;
#	if degree(args[1],xDF)<=1 then RETURN([args[1]]) fi;
# That if caused a bug because the output is in the wrong format for
# the use in formal_sol. So I moved this line to factor_newton.
	[seq(op(factor_riccati(dummy,args[2..nargs])),
	 dummy=factor_newton(args))]
end:

#################################################
#       Factorization using the Newton method   #
#################################################

# Input: an operator f, and (optional) a second argument
# Output (if called with 2 arguments): a list giving 3 things: 
#  1) the coordinates of the extreme points of the Newton polygon
#  2) the slope of a point to the next point
#  3) the Newton polynomial of this slope
# If called with only 1 argument then only the extreme points will be given.
Newtonpolygon:=proc()
	global Newtonpolygon,x,xDF;
	local f,n,vals,dummy,res,m,i,powd,
	 val_powd,npg,slope;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	userinfo(8,'diffop',`Computing Newton polygon of`,infosubs(args));
	f:=args[1];
	n:=degree(f,xDF);
	vals:=[seq(valuation_laurent(coeff(f,xDF,dummy),0),dummy=0..n-1),0];
	powd:=0;
	val_powd:=min(op(vals));
	npg:=[[powd,val_powd]];
	while powd<n do
		m:=min(seq((vals[powd+dummy+1]-val_powd)/dummy
		 ,dummy=1..n-powd));
		i:=n;
		while vals[i+1]-val_powd<>(i-powd)*m do i:=i-1 od;
		powd:=i;val_powd:=vals[i+1];
		npg:=[op(npg),[i,vals[i+1]]]
	od;
	if nargs=1 then RETURN(npg) else Newtonpolygon(args[1]):=npg fi;
	res:=NULL;
	for i from 1 to nops(npg)-1 do
		# Now we compute the Newton polynomial of slope number i
		slope:=(npg[i+1][2]-npg[i][2])/(npg[i+1][1]-npg[i][1]);
		res:=res,[op(npg[i]),slope,convert([seq(
		 nthterm_laurent(coeff(f,xDF,denom(slope)*dummy+npg[i][1]),
		  dummy*numer(slope)+npg[i][2])
		 *x^dummy,dummy=0..(npg[i+1][1]-npg[i][1])/denom(slope)
		)],`+`)]
	od;
	[res,npg[nops(npg)]]
end:

# Output: coprime index 1 factorizations
# Input: operator f, monic in xDF
# a string what:
# `split over k((x))` results in a complete Newton factorization, i.e. the
#  broken Newton polygon and the gcd 1 reducible Newton polynomial cases
#  will be factored
# `all right factors` gives all possible safe right factors according to
#  the Newton method
# `all alg factors` idem
# `alg factor` gives only one right factor, the one where the slope is the
#  least steep
# 'semireg' gives a coprime index 1 LCLM factorization, to be used for
#  computing the semi-regular parts
# A Newton polynomial like (a-1)^2 will only be factored when slope=0
# (regular singular case). No algebraic extensions will be made here, they
# will be made in factor_riccati
factor_newton:=proc(f,what,ext)
	global skipped_factors,x,xDF;
	local np,v,dummy,i,j,k,d,e,res,unsafe,n_unsafe,semi;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if degree(f,xDF)<=1 then RETURN([f]) fi;
	np:=Newtonpolygon(f,`include the Newton polynomials`);
	userinfo(7,'diffop',`Newton method factorizing`,infosubs(args));
	res:=NULL;
	unsafe:={};
for k from 1 to nops(np)-1 do
	v:=g_factors(np[k][4],x,ext);
	if np[k][3]<>0 then
		v:=[seq([g_expand(dummy[1]^dummy[2],ext,x),0],dummy=v)]
	else
		# regular singular case, we compute the unsafe factors
		unsafe:={};
		for i in v do n_unsafe[i]:=0;semi[i]:=i[1]^i[2] od;
		# n_unsafe is the number of factors with an exponent that is
		# the exponent of this factor i minus a positive integer.
		# The "unsafe" means that we can not use those factors as
		# a right hand factor because it is not a priori known whether
		# such a factor can be lifted or not (hence: placing such a
		# factor on the right hand side would cause the risk of an
		# error during the lifting process). If such an "unsafe"
		# factor can be lifted or not depends on if there is a log in
		# the formal solutions.
		# The n_unsafe number is used in the global factorization, in
		# the procedure skipped_factors.
		for i from 1 to nops(v)-1 do for j from i+1 to nops(v) do if
		 degree(v[i][1],x)=degree(v[j][1],x) then
			d:=degree(v[i][1],x)-1;
			e:=g_normal(coeff(v[i][1],x,d)-coeff(v[j][1],x,d));
			if type(e,integer) and e<>0 and irem(e,d+1)=0 and
			 g_expand(v[i][1]-subs(x=x+e/(d+1),v[j][1]),ext,x)=0 then if
				e>0 then
					unsafe:=unsafe union {v[i]};
					n_unsafe[v[j]]:=n_unsafe[v[j]]+v[i][2];
					semi[v[j]]:=semi[v[j]]*v[i][1]^v[i][2]
				else
					unsafe:=unsafe union {v[j]};
					n_unsafe[v[i]]:=n_unsafe[v[i]]+v[j][2];
					semi[v[i]]:=semi[v[i]]*v[j][1]^v[j][2]
			fi fi
		fi od od;
		v:=[op({op(v)} minus unsafe)];
		if what='semireg' then
			v:=[seq([g_expand(semi[i],ext,x),1],i=v)]
		fi
	fi;
	# Now v contains the safe right factors of slope number k
	for i in v do
		if degree(v[1][1],x)*denom(np[k][3])=degree(f,xDF) then
			# f allows no coprime index 1 factorization
			RETURN([f])
		fi;
		j:=factor_newton2(f,i[1],np[k][3],ext);
		if np[k][3]=0 and what<>'semireg' then
			skipped_factors(j[2]):=n_unsafe[i]
		fi;
		if what=`alg factor` then
			RETURN([j[2]])
		elif what=`split over k((x))` then
			RETURN([op(factor_newton(j[1],what,ext)),j[2]])
		fi;
		res:=res,j[2]
	od
od; [res]
end:

# For use with: factor_op( .. , `all right factors`, []) in the global
# factorization.
# It gives for each factor the number of factors we skipped, because of
# the integer difference of the exponents (called the unsafe factors).
# If minimum multiplicity = 1 then skipped_factors=0. If skipped_factors=0
# then not necessarily minmult=1, however, we may still treat this case in
# the global factorization as if it were minmult=1, because then in this
# exponential part there is only 1 unique right hand factor anyway (the others
# give rise to a log in the formal solutions).
#
# This procedure will probably get obsolete soon.
skipped_factors:=proc(f)
	local v;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	v:=[op(indets(f) intersect set_laurents)];
	v:=description_laurent[v[1]];
	if v[1]=lift_rightfactor then
		skipped_factors(v[5])
	elif v[1]=lift_rsplit then
		skipped_factors(v[8])
	elif v[1]=lift_substitute then
		skipped_factors(v[6])
	elif member(v[1],{lift_newton,nterms_expression,truncate_x}) then
	# the options remember must treat the regular singular case
		0
	else
		ERROR(`?`)
	fi
end:

# Factor monic operator f, such that r is the Newton polynomial of
# the right factor
# Output: a list of 2 factors
# r must be monic.
factor_newton2:=proc(f,r,slope,ext)
	global rem_lift,x,xDF;
	local np,l,i,res,shift;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	userinfo(7,'diffop',`Computing factor with Newton polynomial:`,infosubs(r));
	np:=Newtonpolygon(f,`include the Newton polynomials`);
	for i from 1 to nops(np)-1 do if np[i][3]=slope then l:=i fi od;
	if l=evaln(l) then ERROR(`slope not found`) fi;
	shift:=min(0,np[l][2]+(degree(f,xDF)-np[l][1])*slope);
	# We will compute the Newton polynomial of the left factor
	np:=g_expand(x^np[l][1]*subs(x=x^denom(slope),g_quo(np[l][4],r,x,ext)),ext,x);
	res:=op_with_slope(degree(f,xDF)-degree(r,x)*denom(slope)
	 ,slope,shift,[lift_newton,degree(r,x)*denom(slope),f,np,
	 subs(x=x^denom(slope),r),shift,ext]);
	rem_lift[res]:=0;
	res
end:

# This procedure lifts "coprime index 1 factorizations"
# Meaning of the variables:
# llaur: the Laurent series we lift
# acc: desired accuracy of LL. Note that the other series in left, right
# are lifted too.
# n_lifts: number of lifts to do
# n_known: number of computed coefficients
# slope
# mult_args: arguments for nm_mult to determine which terms are needed
# ff, left, right: exact operators f=left*right
# l, r: the lowest parts of left and right
# f: a middle part of ff, precisely the terms we need
# lr: the product l*r
# rem_lift: Those terms of the previous lift of lr, that we can use now.
# l_low, r_low: the lowest line of left and right. These have gcd 1.
# If they have gcd<>1 than it is not uniquely liftable.
# l_extra, r_extra: the new terms we computed of l and r.
# ext: the algebraic extension
lift_newton:=proc(llaur,order_l,slope,order_r,ff,l_low,r_low,shift,ext,v,acc)
	global rem_lift,value_laurent,accuracy_laurent,x,xDF;
	local f,lr,l,r,l_extra,r_extra,s,t,n_known,i,lau
	 ,n_lifts,mult_args,left,right,ll_low,rr_low,fe,le,re;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	n_lifts:=acc-accuracy_laurent[llaur];
	if n_lifts<=0 then RETURN() fi;
	left:=v[1];
	right:=v[2];
	n_known:=accuracy_laurent[coeff(right,xDF,0)]+degree(right,xDF)*slope;
	mult_args:=ceil(n_known+shift-slope*degree(ff,xDF)),
	 ceil(n_known+shift+n_lifts)-1,ext;
	f:=nm_block(ff,n_known+shift,n_known+n_lifts+shift,slope,ceil);
	l:=expand(subsvalueslaurents(left));
	if n_known+shift<=0 then l:=l-xDF^degree(l,xDF) fi;
	if n_known=0 then
		r:=0
	else
		r:=expand(subsvalueslaurents(right))
	fi;
	lr:=rem_lift[v]-f+nm_mult(l,r,mult_args[2]+1-n_lifts
	 ,mult_args[2..nops([mult_args])]);
if slope=0 then # regular singular case
	rr_low:=subs(x=xDF,r_low);
	for i from n_known to n_known+n_lifts-1 do
		if i=0 then
			r_extra:=subs(x=xDF,r_low);
			l_extra:=expand(x^shift*subs(x=xDF,l_low))
		else
			ll_low:=subs(x=xDF,g_expand(subs(x=x+i,l_low),ext));
			s:=g_gcdex(rr_low,ll_low,1,xDF,ext);
			if s[1]<>1 then
				# ERROR(`unsafe factor`)
				bug()
			fi;
			t:=s[3];
			s:=s[2];
			r_extra:=-g_rem(t*coeff(lr,x,i+shift),rr_low,xDF,ext);
			l_extra:=modm(collect(-x^(i+shift)*g_quo(coeff(lr,x
			 ,i+shift)+r_extra*ll_low,rr_low,xDF,ext),x));
			r_extra:=modm(collect(x^i*r_extra,x))
		fi;
		l:=l+l_extra;
		lr:=modm(lr+nm_mult(l,r_extra,mult_args)+nm_mult(l_extra,r
		 ,mult_args));
		r:=r+r_extra
	od
else
	s:=g_gcdex(r_low,l_low,1,x,ext);
	t:=s[3];
	s:=s[2];
	le:=denom(slope)*shift-numer(slope)*degree(left,xDF);
	re:=-numer(slope)*degree(right,xDF);
	fe:=denom(slope)*shift-numer(slope)*degree(ff,xDF);
	for i from n_known*denom(slope) to (n_known+n_lifts)*denom(slope)-1 do
		if i=0 then
			r_extra:=coefs_operator(r_low,slope,re,1);
			l_extra:=coefs_operator(l_low,slope,le,1)
		else
			r_extra:=g_rem(t*coefs_operator(lr,slope,
			 i+fe,0),r_low,x,ext);
			l_extra:=coefs_operator(modm(g_expand(-g_quo((
			 coefs_operator(lr,slope,i+fe,0)-r_extra*l_low),r_low,x,ext)
			 ,ext)),slope,i+le,1);
			r_extra:=coefs_operator(modm(-r_extra),slope,i+re,1)
		fi;
		l:=l+l_extra;
		lr:=lr+nm_mult(l_extra,r,mult_args)+nm_mult(l,r_extra,mult_args);
		r:=r+r_extra
	od;
fi;
	rem_lift[v]:=lr;
	for f in [[left,collect(l,xDF)],[right,collect(r,xDF)]] do
		for i from 0 to degree(f[1],xDF)-1 do
			lau:=coeff(f[1],xDF,i);
			value_laurent[lau]:=coeff(f[2],xDF,i);
			accuracy_laurent[lau]:=accuracy_laurent[lau]+n_lifts
		od
	od;
	NULL
end:

# Input: an operator or a polynomial f
# Output: a polynomial if f is an operator, and an operator if f is
# a polynomial. 
# This procedure either takes coefficients from an operator and gives them as
# a polynomial, or constructs an operator from a given set of coefficients
# (given as a polynomial). Something like a Newton polynomial.
# slope must be > 0
coefs_operator:=proc(f,slope,i,what)
	global x,xDF;
	local dummy,start_x,start_D;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	start_D:=modp(-i/numer(slope),denom(slope));
	start_x:=start_D*slope+i/denom(slope);
	if what=0 then
		convert([seq(coeff(coeff(f,x,start_x+numer(slope)*dummy),xDF,
		 start_D+dummy*denom(slope))*x^(start_D+dummy*denom(slope)),
		 dummy=0..floor(degree(f,xDF)/denom(slope)))],`+`)
	else
		convert([seq(coeff(f,x,start_D+dummy*denom(slope))*x^(start_x
		 +numer(slope)*dummy)*xDF^(start_D+dummy*denom(slope))
		 ,dummy=0..ceil(degree(f,x)/denom(slope)))],`+`)
	fi
end:

#################################################
#       Factorization using a Riccati solution  #
#################################################

# Adapted to Maple 5.4 syntax for []
# f should be monic, have only 1 slope, and the Newton polynomial is the
# power of an irreducible polynomial. The dummy argument is because of the
# options remember in case of a different groundfield.
# output: same as factor_op
# If what='semireg' then the output is a list of things like:
# [semi-regular operator,ext,ram,exp part,'semireg']
# the exp part e is of the form: ramification_of(xDF-real exp part)=xDF-e
factor_riccati:=proc(f,what,ext)
	global x,xDF;
	local np,slope,gr,res,r,l,lr,i,v,n,dummy,ex;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	np:=Newtonpolygon(f,`include the Newton polynomials`)[1];
	slope:=np[3];
if what='semireg' and slope=0 then
	v:=op(factor_op(f,`alg factor`,ext));
	if not type(v,list) then v:=[v,ext,x] fi;
	i:=xDF-nt(v[1],1);
	n:=degree(f,xDF)/degree_ext(v[2],ext);
	if n=1 then
		RETURN([[eval_laurent(v[1]+i,v[2],xDF),op(v[2..3]),i]])
	fi;
	# Only search for integer roots:
	np:=factors(numer(g_expand(subs(x=x+i,np[4]),v[2])))[2];
	r:=1;
	for l in np do if type(x-l[1],integer) then
		r:=r*l[1]^l[2]
	fi od;
	if degree(r,x)<>n then bug() fi;
	n:=min(solve(r));
	r:=expand(subs(x=x+n,r));
	RETURN([[factor_newton2(substitute(xDF+i+n,f,0,0,v[2]),r,0,v[2])[2]
	 ,op(v[2..3]),i+n]])
fi;
	np:=g_factors(np[4],x,ext)[1];
if (degree(f,xDF)<=1) or (np[2]=1 and member(what
 ,{`split over k((x))`,`all right factors`})) then
	if what='semireg' then
		# degree(f,xDF)=1 and slope<>0
		i:=xDF-nt(f,1);
		RETURN([[eval_laurent(f+i,ext,xDF),ext,x,i]])
	else
		# f is irreducible
		RETURN([f])
	fi
elif degree(np[1],x)=1 and denom(slope)=1 then
	# now slope<>0
	np:=(x-np[1])*x^(-slope);
	userinfo(8,'diffop','substituting',infosubs(xDF=xDF+np));
	# Apply a homomorphism xDF -> xDF+np to simplify the problem.
	v:=factor_op(substitute(xDF+np,f,slope,0,ext),what,ext);
	if what='semireg' then
		RETURN([seq([op(i[1..3]),i[4]+
		 solve(ramification_of(xDF-np,i[3],i[2]),xDF)],i=v)])
	fi;
	# Apply the inverse homomorphism. Note: extension could have changed
	RETURN(substitute(xDF-np,v,slope,0,g_ext(v)))
elif what=`split over k((x))` then
	r:=factor_riccati(f,`alg factor`,ext);
	# An order 1 factor over the alg. closure of k((x))
	r:=make_rightfactor(f,r[1],ext); # a righthand factor over k((x))
	if r=f then RETURN([f]) fi;
	l:=rightdivision(f,r,ext,slope)[1];
	# We factorized f here in l and r, however, these l and r lift very
	# slowly, so we try
	userinfo(7,'diffop',`Computing factor using Riccati solution`,infosubs(r));
	lr:=faster_riccati_split(f,l,r,ext,slope);
	RETURN([op(factor_riccati(lr[1],what,ext)),lr[2]])
elif what=`all right factors` then
	v:=factor_riccati(f,`all alg factors`,ext);
	res:=NULL;
	for i in v do
		r:=make_rightfactor(f,i,ext);
		if r=f then res:=f else
			l:=rightdivision(f,r,ext,slope)[1];
			userinfo(7,'diffop',
			 `Computing factor using Riccati solution`,infosubs(r));
			lr:=faster_riccati_split(f,l,r,ext,slope);
			res:=res,lr[2]
		fi
	od;
	RETURN([res])
elif degree(np[1],x)>1 then
	# Now we need an algebraic extension on k
	gr:=g_zero_of(np[1],x,dummy);
	userinfo(7,'diffop',`Making an algebraic extension:`,infosubs(gr));
	ex:=[gr,op(ext)];
	r:=factor_newton2(f,g_expand((x-gr)^np[2],ex)
	 ,slope,ex);
	r:=factor_riccati(r[2],what,ex);
	# Now we denote the algebraic extension in a list:
	res:=NULL;
	for i in r do if type(i,list) then res:=res,i
		else res:=res,[i,ex,x,`alg over k((x))`]
	fi od;
	RETURN([res])
else
	# Now we need a ramification
	userinfo(7,'diffop',`Applying a ramification`);
	n:=mods(1/numer(slope),denom(slope));
	r:=g_normal((x-np[1])^n)*x^denom(slope);
	v:=factor_newton2(ramification_of(f,r,ext),g_expand((x-denom(slope)
	 *(x-np[1])^((1-n*numer(slope))/denom(slope)))^np[2],ext)
	 ,numer(slope),ext);
	v:=factor_riccati(v[2],what,ext);
	res:=NULL;
	for i in v do if type(i,list) then
		res:=res,[i[1],i[2],ramification_of(r,i[3],ext),i[4]]
	else
		res:=res,[i,ext,r,`alg over k((x))`]
	fi od;
	RETURN([res])
fi
end:

# gives a ramification, it maps xDF to xDF* 1/degree(r), and maps x to r
# r is a power of x (then we have a pure ramification), or a power of x
# times a constant. We allow these latter kind of ramifications, because
# then we need less algebraic extensions.
ramification_of:=proc(f,r,ext)
	global x,xDF;
	local s;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if r=x then RETURN(f)
	elif indets(f) intersect set_laurents = {} then
		RETURN(g_expand(subs({x=r,xDF=xDF/degree(r,x)},f),ext)*
		 degree(r,x)^degree(f,xDF))
	fi;
	s:=Newtonpolygon(f);
	s:=s[nops(s)-1];
	s:=s[2]/(s[1]-degree(f,xDF));
	op_with_slope(degree(f,xDF),degree(r,x)*s
	 ,0,[lift_ramification,f,r,ext])
end:

lift_ramification:=proc(laur,order,slope,ff,r,ext,f,acc)
global value_laurent,accuracy_laurent,x,xDF;
local i,res,n_lifts;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	n_lifts:=acc-accuracy_laurent[laur];
	if n_lifts<=0 then RETURN() fi;
	n_lifts:=ceil(n_lifts/degree(r,x))*degree(r,x);
	i:=(accuracy_laurent[coeff(f,xDF,0)]+order*slope)/degree(r,x);
	res:=g_expand(degree(r,x)^order*subs({x=r,xDF=xDF/degree(r,x)},
	 nm_block(ff,i,i+n_lifts/degree(r,x),slope/degree(r,x),ceil)),ext);
	for i from 0 to degree(f,xDF)-1 do
		value_laurent[coeff(f,xDF,i)]:=
		 value_laurent[coeff(f,xDF,i)]+coeff(res,xDF,i);
		accuracy_laurent[coeff(f,xDF,i)]:=
		 accuracy_laurent[coeff(f,xDF,i)]+n_lifts
	od;
	NULL
end:

#savelib('nm_mult','op_with_slope',\
  'nm_block','factor_op','Newtonpolygon',\
  'factor_newton','skipped_factors','factor_newton2','lift_newton',\
  'coefs_operator','factor_riccati','ramification_of','lift_ramification',\
  '`DEtools/diffop/factor_op.m`'):

###########################################################
#  make right-factor over k((x)) using a Riccati solution #
###########################################################

# The contents of this section consists of section 2.7 in my Ph.D thesis,
# and the lift algorithm for factorizations with coprime index > 1.

# I intend to change the global factorization in such a way that this code
# won't be used in there anymore, but for now this code is still used.

macro(
	make_rightfactor=`DEtools/diffop/make_rightfactor`
):

# This section computes LCLM(xDF-r,`and conjugates over k((x))`) where r is a
# Riccati solution. In this way a local irreducible factor is obtained. This
# is only necessary if the coprime index is > 1, otherwise one can factor
# via the Newton method.

# f is an operator
# ric is a right-factor of order 1 over the algebraic closure of k((x))
# output: a right factor over k((x))
make_rightfactor:=proc(f,ric,ext)
	global x,xDF;
	local d,v;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	userinfo(3,'diffop',`Constructing a local right hand factor`);
	d:=degree(ric[3],x);
	v:=-valuation_laurent(coeff(ric[1],xDF,0),0)/d;
	# Now we multiply by the degree of the algebraic extension
	d:=d*degree_ext(ric[2],ext);
	if d=degree(f,xDF) then RETURN(f) fi;
	op_with_slope(d,v,0,[lift_rightfactor,ric,ext,2])
end:

# This procedure lifts the factor generated by make_rightfactor. It is rather
# slow, therefor we also have a faster lift method. We still need this
# procedure because the faster method does not work for the first few lifts.
# Meaning of the variables:
# ric: solution of the Riccati, i.e. 1st order factor over alg. clos. k((x))
# ram: ramification is a map x -> c * x^(ram) for a constant c.
# b[0], b[1], ..: these are used as an Ansatz
# This procedure generates linear equations stating that ric is a righthand
# factor. This righthand factor gets computed by solving these equations.
lift_rightfactor:=proc(dummy_laur,order,slope,ric,ext,acc,som,dummy_ac)
	global value_laurent,accuracy_laurent,description_laurent,x,xDF;
	local s,r,ram,rp,i,b,dummy,fout,l,ld,extl;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	ram:=degree(ric[3],x);
	l:=coeff(ric[1],xDF,0);
	r:=-nmterms_laurent(coeff(ric[1],xDF,0),acc);
	ld:=min(0,ldegree(r,x));
	rp:=1;
	s:=b[0];
	for i from 1 to order do
		l:=(i-1)*ld+acc;
		rp:=truncate(x*diff(rp,x)+r*rp,l,x,ext)/ram+fout[i]*x^l;
		s:=s+b[i]*rp
	od;
	s:=g_expand(subs(b[order]=1,s),ric[2]);
	extl:=NULL;
	if ram>1 then
		s:=subs(x=g_zero_of(subs(x=xDF,ric[3])-x,xDF,'extl'),s);
		extl:=eval(extl);
s:=expand(s/extl^ldegree(s,extl));
		s:=g_expand(s,[extl,op(ric[2])])
	fi;
	s:={ext_to_coeffs(s,ext)};
	assign(g_solve(s,{seq(b[dummy],dummy=0..order-1)},ext,0));
	for i from 0 to order-1 do
		l:=coeff(som,xDF,i);
		b[i]:=nterms_expression(0,g_normal(b[i]),'maple',ext
		 ,accuracy_laurent[l]+3);
		b[i]:=collect(b[i],x,g_normal);
		s:=ldegree(b[i],x)-1;
# Needed because s=ldegree( .. ) can be infinity in the new version of maple.
s:=min(s,accuracy_laurent[l]+3);
		while indets(coeff(b[i],x,s)) intersect {seq(fout[dummy]
		 ,dummy=1..order)} = {} and s<accuracy_laurent[l]+3 do
			s:=s+1
		od;
		value_laurent[l]:=truncate(b[i],s,x,ext);
		accuracy_laurent[l]:=s;
		description_laurent[l]:=[lift_rightfactor,l,args[2..5]
		 ,acc+2*order,som]
	od;
	NULL
end:

############################################################################
#  Faster split method  (is now called: coprime index > 1 lift algorithm)  #
############################################################################

# I think I'll change the global factorizer such that it does not use any
# coprime index>1 factorizations anymore. When that is done, this section
# and the previous section won't be important anymore (only for local
# factorizations).

# Input: f=left*right
# Output: l and r such that l=left and r=right, and such that l and r lift
# faster then left and right.
faster_riccati_split:=proc(f,left,right,ext,slope)
	global xDF;
#       RETURN([left,right]);
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	op_with_slope(degree(left,xDF),slope,0,[lift_rsplit,degree(right,xDF)
	 ,f,left,right,ext,2])
end:

lift_rsplit:=proc(laur,order_l,slope,order_r,f,left,right,ext,n_ind,v,acc)
	global description_laurent,value_laurent,accuracy_laurent,xDF;
	local try,i,lau,aclau;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	try:=lift_rsplit2(laur,f,slope,op(v),n_ind,ext,acc);
	if try=`n_ind is too small` then # I must use more indeterminates
		userinfo(6,'diffop',`Increasing the coprime index of`,infosubs(right),`to`,n_ind+1);
		for i in indets(v) minus {xDF} do
			description_laurent[i]:=
			 [lift_rsplit,i,args[2..8],n_ind+1,v]
		od;
	elif try=`too few terms are known` then
		# the slow procedure lift_rightfactor will compute more terms
		for i from 0 to degree(v[1],xDF)-1 do
			lau:=coeff(v[1],xDF,i);
			aclau:=accuracy_laurent[lau]+1;
			value_laurent[lau]:=nmterms_laurent(coeff(left,xDF,i)
			 ,aclau);
			accuracy_laurent[lau]:=aclau
		od;
		for i from 0 to degree(v[2],xDF)-1 do
			lau:=coeff(v[2],xDF,i);
			aclau:=accuracy_laurent[lau]+1;
			value_laurent[lau]:=nmterms_laurent(coeff(right,xDF,i)
			 ,aclau);
			accuracy_laurent[lau]:=aclau
		od
	fi;
	NULL
end:

lift_rsplit2:=proc(llaur,ff,slope,left,right,n_ind,ext,acc)
	global value_laurent,accuracy_laurent,x,xDF;
	local f,lr,l,r,l_extra,r_extra,l_low,r_low,
		left_ind,right_ind,
		n_known,i,v,all_ind,lau,n_lifts,mult_args;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	n_lifts:=acc-accuracy_laurent[llaur]; # number of lifts
	if n_lifts<=0 then RETURN() fi;
	n_known:=accuracy_laurent[coeff(right,xDF,0)]+degree(right,xDF)*slope;
	# number of computed coefficients
	if n_known<n_ind then
		RETURN(`too few terms are known`)
	fi;
	l:=expand(subsvalueslaurents(left));
	r:=expand(subsvalueslaurents(right));
	mult_args:=ceil(n_known-slope*degree(ff,xDF))
	 ,ceil(n_known+n_lifts+n_ind-2),ext;
	f:=nm_block(ff,n_known,n_known+n_lifts+n_ind,slope,ceil);
	lr:=nm_mult(l,r,mult_args)-f;
	left_ind:=indeterminates_op(slope,n_ind,left);
	right_ind:=indeterminates_op(slope,n_ind,right);
	all_ind:=indets([left_ind,right_ind]) minus {xDF,x};
	l_low:=nm_block(left,0,n_ind,slope,ceil);
	r_low:=nm_block(right,0,n_ind,slope,ceil);
	for i from 1 to n_lifts do
		l_extra:=expand(x^n_known*left_ind);
		r_extra:=expand(x^n_known*right_ind);
		v:=lr+nm_mult(l_extra,r_low,mult_args)
		 +nm_mult(l_low,r_extra,mult_args);
		v:=nm_block2(v,n_known,n_known+n_ind,slope,ceil,degree(ff,xDF));
		v:=g_solve({coeffs(v,[x,xDF])},all_ind);
		l_extra:=g_expand(subs(v,nm_block2(l_extra,n_known,
		 1+n_known,slope,ceil,degree(left,xDF))),ext);
		r_extra:=g_expand(subs(v,nm_block2(r_extra,n_known,
		 1+n_known,slope,ceil,degree(right,xDF))),ext);
		if indets([l_extra,r_extra]) intersect all_ind <> {} then
			RETURN(`n_ind is too small`)
		fi;
		lr:=lr+nm_mult(l+l_extra,r_extra,mult_args)+nm_mult(l_extra,r,mult_args);
		l:=l+l_extra;
		r:=r+r_extra;
		n_known:=n_known+1
	od;
	for i from 0 to degree(left,xDF)-1 do
		lau:=coeff(left,xDF,i);
		value_laurent[lau]:=coeff(l,xDF,i);
		accuracy_laurent[lau]:=accuracy_laurent[lau]+n_lifts
	od;
	for i from 0 to degree(right,xDF)-1 do
		lau:=coeff(right,xDF,i);
		value_laurent[lau]:=coeff(r,xDF,i);
		accuracy_laurent[lau]:=accuracy_laurent[lau]+n_lifts
	od;
	NULL
end:

indeterminates_op:=proc(slope,number,f)
	global x,xDF;
	local i,j,res;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	res:=0;
	for i from 0 to degree(f,xDF)-1 do for j from 0 to number-1 do
		res:=res+indeterminate()*x^(j+ceil(
		 (i-degree(f,xDF))*slope))*xDF^i
	od od;
	res
end:

indeterminate:=proc() local r;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
r end:

nm_block2:=proc(f,n,m,slope,round_off,deg)
	global x,xDF;
	local i,j;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	modm(convert([seq(seq(coeff(coeff(f,xDF,i),x,j)*x^j*xDF^i,j=
	 round_off(n+slope*(i-deg))..round_off(m+slope*(i-deg))-1
	),i=0..degree(f,xDF))],`+`))
end:

#savelib('make_rightfactor','lift_rightfactor','faster_riccati_split',\
  'lift_rsplit','lift_rsplit2','indeterminates_op','indeterminate',\
  'nm_block2','`DEtools/diffop/make_rightfactor.m`'):

macro(
	mult=readlib(`DEtools/diffop/mult`),
	subs_local=`DEtools/diffop/subs_local`,
	substitute=`DEtools/diffop/substitute`
):

# This procedure substitutes a for DF in a global (no laurent series)
# operator f.
# a must be DF + something
# Arguments: a,f,ext (rootof syntax)
#  or  a,f (RootOf syntax)
substitute:=proc()
	global DF,x,xDF;
	local i,ide,res,D;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if nargs=3 then
	if has(args[2],DF) then D:=DF else D:=xDF fi;
	res:=coeff(args[2],D,0);
	for i from 1 to degree(args[2],D) do
		if i=1 then ide:=args[1] else ide:=mult(args[1],ide,args[3]) fi;
		res:=res+coeff(args[2],D,i)*ide
	od;
	if type(res,polynom(anything,x)) or D=xDF then
		g_expand(res,args[3])
	else
		collect(res,D,g_normal)
	fi
elif nargs=2 then
	# RootOf syntax:
	i:=g_ext([args]);
	subs(g_conversion2,
	 procname(op(subs(g_conversion1,[args])),i))
else
	# the syntax with 5 arguments is no longer supported, it's moved
	# to the procedure below.
	ERROR(`wrong number of arguments`)
fi
end:

# Adapted to Maple 5.4 syntax for []
# This procedure substitutes a for xDF in an operator f in internal form.
# a must be xDF + something of valuation at least -slope
# arguments: a,f,slope,shift,ext
subs_local:=proc()
	local f,i;
	global xDF;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if type(args[2],list) then
	f:=args[2];
	if nops(f)=4 and f[4]=`alg over k((x))` then
		[procname(ramification_of(
		 args[1],f[3],args[nargs]),f[1],args[3..nargs]),op(f[2..4])]
	else
		[seq(procname(args[1],i,args[3..nargs]),i=f)]
	fi
else
	op_with_slope(degree(args[2],xDF),args[3..4],
	 [lift_substitute,args[1..2],args[5]])
fi
end:

# ff is f with xDF + a for xDF substituted. This procedure lifts ff.
# Note: the resulting Laurents may have a higher degree then their
# accuracy. The terms higher than the accuracy are not yet correct.
lift_substitute:=proc(l,d,slope,a,f,ext,ff,acc)
global value_laurent,accuracy_laurent,xDF;
local i,res,n_lifts;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	n_lifts:=acc-accuracy_laurent[l];
	if n_lifts<=0 then RETURN() fi;
	i:=accuracy_laurent[coeff(ff,xDF,0)]+d*slope;
	res:=modm(substitute(a,nm_block(f,i,i+n_lifts,slope,ceil),ext));
	for i from 0 to degree(f,xDF)-1 do
		value_laurent[coeff(ff,xDF,i)]:=coeff(res,xDF,i)
		 +value_laurent[coeff(ff,xDF,i)];
		accuracy_laurent[coeff(ff,xDF,i)]:=
		 accuracy_laurent[coeff(ff,xDF,i)]+n_lifts
	od;
	NULL
end:

#savelib('substitute','`DEtools/diffop/substitute.m`'):
#savelib('subs_local','lift_substitute','`DEtools/diffop/subs_local.m`'):