# $Source: /u/maple/research/lib/LREtools/src/RCS/p_curv,v $
# $Notify: mvanhoei@daisy.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:=Primpart(subs(tau=n*tau,Rem(
	 LinearAlgebra:-LA_Main:-CharacteristicPolynomial(M,tau)
	 mod p,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:

#savelib('order1','AFactor','p_curv','co','`LREtools/p_curv.m`'):