TauLocalDatas:=proc(M,x,lambda)
	local B, L, p,lj;
	B := lcm(op(convert(map(denom,eval(M)),set)));
	p := degree(B,x) + 1;
	L := g_super_reduce(evalm(x^p*M-1), x, 1/x, lambda);
	[seq(`if`(degree(lj[2],lambda) > 0,  [lj[1]-1-p, subs(lambda=-lambda,lj[2])], NULL), lj = L[3])]
end:


# scIndicialEquation(M,x, s,c,lambda, T, invT)
# In : An invertible difference system Y(x+1)=M(x)Y(x)   (1)
# M is an invertible nxn matrix  with rational functions (in x) entries, 
# x is a name
# lambda is a name
# s is an integer (a slope of the tau-Newton polygon of (1))
# c is a constant  (c a root of the Newton polynomial associated with n). This means that (1) has a formal solution of the form : c^x*Gamma(x)^n*Z(x)
# Z(x) satisfies then the system : Z(x+1)= (x^(-n)/c)*M(x)*Z(x)    (2)
# Here we are interested in those Z(x) of the form Z(x) = x^d*(a nonzero  formal power series in 1/x) where d is a constant
# If such a Z exists then  the delta-polygon of (2) has a horizontal side (i.e side with slope 0)  and d is a root of the corresponding indicial equation.
# Out : the  indicial equation of (2), that is a polynomial in lambda
#
scIndicialEquation:=proc(M,x,s,c,lambda)
	local L, T, invT,nu,eq;
	L := g_super_reduce(evalm((x^(-s)/c)*M-1), x, 1/x, lambda);
	L[1], subs(lambda=-lambda,L[3][1][2])
end:

LocalDatas:=proc(M,x,lambda)
	local l, j, lj, s, lc, c, ll, ld, eqsc, degT;
	ll := NULL;
	if _Env_ratsols_only = true then
		# Uncomment the next line "return .." if you want to return degree-information
		return [ [0, [1,1], [[0,0]], degT] ];
		l := [ [0, lambda-1] ]
	else
		l := TauLocalDatas(M, x, lambda)
	fi;
	for j from 1 to nops(l) do 
		lj := l[j];
		s := lj[1];
		lc := roots(lj[2],lambda);
		for c in lc do
			degT, eqsc := scIndicialEquation(M,x,s,c[1],lambda);
			ld := roots(eqsc,lambda);
			if _Env_ratsols_only = true then
				ld := [seq(`if`(type(i[1],integer),i,NULL), i=ld)]
			fi;
			if nops(ld) > 0 then
				ll := ll, [s,c, ld, degT]
			fi:
		od:
	od:
	[ll]
end:


PolySols := proc(M,x) local n,l,i,Z,y,Y,l1;
	n := linalg['rowdim'](M);
	l := map(factor,map(primpart, MatGcd(M,n), x));
	l1 := map(factor,subs(x=x-1, l));
	if has(l,x) and nargs < 10 then
		# Without the "nargs < 10" this may not terminate if there are no polynomial solutions.
 		return [seq([seq(l1[i]*Y[i], i=1..n)], Y =
		procname(map(normal,evalm(
			linalg['diag'](seq(1/i,i=l)) &* M &*
			linalg['diag'](op(l1))
		)),args[2..-1],x))]
	fi;
	Y := [seq( (y||i)(x), i=1..n) ];
	# lprint("time before pol-sol", time());
	l := LinearFunctionalSystems:-PolynomialSolution(convert(
		evalm(evalm( M &* Y ) - subs(x=x+1,Y)),list), Y);
	# lprint("time after pol-sol", time());
	Y := indets(convert(l,list)) minus indets(convert(eval(M),set));
	l := subs(RowReduceR(l, n, Y), l);
	l := [seq(map(factor,map(coeff, l, i)), i=Y)];
	# lprint("Found degrees", seq(max(map(degree,i,x)),i=l));
	# lprint("Found", seq(map(degree,i,x),i=l)); l;
end:


