# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/DFactor,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca $

# 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(
	infosubs=`DEtools/diffop/infosubs`,
	modulus=`DEtools/diffop/modulus`,
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	xDF=`DEtools/diffop/xDF`,

	g_conversion1=`algcurves/g_conversion1`,
	g_conversion2=`algcurves/g_conversion2`,
	degree_ext=`algcurves/degree_ext`,
	g_zero_of=`algcurves/g_zero_of`,
	v_ext_m=`algcurves/v_ext_m`,

	g_ext=`DEtools/diffop/g_ext`,
	g_factors=`DEtools/diffop/g_factors`,
	g_normal=`DEtools/diffop/g_normal`,
	g_expand=`DEtools/diffop/g_expand`,

	valuation_laurent=`DEtools/diffop/valuation_laurent`,
	nt=`DEtools/diffop/nt`,

	Newtonpolygon=`DEtools/diffop/Newtonpolygon`,
	factor_op=`DEtools/diffop/factor_op`,
	ramification_of=`DEtools/diffop/ramification_of`,
	skipped_factors=`DEtools/diffop/skipped_factors`,
	substitute=`DEtools/diffop/substitute`,

	faster_riccati_split=`DEtools/diffop/faster_riccati_split`,

	GCRD=`DEtools/diffop/GCRD`,
	adjoint=`DEtools/diffop/adjoint`,
	eigenring=`DEtools/diffop/eigenring`,
	endomorphism_charpoly=`DEtools/diffop/endomorphism_charpoly`,
	l_p=`DEtools/diffop/l_p`,
	leftdivision=`DEtools/diffop/leftdivision`,
	make_rightfactor=`DEtools/diffop/make_rightfactor`,
	mult=`DEtools/diffop/mult`,
	ratbeke=`DEtools/diffop/ratbeke`,
	rightdivision=`DEtools/diffop/rightdivision`,
	subs_local=`DEtools/diffop/subs_local`,
	substitute=`DEtools/diffop/substitute`,
	try_factorization=`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,ext;
  option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. \
  Author: Mark van Hoeij`;

	l_p();
	ext:=g_ext([args]);
	f:=collect(ff,DF,g_normal);
	if degree(f,DF)<=1 then
		return [f]
	elif not has(ff,x) then
		# constant coeffs.
		b:=polytools['splits'](ff,DF,'i',{args});
		b:=[b[1],seq(i[1]$i[2],i=b[2])];
		return [b[1]*b[2],op(b[3..-1])]
	fi;
	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:=evala(Primpart(f,DF));
	b:=gcd(subs(x=0,b), subs(x=1,b));
	if has(b,DF) then
		b:=GCRD(b,f)
	fi;
    if has(b,DF) then
	# Found a right-hand factor with constant coeffs.
	i:=rightdivision(f,b);
	if i[2]<>0 then
		error "failed"
	fi;
	b:=[i[1],b]
    else
	b:=factor_global(subs(g_conversion1,f),ext,eb,
	 {args} intersect {`skip left`} );
	b:=subs(g_conversion2,b);
    fi;
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;
 options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. \
 Author: Mark van Hoeij`;

	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 or not type(extra_bound,integer) 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
	else
		eb:=extra_bound
	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,eb,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`],eb,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;
	options remember,
 `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. \
 Author: Mark van Hoeij`;

	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 error "should not happen"  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,eb,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:=roots(collect(ff[1][4]/lcoeff(ff[1][4],x),x,g_normal),x);
	for r in ff do if type(r[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 `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
			    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]]
		elif _Env_diffop_skipleft=true then
			v:=[0,v]
		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;
 options remember,`Copyright (c) 1997 Waterloo Maple Inc. All rights reserved. \
 Author: Mark van Hoeij`;
	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(`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:=`DEtools/diffop/symmetric_power`(f,2);
	j:=DEtools[ratsols]([seq(coeff(j,DF,i),i=0..3)],0,x,'skiptest');
	if nops(j)>1 then
		error "should not happen" # 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;
	option `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);
		[`if`(_Env_diffop_skipleft=true,0,rightdivision(f,k)[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'):