# $Source: /u/maple/research/lib/LREtools/src/RCS/p_curv,v $
# $Notify: mvanhoei@scg.math.uwaterloo.ca

macro(
	p_curv=`LREtools/p_curv`,
	co=`LREtools/p_curv/co`,
	AFactor=`LREtools/p_curv/AFactor`,
	order1=`LREtools/p_curv/order1`
):

# Compute the factors of the characteristic polynomial of
# the p-curvature of a difference operator.
#
# Mark van Hoeij, 1998.

# L in C[x,tau]
# args[5..-1] can be used to specify additional field extensions.
#
# If p=list of primes then just check if an order1 factor could
# exist, then the output is true/false.
#
# If p=prime then compute the characteristic polynomial of the
# p-curvature and factor it over the algebraic closure.
p_curv:=proc(L,x,tau,p)
	local M,N,n,i,c,b,F;
	option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
	if type(p,list) then
		for i in p do
			M:=traperror(procname(L,x,tau,i,'test'));
			if M=false then return M fi;
		od;
		return true
	fi;
	M:=ReduceField([L,indets([args],{'RootOf','radical','nonreal'})]
	 ,indets([args[2..-1]],'name')) mod p;
	M, F:=Primpart(M[1],tau) mod p, indets(M,'RootOf');
	if args[-1]='test' and degree(M,tau)<degree(L,tau) then
		return true
	fi;
	b:=x^p-x-c,x;
	n:=Factors(lcoeff(M,tau)) mod p;
	n:=mul((Rem(mul(subs(x=x+i,n[1]),i=0..p-1),b) mod p)^n[2]
	 ,n=[[n[1],1],op(n[2])]);
	N:=degree(M,tau);
	M:=Rem(M,b) mod p;
	N:=[seq(co(Expand(subs(x=x+i,M)) mod p,tau,N),i=0..p-1)];
	for i from 1 to p do
		if i=1 then
		   M := N[1];
		else
		   M := LinearAlgebra:-MatrixMatrixMultiply(M,N[i],
		      'inplace'=false,'outputoptions'=[]);
		   # the Matrix has no special options or indexing
		   # functions, so this is okay
		   M := map(`mod/Rem`,M,b,p);
		fi;
	od;
	M := LinearAlgebra:-LA_Main:-CharacteristicPolynomial(M,tau) mod p;
	M:=Primpart(subs(tau=n*tau,Rem( M ,b) mod p),tau) mod p;
	if has(M,x) then
		userinfo(1,'LREtools',`Found wrong p-curvature`);
		error ""
	fi;
	userinfo(4,'LREtools',cat(`Characteristic polynomial of the `,
	 p,`-curvature is`),M);
	M:=subs(c=x,M);
	if args[-1]=`dont factor` then
		return M
	fi;
	M:={seq(`if`(has(i,tau),i,NULL)
	,i=(Factors(M,op(F)) mod p)[2])};
	if args[-1]='test' and member(1,[seq(degree(i[1],tau),i=M)]) then
		true
	else
		M:=order1(L,x,tau,p,M,args[-1]);
		if args[-1]='test' then
			member(1,[seq(degree(i[1],tau),i=M)])
		else
			M
		fi
	fi
end:

# Companion matrix for a non-monic polynomial.
co:=proc(M,tau,n)
	local i,j;
	option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
	Matrix(n,n,[seq([seq(`if`(i=j+1,lcoeff(M,tau),
	 `if`(j=n,-coeff(M,tau,i-1),0)),j=1..n)],i=1..n)])
end:


# The following procedure can be used to give
# upper bound the number of first order factors, and
# lower bounds for the degrees of the field extensions:
#
# Give the order 1 factors mod p, with the degrees of the
# necessary field extensions. If factors have a multiplicity
# just give them multiple times.
order1:=proc(L,x,tau,p,D)
	local n,C,C1,i,j;
	# The n is how much the order dropped mod p.
	option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
	n:=degree(L,tau)-add(degree(i[1],tau)*i[2],i=D);
	C1:={seq(`if`(degree(i[1],tau)=1,i,NULL),i=D)};
	C:=[seq([op(i),AFactor(i[1],x,tau,p,args[6..-1])],i=D minus C1)];
	C:=[seq(`if`(nops(i[3])=degree(i[1],tau),i,NULL),i=C)];
	[`if`(n>0,[tau-infinity,1,n],NULL),seq([i[1],1,i[2]],i=C1),
	 seq(seq([j,degree(i[1],tau),i[2],i[1]],j=i[3]),i=C)]
end:


# If a polynomial in 2 variables is a product of factors of deg 1 in tau,
# and irreducible over the base field, then factor it over closure of
# the base field.
AFactor:=proc(f,x,tau,p)
	local i,n,R,S;
	option `Copyright (c) 1998 Waterloo Maple Inc. All rights reserved.`;
	n:=degree(f,tau);
if n=1 or {seq(irem(degree(f,i),n),i=[x,{x,tau}])}<>{0} or member(1
 ,{seq(irem(i[2],2),i=(Sqrfree(Discrim(f,tau) mod p) mod p)[2])}) then
	R:=f # f is not a product of deg 1 factors.
elif args[-1]='test' then
	# Just check if a 1st order factor could exist, which is
	# the case.
	R:=tau$n
else
	R:=op(indets(f,'RootOf'));
	if R=NULL then
		R:=RootOf(Randprime(n,x) mod p) mod p;
		R:=seq(i[1],i=(Factors(f,R) mod p)[2])
	else
		S:=RootOf(Randprime(n*degree(op(1,R),_Z),x) mod p) mod p;
		R:=seq(i[1],i=(Factors(subs(R=(Roots(subs(_Z=x,op(1,R))
		,S) mod p)[1][1],f),S) mod p)[2])
	fi
fi;
	[seq(collect(i/lcoeff(i,tau),tau,Normal) mod p,i=[R])]
end:

#########################################################
#  Extract local information at x=infinity from p_curv  #
#########################################################

macro(
	expansion=`LREtools/cnd_count/expansion`,
	cnd_count=`LREtools/cnd_count`,
	x   = _Env_LRE_x,
	tau = _Env_LRE_tau,
	NO_RAT = 3  # See hsolsR
):

# F is the p_curvature polynomial in x and tau.
# R is the c*x^n*(1+d*(1/x))
# We want the number of roots with c*x^n+lower terms
expansion:=proc(F,R,p,deg_R)
	local l,n,i;
	l:=lcoeff(F,tau);
	if l<>1 then
		n:=degree(F,tau);
		procname(Expand(tau^n+
		 add(coeff(F,tau,i)*l^(n-1-i)*tau^i,i=0..n-1)) mod p
		 ,collect(R*l,x,Normal) mod p,p)
	elif ldegree(R,x)<degree(R,x)-1 then
		n:=degree(R,x);
		procname(F,add(coeff(R,x,l)*x^l,l=n-1..n),p)
	elif ldegree(R,x)<0 then
		0
	elif nargs=3 then
		procname(Expand(subs(tau=tau+R,F)) mod p
		 ,0,p,degree(R,x)-2,1)
	elif deg_R<0 then
		ldegree(F,tau)
	else
		l:=Expand(subs(tau=n*x^deg_R,F)) mod p;
		l:=subs(n=x,lcoeff(l,x));
		if not has(l,x) then
			return 0
		fi;
		n:=Factors(l,op(indets(F,RootOf))) mod p;
		add(degree(l[1],x)*procname(Expand(`mod/ReduceField`(
		 subs(tau=tau+RootOf(l[1],x)*x^deg_R,F),{x,tau},p)) mod p
		,0,p,deg_R-1,1),l=n[2])
	fi
end:


# cnd[-1] = max number of hypergeomsols with that local type.
cnd_count:=proc(L,p,cnds,TO_DO)
	local v,w,cnd,delta,pc,res,i,f;
	# return cnds;

	v:=indets(L,'nonreal');
	if v<>{} then v:=v union {I} fi;
	v:=[L,op(indets(L,{'RootOf','radical','name'})),op(v)];
	# Could generate an error:
	w:=traperror(`mod/ReduceField`(v,{x,tau},p));
	if w=lasterror then
		return cnds
	fi;
	delta:=degree(L,tau)-degree(w[1],tau) + ldegree(w[1],tau)-ldegree(L,tau);
	if delta >= max(seq(cnd[-1],cnd=cnds)) then
		return cnds
	fi;
	pc:=p_curv(w[1],x,tau,p,`dont factor`);
	res:=NULL;
	for cnd in cnds do
		if TO_DO=NO_RAT and cnd[-1]=1 then
			# are not searching for rational hypergeomsols.
			next
		fi;
		if delta >= cnd[-1] then
			res:=res,cnd;
			next
		fi;
		f:=Normal(cnd[1]^p*x^(-cnd[2])*(1-(cnd[3]^p-cnd[3])/x));
		f:=traperror(`mod/ReduceField`([subs(
		 {seq(v[i]=w[i],i=2..nops(v))},f),pc],{x,tau},p));
		if f=lasterror or f[1]=0 then
			res:=res,cnd;
			next
		fi;
		i, f := op(f);
		i:=expansion(f,i,p);
		if i=0 or (i=1 and TO_DO = NO_RAT) then next fi;
		res:=res, subsop(-1=min(cnd[-1],i+delta),cnd)
	od;
	[res]
end:

#savelib('order1','AFactor','p_curv','co'):
#savelib('expansion','cnd_count'):