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

# Procedures for users
#    none
# Procedures for use only in the rest of the diffop package:
#    try_factorization (for in DFactor)
# Procedures for use only in this file:
#    GCLD, xDFn_modr, flist, try_factorization2
#

macro(
	adjoint=`DEtools/diffop/adjoint`,
	GCRD=`DEtools/diffop/GCRD`,
	rightdivision=`DEtools/diffop/rightdivision`,
	leftdivision=`DEtools/diffop/leftdivision`,
	mult=`DEtools/diffop/mult`,
	substitute=`DEtools/diffop/substitute`,
	pade2=`DEtools/diffop/pade2`,
	lift_pade2=`DEtools/diffop/lift_pade2`,
	smallest_pade2=`DEtools/diffop/smallest_pade2`,

	compute_modm=`DEtools/diffop/compute_modm`,
	reconstruct=`DEtools/diffop/reconstruct`,
	modulus2=`DEtools/diffop/modulus2`
):

macro(
	flist=`DEtools/diffop/flist`,
	try_factorization=`DEtools/diffop/try_factorization`,
	try_factorization2=`DEtools/diffop/try_factorization2`,
	xDFn_modr=`DEtools/diffop/xDFn_modr`,
	GCLD=`DEtools/diffop/GCLD`
):

macro(
	infosubs=`DEtools/diffop/infosubs`,
	x=`DEtools/diffop/x`,
	DF=`DEtools/diffop/DF`,
	xDF=`DEtools/diffop/xDF`,
	differentiate=`DEtools/diffop/differentiate`,
	eval_laurent=`DEtools/diffop/eval_laurent`
):


# Computes xDF^n mod r.
# Only used in flist.
xDFn_modr:=proc(n,r,ext)
	global x,xDF;
	local a;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	if n<degree(r,xDF) then return xDF^n fi;
	a:=procname(n-1,r,ext);
	eval_laurent(expand(x*differentiate(a)
	 +a*xDF-coeff(a,xDF,degree(r,xDF)-1)*r),ext,xDF)
end:

# Gives the function list for pade2
# r is a local factor
# Only used in try_factorization
flist:=proc(r,order,ext)
	global xDF;
	local res,i,p,j;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	res:=NULL;
	for i from 0 to order do
		p:=xDFn_modr(i,r,ext);
		res:=res,[seq(eval_laurent(coeff(p,xDF,j)),j=0..degree(r,xDF)-1)]
	od;
	[res]
end:

# Input: local factor,maximal order to look for,Newton polygon bound,
#  point,global operator,extra singularities bound,algebraic extension
# Output: a string 'failed', or a global factorization.
# Sep 1994: Added the option that f=[operator,`use adjoint`]. This means that
# r is a right factor of the adjoint of the operator.
try_factorization:=proc(r,max_order,bound,point,f,eb,ext,min_order)
	global x,DF,xDF;
	local i,flm,sr,j,re,fl,exponent,nstep;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	nstep:=4;
	userinfo(3,'diffop',`trying local factor`,infosubs(r),`in the point`,point,
	 `for global factors of order:`,seq(i,i=degree(r,xDF)..max_order));
for i from max(degree(r,xDF),args[8..nargs]) to max_order do
fl:=pade2(flist(r,i,ext),infinity,ext);
if compute_modm(1,try_factorization2,fl,bound,eb,nstep)<>'failed' then
	userinfo(3,'diffop',`possible factor found modulo prime`,modulus2
	 ,`of order`,i);
	# Computing with higher accuracy:
	exponent:=ceil(max(12,op(map(length,indets(f,rational))))/3);
	flm:=compute_modm(exponent,try_factorization2,fl,bound,eb,nstep);
	while flm<>'failed' do
		sr:=expand(convert([seq(
		 flm[1][j+1]*xDF^j,j=0..nops(flm[1])-1)],`+`));
		sr:=expand(reconstruct(sr/lcoeff(sr,indets([x,sr])) mod
		 modulus2^exponent,modulus2^exponent));
		if sr<>'failed' then
			if type(f,list) and f[2]=`use adjoint` then
			# Note: adjoint and substitute(x*DF, ) do not commute
				sr:=adjoint(sr)
			fi;
			sr:=substitute(x*DF,sr,ext);
			if point=infinity then
				sr:=substitute(-x^2*DF,subs(x=1/x,sr))
			else
				sr:=subs(x=x-point,sr)
			fi;
			if not(type(f,list)) then
				re:=GCRD(f,sr);
				if has(re,DF) then
					userinfo(3,'diffop',`found factor`,infosubs(sr));
					return [rightdivision(f,re)[1],re]
				fi;
				userinfo(3,'diffop',`unlikely failed attempt`);
				# to prevent always getting failed attempts:
				nstep:=floor(3/2*nstep)+2
			else
				# Because of f was made monic after localizing,
				# the adjoint of this local operator is not
				# the adjoint of f. To fix this sr must be
				# multiplied with the leading coefficient of
				# the localized f.
				if point=infinity then
					re:=x^(-degree(f[1],DF))
				else
					re:=(x-point)^(-degree(f[1],DF))
				fi;
				sr:=collect(re*lcoeff(f[1],DF)*sr,DF);
				sr:=mult(sr,1/lcoeff(sr,DF),ext);
				re:=GCLD(f[1],sr);
				if has(re,DF) then
					userinfo(3,'diffop',`found factor`,infosubs(re));
					return [re,leftdivision(f[1],re)[1]]
				fi;
				userinfo(3,'diffop',`unlikely failed attempt`);
				nstep:=floor(3/2*nstep)+2
			fi
		fi;
		userinfo(4,'diffop',`higher prime power needed`);
		exponent:=ceil(exponent*1.7);
		flm:=compute_modm(exponent,try_factorization2,fl,bound,eb
		,nstep)
	od
fi od;
	'failed'
end:

# Used by try_factorization
try_factorization2:=proc(fl,bound,eb,nstep)
	global x;
	local i,v,sm,sm_old;
	option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
	v:=lift_pade2(op(fl),nstep);
	sm:=smallest_pade2(op(v));
	while sm<>sm_old or sm[nops(sm)]=0 do
		sm_old:=sm;
		v:=lift_pade2(op(v),nstep);
		sm:=smallest_pade2(op(v));
		for i in sm do if degree(i,x)>bound[nops(sm)-1]+eb then
			return 'failed'
		fi od
	od;
	[sm,v]
end:

GCLD:=proc(f,g)
	global DF;
	local a;
option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`;
if f=0 or g=0 then
	f+g
else
	a:=mult(g,1/lcoeff(g,DF));
	procname(a,leftdivision(f,a)[2])
fi
end:

#savelib('try_factorization','try_factorization2','xDFn_modr','flist','GCLD'):