MatGcd := proc(M,n) local Z,i,J, P, j, v1;
	global x, cct;
	for i to n do v1[i] := content(add(numer(M[i,j])*Z^j, j=1..n),Z) od:
	for i to n do
		J := 0;
		# Could do better here and search for some irreducible P
		# that appears only once in the denominator with max-pole-order.
		for j to n do if J<>-1 and has(denom(M[i,j]), x) then
			if J=0 then
				J := j; P := denom(M[i,j])
			else
				J := -1
			fi
		fi od;
		if J > 0 then
			v1[J] := lcm( v1[J], subs(x=x+1,P) )
		fi
	od:
	[seq( factor(primpart(v1[i],x)), i=1..n) ]
end:

infolevel[HyperSols] := 2:  # Used for printing info, higher number prints more data.

# Check if two monic irreducible poly's are an integer shift from each other.
IntShift := proc(A, P, x, d)
	local n,e;
	n := degree(A,x);
	if n=degree(P,x) then
		e := Normalizer(coeff(A,x,n-1)-coeff(P,x,n-1))/n;
		if type(e,integer) and Normalizer(A-subs(x=x+e,P))=0 then
			if nargs=4 then
				d := e
			fi;
			return true
		fi
	fi;
	false
end:

# Example:
#	A = [[x+1/2,1], [x,2], [x-1,1]];
#	P = x
# Then Output = [[x + 1/2, 1]],
#		[[x, 0], [x, 0], [x - 1, 1]]
#
FindIntShifts := proc(A, P, x)
	local res, resA, i, d;
	resA := NULL;
	res := NULL;
	for i in A do
		if IntShift(P, i[1], x, 'd') then
			res := res, [i[1],d]$i[2]
		else
			resA := resA, i
		fi
	od;
	[resA], sort([res],proc(i,j) evalb(i[2]<j[2]) end)
end:

AddSet := proc(L)
	local i,j;
	if nargs=0 then
		{0}
	elif nargs=1 then
		L
	else
		{seq(seq(i+j,j=L),i=procname(args[2..-1]))}
	fi
end:

CombinationsAB := proc(L)
	local i,j;
	if nargs=0 then
		[ [1,1,0,0] ]
	elif nargs=1 then
		L
	else
		[seq(seq([i[1]*j[1], i[2]*j[2], i[3]+j[3], i[4]+j[4]],
		 j=L),i=procname(args[2..-1]))]
	fi
end:

Count := proc(i, L)
	local j;
	add(`if`(j[2]=i,1,0),j=L)
end:

Remove := proc(i, v)
	local j;
	for j from 1 to nops(v) do if v[j,2]=i then
		return [op(v[1..j-1]),op(v[j+1..-1])], v[j]
	fi od;
	v, FAIL
end:
		

# DenomMatrix := proc(M)
#	lcm(op({op(map(denom,map(op,convert(M,listlist))))}))
# end:
#
# For some reason this is slow, so I replaced it by:
#
#
DenomMatrix := proc(M)
	local s,d,i;
	s := sort([op({op(map(denom,map(op,convert(M,listlist))))})],length);
	d := s[-1];
	for i from 2 to nops(s) do
		d := d*denom(Normalizer(d/s[-i]))
	od;
	d
end:

