#######################################################################
## Research: Save this file as PositiveCoefficients.txt. To use it, stay in the#
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read "PositiveCoefficients.txt"; <Enter>           #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Bryan Ek, Rutgers University,                           #
##  bryan.t.ek@math.rutgers.edu.                                      # 
#######################################################################

print(`This is PositiveCoefficients`):
print(`by Bryan Ek, advisee to Dr. Zeilberger, Rutgers University.`):
print(``):
print(`If you have comments or notice any errors, please email`):
print(`bryan [dot] t [dot] ek [at] math [dot] rutgers [dot] edu`):
print(`bryan.t.ek@math.rutgers.edu`):
print(``):
print(`The latest version (10/4/2017) of this package is available from`):
print(`http://www.math.rutgers.edu/~bte14/Code/PositiveCoefficients/PositiveCoefficients.txt`):
print(``):
print(`For a list of functions this package provides type Help().`):
print(``):

with(combinat):

######################################################################################################
#Help function
Help:=proc():
	if args=`guessConstantSymmetricRecurrence ` then
		print(`guessConstantSymmetricRecurrence(A,n,m): A is a table of values. n is the dimension of A. m is the maximum total depth to search. Individual max depth is currently set to 1.`):
		print(`Guessing a constant coefficient recurrence.`):
		print(`Given an nD organized "table" of values (A(0,1,0,0,1,1) etc).`):
		print(`A(x1,x2,x3,...)=alpha[1,0,0,...]*A(x1-1,x2,x3,...)+alpha[0,1,0,...]*A(x1,x2-1,x3,...)+...`):
		print(``):
		print(``):
		print(``):
		print(``):
		print(``):
		print(``):
		print(``):
	else
		print(`Useful functions: guessConstantSymmetricRecurrence, subsetvectorsm, subsetvectorsM, gDRN`):
	fi:
end:
######################################################################################################


######################################################################################################
#Guessing a constant coefficient recurrence.
#Given an nD organized "table" of values (A(0,1,0,0,1,1) etc).
#A(x1,x2,x3,...)=alpha[1,0,0,...]*A(x1-1,x2,x3,...)+alpha[0,1,0,...]*A(x1,x2-1,x3,...)+...
#Given max total depth, m, (how much all coordinates go back) and max individual depth (how much each coordinate can go back). Starting with max ind. depth = 1.
#for [1,2,...,xn] choose 1,2,3,...,m.
#create add(n choose k,k=1..m) equations of the form: A(x1,...,xn)=add(alpha[0,...,1,...,0]*A(x1,...,xi-1,...,xn),over all desired subsets)
#Start with A(1,1,...,1)=...
#Build up from there
#Again, go over all desired subsets and add 1 to each coordinate to create a new equation.
#Solve for the alpha[0,...,1,...,0].
#Output (if it exists) the solution.
#guessConstantSymmetricRecurrence(A,n,m): A is a table of values. n is the dimension of A. m is the maximum total depth to search. Individual max depth is currently set to 1.
#We are assuming that A is symmetric in its parameters.
guessConstantSymmetricRecurrence:=proc(A,n,m) local kappa,B,x,vec,vectors,vars,eqns,sols:
	vectors:=subsetvectorsM(n,m):
	vars:={seq(kappa[op(sort(vect,`>`))],vect in vectors)}:
	eqns:={seq(A[op([1$n]+vect2)]=add(kappa[op(sort(vect,`>`))]*A[op([1$n]+vect2-vect)],vect in vectors),vect2 in vectors union {[0$n]}),seq(A[op([2$n]+vect2)]=add(kappa[op(sort(vect,`>`))]*A[op([2$n]+vect2-vect)],vect in vectors),vect2 in vectors union {[0$n]})}:#,seq(A[op([3$n]+vect2)]=add(kappa[op(sort(vect,`>`))]*A[op([3$n]+vect2-vect)],vect in vectors),vect2 in vectors union {[0$n]})}:
#Not sure if the second set of equations is necessary. I guess a nice check. Helps to eliminate free variables.
#Third set of equations should eliminate even more false positives. Really only useful in 2D case.
#	vec:=[seq(x[i],i=1..n)]:
	solve(eqns,vars):
#	sols:=solve(eqns,vars):
#	if sols<>NULL then
#		[subs(sols,B(op(vec))=add(kappa[op(sort(vect,`>`))]*B(op(vec-vect)),vect in vectors)),B,kappa,x]:
#	else
#		"No solution found":
#	fi:
end:


######################################################################################################
#subsetvectorsm(n,m): returns a set of all binary vectors of length n with # of ones=m.
#Must have n>=m.
subsetvectorsm:=proc(n,m) option remember:
	if m=n then
		{[1$n]}:
	elif m=0 then
		{[0$n]}:
	else
		{seq([1,op(psub)],psub in subsetvectorsm(n-1,m-1)),seq([0,op(psub)],psub in subsetvectorsm(n-1,m))}:
	fi:
end:
###################################################
#subsetvectorsM(n,M): returns a set of all binary vectors of length n with 1<=# of ones<=M.
subsetvectorsM:=proc(n,M):
	{seq(op(subsetvectorsm(n,m)),m=1..M)}:
end:


######################################################################################################
#findA(m,n,k): finds the coefficient of x^m*y^n*z^k in the taylor expansion of 1/(1-x-y-z+4*x*y*z)
#Uses a special recurrence found by Gillis and Kleeman (1979).
findA:=proc(m,n,k) local newM,newK,newN: option remember:
	newM:=max(m,n,k):
	newK:=min(m,n,k):
	newN:=m+n+k-newM-newK:
	if newK<0 then
		return(0):
	elif newM=0 then
		return(1):
	else
		((newM+newN-newK)*findA(newM-1,newN,newK)+2*(newM-newN+newK-1)*findA(newM-1,newN,newK-1))/newM:
	fi:
end:
###################################################
#gDR3(): genericDerivativeRecurrence. 3D case.
#gDR3:=proc(M,N,K,alpha,beta,gamma,delta) local m,n,k:
#	m:=max(M,N,K):
#	k:=min(M,N,K):
#	n:=M+N+K-m-k:
#	if k<0 then
#		return(0):
#	elif m=0 then
#		return(1):
#	else
#		#How to set anything that does not have a value to 0...
#		((alpha[1,0,0]*(m-1-n)+beta[1,0,0]*(n-k)+gamma[1,0,0]*k+delta[1,0,0])*A(m-1,n,k)+(alpha[0,1,0]*(m-n+1)+beta[0,1,0]*(n-1-k)+gamma[0,1,0]*k+delta[0,1,0])*A(m,n-1,k)+(alpha[0,0,1]*(m-n)+beta[0,0,1]*(n-k+1)+gamma[0,0,1]*(k-1)+delta[0,0,1])*A(m,n,k-1)+(alpha[1,1,0]*(m-n)+beta[1,1,0]*(n-1-k)+gamma[1,1,0]*k+delta[1,1,0])*A(m-1,n-1,k)+(alpha[1,0,1]*(m-1-n)+beta[1,0,1]*(n-k+1)+gamma[1,0,1]*(k-1)+delta[1,0,1])*A(m-1,n,k-1)+(alpha[0,1,1]*(m-n+1)+beta[0,1,1]*(n-k)+gamma[0,1,1]*(k-1)+delta[0,1,1])*A(m,n-1,k-1)+(alpha[1,1,1]*(m-n)+beta[1,1,1]*(n-k)+gamma[1,1,1]*(k-1)+delta[1,1,1])*A(m-1,n-1,k-1))/(alpha[0,0,0]*(m-n)+beta[0,0,0]*(n-k)+gamma[0,0,0]*k+delta[0,0,0]):
#	fi:
#end:
###################################################
#gDRN(): values is a list of [x1,x2,x3,...,xn]. red (reductions) is a list of tables [alpha,beta,gamma,...,delta] (length n+1). M is the deepest depth of reductions.
#Requires pre-parsing of red so that every possible reduction vector has a value. 0 if necessary.
gDRN:=proc(values,M,red) local sV,n,newValue,vect:
	sV:=sort(values,`>`):
	if sV[-1]<0 then
		return(0):
	elif sV[1]=0 then
		return(1):
	fi:
	n:=nops(values):
	newValue:=0:
	for vect in subsetvectorsM(n,M) do
		newValue:=newValue+(add(red[i][op(vect)]*(sV[i]-sV[i+1]),i=1..n-1)+red[n][op(vect)]*sV[n]+red[n+1][op(vect)])*gDRN(sV-vect,M,red):
	od:
	newValue/(add(red[i][0$n]*(sV[i]-sV[i+1]),i=1..n-1)+red[n][0$n]*sV[n]+red[n+1][0$n]):
end:
###################################################
#gDRN2(): values is a list of [x1,x2,x3,...,xn]. red (reductions) is a list of matrices [alpha,beta,gamma,...,delta] (length n+1). M is the deepest depth of reductions.
#Requires pre-parsing of red so that every possible reduction vector has a value. 0 if necessary. 0 is default for matrices.
gDRN2:=proc(values,M,red) local sV,n,newValue,vect: option remember:
	sV:=sort(values,`>`):
	if sV[-1]<0 then
		return(0):
	elif sV[1]=0 then
		return(1):
	fi:
	n:=nops(values):
	newValue:=0:
	for vect in subsetvectorsM(n,M) do
		newValue:=newValue+(add(red[i][op(vect)]*(sV[i]-sV[i+1]),i=1..n-1)+red[n][op(vect)]*sV[n]+red[n+1][op(vect)])*gDRN(sV-vect,M,red):
	od:
	newValue/(add(red[i][0$n]*(sV[i]-sV[i+1]),i=1..n-1)+red[n][0$n]*sV[n]+red[n+1][0$n]):
end:


######################################################################################################
#randomTables3D(upperBound): produces random tables of length 3.
randomTables3D:=proc(upperBound,alpha,beta,gam,delta) local ra,ra01,i,vec:
	ra:=rand(0..upperBound):
	ra01:=rand(0..1):
	for i from 0 to 2^3-1 do
		vec:=convertToBinaryList(i,3):
		alpha[op(vec)]:=ra01()*ra():
		beta[op(vec)]:=ra01()*ra():
		gam[op(vec)]:=ra01()*ra():
		if vec[-1]=1 then
			delta[op(vec)]:=ra01()*rand(-gam[op(vec)]..upperBound)():
		else
			delta[op(vec)]:=ra01()*ra():
		fi:
	od:
	[alpha,beta,gam,delta]:
end:
###################################################
#randomMatrices3D(upperBound): produces random 3D matrices of length 2X2X2.
randomMatrices3D:=proc(upperBound) local ra,ra01,i,vec,a,b,g,d:
	ra:=rand(0..upperBound):
	ra01:=rand(0..1):
	a:=Array(0..1,0..1,0..1):
	b:=Array(0..1,0..1,0..1):
	g:=Array(0..1,0..1,0..1):
	d:=Array(0..1,0..1,0..1):
	for i from 0 to 2^3-1 do
		vec:=convertToBinaryList(i,3):
		a[op(vec)]:=ra01()*ra():
		b[op(vec)]:=ra01()*ra():
		g[op(vec)]:=ra01()*ra():
		if vec[-1]=1 then
			d[op(vec)]:=ra01()*rand(-g[op(vec)]..upperBound)():
		else
			d[op(vec)]:=ra01()*ra():
		fi:
	od:
	[a,b,g,d]:
end:


######################################################################################################
#convertToBinaryList(num,n): num=number to convert, n=number of binary digits. Should have 2^n>num.
convertToBinaryList:=proc(num,n):
	if num=0 then
		return([0$n]):
	fi:
	[0$(n-1-floor(log[2](num))),op(ListTools[Reverse](convert(num,base,2)))]:
end: