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

# Procedures for users
#    DFactor
# Procedures for use only in the rest of the diffop package:
#    none
# Procedures for use only in this file:
#    compute_bound, factor_global, factor_minmult1, 
#    same_charclass, other_methods, order2irrat, try_eigenring

# 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`
);
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`,
	nm_block2=`DEtools/diffop/nm_block2`
):

# readlib'ed procedures.
macro(
	GCRD=readlib(`DEtools/diffop/GCRD`),
	LCLM=readlib(`DEtools/diffop/LCLM`),
	adjoint=readlib(`DEtools/diffop/adjoint`),
	eigenring=readlib(`DEtools/diffop/eigenring`),
	endomorphism_charpoly=readlib(`DEtools/diffop/endomorphism_charpoly`),
	l_p=readlib(`DEtools/diffop/l_p`),
	leftdivision=readlib(`DEtools/diffop/leftdivision`),
	make_rightfactor=readlib(`DEtools/diffop/make_rightfactor`),
	mult=readlib(`DEtools/diffop/mult`),
	ratbeke=readlib(`DEtools/diffop/ratbeke`),
	rightdivision=readlib(`DEtools/diffop/rightdivision`),
	subs_local=readlib(`DEtools/diffop/subs_local`),
	substitute=readlib(`DEtools/diffop/substitute`),
	try_factorization=readlib(`DEtools/diffop/try_factorization`),

	modulus2=`DEtools/diffop/modulus2`
):

####################################
#      Global factorizations       #------------------------------------------
####################################			diffop_global

macro(
	DFactor=`DEtools/diffop/DFactor`,
	compute_bound=`DEtools/diffop/compute_bound`,
	factor_global=`DEtools/diffop/factor_global`,
	factor_minmult1=`DEtools/diffop/factor_minmult1`,
	same_charclass=`DEtools/diffop/same_charclass`,
	try_eigenring=`DEtools/diffop/try_eigenring`,
	other_methods=`DEtools/diffop/other_methods`,
	order2irrat=`DEtools/diffop/order2irrat`
):

DFactor:=proc(ff)
	global modulus2,DF,xDF;
	local i,b,f,eb;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
# Hmm, that radical -> RootOf conversion will have to be done somewhere else.
#	f:=convert(ff,RootOf);
#	if f<>ff then RETURN(procname(f,args[2..nargs])) fi;
	f:=collect(subs(g_conversion1,ff),DF,g_normal);
	if degree(f,DF)<=1 then
		RETURN([subs(g_conversion2,ff)])
	fi;
#	modulus:=0; # Should already have been done.
	modulus2:=503;
	eb:=NULL;
	for i in [args[2..nargs]] do if type(i,`=`)
		and op(i)[1]=`bound apparant singularities` then
		eb:=op(i)[2]
	fi od;
	b:=factor_global(f,g_ext([args]),eb);
	b:=subs(g_conversion2,b);
if nops(b)=1 or member(`one step`,{args}) then
	b
elif member(`right factor`,{args}) then
	procname(b[2],args)
else
	[op(procname(b[1],args)),op(procname(b[2],args))]
fi
end:

# Input: a global operator f, the coefficients field and a bound for the
# number of extra singularities.
# Output: either [f] or f factored in 2 factors.
factor_global:=proc(f,ext,extra_bound,what)
	global DF,x,xDF;
	local t,v,i,j,k,l,bound,eb,singularities,min_deg,done_s,all_one,R_left;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if degree(f,DF)<=1 then RETURN([f]) fi;
	if lcoeff(f,DF)<>1 then
		v:=procname(collect(f/lcoeff(f,DF),DF,g_normal),args[2..nargs]);
		RETURN([lcoeff(f,DF)*v[1],op(v[2..nops(v)])])
	fi;
	if degree(f,DF)=2 or (indets(f,name) minus {x,DF} <> {}) then
		# the case with parameters is not yet implemented in the
		# main algorith, so for that case I can only do:
		RETURN(other_methods(f))
		# and for order 2 this is faster.
	fi;
	# point, localized f, algebraic extensions:
	singularities:={[infinity,l_p(f,infinity),ext],seq([i[1],
	 l_p(f,i[1]),[op(i[3]),op(ext)]],i=v_ext_m(g_factors(
	 denom(g_normal(f/lcoeff(f,DF))),x,ext),x))};
  userinfo(5,'diffop',`done computing localizations`);
	for i from 1 to degree(f,DF)-1 do bound[i]:=0 od;
	eb:=-1;
	for k in singularities do
		j:=Newtonpolygon(k[2]);
		v:=[seq(seq(j[i][2]+(l-j[i][1])*(j[i+1][2]-j[i][2])
		/(j[i+1][1]-j[i][1]),l=j[i][1]..j[i+1][1]-1),i=1..nops(j)-1)];
		for i from 1 to degree(f,DF)-1 do
			bound[i]:=bound[i]+degree_ext(k[3],ext)*v[i+1]
		od;
		eb:=eb+degree_ext(k[3],ext)
	od;
	# list of bounds for the coefficients of possible global factors
	bound:=[seq(-bound[degree(f,DF)-j]+j*eb,j=1..degree(f,DF)-1)];
	if nargs<3 then
		eb:=compute_bound(singularities,ext,f);
		if eb>20 then
			# userinfo(1,'diffop',
			# `Computed bound number of extra singularities is`,
			# eb,`Use option bound=computed to use this bound.`
			# ,`Now using a lower bound`);
			# eb:=15
			v:=other_methods(f);
			if nops(v)<>1 then
				RETURN(v)
			fi
		fi
#	elif extra_bound=computed then
#		eb:=compute_bound(singularities,ext,f)
	elif type(extra_bound,integer) then
		eb:=extra_bound
	else
		ERROR(`wrong arguments`)
	fi;
	userinfo(5,'diffop',`bounds are`,bound,eb);
	all_one:=true;
	done_s:=[]; # The singularities we have considered
while singularities<>{} do
	min_deg:=min(seq(degree_ext(i[3],ext),i=singularities));
	for i in singularities do if degree_ext(i[3],ext)=min_deg then
		singularities:=singularities minus {i};
		v:=factor_op(i[2],`all right factors`,i[3]);
		done_s:=[op(done_s),[op(i),v]];
		for j in v do if skipped_factors(j)=0 then
			userinfo(3,'diffop',cat(`Minimum multiplicity 1 in x=`,i[1]));
			RETURN(factor_minmult1(bound,i[1],f,i[2],eb,j,i[3]))
		fi od;
		all_one:=evalb(all_one and nops(v)=1 and degree(v[1],xDF)=1)
	fi od;
od;
	userinfo(3,'diffop',`Minimum multiplicity >1`);
if all_one then # f has only 1 exponential part
	RETURN(other_methods(f))
fi;
	# if degree(f,DF)<=3 then bug() fi;
	v:=try_eigenring(subs(g_conversion2,f));
	if nops(v)>1 then
		RETURN(v)
	fi;
	all_one:=true;
	for i in done_s do
		all_one:=evalb(all_one and nops(i[4])=1);
		for j in i[4] do
			t:=try_factorization(j,floor(degree(f,DF)/2),bound,
			 i[1],f,eb,i[3]);
			if t<>'failed' then RETURN(t) fi
		od;
		v:=factor_op(adjoint(i[2],i[3]),`all right factors`,i[3]);
		for j in v
		do
			t:=try_factorization(j,degree(f,DF)-1,bound,
			 i[1],[f,`use adjoint`],eb,i[3]);
			if t<>'failed' then RETURN(t) fi
		od;
		# Now try the same, but searching for higher order factors:
		for j in i[4] do
			t:=try_factorization(j,degree(f,DF)-1,bound,
			 i[1],f,eb,i[3],floor(degree(f,DF)/2)+1);
			if t<>'failed' then RETURN(t) fi
		od;
	od;
if all_one then # up to conjugation only 1 exp part in every singularity
	userinfo(3,'diffop',`Trying with algebraic extensions`);
	for i in done_s do if degree(i[4][1],xDF)>1 then
		j:=factor_op(i[4][1],`alg factor`,i[3])[1];
		if degree(j[3],x)=1 then
			k:=j[1]
		elif degree_ext(j,i[3])=1 then
			next
		else
			k:=make_rightfactor(i[4][1],j,g_ext(j));
			# the slope:
			l:=-valuation_laurent(coeff(i[4][1],xDF,0),0)
			 /degree(i[4][1],xDF);
			k:=faster_riccati_split(i[4][1],rightdivision(i[4][1]
			,k,g_ext(j),l),k,g_ext(j),l)[2]
		fi;
		t:=try_factorization(k,degree(f,DF)-1,bound,i[1],f,0,g_ext(j));
		if t<>'failed' then RETURN(t) fi;
		v:=factor_op(adjoint(i[2],g_ext(j)),`all right factors`,ext);
		R_left:=0;
		for l in v while R_left=0 do if same_charclass(l
		,adjoint(k,g_ext(j)),g_ext(j))
			then R_left:=l
		fi od;
		# if R_left=0 then bug() fi;
		# Now try if R gives a left hand global factor
		t:=try_factorization(R_left,degree(f,DF)-1,bound,i[1]
		 ,[f,`use adjoint`],0,g_ext(j));
		if t<>'failed' then RETURN(t) fi
	fi od;
fi;
	# userinfo(1,'diffop',`factorization of`,infosubs(f),`may be incomplete`);
	# [f]
	other_methods(f)
end:


# R is a local factor with multiplicity 1.
# The result of the procedure is a factorization. If f is returned unfactored
# it is irreducible (if the specified bound is correct).
factor_minmult1:=proc(bound,point,f,f_local,eb,R,ext)
	global DF,x,xDF;
	local t,w,fl,slope,k,l,R_left,i;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if g_ext(f)<>g_ext([args]) then
		t:=try_eigenring(subs(g_conversion2,f));
		if nops(t)>1 then RETURN(t) fi
	fi;
	fl:=floor(degree(f,DF)/2);
	t:=try_factorization(R,fl,bound,point,f,eb,ext);
	if t<>'failed' then RETURN(t) fi;
	w:=factor_op(adjoint(f_local,ext),`all right factors`,ext);
	R_left:=0;
	for i in w while R_left=0 do if same_charclass(i,adjoint(R,ext),ext)
		then R_left:=i
	fi od;
	if R_left=0 then bug() fi;
	# Now try if R gives a left hand global factor
	t:=try_factorization(R_left,fl,bound,point,[f,`use adjoint`],eb,ext);
	if t<>'failed' then RETURN(t) fi;
	# Now try more terms:
	t:=try_factorization(R,degree(f,DF)-1,bound,point,f,eb,ext,fl);
	if t<>'failed' then RETURN(t) fi;
	t:=try_factorization(R_left,degree(f,DF)-1,bound,point,[f,`use adjoint`],eb
	 ,ext,fl);
	if t<>'failed' then RETURN(t) fi;
	userinfo(5,'diffop',`trying if algebraic extensions are needed`);
#	if igcd(op(map(degree,w,xDF)))=1 then RETURN([f]) fi;
#	if igcd(seq(`if`(skipped_factors(i)=0,degree(i,xDF),0),i=w))=1 then
#		RETURN([f])
#	fi;
	if igcd(degree(f,DF),degree(R,xDF))=1 or g_ext(f)<>g_ext([args]) then
		RETURN([f])
	fi;
	slope:=-valuation_laurent(coeff(R,xDF,0),0)/degree(R,xDF);
	k:=op(factor_op(R,`alg factor`,ext));
	if degree(k[3],x)=1 then
		l:=k[1]
	elif degree_ext(k,ext)=1 then
		RETURN([f])
	else
		l:=make_rightfactor(R,k,g_ext(k));
		l:=faster_riccati_split(R,rightdivision(R,l,g_ext(k),slope)[1]
		 ,l,g_ext(k),slope)[2]
	fi;
	t:=try_factorization(l,fl,bound,point,f,0,g_ext(k));
	if t<>'failed' then RETURN(t) fi;
	[f]
end:

# f and g irreducible.
# Output: true if f and g have the same exponential parts.
# Adapted to Maple 5.4 syntax for []
same_charclass:=proc(f,g,ext)
	global x,xDF;
	local gg,ff,r,c,d;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if degree(f,xDF)<>degree(g,xDF) then RETURN(false) fi;
	if degree(f,xDF)=1 then RETURN(type(nt(f-g,1),integer)) fi;
	r:=map(Newtonpolygon,[f,g],`include the Newton polynomials`);
	if r[1][1][3]=0 then
		c:=r[1][1][4];
		d:=degree(c,x);
		c:=(coeff(c,x,d-1)-coeff(r[2][1][4],x,d-1))/d;
		if not type(c,integer) then RETURN(false) fi;
		r:=[[[op(r[1][1][1..3]),g_expand(subs(x=x-c,r[1][1][4]),ext)]
		 ,op(r[1][2..nops(r[1])])],r[2]]
	fi;
	if r[1]<>r[2] then RETURN(false) fi;
	gg:=factor_op(g,`alg factor`,ext)[1];
	r:=nt(gg[1],1);
	ff:=subs_local(2*xDF-r,ramification_of(f,gg[3],gg[2]),ldegree(r,x),0,gg[2]);
	ff:=Newtonpolygon(ff,`include the Newton polynomials`);
	if ff[1][3]<>0 then RETURN(false) fi;
	ff:=factors(collect(ff[1][4]/lcoeff(ff[1][4],x),x,g_normal));
	for r in ff[2] do if degree(r[1],x)=1 and
	 type(coeff(r[1],x,0)/coeff(r[1],x,1),integer) then RETURN(true)
	fi od;
	false
end:

# try ratbeke, try_eigenring and order2irrat
other_methods:=proc(f)
	global DF,g_conversion2,infolevel;
	local v;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	v:=subs(g_conversion2,f);
	if v<>f then RETURN(procname(v,args[2..nargs])) fi;
	userinfo(3,'diffop',`Looking for first order factors`);
	if not has({args},`second attempt`) then
		v:=ratbeke(f)
	elif degree(f,DF)=2 then
		RETURN(order2irrat(f))
	else
		v:=ratbeke(adjoint(f))
	fi;
	if v=[] and has({args},`second attempt`) then
		# Play our last card:
		v:=try_eigenring(f);
		if nops(v)=1 and degree(f,DF)>3 then
			if _Envdiffopwarnings=true and not(type(infolevel[diffop],integer)) then
#			# Make sure that this message reaches to the user:
#				lprint(`factorization of`,infosubs(f),
#				`may be incomplete`)
                        
			    userinfo(1,'hints',`factorization of`
			    ,infosubs(f),`may be incomplete`)
			fi;
			userinfo(1,'diffop',`factorization of`
			,infosubs(f),`may be incomplete`)
		fi
	elif v=[] then
		v:=procname(f,`second attempt`)
	else
		v:=collect(v[1][2]-diff(v[1][1][1],x)/v[1][1][1],DF,g_normal);
		if has({args},`second attempt`) then
			v:=adjoint(v);
			v:=[v,leftdivision(f,v)[1]]
		else
			v:=[rightdivision(f,v)[1],v]
		fi
	fi;
	v
end:

# this is something similar to the try_eigenring approach
# but it seems to be faster in case the groundfield contains
# trancendental numbers.
# Assumptions: f is of order 2 and is irreducible over the base field.
order2irrat:=proc(f)
global DF;
local i,j,c2;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	i:=coeff(f,DF,1)/coeff(f,DF,2)/2;
	# Make sure coeff(f,DF,1)=0:
	if i<>0 then RETURN([seq(substitute(DF+i,j),j=
		procname(substitute(DF-i,f),i))])
	fi;
	# rational solutions of the second symmetric power:
	# j:=ratbeke(readlib(`DEtools/diffop/symmetric_power`)(f,2),rational);
	# Could also have used ratbeke( ... , type=... ) instead of using substitute
	userinfo(3,'diffop',`trying to factor over a quadratic extension`);
	j:=readlib(`DEtools/diffop/symmetric_power`)(f,2);
	j:=DEtools[ratsols]([seq(coeff(j,DF,i),i=0..3)],0,x);
	if nops(j)>1 then
		ERROR() # bug
	elif nops(j)=0 then
		userinfo(3,'diffop',`verified irreducibility of`,infosubs(f));
		RETURN([f])
	fi;
	j:=j[1];
	c2:=g_normal(diff(j,x)^2-2*diff(j,x,x)*j-4*coeff(f,DF,0)/coeff(f,DF,2)*j^2);
	if has(c2,x) then
		userinfo(3,'diffop',`verified irreducibility of`,infosubs(f));
		[f]
	else
		userinfo(3,'diffop',`factorized `,infosubs(f));
		j:=g_normal(1/2*(2*RootOf(_Z^2-c2/4)+diff(j,x))/j);
		[collect(coeff(f,DF,2)*(DF+j),DF),DF-j]
	fi
end:

# Quick hack, this could be done somewhat more efficient.
# RootOf syntax

# Idea: set eigenring(of the factors) := [1] (first make sure that it's true).
try_eigenring:=proc(f)
	global DF,g_conversion1,x;
	local k,i,m,l,ext;
	options remember, 
	    `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	ext:=g_ext(f);
	k:=eigenring(f);
	if nops(k)>1 then
		m:=min(op({seq(degree(i,DF),i=k)} minus {0}));
		k:=[seq(`if`(degree(i,DF)=m,i,NULL),i=k)];
		m:=min(seq(length(i),i=k));
		# k := the minimal element from k:
		for i in k do if length(i)=m then k:=i; break fi od;
		m:=g_factors(subs(g_conversion1,endomorphism_charpoly(f,k))
		 ,x,ext);
		l:=min(seq(degree(i[1],x),i=m));
		for i in m do if degree(i[1],x)=l then l:=i[1]; break fi od;
		k:=GCRD(k-g_zero_of(l,x,'ext'),f);
# shouldn't always do that:
		l:=rightdivision(f,k);
		# if l[2]<>0 then bug() fi;
		[l[1],k]
	else
		[f]
	fi
end:



# Bound for the number of extra singularities in a 1st order right factor
compute_bound:=proc(singularities,ext,f)
	local v,i,res,j,ma,c;
	global g_conversion2,x,xDF,DF;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	userinfo(5,'diffop',`computing number of extra singularities bound`);
	res:=0;
for i in singularities do
	# compute the generalized exponents:
	v:=factor_op(i[2],'semireg',i[3]);
	ma:=-infinity;
	for j in v do if degree(j[3],x)=1 then
		# c = - an exponent. Bound= - sum exponents
		c:=coeff(min(solve(nt(j[1],1)))+j[4],x,0);
		# Now compute the trace of this algebraic number divided by
		# the degree of this number over Q:
		c:=subs(g_conversion2,[c,{op(g_ext(c))}]);
		c:=evala(Trace(op(c)))/degree_ext(c,0);
		# if not type(c,rational) then bug() fi;
		ma:=max(-c,ma)
	fi od;
	if ma=-infinity then # there are no 1st order factors
		res:=0;
		break
	fi;
	res:=res+degree_ext(i[3],ext)*ma
od;
	res:=floor(res);
	userinfo(5,'diffop',`bound for extra singularities is`,res);
	max(0,res)
end:



#savelib('DFactor','other_methods','order2irrat', \
  'factor_global','try_eigenring','factor_minmult1','same_charclass', \
  'compute_bound','`DEtools/diffop/DFactor.m`'):