# Look for pairs A,B such that we have to find solutions of the form
#  	Sol(tau - c * A/B) * Polynomial
#
#
PairsAB := proc(Minp, sclist, x, ext)
	local M, tM, MM, Mi, Mti, MMi, dM, dMi, dMM, dMMi, 
	 Af, Bf, dMMf, dMMif, i, DBoundA, DBoundB, Pseq, P,
	 vA, vB, w, q1, L, i_min, i_max, removed, q2, dP, db,
	 resP, sSet, SClist, dMf, dMif;
	if nops(sclist)=0 then
		return []
	elif not type(args[-1],set) then
		return procname(args, indets([args],{'RootOf','radical','nonreal'}))
	elif ext = {} then
		# These should be moved to the top procedure
		Normalizer := normal
	else
		Normalizer := evala
	fi;
	# lprint(`enter PairsAB`,time());
	M  := map(Normalizer, Minp);
	tM := map(Normalizer, eval(M, x=x+1));
	MM := map(Normalizer, evalm(tM &* M));
	Mi := map(Normalizer,linalg['inverse'](M));
	Mti:= map(Normalizer, eval(Mi, x=x-1));
	MMi :=map(Normalizer, evalm(Mti &* Mi));
	dM, dMi, dMM, dMMi:= seq(DenomMatrix(i), i = [M, Mi, MM, MMi]);
	dMf, dMif := seq(factor(i), i=[dM, dMi]);
	dMM := gcd(dMM,  dMf*eval(dMf, x=x+1)); # This simply factors dMM
	dMMi:= gcd(dMMi,dMif*eval(dMif,x=x-1)); # and dMMi

	# This makes the factors monic and Normalizes the coefficients:
	Af, Bf, dMMf, dMMif := seq(`LREtools/g_p/factors`(i,x,ext),
		i=[dMif, dMf, dMM, dMMi]);

	# Example:
	#	Af = [[x+1/2,1], [x,2], [x-1,1]]
	#	Bf = [[x-5/2,2], [x,1], [x-2,1], [x+2,1]]

	# Denominator bounds and their shifts
	DBoundA := 1;
	DBoundB := 1;

	Pseq := NULL;
	do
		# Select the next factor P:
		i := [op(Af), op(Bf)];
		if i = [] then break fi;
		P := i[1,1];
		Af, vA := FindIntShifts(Af, P, x);
		Bf, vB := FindIntShifts(Bf, P, x);
		#
		# Example: If P was x, then we now have
		# 	Af = [[x + 1/2, 1]]
		#	Bf = [[x - 5/2, 2]]
		#	vA = [[x, 0], [x, 0], [x - 1, 1]]
		#	vB = [[x + 2, -2], [x, 0], [x - 2, 2]]
		#
		userinfo(3,'HyperSols',`number of combinations at`,P,
			`is at most`,nops(vA)+nops(vB)+1);

		if vB<>[] then
			dMMf, q1 := FindIntShifts(dMMf, P, x);
			L[vB[1,2]-1] := 0;
			L[vB[1,2]-2] := 0;
			i_min,i_max := vB[1,2], vB[-1,2];
			for i from i_min to i_max + 1 do
				L[i] := L[i-1] + Count(i, vB);
				w := L[i-2] + Count(i-1, q1);
				while w < L[i] do
					vB, removed := Remove(i, vB);
					if removed=FAIL then
						vB, removed := Remove(i-1, vB);
						if removed=FAIL then
							error "counted wrong"
						fi;
						# This will later be lcm'd with DBoundA
						DBoundB := DBoundB * eval(removed[1],x=x-1)
					fi;
					L[i] := L[i]-1
				od
			od
		fi;

		if vA<>[] then
			dMMif, q2 := FindIntShifts(dMMif, P, x);
			L[vA[-1,2]+1] := 0;
			L[vA[-1,2]+2] := 0;
			i_min,i_max := vA[1,2], vA[-1,2];
			for i from i_max to i_min - 1 by -1 do
				L[i] := L[i+1] + Count(i, vA);
				w := L[i+2] + Count(i+1, q2);
				while w < L[i] do
					vA, removed := Remove(i, vA);
					if removed=FAIL then
						vA, removed := Remove(i+1, vA);
						if removed=FAIL then
							error "counted wrong"
						fi;
						DBoundA := DBoundA * removed[1]
					fi;
					L[i] := L[i]-1
				od
			od
		fi;

		userinfo(3,'HyperSols',`number of combinations at`,P,`is now`,nops(vA)+nops(vB)+1);

		# For A, use the last ones from vA in resP.
		# For B, use the first ones from vB in resP (these are moved from vB to w)
		#
		w := []; resP := [vA,w];

		# Now loop over the nops(vA)+nops(vB)+1 combinations:
		do
			if [vA,vB]=[[],[]] then
				break
			elif vA = [] or (vB<>[] and vA[1,2] > vB[1,2]) then
				w := [op(w), vB[1]];
				vB := vB[2..-1]
			else
				vA := vA[2..-1]
			fi;
			if _Env_ratsols_only <> true or nops(vA)=nops(w) then
				resP := resP, [vA,w]
			fi
		od;

		dP := degree(P,x);
		# Give list of all [A,B,s,d] and add this list to Pseq
		# Here s corresponds to a slope of the NewtonPolygon
		# So s = -v where v is as in AAECC, Vol. 17, 83-115 (2006).
		# The d is the same as in that paper (e.g. a polynomial solution
		# of degree 5 would correspond to some [c,s,d] with c=1, s=0, d=5.
		# Here [c,s,d] corresponds to a factor tau - c*x^s*(1+d*t+O(t^2))
		# where t = 1/x.  In other words:   tau - c*t^(-s)*(1+d*t+O(t^2))
		Pseq := Pseq, [seq([mul(i[1],i=w[1]), mul(i[1],i=w[2]),
			(nops(w[1])-nops(w[2])) * dP, Normalizer(
				 add(coeff(i[1], x, dP-1), i=w[1])
				-add(coeff(i[1], x, dP-1), i=w[2]))], w=[resP])]
	od;
	db := lcm(DBoundA, DBoundB);
	if has(db,x) then
		Pseq := Pseq, factor( [[db, eval(db,x=x+1), 0, -degree(db,x)]] )
	fi;

	# Select those elements from sclist for which s is consistent with some A,B.
	sSet := AddSet(seq({seq(i[3],i=w)},w=[Pseq]));
	SClist := [seq(`if`(member(i[1],sSet),i,NULL),i=sclist)];
	i := nops(sclist)-nops(SClist);
	if i>0 then userinfo(2, 'HyperSols', i,
		`pairs s,c removed because s is incompatible with degrees A,B`)
	fi;
	if add(nops(i[3]),i=SClist)=0 then
		return []
	fi;

	# SClist and sclist are of the form [ [s, c, dlist, degT], ...
	# Slope s is in Z,  c = [c, mult] and dlist = [[d,mult], ... ].

	# lprint(`Now forming all combinations and testing which ones match c,s,d`,time());
	# This is [ [A, B, s, c, [integer d's], degT], ... ]
	P := [seq(seq(MatchAB(op(i), w),w=SClist), i=CombinationsAB(Pseq))];
	userinfo(2,'HyperSols', `number of combinations left after comparing with local data c,s,d`, nops(P));
	Check_p_curvature(P,x,2,ext,MM,MMi)
end:

# Check if an A/B can be matched.
# Still need to add: p-curvature tests.
MatchAB := proc(A,B,s,d,L)
	local dd,i;
	if s=L[1] then
		dd := select(type,[seq(Normalizer(i[1]-d), i=L[3])], integer);
		# Can replace "integer" by "nonnegint" if R = dlist in deltaPS
		if dd<>[] then
			return [A,B,s,L[2],dd,L[4]]
		fi
	fi;
	NULL
end:

Check_p_curvature := proc(PP, x, p, ext, MM, MMi)
	local P, w, pc, pc2, i;
	P := PP;
	w := traperror((Primfield(ext union {i}) mod p)[2]);
	if w = lasterror then
		userinfo(3, 'HyperSols', `Skipping prime`,p, `due to difficulties reducing coefficients mod p`);
	elif p=2 then
		pc := traperror(p_val(convert(MM,listlist),p,w));
		if pc = lasterror then
			userinfo(3, 'HyperSols', `Skipping prime`,p, `due to difficulties reducing coefficients mod p`);
			return P
		fi;
		pc2 := traperror(p_val(convert(MMi,listlist),p,w));
		if pc2 = lasterror then
			pc2 := 0, -infinity
		fi;
		P := [seq(`if`(Check_2_curvature(i,w,x,pc,pc2),i,NULL), i=P)];
		forget(DetMod2)
	else
		lprint( "decided to do p-curvature only for p=2" )
	fi;
	userinfo(3, 'HyperSols',cat( `number of combinations left after checking `,p,`-curvature is `, nops(P)),`time is:`,time());
	P
end:

# v is a list [A, B, s, c, ...]
# Check if it is consistent with 2-adic data.
#
Check_2_curvature := proc(v,w,x,M,vM,Mi,vMi)
	local A,a,d;
	A := v[4,1]*v[1]/v[2];
	A := Normalizer(A * eval(A,x=x+1));
	A := traperror(p_val(numer(A),2,w), p_val(denom(A),2,w));
	if A=lasterror then return true fi;
	# This line should help speed up the computation of the determinant:
	Normalizer := proc(a) Normal(a) mod 2 end:
	A, a := Normalizer(A[1]/A[3]), A[2]-A[4];
	if a < vM or -a < vMi then
		return false
	fi;
	d := DetMod2(M, `if`(a=vM, A, NULL));
	if d <> 0 then
		false
	elif a=vM or (-a > vMi) then
		true
	else
		d := DetMod2(Mi, 1/A);
		evalb( d = 0 )
	fi
end:

# Determine the p-valuation of an expression a, if this
# p-valuation is an integer.  Determine also the reduction
# of a/p^v mod p.
p_val := proc(a,p,w)
	local b,c,v;
	if nargs=2 then
		b := traperror(a mod p);
		if b = lasterror then
			b, v := procname(p*a, p);
			return b, v-1
		fi;
		c := {b};
		# To handle type list and listlist:
		while type(c[1],list) do c:=map(op,c) od;
		if c={0} then
			b, v := procname(a/p, p);
			b, v+1
		else
			b, 0
		fi
	else
		b, v := frontend(procname,[a,p],[{`+`,`*`,list,RootOf,radical,nonreal},{}]);
		b := Normal(eval(b,w)) mod p;
		c := {b};
		while type(c[1],list) do c:=map(op,c) od;
		if c={0} then
			# The two-valuation is not an integer, and we only have
			# code for the integer valuation case. All we can say
			# is that the valuation is greater than v, which we
			# indicate by having v not be an integer.
			# v := v + 1/2
			# Actually, lets just give up in this case:
			error "failed reduction"
		fi;
		b, v
	fi
end:


DetMod2 := proc(l,e)
	options remember;
	local i;
	if nargs=2 then
		procname([seq(subsop(i=Normalizer(l[i,i]-e),l[i]),i=1..nops(l))])
	else
		Normalizer(LinearAlgebra:-Determinant(Matrix(l)))
	fi
end:




# deltaHS(M,x)
#
#  INPUT:  M   - a nonsingular matrix with rational function entries (over Q)
#          x   - a variable
#          
#  OUTPUT: the general hypergeometric 
#          solutions of the the first order difference system 
#		y(x+1) = A(x)y(x)
#	   The arbitrary constants are named c_1, c_2, etc...
#
deltaHS:= proc(M,x)
	local ABlist, sclist, j, lpolysol,Nsc,HSseq,lambda;
	# lprint("Entrance",time());
	sclist := LocalDatas(M, x, lambda);
	userinfo(5,'HyperSols', "Computed Local Datas",time());
	ABlist:=PairsAB(M, sclist, x, {});
	userinfo(1,'HyperSols', "Computed",nops(ABlist), "ABPairs",time());
	HSseq:=NULL; 
	for j in ABlist do
		# j = [A,B,s,[c,m],dlist]
		Nsc:=map(normal,evalm(j[2]/(j[1] * j[4][1]) * M));
		# lprint("Indicial equation non-negative Roots:", j[5]);
		# lprint("degT", j[6]);
		lpolysol:= PolySols(Nsc,x);
		# lprint("Found degrees:", seq(map(degree,i,x),i=lpolysol));
		if type(degT,list) then for i in lpolysol do for ii to nops(i) do if degree(i[ii],x) > max(j[5]) + j[6][ii] then
			error "degrees do not match the bound"
		fi od od fi;
		# lprint("Computed", nops(lpolysol),"polynomial solutions",time());
		if lpolysol<>[] then 
			HSseq:= HSseq, [j[1],j[2],j[4][1], lpolysol]
		fi;
	od:
	# lprint("end time",time());
	HSseq
end: