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

# 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=readlib(`DEtools/diffop/GCRD`),
	rightdivision=readlib(`DEtools/diffop/rightdivision`),
	leftdivision=readlib(`DEtools/diffop/leftdivision`),
	mult=readlib(`DEtools/diffop/mult`),
	substitute=readlib(`DEtools/diffop/substitute`),
	pade2=readlib(`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'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`
):


# 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.`;
readlib(`DEtools/diffop/compute_modp`):
	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:=readlib(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','`DEtools/diffop/try_factorization.m`'):