# $Source: /u/maple/research/lib/LREtools/src/RCS/g_p,v $
# $Notify: mvanhoei@daisy.uwaterloo.ca
#
# The code in this file computes local information at a singularity,
# and contains a few other tools
#
# Calling sequence:
#   g_p(L,x,tau,point)
#
# Here point is in C union {infinity}.
#
# L is a difference operator,  L=sum(a.i*tau^i,i=0..n), where the a.i
# are rational functions in x. The automorpism tau acts as tau(f(x))=f(x+1).
#
# See "Finite Singularities and Hypergeometric Solutions of Linear Recurrence
#      Equations", http://math.fsu.edu/~hoeij/papers.html
# for a definition of g_p.

macro(	g_p=`LREtools/g_p`,
	do_it=`LREtools/g_p/do_it`,
	extend_right=`LREtools/g_p/ER`,
	normal_poly=`LREtools/g_p/normal_p`,
	Newtonpolygon=`LREtools/g_p/NP`,
	g_factors=`LREtools/g_p/factors`,
	symprod1=`LREtools/g_p/symprod1`,
	cnd=`LREtools/g_p/cnd`,
	int_difference=`LREtools/g_p/int_dif`,
	singularities=`LREtools/g_p/sing`
): # :)

# L in C[x,tau]
g_p:=proc(L,x,tau,p)
	local i,M,n,a0,an,r0,rn,g_min,g_max,left,right,lr,rl,LL,j,acc,rang,res,ext;
	ext:=indets([args],{'RootOf','radical','nonreal'});
	if ext<>{} then
		Normalizer:=`evala/Normal`
	else
		Normalizer:=normal
	fi;
	_Env_LRE_x:=x; # used by normal_poly
	if not has(L,tau) then
		error `wrong input`
	elif Normalizer(lcoeff(L,tau))*Normalizer(coeff(L,tau,0))=0 then
		i:=ldegree(L,tau);
		if i<>0 then
			userinfo(1,'LREtools',L,`has factor`,tau^i)
		fi;
		return procname(collect(subs(x=x+i,L/tau^i)
		 ,tau,normal_poly),args[2..-1])
	elif p=infinity then
		return cnd(L,x,tau,ext)
	fi;
	M:=collect(subs(x=x+p,L),tau,normal_poly);
	if ldegree(M,x)>0 then
		M:=collect(M/x^ldegree(M,x),tau,normal_poly)
	fi;
	n:=degree(M,tau);
	a0:=coeff(M,tau,0);
	an:=lcoeff(M,tau);
	r0,rn:=seq([seq(`if`(type(j[1],integer),j,NULL)
	 ,j=frontend(roots,[i,x]))],i=[a0,an]);
	g_min,g_max:=-add(i[2],i=rn),add(i[2],i=r0);
	userinfo(4,'LREtools',`lower and upper bound for g_p`,g_min,g_max);
	# left basis: left-n..left-1
	left:=min(seq(i[1],i=rn),seq(i[1],i=r0))+n;
	# right basis: right+1..right+n
	right:=max(seq(i[1],i=rn),seq(i[1],i=r0));
	if {g_min,g_max}={0} then
		userinfo(1,'LREtools',p,`is not a singularity`);
		return {0},n,n,[]$2,left,right
	elif n=1 then # g_p and bound can easily be computed
		lr:=[0];
		for i from left-1 to right do
			lr:=[op(lr),lr[-1]+add(`if`(j[1]=i,j[2],0),j=r0)
					  -add(`if`(j[1]=i,j[2],0),j=rn)]
		od; # lr[-1] = g_min+g_max = the only valuation growth.
		return {lr[-1]},1,1,lr$2,left,right
	fi;
	# Now determine gp and the bounds.
	for i from 1 to n do
		# starting bases left and right.
		lr[i]:=[0$(i-1),1,0$(n-i)]; # must be extended to the right.
		rl[i]:=lr[i]    # rl (r to l) must be extended to the left.
	od;

for j from 1 to 3 do
    for i from 1 to 3 do
	# Don't go higher than 46327 because that would be too high for modp1.
	_Env_LRE_p:=nextprime(rand(20000..46326)());

	# Remove parameters by random substitution in Fp,
	# and reduce RootOf's to one single RootOf by factorization:
	userinfo(2,'LREtools',`reducing mod`,_Env_LRE_p);
	Normalizer:=proc(a) `mod/Expand`(a,_Env_LRE_p) end:
	rang:=proc(n,M) local r; Gausselim( Matrix(n,n,M)
	 ,'r') mod _Env_LRE_p; r end:
	res[i]:=traperror(do_it(ReduceField(M,{tau,x}) mod _Env_LRE_p
		,x,tau,p,n,g_min,g_max,left,right,rang,i))
	od;
    if res[1]=res[2] and res[1]=res[3] then
	if res[1]=lasterror then error res[1] else return res[1] fi
    fi
od;
	error "failed"

	# Characteristic 0:
	# rang:=proc(n,M) local r; linalg[gausselim](matrix(n,n,M),'r'); r end:
	# do_it(M,x,tau,p,n,g_min,g_max,left,right,rang)
end:


do_it:=proc(M,x,tau,p,n,gmi,gma,left,right,rang,ii)
	local i,j,LL,a0,an,r0,rn,lr,rl,g_min,g_max,acc;
	g_min, g_max := gmi, gma;
	# Now determine gp and the bounds.
	for i from 1 to n do
		# starting bases left and right.
		lr[i]:=[0$(i-1),1,0$(n-i)]; # must be extended to the right.
		rl[i]:=lr[i]    # rl (r to l) must be extended to the left.
	od;
	LL:=add(coeff(M,tau,i)*tau^(n-i),i=0..n);
	acc:=g_max-g_min+1;
	seq(extend_right(lr,acc,subs(x=x+i+left-n,M),n,x,tau),i=0..right-left+n):
	# Since we computed with accuracy acc, if a valuation (= ldegree) is >acc
	# we can at best take acc as a lower bound for the valuation, but not higher.
	a0:=[seq(min(acc,seq(ldegree(lr[i][j],x),i=1..n)),j=1..nops(lr[1]))];
	r0:=min(op(a0[-n..-1]));
	# Let lr be the rank of the left-to-right map E_{p,r}:
	lr:=rang(n,[seq([seq(coeff(j,x,r0),j=lr[i][-n..-1])],i=1..n)]):
	a0:=[seq(i+g_min,i=a0)]; # set the first n entries to 0.
	g_min:=g_min+r0; # Now g_min is the {minimum valuation growth} g_{p,r}
if lr<n then
	seq(extend_right(rl,acc,subs(x=x-i+right,LL),n,x,tau),i=0..right-left+n):
	an:=[seq(min(acc,seq(ldegree(rl[i][-j],x),i=1..n)),j=-nops(rl[1])..-1)];
	rn:=min(op(an[1..n]));
	# Let rl be the rank of the right-to-left map E_{p,l}:
	rl:=rang(n,[seq([seq(coeff(j,x,rn),j=rl[i][-n..-1])],i=1..n)]):
	if rl+lr>n then
		# At an essential singularity, the composition of the linear maps
		# E_{p,r} and E_{p,l} is zero (last lemma in section 4 in my
		# paper). Hence the sum of the ranks is at most n.
		error "sum of ranks too high"
	fi;
	an:=[seq(i-g_max,i=an)]; # set the first n entries to 0.
	g_max:=g_max-rn # Now g_max is the {\em maximum valuation growth} g_{p,l}
else
	rl,an,g_max := n,[seq(i-g_min,i=a0)],g_min
fi;
if ii=3 then
	userinfo(3,'LREtools',`p =`,p,
	 `if`(lr<n,`is an essential singularity`,
	 cat(`if`(g_min=0,`is an `,`is a semi-`),`apparant singularity`)),
	  `[g_{p,r}, g_{p,l}, rank(E_{p,r}), rank(E_{p,l}] =`
	 , [g_min  , g_max  ,    lr        ,     rl      ]);
fi:
	{seq(i,i=g_min..g_max)},lr,rl,a0,an,left,right
end:



# lr[1..n] is a basis
# acc is: compute modulo x^acc
# extend basis one step to the right.
extend_right:=proc(lr,acc,LL,n,x,tau)
	local ac,an,ld,co,i,j,L;
	ac:=acc-degree(lr[1][1],x);
	# Normalizer=Expand mod p, so this should make ldegree correct,
	# because the parameters have been removed.
	L:=collect(LL,tau,Normalizer);
	an:=coeff(L,tau,n);
	ld:=ldegree(an,x);
	if ld>0 then
		an:=collect(an/x^ld,x)
	fi;
	if ac-ld < 1 then error "took bad prime" fi;
	# Note: for series it is important that Normalizer
	# has the correct value.
	an:=Normalizer(convert(series(1/an,x=0,ac),polynom));
	# coefficients of the operator made monic.
	co:=[seq(Normalizer(convert(series(an*coeff(L,tau,i),x=0,ac),polynom)),i=0..n-1)];
	for i from 1 to n do
		lr[i]:=[`if`(ld=0,op(lr[i]),seq(j*x^ld,j=lr[i])),Normalizer(
		 convert(series(-add(lr[i][-n+j-1]*co[j],j=1..n),x=0,acc),polynom))]
	od;
	NULL
end:

normal_poly:=proc(A)
	collect(A,_Env_LRE_x,
	 `if`(indets(A,{'RootOf','name','radical','nonreal'}) minus {x}={},NULL,Normalizer))
end:


# The code below computes local information at infinity of a
# linear recurrence operator.
# The procedure cnd returns the local types at infinity of all
# possible hypergeometric solution, and also a set of which the
# number of elements bounds the dimension the hypergeometric
# solutions with that local type. It also outputs a number (the
# d in the (c,n,d)) that can be used for bounding the valuation of
# polynomial or rational solutions at infinity. The c in (c,n,d) is
# called Z in the Maple's previous hypergeomsols.
#
# Calling sequence:
#   g_p(L,x,tau,infinity)
# Then g_p calls the procedure cnd.


# Input: L=sum(a.i*tau^i,i=0..n), a.0<>0
#  a.i=collect(a.i,x,Normalizer)
# Output:  a list giving 3 things: 
#  1) the coordinates of the extreme points of the Newton polygon
#  2) the slope of a point to the next point
#  3) the Newton polynomial of this slope
#
# Similar to Newtonpolygon in algcurves/puiseux and
# in DEtools/diffop/factor_op
# There are two differences with algcurves:
# 1) the line val_powd:=min(op(vals)) in algcurves, which causes
#    it to search for slopes >= 0 only. This procedure should
#    also compute negative slopes.
# 2) the ldegree instead of degree.
Newtonpolygon:=proc(f,x,tau)
	local n,vals,j,res,m,i,powd,npg,slope;
	n:=degree(f,tau);
	vals:=[seq(-degree(coeff(f,tau,i),x),i=0..n)];
	powd:=0;
	npg:=[[powd,vals[powd+1]]];
	while powd<n do
		m:=min(seq((vals[powd+j+1]-vals[powd+1])/j
		 ,j=1..n-powd));
		i:=n;
		while vals[i+1]-vals[powd+1]<>(i-powd)*m do i:=i-1 od;
		powd:=i;
		npg:=[op(npg),[powd,vals[i+1]]]
	od;
	res:=NULL;
	for i from 1 to nops(npg)-1 do
		# Now we compute the Newton polynomial of slope number i
		slope:=(npg[i+1][2]-npg[i][2])/(npg[i+1][1]-npg[i][1]);
		res:=res,[op(npg[i]),slope,add(
		 coeff(coeff(f,tau,denom(slope)*j+npg[i][1]),x,
		  -j*numer(slope)-npg[i][2])
		 *x^j,j=0..(npg[i+1][1]-npg[i][1])/denom(slope))]
	od;
	[res,npg[-1]]
end:

# Compute monic factors over field F, with multiplicity.
# args[-1]='roots' means take RootOf of it.
g_factors:=proc(f,x,F)
	local a,i;
	a:=`DEtools/kovacicsols/factors`(f,x,F);
	[seq(`if`(args[-1]='roots',RootOf(i[1],x),
	 [collect(i[1]/lcoeff(i[1],x),x,Normalizer),i[2]]
	 ),i=a)]
end:

# L is an operator.
# Output: the symmetric product of L with tau-1/r
symprod1:=proc(L,r,x,tau)
	local R,S,i;
	R:=r;
	S:=coeff(L,tau,0)+R*coeff(L,tau,1)*tau;
	for i from 2 to degree(L,tau) do
		R:=R*subs(x=x+i-1,r);
		S:=S+R*coeff(L,tau,i)*tau^i
	od;
	collect(primpart(S,tau),tau,normal_poly)
end:

# Input: L=sum(a.i*tau^i,i=0..n), a.0<>0,
#  a.i=collect(a.i,x,Normalizer) in C[x]
# To be called by g_p(L,x,tau,infinity)
#
cnd:=proc(L,x,tau,ext)
	local v,i,j,Ln,Lnc,F,res,np,n;
	# Remove non-integer slopes from Newtonpolygon.
	v:=[seq(`if`(nops(i)>3 and type(i[3],integer)
	  ,seq([j,i[3]],j=g_factors(i[4],x,ext,'roots'))
	,NULL),i=Newtonpolygon(L,x,tau))];
	if v=[] then
		return NULL
	fi;
	# Now v is a list  [[c1,n1],[c2,n2],...]
	# Note: the c.i were called Z in the previous hypergeomsols
	# algorithm.
	res:=NULL;
	if ext={} and has(v,'RootOf') then
		Normalizer:=`evala/Normal`
	fi;
	n:=degree(L,tau);
	for i in v do
		if not assigned(Ln[i[2]]) then
			Ln[i[2]]:=symprod1(L,x^i[2],x,tau)
		fi;
		Lnc:=collect(subs(tau=i[1]*(Delta+1),Ln[i[2]]),Delta);
		F:=max(seq(degree(coeff(Lnc,Delta,j),x)-j,j=0..n));
		do	# do find the Newton polynomial with slope 1
			# in the Delta polygon. If it is zero (which can happen
			# because we didn't normalize, so degree could have given
			# a result F which was too high) then we'll decrease F.
			np:=add(Normalizer(coeff(coeff(Lnc,Delta,j),x,F+j))*x^j
			 ,j=0..n);
			if np<>0 then break
			elif F<-n then error
			fi;
			F:=F-1
		od;
		np:=add(coeff(np,x,F)*mul(-(x+j),j=0..F-1),F=0..degree(np,x));
		# The factors are monic.
		np:=[seq(j[1],j=g_factors(np,x,indets([args,i]
		 ,{'RootOf','radical','nonreal'})))];
		np:=int_difference(np,x);
		res:=res,seq([i[1],-i[2],RootOf(j[1],x),j[2]],j=np)
	od;
	userinfo(3,'LREtools',`c,n,d,{integers}`,res);
	res
end:

# v is a list of irreducible polynomials. Compute which of them
# are integer shifts of each other, determine the "minimal" ones, and
# the non-negative integers you have to shift them with to get the
# other ones.
int_difference:=proc(v,x)
	local i,j,d,e,large,dif;
	large:={};
	for i in v do dif[i]:={0} od;
	for i from 1 to nops(v)-1 do for j from i+1 to nops(v) do if
	 degree(v[i],x)=degree(v[j],x) then
		d:=degree(v[i],x)-1;
		e:=Normalizer(coeff(v[i],x,d)-coeff(v[j],x,d));
		if type(e,integer) and e<>0 and irem(e,d+1)=0 and
		 Normalizer(v[i]-subs(x=x+e/(d+1),v[j]))=0 then if
			e>0 then
				large:=large union {v[j]};
				dif[v[i]]:=dif[v[i]] union {e}
			else
				large:=large union {v[i]};
				dif[v[j]]:=dif[v[j]] union {-e}
		fi fi
	fi od od;
	d:={op(v)} minus large;
	[seq([i,dif[i]],i=d)]
end:

# L=operator
# LP = collect(primpart(L,tau,'co'),tau,normal_poly)
#
singularities:=proc(L,x::name,tau::name,ext::set,LP,co)
local F,i;
  if LP=0 then
	error "wrong input"
  fi;
  F:=evala(Normal(lcoeff(L,tau)*tcoeff(L,tau)/co^2));
  if F=0 then
	return procname(LP,x,tau,ext,LP,1)
  fi;
  F:={seq(i[1],i=`DEtools/kovacicsols/factors`(F,x,ext))};
  F:={seq(`polytools/shorten_int`(i[1],x),i=int_difference(F,x))};
  `polytools/sort_poly`([op(F)],x)
end:

#savelib('g_p','do_it','extend_right','normal_poly','singularities',\
 'Newtonpolygon','g_factors','symprod1','cnd','int_difference'):