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

print(`This is qWZidentities`):
print(`by Bryan Ek, advisee to Dr. Zeilberger`):
print(``):
print(`The latest version of this package is available from`):
print(`http://www.math.rutgers.edu/~bte14/Code/qIdentities/qWZidentities.txt`):
print(``):
print(`For a list of functions this package provides type Help().`):
print(`For a list of functions taken (with permission) from qEKHAD.txt, type ezra().`):
print(``):
print(`WARNING! q is a global variable. Do not assign it a value.`):
print(`We treat qfac(n)=(q)_n=mul(1-q^i,i=1..n)=(1-q)*(1-q^2)*...*(1-q^(n-1)).`):
print(`An alternate convention is qfac(n)=(q)_n/(1-q)^n=mul(1-q^i,i=1..n)/(1-q)^n.`):
print(`Under either convention, qbin(n,k)=qfac(n)/(qfac(n-k)*qfac(k)) so many identities remain the same.`):
print(`WARNING! The arguments of qfac should be linear. It may not work correctly if the argument is anything else.`):
print(`Normal functions use the WZ-pair convention F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k). If using =G(n,k)-G(n,k-1) input G(n,k-1).`):
print(`q-Functions use the WZ-pair convention F(n+1,k)-F(n,k)=G(n,k)-G(n,k-1). If using =G(n,k+1)-G(n,k) input G(n,k+1).`):
print(`Note that qEKHAD uses the convention =G(n,k)-G(n,k-1).`):

#interface(showassumed=0):


Help:=proc():
	if args=NULL then
		print(`Current useable functions:`):
		print(`SymFunc, CheckWZpair, CheckWZpairEmp, Shadow, Dual, DualPair, qfactor, qfac1, evalqfac, simplifyqfac, simplifyqfacsum, simpqfac, normqfac, CheckqWZpair, CheckqWZpairEmp, qShadow, qDual, qDualPair, qIdentities, qWZpairIdentitites, CheckqIdentity, seqCheckqIdentity.`):
		print(`The most important function is qIdentities. Given a q-hypergeometric function, it will attempt to find many q-identities.`):
		print(`Any identities found can be empirically checked with CheckqIdentity or seqCheckqIdentity`):
		print(`For help with a specific function type Help(function).`):
	elif args=`SymFunc` then
		print(`SymFunc(F,x,p): symmetrizes the function F in the variable x around the pivot point p.`):
		print(`Try SymFunc(binomial(n,k),k,n).`):
	elif args=`CheckWZpair` then
		print(`CheckWZpair(F,G,n,k): inputs F,G as expressions of n,k. Certifies that F,G are a WZ pair.`):
		print(`It tests whether F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k). Another convention is F(n+1,k)-F(n,k)=G(n,k)-G(n,k-1). If using that, input G(n,k-1).`):
		print(`Try CheckWZpair(binomial(n,k)/2^(n),-binomial(n,k-1)/2^(n+1),n,k).`):
	elif args=`CheckWZpairEmp` then
		print(`CheckWZpairEmp(F,G,n,n1,n2,k,k1,k2): inputs F,G as expressions of n,k. Empirically certifies that F,G are a WZ pair.`):
		print(`It tests whether F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k) for all max(0,n1)<=n<=n2 and k1<=k<=k2.`):
		print(`Try CheckWZpairEmp(binomial(n,k)^2/binomial(2*n,n),-(3*n-2*k+3)/(4*n+2)*binomial(n,k-1)^2/binomial(2*n,n),n,0,5,k,0,5).`):
	elif args=`Shadow` then
		print(`Shadow(F,n,k): takes a function F(n,k) and outputs the shadow form.`):
		print(`Each factor (a*n+b*k+c)! is replaced by (-1)^(a*n+b*k+c)/(-a*n-b*k-c-1)!`):
		print(`It is necessary to write F expanded into factorials and polynomials.`):
		print(`I.e. Shadow(binomial(n,k),n,k)=binomial(n,k) but Shadow(n!/(k!*(n-k)!),n,k)=(-1)^(n+k)*factorial(-k-1)/(factorial(n-k)*factorial(-1-n)).`):
	elif args=`Dual` then
		print(`Dual(F,n,nShift,k,kShift): inputs a hypergeometric function. Outputs a dual function obtained by shadowing and shifting.`):
		print(`Mainly used as a helper function for DualPair(F,G,n,k).`):
		print(`Try Dual(-n!/((k-1)!*(n-k+1)!*2^(n+1)),n,0,k,1),Dual(n!/(k!*(n-k)!*2^n),n,1,k,0).`):
	elif args=`DualPair` then
		print(`DualPair(F,G,n,k): inputs a WZ-pair. Outputs the dual pair F’,G’. Output can typically be cleaned up.`):
		print(`Requires that G(n,k)/F(n,k-1) is a rational function of n,k. Essentially that G is obtained from a certificate of F.`):
		print(`Uses Theorem B from Wilf and Zeilberger “Rational Functions Certify Combinatorial Identities” (1990).`):
		print(`Try DualPair(n!^4/(k!^2*(n-k)!^2*(2*n)!),-(3*n-2*k+3)*n!^4/((k-1)!^2*(n-k+1)!^2*(2*n)!*2*(2*n+1)),n,k).`):
	elif args=`qfactor` then
		print(`qfactor(F): attempts to factor F which can have powers of q.`):
		print(`Essentially shorthand for factor(expand(F)) with some clean up.`):
		print(`Try qfactor(1-q^k-q^n+q^(n+k)) versus factor(1-q^k-q^n+q^(n+k)).`):
	elif args=`qfac1` then
		print(`Shorthand for qfac(n) in product form. I.e. qfac1(n)=Product(1-q^i,i=1..n).`):
	elif args=`evalqfac` then
		print(`evalqfac(n): Inputs n. Outputs qfac(n) evaluated if n is an integer or -infinity. Otherwise outputs qfac(n).`):
		print(`Try evalqfac(2*n), evalqfac(3), evalqfac(-1).`):
	elif args=`simplifyqfac` then
		print(`simplifyqfac(F,nk,nkVal): Takes limit of F as nk->nkVal. Subject to |q|<1 so q^infinity=0.`):
		print(`Try simplifyqfac(-qfac(n)^4*(q^(k^2+3*n+2)-2*q^(k^2+k+n+1)+q^(k^2+2*n+1))/(qfac(k)^2*qfac(n-k)^2*qfac(2*n)*(q^(n+1)+1)*(q^(2*n+1)-1)),k,infinity).`):
		print(`Try F:=q^(k^2)*qfac(n)^4/(qfac(k)^2*qfac(n-k)^2*qfac(2*n)): simplifyqfac(F,n,infinity)-simplifyqfac(F,n,-infinity).`):
	elif args=`simplifyqfacsum` then
		print(`simplifyqfacsum(F,nkSumvar,lowerBound,upperBound): Sums F from lowerBound to upperBound.`):
		print(`It checks whether the sum is actually finite. qfac(<0)=infinity so large portions of the sum can potentially be ignored.`):
		print(`Try simplifyqfacsum(simplifyqfac(q^(k^2)*qbin(n,k)^2/qbin(2*n,n),n,0),k,-infinity,infinity)`):
	elif args=`simpqfac` then
		print(`simpqfac(F): Simplifies F (as a function of powers of q and qfac()) as much as possible.`):
		print(`Combines (1-q^a)*qfac(a-1)=qfac(a), qfac(b)/qfac(b-1)=1-q^b, and qfac(c)/(1-q^c)=qfac(c-1).`):
		print(`In the end there will be no qfac() with the opposite power sign that differs only by a constant.`):
		print(`Any factor of (1-q^i) that could be multiplied into (or divided from) a qfac will have that done.`):
		print(`The function runs by looking at all factors. If a factor is of the form qfac(), it compares it to other factors.`):
		print(`Check all other qfac factors. If one differs by only a constant (and opposite power), cancel them, leaving remainder.`):
		print(`Check all other non-qfac factors. If possible, combine with current qfac (either increasing or decreasing argument). Then try again.`):
		print(`Try simpqfac((1-q^n)*qfac(n-1)/qfac(n+2)*qfac(n+2*k+2)/(1-q^(n+2*k+2))).`):
	elif args=`normqfac` then
		print(`normqfac(F): Replaces qfac(n+2)=qfac(n)*(1-q^(n+1))*(1-q^(n+2)) so that you can potentially add/subtract products of qfac().`):
		print(`Need that F does not have any factors of the form qfac(a)-qfac(b).`):
		print(`Try simpqfac(normqfac(qfac(n)*qfac(n+k+1)*qfac(k-1))-normqfac(qfac(n+1)*qfac(n+k)*qfac(k-1))).`):
	elif args=`CheckqWZpair` then
		print(`CheckqWZpair(F,G,n,k): inputs F,G as hypergeometric expressions of q^n,q^k. Certifies that F,G are a q-WZ pair.`):
		print(`It tests whether F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k). Another convention is F(n+1,k)-F(n,k)=G(n,k)-G(n,k-1). If using that, input G(n,k-1).`):
		print(`Try CheckqWZpair(q^(k^2-m^2)*qfac(n+b+c-k)*qfac(n+m)*qfac(b-m)*qfac(c+m)*qfac(n-m)*qfac(b+m)*qfac(c-m)/(qfac(n-k)*qfac(b-k)*qfac(c-k)*qfac(k+m)*qfac(k-m)*qfac(n+b)*qfac(n+c)*qfac(b+c)),-q^(n-(k-1)+1)*(q^(k-1)-q^c)*(q^(k-1)-q^b)*q^((k-1)^2-m^2)*qfac(n+b+c-(k-1))*qfac(n+m)*qfac(b-m)*qfac(c+m)*qfac(n-m)*qfac(b+m)*qfac(c-m)/(qfac(n-(k-1))*qfac(b-(k-1))*qfac(c-(k-1))*qfac((k-1)+m)*qfac((k-1)-m)*qfac(n+b)*qfac(n+c)*qfac(b+c))/((1-q^(n+c+1))*(1-q^(n+b+1))),n,k).`):
	elif args=`CheckqWZpairEmp` then
		print(`CheckqWZpairEmp(F,G,n,n1,n2,k,k1,k2): inputs F,G as hypergeometric expressions of q^n,q^k. Empirically certifies that F,G are a q-WZ pair.`):
		print(`It tests whether F(n+1,k)-F(n,k)=G(n,k+1)-G(n,k) for all max(0,n1)<=n<=n2 and k1<=k<=k2.`):
		print(`Try CheckqWZpairEmp(q^(k^2)*qbin(n,k)^2/qbin(2*n,n),q^(n+1)*(q^(2*n+1)+q^n-2*q^(k-1))/((1-q^(2*n+1))*(1+q^(n+1)))*q^((k-1)^2)*qbin(n,(k-1))^2/qbin(2*n,n),n,0,5,k,0,5).`):
	elif args=`qShadow` then
		print(`qShadow(F,n,k): takes a q-function F(n,k) and outputs the shadow form.`):
		print(`Each factor qfac(a*n+b*k+c) is replaced by (-1)^(a*n+b*k+c)*q^((a*n+b*k+c)*(a*n+b*k+c+1)/2)/qfac(-a*n-b*k-c-1)`):
		print(`It is necessary to write F expanded into qfactorials and polynomials.`):
		print(`Try qShadow(q^(k^2)*qfac(n)^4/(qfac(k)^2*qfac(n-k)^2*qfac(2*n)),n,k).`):
	elif args=`qDual` then
		print(`qDual(F,n,nShift,k,kShift): inputs a q-hypergeometric function. Outputs a dual function obtained by shadowing and shifting.`):
		print(`Mainly used as a helper function for qDualPair(F,G,n,k).`):
		print(`Try qDual(q^(n+1)*(q^(2*n+1)+q^n-2*q^(k-1))/((1-q^(2*n+1))*(1+q^(n+1)))*q^((k-1)^2)*qfac(n)^4/(qfac((k-1))^2*qfac(n-(k-1))^2*qfac(2*n)),n,0,k,1),qDual(q^(k^2)*qfac(n)^4/(qfac(k)^2*qfac(n-k)^2*qfac(2*n)),n,1,k,0).`):
	elif args=`qDualPair` then
		print(`qDualPair(F,G,n,k): inputs a qWZ-pair. Outputs the dual pair F’,G’. Output can typically be cleaned up.`):
		print(`Requires that G(n,k)/F(n,k-1) is a rational function of q^n,q^k. Essentially that G is obtained from a certificate of F.`):
		print(`Uses Theorem B from Wilf and Zeilberger “Rational Functions Certify Combinatorial Identities” (1990).`):
		print(`Try qDualPair(q^(k^2)*qfac(n)^4/(qfac(k)^2*qfac(n-k)^2*qfac(2*n)),q^(n+1)*(q^(2*n+1)+q^n-2*q^(k-1))/((1-q^(2*n+1))*(1+q^(n+1)))*q^((k-1)^2)*qfac(n)^4/(qfac((k-1))^2*qfac(n-(k-1))^2*qfac(2*n)),n,k).`):
	elif args=`qIdentities` then
		print(`qIdentities(F,n,k,constants,toprint:=false,checkDual:=false): inputs a hypergeometric function in q^n,q^k (and a list of constants).`):
		print(`Uses qzeil to find pair. If a WZ-pair was found, it continues.`):
		print(`Outputs all of the identities that hold given certain conditions: identity, companion, dual, dual companion.`):
		print(`Takes the limiting case of all of the identities as well.`):
		print(`Also outputs F,G and the dual pair. These can be used to certify the identities.`):
		print(`For identities while the program is running, set toprint to true.`):
		print(`To ensure the dual pair found is a WZ pair, set checkDual to true.`):
		print(`Output is in the form [[Original Identity,Limit,Companion,Limit],[F,G]],[[Dual,Limit,Dual Companion,Limit],[dualF,dualG]].`):
		print(`Note that in many cases, the identities can likely be cleaned up w.r.t. summation bounds and overall appearance.`):
		print(`An identity of 0=0 indicates a certain condition was not met, or an error occurred, or you found a trivial identity.`):
		print(`Try q-Chu-Vandermonde qIdentities(q^(k^2)*qbin(n,k)^2/qbin(2*n,n),n,k,[]).`):
		print(`q-Saalschutz qIdentities(q^(k^2-m^2)*qfac(n+b+c-k)*qfac(n+m)*qfac(b-m)*qfac(c+m)*qfac(n-m)*qfac(b+m)*qfac(c-m)/(qfac(n-k)*qfac(b-k)*qfac(c-k)*qfac(k+m)*qfac(k-m)*qfac(n+b)*qfac(n+c)*qfac(b+c)),n,k,[b,c,m]).`):
		print(`q-Dixon qIdentities((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(2*n,n+k)^3/qfac(3*n)*qfac(n)^3,n,k,[]).`):
	elif args=`qWZpairIdentities` then
		print(`inputs a qWZ-pair and a list of constants.`):
		print(`Outputs original and companion identities as well as the limiting cases.`):
		print(`Try q-Saalschutz qWZpairIdentities(q^(k^2-m^2)*qfac(n+b+c-k)*qfac(n+m)*qfac(b-m)*qfac(c+m)*qfac(n-m)*qfac(b+m)*qfac(c-m)/(qfac(n-k)*qfac(b-k)*qfac(c-k)*qfac(k+m)*qfac(k-m)*qfac(n+b)*qfac(n+c)*qfac(b+c)),q^(k^2-m^2)*qfac(n+b+c-k)*qfac(n+m)*qfac(b-m)*qfac(c+m)*qfac(n-m)*qfac(b+m)*qfac(c-m)/(qfac(n-k)*qfac(b-k)*qfac(c-k)*qfac(k+m)*qfac(k-m)*qfac(n+b)*qfac(n+c)*qfac(b+c))*-(-q^c+q^k)*(-q^b+q^k)*q^n*q/((q^n*q^c*q-1)*(q^n*q^b*q-1)*q^k),n,k,[b,c,m],true).`):
	elif args=`CheckqIdentity` then
		print(`CheckqIdentity(ID,nkvar,nkval): Empirically checks the q-identity ID in the substitution nkvar=nkval.`):
		print(`It is highly recommended to trim any bounds of summation to ignore zero terms, i.e. anything with 1/qfac(<0).`):
		print(`ID should be an equation. If you have an expression, input ID=0.`):
		print(`Outputs the left side minus the right side and truncated to q^(n+1)`):
		print(`Try CheckqIdentity(Sum(2*qfac(2*k)*q^k/(qfac(k)^4*(q^k+1)), k = 0 .. infinity) = Sum(q^k/(qfac(infinity)*qfac(k)^2), k = 0 .. infinity),n,5).`):
	elif args=`seqCheckqIdentity` then
		print(`seqCheckqIdentity(ID,nkvar,nkvalmin,nkvalmax): Empirically checks the q-identity ID in the substitution nkvar=nkval for nkval=nkvalmin…nkvalmax.`):
		print(`Outputs true/false.`):
		print(`Try seqCheckqIdentity(Sum(q^(k^2)*qbin(n,k)^2/qbin(2*n,n),k=0..n)=1,n,0,20).`):
	fi:
end:


####################################################################################################
#SymFunc(f,x,p): symmetrizes the function f in the variable x around the pivot point p.
SymFunc:=proc(f,x,p):
	(f+subs(x=p-x,f))/2:
end:


####################################################################################################
#CheckWZpair(F,G,n,k): inputs F,G as hypergeometric expressions of n,k. Certifies that F,G are a WZ pair.
CheckWZpair:=proc(F,G,n,k):
	evalb(simplify(expand((subs(n=n+1,F)-subs(k=k+1,G)+G))/F-1)=0):
end:
#CheckWZpairEmp(F,G,n,n1,n2,k,k1,k2): inputs F,G as expressions of n,k.
#Empirically checks that F,G are a WZ pair for all n1<=n<=n2 and k1<=k<=k2.
CheckWZpairEmp:=proc(F,G,n,n1,n2,k,k1,k2) local expr:
	expr:=simplify(expand((subs(n=n+1,F)-subs(k=k+1,G)+G))/F-1):
	evalb({seq(seq(simplify(expand(subs({n=nn,k=kk},expr))),nn=max(0,n1)..n2),kk=k1..k2)}={0});
end:


####################################################################################################
#Shadow(F,n,k): converts the “normal” function f to its shadow form.
#Each factor (a*n+b*k+c)! is replaced by (-1)^(a*n+b*k+c)/(-a*n-b*k-c-1)!
Shadow:=proc(F,n,k) local facts,newF,fact,fac:
	facts:=factors(F):
	newF:=facts[1]:
	for fact in facts[2] do
		if type(fact[1],`!`) then
			fac:=op(fact[1]):
			if coeff(fac,n,1)+coeff(fac,k,1)=0 then
				newF:=newF*fact[1]^fact[2]:
			else
				newF:=newF*((-1)^fact[2])^(fac)/((-fac-1)!)^fact[2]:
			fi:
		else
			newF:=newF*fact[1]^fact[2]:
		fi:
	od:
	convert(simplify(newF),factorial):
end:

####################################################################################################
#Dual(F,n,nShift,k,kShift): inputs a hypergeometric function. Outputs a dual function obtained by shadowing and shifting.
#Mainly used as a helper function for DualPair(F,G,n,k).
Dual:=proc(F,n,nShift,k,kShift):
	convert(simplify(subs({n=-k-kShift,k=-n-nShift},Shadow(F,n,k))),factorial):
end:
#DualPair(F,G,n,k): inputs a WZ-pair. Outputs the dual pair F’,G’.
#Can it check that G(n,k)/F(n,k-1) is a rational function of n,k? This is necessary to ensure the dual is a WZ pair.
DualPair:=proc(F,G,n,k):
	Dual(G,n,0,k,1),Dual(F,n,1,k,0):
end:


#####################################################Extra q-functions.
#####################################################
#qfactor(F): Takes a polynomial, F, in powers of q, and attempts to factor.
qfactor:=proc(F) local newF:
	newF:=factors(expand(simplify(F))):
	simplify(newF[1])*mul(simplify(fact[1])^fact[2],fact in newF[2]):
end:
#####################################################
#qfac1(n): Shorthand for qfac(n) in product form.
qfac1:=proc(n):
	Product(1-q^i,i=1..n):
end:
###############################
#evalqfac(n): Inputs n. Outputs qfac(n) evaluated if n is an integer. Otherwise outputs qfac(n).
evalqfac:=proc(n):
	if type(n,integer) or n=-infinity then
		if n<0 then
			infinity:
		else
			mul(1-q^i,i=1..n):
		fi:
	else
		qfac(n):
	fi:
end:
###############################
#ConvertToExpr(k,A,B,z1,q1:=q): Converts a hypergeometric series given as phi(A;B;q1,z1) into a hypergeometric term.
#A,B are lists. k is the summation variable.
ConvertToExpr:=proc(k,A,B,z1,q1:=q):
	mul(qrf(a,k),a in A)*((-1)^k*q1^(k*(k-1)/2))^(1+nops(B)-nops(A))*z1^k/(mul(qrf(b,k),b in B)*qrf(q1,k)):
end:


#q-Functions for simplifying expressions.
########################################################
#simplifyqfac(F,nk,nkVal): Takes limit of F as nk->nkVal. Subject to |q|<1 so q^infinity=0.
simplifyqfac:=proc(F,nk,nkVal) local facts,newF1,newF2,numqfacinfinity,fact,newValue,L,dummy:
	facts:=factors(F):
	newF1:=1:
	newF2:=facts[1]:
	numqfacinfinity:=0:
	for fact in facts[2] do
		if qfac(op(fact[1]))=fact[1] then #Then fact[1] is qfac(something).
			newValue:=subs(nk=nkVal,op(fact[1])):
			L:=convert(newValue+dummy,set) minus {dummy}:#dummy is added so that op(fact[1]) is broken up by sum and not multiplication.
			if member(-infinity,L) then
				numqfacinfinity:=numqfacinfinity+fact[2]:
			elif type(newValue,integer) then
				if newValue<0 then
					newF2:=newF2/mul(1-q^(op(fact[1])+i),i=1..-newValue)^fact[2]:
				else
					newF1:=newF1*mul(1-q^i,i=1..newValue)^fact[2]:
				fi:
			elif member(infinity,L) then
				newF1:=newF1*qfac(infinity)^fact[2]:
			else
				newF1:=newF1*qfac(newValue)^fact[2]:
			fi:
		else
			newF2:=newF2*fact[1]^fact[2]:
		fi:
	od:
	newF2:=limit(simplify(newF2),nk=nkVal) assuming abs(q)<1:
	if op(0,newF2)=limit and numqfacinfinity>=0 then
		undefined:
	else
		infinity^numqfacinfinity*simpqfac(newF1)*newF2:
	fi:
end:

########################################################
#simplifyqfacsum(F,nkSumvar,lowerBound,upperBound): Sums F from lowerBound to upperBound.
#It first checks whether the sum is actually finite.
simplifyqfacsum:=proc(F,nkSumvar,lowerBound,upperBound) local newF,kValmin,kValmax,fact,coe,newBound:
	kValmin:=ceil(lowerBound):
	kValmax:=floor(upperBound):
	for fact in factors(F)[2] do
		if fact[2]<0 and qfac(op(fact[1]))=fact[1]  then
			coe:=coeff(op(fact[1]),nkSumvar):
			if coe<>0 then
				newBound:=solve(op(fact[1])=0,nkSumvar):
				if coe<0 then
					kValmax:=min(kValmax,floor(newBound)):
				else
					kValmin:=max(kValmin,ceil(newBound)):
				fi:
			fi:
		fi:
	od:
	if type(kValmin,constant) and abs(kValmin)<>infinity and type(kValmax,constant) and abs(kValmax)<>infinity then
		simpqfac(add(normqfac(simplifyqfac(F,nkSumvar,kVal)),kVal=kValmin..kValmax)):
	else
		Sum(F,nkSumvar=kValmin..kValmax):
	fi:
end:

########################################################
#simpqfac(F): Simplifies F as much as possible.
#Combines (1-q^k)*qfac(k-1)=qfac(k), qfac(n)/qfac(n-1)=1-q^n, and qfac(k)/(1-q^k)=qfac(k-1).
#In the end there will be no qfac() with the opposite power sign that differs only by a constant.
#Any factor of (1-q^i) that could be multiplied into (or divided from) a qfac will have that done.
#Run through factors. If qfac, try looking at other factors.
#Check all other qfacs. If one differs by only a constant (and opposite power), cancel them, leaving remainder.
#Check all other factors. If possible, combine with current qfac (either increasing or decreasing argument). Then try again.
#Everything should be in a qfac.
simpqfac:=proc(F) local qfacSet,polySet,facts,fact,newF,newF2,furtherTest,factqfac,difference,newPower,polyToAdd,polytoadd,factPoly,dummy:
	qfacSet:={}:
	polySet:={}:
	facts:=factors(F):
	newF:=facts[1]:
	newF2:=1:
	for fact in facts[2] do #Initial run through.
		if qfac(op(fact[1]))=fact[1] then #Then fact[1]=qfac()
			qfacSet:=qfacSet union {fact}:
		else
			factPoly:=convert(fact[1]+dummy,set) minus {dummy}:
			if nops(factPoly)=1 then #This factor is not of the form 1-q^a.
				newF:=newF*fact[1]^fact[2]:
			elif nops(factPoly)>2 then #This factor is not of the form 1-q^a. Let’s try to factor it.
				newF2:=newF2*qfactor(fact[1])^fact[2]:
			else
				newF2:=newF2*fact[1]^fact[2]:
			fi:
		fi:
	od:
	facts:=factors(newF2): #The reason this is done twice is so that polynomials can be factored and combined properly.
	newF:=newF*facts[1]:
	newF2:=1:
	for fact in facts[2] do
		factPoly:=convert(fact[1]+dummy,set) minus {dummy}:
		if nops(factPoly)=1 then #This factor is not of the form 1-q^a.
			newF:=newF*fact[1]^fact[2]:
		elif nops(factPoly)>2 then #This factor is not of the form 1-q^a. The program has already attempted to factor it.
			newF:=newF*factPoly[1]^fact[2]:
			newF2:=newF2*simplify(fact[1]/factPoly[1])^fact[2]:
		else
			if member(1,factPoly) then
				polytoadd:=fact[1]:
			elif member(-1,factPoly) then
				polytoadd:=-fact[1]:
				newF:=newF*(-1)^fact[2]:
			else
				newF:=newF*factPoly[1]^fact[2]:
				polytoadd:=simplify(fact[1]/factPoly[1]):
			fi:
			if polytoadd=1-q or (nops(1-polytoadd)=2 and polytoadd=1-q^(op(2,1-polytoadd))) then
				polySet:=polySet union {[polytoadd,fact[2]]}:
#			elif polytoadd=-1+q or (nops(1+polytoadd)=2 and polytoadd=-1+q^(op(2,polytoadd+1))) then #This shouldn’t happen anymore.
#				polySet:=polySet union {[-polytoadd,fact[2]]}:#So that polytoadd=1-q^n in the future. Maple syntax makes it -1+q^n when factoring.
#				newF:=newF*(-1)^fact[2]:
			else
				newF2:=newF2*polytoadd^fact[2]:
			fi:
		fi:
	od:
	while qfacSet<>{} do
		fact:=qfacSet[1]:
		qfacSet:=qfacSet minus {fact}:
		furtherTest:=true:
		for factqfac in qfacSet do
			difference:=op(fact[1])-op(factqfac[1]):
			if type(difference,integer) and sign(fact[2])<>sign(factqfac[2]) then
				newPower:=min(abs(fact[2]),abs(factqfac[2])):
				if difference>0 then #Should not have difference=0 because that would have been simplified.
					polyToAdd:={seq([1-q^(op(factqfac[1])+i),sign(fact[2])*newPower],i=1..difference)}:
				else
					polyToAdd:={seq([1-q^(op(fact[1])+i),sign(factqfac[2])*newPower],i=1..-difference)}:
				fi:
				for factPoly in polySet do
					for polytoadd in polyToAdd do
						if polytoadd[1]=factPoly[1] then
							polySet:=polySet minus {factPoly}:
							if factPoly[2]+polytoadd[2]<>0 then
								polySet:=polySet union {[factPoly[1],factPoly[2]+polytoadd[2]]}:
							fi:
							polyToAdd:=polyToAdd minus {polytoadd}:
							break:
						fi:
					od:
				od:
				polySet:=polySet union polyToAdd:
				qfacSet:=qfacSet minus {factqfac}:
				if newPower<abs(factqfac[2]) then
					qfacSet:=qfacSet union {[factqfac[1],factqfac[2]-sign(factqfac[2])*newPower]}:
				fi:
				fact:=[fact[1],fact[2]-sign(fact[2])*newPower]:
				if fact[2]=0 then
					furtherTest:=false:
					break:
				fi:
			fi:
		od:
		while furtherTest do
			furtherTest:=false:
			for factPoly in polySet do
				if 1-q^(op(fact[1])+1)=factPoly[1] and sign(fact[2])=sign(factPoly[2]) then
					newPower:=sign(fact[2])*min(abs(fact[2]),abs(factPoly[2])):
					newF2:=newF2*fact[1]^(fact[2]-newPower):
					fact[1]:=qfac(op(fact[1])+1):
					fact[2]:=newPower:
					polySet:=polySet minus {factPoly}:
					if factPoly[2]<>newPower then
						polySet:=polySet union {[factPoly[1],factPoly[2]-newPower]}:
					fi:
					furtherTest:=true:
					break:
				elif 1-q^(-op(fact[1])-1)=factPoly[1] and sign(fact[2])=sign(factPoly[2]) then
					newPower:=sign(fact[2])*min(abs(fact[2]),abs(factPoly[2])):
					newF:=newF*(-q^(op(fact[1])+1))^(-newPower):
					newF2:=newF2*fact[1]^(fact[2]-newPower):
					fact[1]:=qfac(op(fact[1])+1):
					fact[2]:=newPower:
					polySet:=polySet minus {factPoly}:
					if factPoly[2]<>newPower then
						polySet:=polySet union {[factPoly[1],factPoly[2]-newPower]}:
					fi:
					furtherTest:=true:
					break:
				elif 1-q^(op(fact[1]))=factPoly[1] and sign(fact[2])<>sign(factPoly[2]) then #Both cannot happen as those polynomials would have canceled.
					newPower:=sign(fact[2])*min(abs(fact[2]),abs(factPoly[2])):
					newF2:=newF2*fact[1]^(fact[2]-newPower):
					fact[1]:=qfac(op(fact[1])-1):
					fact[2]:=newPower:
					polySet:=polySet minus {factPoly}:
					if factPoly[2]+newPower<>0 then
						polySet:=polySet union {[factPoly[1],factPoly[2]+newPower]}:
					fi:
					furtherTest:=true:
					break:
				elif 1-q^(-op(fact[1]))=factPoly[1] and sign(fact[2])<>sign(factPoly[2]) then #Both cannot happen as those polynomials would have canceled.
					newPower:=sign(fact[2])*min(abs(fact[2]),abs(factPoly[2])):
					newF:=newF*(-q^(op(fact[1])))^newPower:
					newF2:=newF2*fact[1]^(fact[2]-newPower):
					fact[1]:=qfac(op(fact[1])-1):
					fact[2]:=newPower:
					polySet:=polySet minus {factPoly}:
					if factPoly[2]+newPower<>0 then
						polySet:=polySet union {[factPoly[1],factPoly[2]+newPower]}:
					fi:
					furtherTest:=true:
					break:
				fi:
			od:
		od:
		newF:=newF*fact[1]^fact[2]:
	od:
	simplify(newF)*newF2*mul(factPoly[1]^factPoly[2],factPoly in polySet):
end:

########################################################
#normqfac(F): Replaces qfac(n+2)=qfac(n)*(1-q^(n+1))*(1-q^(n+2)) so that you can potentially add/subtract products of qfac().
#Need that F does not have any factors of the form qfac(n-1)-qfac(k+1).
normqfac:=proc(F) local facts,newF,fact,changed,L,l,dummy:
	facts:=factors(F):
	newF:=facts[1]:
	for fact in facts[2] do
		if qfac(op(fact[1]))=fact[1] then #Then fact[1]=qfac()
			changed:=false:
			L:=convert(op(fact[1])+dummy,set) minus {dummy}:#dummy is added so that op(fact[1]) is broken up by sum and not multiplication.
			for l in L do
				if type(l,integer) then
					newF:=newF*qfac(op(fact[1])-l)^fact[2]:
					if l>0 then
						newF:=newF*mul(1-q^(op(fact[1])-i),i=0..l-1)^fact[2]:
					else
						newF:=newF/mul(1-q^(op(fact[1])+i),i=1..-l)^fact[2]:
					fi:
					changed:=true:
					break:
				fi:
			od:
			if not changed then
				newF:=newF*fact[1]^fact[2]:
			fi:
		else
			newF:=newF*fact[1]^fact[2]:
		fi:
	od:
	newF:
end:


####################################################################################################
#CheckqWZpair(F,G,n,k): inputs F,G as hypergeometric expressions of q^n,q^k. Certifies that F,G are a q-WZ pair.
CheckqWZpair:=proc(F,G,n,k):
	evalb(simplify((normqfac(subs(n=n+1,F))-normqfac(G)+normqfac(subs(k=k-1,G)))/normqfac(F)-1)=0):
end:
#CheckqWZpairEmp(F,G,n,n1,n2,k,k1,k2): inputs F,G as hypergeometric expressions of q^n,q^k.
#Empirically checks that F,G are a q-WZ pair for all n1<=n<=n2 and k1<=k<=k2.
CheckqWZpairEmp:=proc(F,G,n,n1,n2,k,k1,k2) local expr:
	expr:=simpqfac((normqfac(subs(n=n+1,F))-normqfac(G)+normqfac(subs(k=k-1,G)))/normqfac(F)-1):
	evalb({seq(seq(simplifyqfac(simplifyqfac(expr,n,nn),k,kk),nn=max(0,n1)..n2),kk=k1..k2)}={0}):
end:




####################################################################################################
#qShadow(F,n,k): converts the q-function f to its shadow form.
#Each factor qfac(a*n+b*k+c) is replaced by (-1)^(a*n+b*k+c)*q^((a*n+b*k+c)*(a*n+b*k+c+1)/2)/qfac(-a*n-b*k-c-1)
qShadow:=proc(F,n,k) local facts,newF,newFsign,newFqpower,fact,fac:
	facts:=factors(F):
	newF:=facts[1]:#This is broken up because of the simplify procedure.
	newFsign:=0:
	newFqpower:=0:
	for fact in facts[2] do
		if qfac(op(fact[1]))=fact[1] then
			fac:=op(fact[1]):
			if coeff(fac,n,1)+coeff(fac,k,1)=0 then
				newF:=newF*fact[1]^fact[2]:
			else
				newFsign:=newFsign+fact[2]*fac mod 2:
				newFqpower:=newFqpower+fac*(fac+1)*fact[2]/2:
				newF:=newF/(qfac(-fac-1))^fact[2]:
			fi:
		else
			newF:=newF*fact[1]^fact[2]:
		fi:
	od:
	newF*q^simplify(newFqpower)*(-1)^(newFsign):
end:

####################################################################################################
#qDual(F,n,nShift,k,kShift): inputs a hypergeometric function. Outputs a dual function obtained by shadowing and shifting.
#Mainly used as a helper function for DualPair(F,G,n,k).
qDual:=proc(F,n,nShift,k,kShift):
	subs({n=-k-kShift,k=-n-nShift},qShadow(F,n,k)):
end:
#qDualPair(F,G,n,k): inputs a WZ-pair. Outputs the dual pair F’,G’.
#Can it check that G(n,k)/F(n,k-1) is a rational function of q^n,q^k? This is necessary to ensure the dual is a WZ pair.
qDualPair:=proc(F,G,n,k):
	qDual(subs(k=k-1,G),n,0,k,1),subs(k=k+1,qDual(F,n,1,k,0)):
end:


####################################################################################################
#qWZpairIdentities(F,G,n,k,constants,toprint:=false): inputs a qWZ-pair and a list of constants.
#Outputs original and companion identities as well as the limiting cases (if non-trivial).
qWZpairIdentities:=proc(F,G,n,k,constants,toprint:=false) local Fninfinity,Gkinfinity,Identities,constant,con,identity,upperBound,upperLimit,lowerBound,lowerLimit,FGtemp:
	Fninfinity:=simplifyqfac(F,n,infinity):#Computed often enough that I might as well assign it to a local variable.
	Gkinfinity:=simplifyqfac(G,k,infinity):

	#One of the first 2 conditions should be met since F,G are a qWZ pair.
	#Condition that G(n,infinity)-G(n,-infinity)=0.
	#Produces the identity Sum(F,k=-infinity..infinity)=constant.
	if simplify(normqfac(Gkinfinity)-normqfac(simplifyqfac(G,k,-infinity)))=0 then
		constant:=simplifyqfac(F,n,0):
		for con in constants do
			if subs({n=con,con=n},F)=F then
				constant:=simplifyqfac(constant,con,0):
			fi:
		od:
		constant:=simplifyqfacsum(constant,k,-infinity,infinity):
		identity:=simplifyqfacsum(F,k,-infinity,infinity)=constant assuming n::integer,n>0:
	#Condition that G(n,n)-G(n,-1)=-F(n+1,n+1).
	#Produces the identity Sum(F,k=0..n)=constant.
	elif simplify(normqfac(simplifyqfac(G,k,n))-normqfac(simplifyqfac(G,k,-1))+normqfac(simplifyqfac(subs(n=n+1,F),k,n+1)))=0 then
		constant:=simplifyqfac(F,n,0):
		for con in constants do
			if subs({n=con,con=n},F)=F then
				constant:=simplifyqfac(constant,con,0):
			fi:
		od:
		constant:=simplifyqfac(constant,k,0):
		identity:=simplifyqfacsum(F,k,0,n)=constant assuming n::integer,n>0:
	#Can at least print this out which might produce something useful. Generically summing and telescoping.
	#Produces the identity Sum(F(n+1,k)-F(n,k),k=0..n)=G(n,n)-G(n,-1).
	else
		identity:=simplifyqfacsum(simpqfac(normqfac(subs(n=n+1,F))-normqfac(F)),k,0,n)=simpqfac(normqfac(simplifyqfac(G,k,n))-normqfac(simplifyqfac(G,k,-1))) assuming n::integer:
	fi:
	if toprint then
		print(identity):
	fi:
	Identities:=[identity]:

	#Now we take the limit of identity as n->infinity.
	if Fninfinity<>undefined and Fninfinity<>infinity then
		if op(0,op(1,identity))=Sum then
			upperBound:=op(2,op(2,op(2,op(1,identity)))):
			upperLimit:=limit(upperBound,n=infinity):
			lowerBound:=op(1,op(2,op(2,op(1,identity)))):
			lowerLimit:=limit(lowerBound,n=infinity):
			constant:=simplifyqfac(op(2,identity),n,infinity):
			if (upperBound<>infinity and upperLimit=infinity) or (abs(lowerBound)<>infinity and abs(lowerLimit)=infinity) then
				FGtemp:=simplifyqfac(subs(k=upperBound-lowerBound-k,F),n,infinity) assuming k::real:
				if FGtemp<>undefined then
					if FGtemp=0 or subs(k=upperBound-lowerBound-k,FGtemp)=Fninfinity then
						identity:=Sum(Fninfinity,k=lowerLimit..upperLimit)=constant:
					elif constant=0 then
						identity:=Sum(Fninfinity,k=lowerLimit..upperLimit)=Sum(-FGtemp,k=lowerLimit..upperLimit):
					else
						identity:=Sum(Fninfinity+FGtemp,k=lowerLimit..upperLimit)=2*constant:
					fi:
				else
					identity:=0=0:
				fi:
			else
				identity:=Sum(Fninfinity,k=lowerLimit..upperLimit)=constant:
			fi:
		else
			identity:=simplifyqfac(op(1,identity),n,infinity)=simplifyqfac(op(2,identity),n,infinity):
		fi:
	else
		identity:=0=0:
	fi:
	if toprint then
		print(identity):
	fi:
	Identities:=[op(Identities),identity]:
	

	#Condition that F(infinity,k)-F(0,k)=0.
	#Produces the identity Sum(G,n=0..infinity)=constant.
	if simplify(normqfac(Fninfinity)-normqfac(simplifyqfac(F,n,0)))=0 then
		constant:=simplifyqfac(G,k,0):
		for con in constants do
			if subs({k=con,con=k},G)=G then
				constant:=simplifyqfac(constant,con,0):
			fi:
		od:
		constant:=simplifyqfacsum(constant,n,0,infinity):
		identity:=simplifyqfacsum(G,n,0,infinity)=constant assuming k::integer,k>=0:
	#Condition that F(infinity,k)-F(k,k)=-G(k-1,k-1).
	#Produces the identity Sum(G,n=k..infinity)=constant.
	elif simplify(normqfac(Fninfinity)-normqfac(simplifyqfac(F,n,k))+normqfac(simplifyqfac(subs(k=k-1,G),n,k-1)))=0 then
		constant:=simplifyqfac(G,k,0):
		for con in constants do
			if subs({k=con,con=k},G)=G then
				constant:=simplifyqfac(constant,con,0):
			fi:
		od:
		constant:=simplifyqfacsum(constant,n,0,infinity):
		identity:=simplifyqfacsum(G,n,k,infinity)=constant assuming k::integer,k>=0:
	#Condition that G(n,-1)=0. Technically only need condition that Sum(G(n,-1),n=0..infinity)=0. But that is much harder to compute.
	#Produces identity Sum(G,n=0..infinity)=Sum(F(infinity,j)-F(0,j),j=0..k).
	elif simplifyqfac(G,k,-1)=0 and Fninfinity<>undefined then
		constant:=simplifyqfacsum(subs(k=j,Fninfinity),j,0,k)-simplifyqfacsum(subs(k=j,simplifyqfac(F,n,0)),j,0,k) assuming k::integer,k>=0:
		identity:=simplifyqfacsum(G,n,0,infinity)=constant assuming k::integer,k>=0:
	#Can at least print this out which might produce something useful. Generically summing and telescoping.
	#Produces the identity Sum(G(n,k)-G(n,k-1),n=0..infinity)=F(infinity,k)-F(0,k).
	elif Fninfinity<>undefined then
		constant:=simpqfac(normqfac(Fninfinity)-normqfac(simplifyqfac(F,n,0))):
		identity:=Sum(simpqfac(normqfac(G)-normqfac(subs(k=k-1,G))),n=0..infinity)=constant:
	else
		identity:=0=0:
	fi:
	if toprint then
		print(identity):
	fi:
	Identities:=[op(Identities),identity]:

	#Now we take the limit of identity as k->infinity.
	if Gkinfinity<>undefined and Gkinfinity<>infinity then
		if op(0,op(1,identity))=Sum then
			upperBound:=op(2,op(2,op(2,op(1,identity)))):
			upperLimit:=limit(upperBound,k=infinity):
			lowerBound:=op(1,op(2,op(2,op(1,identity)))):
			lowerLimit:=limit(lowerBound,k=infinity):
			constant:=simplifyqfac(op(2,identity),k,infinity):
			if (upperBound<>infinity and upperLimit=infinity) or (abs(lowerBound)<>infinity and abs(lowerLimit)=infinity) then
				FGtemp:=simplifyqfac(subs(n=n+lowerBound,G),k,infinity) assuming n::real:
				if FGtemp<>undefined then
					lowerBound:=0:
					if Gkinfinity<>0 then
						if FGtemp=0 or FGtemp=Gkinfinity then
							identity:=simplifyqfacsum(Gkinfinity,n,lowerBound,upperLimit)=constant:
						elif constant=0 then
							identity:=simplifyqfacsum(Gkinfinity,n,lowerBound,upperLimit)=simplifyqfacsum(-FGtemp,n,lowerBound,upperLimit):
							else
							identity:=simplifyqfacsum(Gkinfinity+FGtemp,n,lowerBound,upperLimit)=2*constant:
						fi:
					else
						if FGtemp=0 then
							identity:=0=constant:
						else
							identity:=simplifyqfacsum(FGtemp,n,lowerBound,upperLimit)=constant:
						fi:
					fi:
				else
					identity:=0=0:
				fi:
			else
				if Gkinfinity<>0 then
					identity:=simplifyqfacsum(Gkinfinity,n,lowerLimit,upperLimit)=constant:
				else
					identity:=0=constant:
				fi:
			fi:
		else
			identity:=simplifyqfac(op(1,identity),k,infinity)=simplifyqfac(op(2,identity),k,infinity):
		fi:
	else
		identity:=0=0:
	fi:
	if toprint then
		print(identity):
	fi:
	[op(Identities),identity]:
end:
####################################################################################################
#qIdentities(F,n,k,constants,toprint:=false): inputs a hypergeometric function in q^n,q^k (and a list of constants). Uses qzeil to find pair. If a WZ-pair was found, it continues.
#Outputs identity, companion, dual, and dual companion as well as the limiting identities (if they are nontrivial).
#Also outputs F,G and the dual pair. These can be used to certify the identities.
#For identities while the program is running, set print to true.
qIdentities:=proc(F,n,k,constants,toprint:=false,checkDual:=false) local N,G,Identities,dualPair:
	G:=qzeil(F,N,k,n,constants):
	if G[1]<>N-1 then
		return(`qWZ pair not found.`):
	fi:
	G:=simpqfac(simplify(F*G[2])):
	#We have F(n+1,k)-F(n,k)=G(n,k)-G(n,k-1) for all integer k and n>=0.
	Identities:=qWZpairIdentities(F,G,n,k,constants,toprint):

	######################################Now repeat for the dual pair.
	dualPair:=qDualPair(F,G,n,k):#qDualPair requires that F(n+1,k)-F(n,k)=G(n,k)-G(n,k-1).
	dualPair:=simpqfac(simplify(dualPair[1])),simpqfac(simplify(dualPair[2])) assuming k::integer,n::integer:
	if checkDual and not CheckqWZpair(dualPair,n,k) then
		print(`Possible error constructing the dual pair.`):
		print(`To force construction of identities using ` ,dualPair):
		print(`qWZpairIdentities(`, dualPair,n,k,constants,toprint ,`)`):
		print(`At least you have the following original identities:`):
		[Identities,[F,G]]:
	else
		#We have dualPair[1](n+1,k)-dualPair[1](n,k)=dualPair[2](n,k)-dualPair[2](n,k-1) for all integer k and n>=0.
		[Identities,[F,G]],[qWZpairIdentities(dualPair[1],dualPair[2],n,k,constants,toprint),[dualPair]]:
	fi:
end:


####################################################################################################
#CheckqIdentity(ID,nkvar,nkval): Empirically checks the q-identity ID in the substitution nkvar=nkval.
#It is highly recommended to trim any bounds of summation to ignore zero terms, i.e. anything with 1/qfac(<0).
#ID should be an equation. If you have an expression, input ID=0.
CheckqIdentity:=proc(ID,nkvar,nkval) local expr:
	expr:=taylor(simplify(eval(subs({infinity=2*nkval,nkvar=nkval,qfac=evalqfac,Sum=add},op(1,ID)-op(2,ID)))),q,nkval+1):
	#Inefficient way of checking whether this is O(q^#).
	if type(expr,series) and nops(expr)=2 and op(1,expr)=op(1,taylor(q^2,q,1)) then
		0:
	else
		expr:
	fi:
end:
####################################################################################################
#seqCheckqIdentity(ID,nkvar,nkvalmin,nkvalmax): Empirically checks the q-identity ID in the substitution nkvar=nkval for nkval=nkvalmin…nkvalmax.
#ID should be an equation.
seqCheckqIdentity:=proc(ID,nkvar,nkvalmin,nkvalmax):
	evalb({seq(CheckqIdentity(ID,nkvar,nkval),nkval=nkvalmin..nkvalmax)}={0}):
end:





###############################################################################################################################################
###############################################################################################################################################
###############################################################################################################################################
###############################################################################################################################################
##########################################################Taken from qEKHAD.txt################################################################
###################################Credit to the book "A=B" by Marko Petkovsek, Herb Wilf, and Doron Zeilberger.###############################
###############################################################################################################################################
###############################################################################################################################################
###############################################################################################################################################
###############################################################################################################################################
ezra:=proc():
	if args=NULL then
		print(`q-EKHAD: A Maple package for proving terminating basic hypergemetric`):
		print(` (alias q-binomial coeff.) and other q-identities`):
		print(`Handles only single-summation`):
		print(`For help with a specific prodecure, type "ezra(procedure_name);"`):
		print(`Contains procedures: `):
		print(`hafokh,qzeil`):
		print(`Warning: q is a global variable`):
		print(`The summand can be expressed in terms of qfac(k):=(q)_k,`):
		print(`or the Gaussian polynomials (alias q-binomial coefficients)`):
		print(`qbin(a,b):=qfac(a)/(qfac(b)*qfac(b-a)), `):
		print(`or qrf(q^(something),k):=(q^(something))_k, `):
		print(`where the "something" must be affine linear in its variables`):
	elif nops([args])=1 and op(1,[args])=`hafokh` then
		print(`hafokh(EXPRESSION,k,n,q,t,s): converts a q-expression in qfac(a*n+b*k+c) and rf(a,k)`):
		print(`to the format required in procedure qctp`):
		print(`hafokh  converts a q-expression involving powers of qfac(a*n+b*k+c)`):
		print(`and qrf(a,k)'s to the format required in procedure qctp`):
		print(`Try:`):
		print(`hafokh(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,q,t,s);`):
	elif nops([args])=1 and op(1,[args])=`qzeil` then
		print(`qzeil(SUMMAND,N,k, n,para)`):
		print(`Inputs SUMMAND in terms of q-factorials of the form`):
		print(`qfac(r):=(1-q)(1-q^2)...(1-q^r)`):
		print(`qbin(a,b):=qfac(a)/(qfac(b)qfac(a-b),`):
		print(` and the familiar q-rising factorial,`):
		print(` (something)_k:=(1-something)*(1-q*something)*...*(1-q^(k-1)*something),`):
		print(` which should be entered as  qrf(something,k).`):
		print(` times q raised to a quadratic in k and n. `):
		print(`It outputs the recurrence operator ope(N,n) followed by`):
		print(`the certificate R(n,k) such that  G(n,k):=F(n,k)*R(n,k)`):
		print(`satisfies ope(N,n)F(n,k)=G(n,k)-G(n,k-1)`):
		print(`All parameters are discrete, so if you have, say (a)_k, you must`):
		print(`change a to q^(alpha), say. k is the summation variable, n is the`):
		print(`the auxilliary variable, with respect to which the recurrence`):
		print(`is given, and para is a list of all the other (necessarily) discrete`):
		print(`variables present. On output, arbitrary constants are denoted by b1, b2, etc. `):
		print(`Examples: The left side of the Finite-Form Rogers-Ramanujan is given`):
		print(`by qzeil(q^(k^2)/qfac(k)/qfac(n-k),N,k,n,[]);,`):
		print(`( while the right side (after Peter Paule's symmetrization) is:`):
		print(` qzeil(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),N,k,n,[]);`):
		print(`Since the second-order recurrences outputted are identical (try it!)`):
		print(`we have a proof of Rogers-Ramanujan.`):
		print(`Another example: To prove q-Chu do:`):
		print(` qzeil(q^(k^2)*qbin(n,k)*qbin(m,k)/qbin(n+m,m),N,k,n,[m]);`):
		print(`Another example: To prove q-Dixon (after the Paule symmetrization) do:`):
		print(` qzeil((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(2*n,n+k)^3/qfac(3*n)*qfac(n)^3,N,k,n,[]);`):
		print(`and with three parameters`):
		print(`qzeil((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(a+b,a+k)*qbin(a+c,c+k)*qbin(b+c,b+k)/qfac(a+b+c)*qfac(a)*qfac(b)*qfac(c),A,k,a,[b,c]);`):
		print(`The syntax, once more, is: qzeil(SUMMAND,N,k, n,para); `):
	else
		print(`That is not an included function.`):
	fi:
end:
###############################################################################################################################################
###############################################################################################################################################
#qbin converts the q-binomial coefficient to q-factorials
qbin:=proc(a,b):
	qfac(a)/qfac(b)/qfac(a-b):
end:
###############################################################################################################################################
#qrf:converts (q^(linear expressions of variables and parameters))_k
#linear expression in parameters to an expression in qfac
qrf:=proc(qal,k) local qal1:
	qal1:=simplify(qal):
	if qal=q then
		RETURN(qfac(k)):
	fi:
	if not type(qal1,`^`) then
		ERROR(`Argument must be a power of q`):
	fi:
	if not op(1,qal1)=q then
		ERROR(`Argument must be a power of q`):
	fi:
	qfac(op(2,qal1)+k-1)/qfac(op(2,qal1)-1):
end:
###############################################################################################################################################
qzeil:=proc(SUMMAND,N,k, n,para) local ORDER1,ope,mu:
	mu:=subs(n=n+1,SUMMAND)/SUMMAND:
	mu:=normal(mu):
	mu:=subs(k=k+1,mu)/mu:
	mu:=normal(mu):
	if mu=1 then
		ERROR(`Summand is seperable`):
	fi:
	ope:=0:
	for ORDER1 from 1 while ope=0 do
		ope:=qzeilsp(SUMMAND,ORDER1,N, k, n,para):
	od:
	ope:
end:
###############################################################################################################################################
qzeilsp:= proc(SUMMAND,ORDER1,N, k, n,para) local s,t,kak,POL,NUM,DEN,ARG,QPOWER1,pip,pip1,pip2,i,coe,i1:
	kak:=hafokh(SUMMAND,k,n,q,t,s):
	POL:=kak[1]:NUM:=kak[2]:DEN:=kak[3]:ARG:=kak[4]:QPOWER1:=kak[5]:
	pip:=qctpnis(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para):
	if pip=0 then
		RETURN(0):
	fi:
	pip:=qctp(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para):
	if pip=0 or [pip]=[0,0] then
		RETURN(0):
	fi:
	pip:=subs({s=q^n,t=q^k},pip[1]),subs({s=q^n,t=q^k},pip[2]):
	pip1:=simplify(pip[1]):
	pip1:=add(coeff(pip1,N,i)*N^i,i=ldegree(pip1,N)..degree(pip1,N)):
	pip2:=expand(normal(simplify(pip[2]))):
	coe:=lcoeff(pip1,N):
	pip1:=add( factor(expand(coeff(pip1,N,i1)/coe))*N^i1,i1=ldegree(pip1,N)..degree(pip1,N)):
	pip2:=factor(expand(pip2/coe)):
	pip1,pip2:
end:
###############################################################################################################################################
#hafokh(EXPRESSION,k,n,q,t,s): converts a q-expression in qfac(a*n+b*k+c) and rf(a,k)
#to the format required in procedure qctp
#hafokh  converts a q-expression involving powers of qfac(a*n+b*k+c)
#and qrf(a,k)'s to the format required in procedure qctp
#Try:
#hafokh(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,q,t,s);
hafokh:=proc(bitui1,k,n,q,t,s) local gu,bitui,POL,NUM,DEN,ARG,QPOWER1,gu1,khezka,gu11,gu11n,gu11k,gu110,i,r:
	bitui:=factor(bitui1):
	NUM:=[]:
	DEN:=[]:
	ARG:=1:
	POL:=1:
	QPOWER1:=0:
	if nops(bitui)>1 then
		for i from 1 to nops(bitui) do
			gu:=op(i,bitui):
			if type(gu,`+`) then
				POL:=POL*gu:
				POL:=hafokhpol(POL,n,k,q,t,s):
				next:
			fi:
			if type(gu,`^`) then 
				gu1:=op(1,gu):
				khezka:=op(2,gu):
				if gu1=q then
					if not  degree(khezka,k)<=2  then
						ERROR(q,`can be only raised to a power a quadratic in`,k):
					fi:
					QPOWER1:=QPOWER1+khezka:
					next:
				fi:
				if not (type(normal(khezka/k),integer) or type(khezka,integer)) then
					ERROR(`Wrong exponent`):
				fi:
				if not type(khezka, integer) then
					if not type(khezka/k ,integer) then
						ERROR(`Exponents must be in`,k):
					else
						ARG:=ARG*gu1^(khezka/k):
					fi:
					next:
				fi:
				if type(khezka,integer) then
					if type(gu1,function) then
						if op(0,gu1)<>'qfac' then
							ERROR(`The only functions allowed are qfac and qrf and qbin`):
						fi:
						if op(0,gu1)='qfac' then
							gu11:=op(1,gu1):
							gu11n:=coeff(gu11,n,1):
							gu11k:=coeff(gu11,k,1):
							gu110:=expand(gu11-gu11n*n-gu11k*k):
							if not type(gu11n,integer) or not type(gu11k,integer) then
								ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
							fi:
							gu11:=t^gu11k*s^gu11n*q^gu110:
						fi:
						if khezka>0 then
							for r from 1 to khezka do
								NUM:=[op(NUM),gu11]:
							od:
						else
							for r from 1 to -khezka do
								DEN:=[op(DEN),gu11]:
							od:
						fi:
					fi:
					next:
				fi:
			fi:
			if type(gu,function) then
				if op(0,gu)<>'qfac'  then
					ERROR(`The only functions allowed are qfac and qrf and qbin`):
				fi:
				if  op(0,gu)='qfac' then
					gu11:=op(1,gu):
					gu11n:=coeff(gu11,n,1):
					gu11k:=coeff(gu11,k,1):
					gu110:=expand(gu11-gu11n*n-gu11k*k):
					if not type(gu11n,integer) or not type(gu11k,integer) then
						ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
					fi:
					gu11:=t^gu11k*s^gu11n*q^gu110:             
					NUM:=[op(NUM),gu11]:
				fi:
			fi:
		od:
	fi:
 
	gu:=bitui:
	if type(gu,`+`) then
		POL:=POL*gu:
		POL:=hafokhpol(POL,n,k,q,t,s):
	fi:
	if type(gu,`^`) then 
		gu1:=op(1,gu):
		khezka:=op(2,gu):
		if gu1=q then
			if not  degree(khezka,k)<=2  then
				ERROR(q,`can be only raised to a power a quadratic in`,k):
			fi:
			QPOWER1:=QPOWER1+khezka:
		fi:
		if not (type(normal(khezka/k),integer) or type(khezka,integer)) then
			ERROR(`Wrong exponent`):
		fi:
		if not type(khezka, integer) then
			if not type(khezka/k ,integer) then
				ERROR(`Exponents must be in`,k):
			else
				ARG:=ARG*gu1^(khezka/k):
			fi:
		fi:
		if type(khezka,integer) then
			if type(gu1,function) then
				if op(0,gu1)<>'qfac' then
					ERROR(`The only function allowed are qfac and qrf`):
				fi:
				if op(0,gu1)='qfac' then
					gu11:=op(1,gu1):
					gu11n:=coeff(gu11,n,1):
					gu11k:=coeff(gu11,k,1):
					gu110:=expand(gu11-gu11n*n-gu11k*k):
					if not type(gu11n,integer) or not type(gu11k,integer)  then
						ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
					fi:
					gu11:=t^gu11k*s^gu11n*q^gu110:             
				fi:
				if khezka>0 then
					for r from 1 to khezka do
						NUM:=[op(NUM),gu11]:
					od:
				else
					for r from 1 to -khezka do
						DEN:=[op(DEN),gu11]:
					od:
				fi:
			fi:
		fi:
	fi:

	if type(gu,function) then
		if op(0,gu)<>'qfac'  then
			ERROR(`The only functions allowed are qfac and qrf and qbin`):
		fi:
		if op(0,gu)='qfac' then
			gu11:=op(1,gu):
			gu11n:=coeff(gu11,n,1):
			gu11k:=coeff(gu11,k,1):
			gu110:=expand(gu11-gu11n*n-gu11k*k):
			if not type(gu11n,integer) or not type(gu11k,integer) then
				ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
			fi:
			gu11:=t^gu11k*s^gu11n*q^gu110:             
		fi:
		NUM:=[op(NUM),gu11]:
	fi:
	POL,NUM,DEN,ARG,QPOWER1:         
end:
###############################################################################################################################################
hafokhpol:=proc(POL,n,k,q,t,s):
	subs({q^k=t,q^n=s},expand(POL)):
end:
###############################################################################################################################################
qctpnis:= proc(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para)
local  ARG1,  DEN1,  KAMA, 
      KHE, L, L1, LU, LU1, NUM1, ORDER, P, Q, QPOWER, R, RAT, 
      RATB,  S, degg, eq, fu, g, gorem, gu, 
      gugu, i, ia1, j1, ja1, kak, khe, l1, l2, lu1a, meka1, mekb1, 
      mumu, va1, va2, y1, yakhas,k1,j2,kapi1,kapi2,misht,a,b,va1q,goremq,fuq:
	KAMA:=5;
	NUM1:=NUM:DEN1:=DEN:ARG1:=ARG:QPOWER:=QPOWER1:
	ORDER:=ORDER1:
	gorem:=0:
	for i from 0 to ORDER do
		gorem:=gorem+(b[i])*N^i:
	od:
	LU:=0:
	KHE:=degree(gorem,N):
	for khe from 0 to KHE do
		RAT:=product('ratn(NUM1[y1],khe,q,s)','y1'=1..nops(NUM1)):
		RATB:=product('ratn(DEN1[y1],khe,q,s)','y1'=1..nops(DEN1)):
		RAT:=RAT/RATB:
		kak:=subs(n=n+khe,QPOWER)-QPOWER:
		kak:=expand(kak):
		gu:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-n*coeff(kak,n,1)):
		RAT:=RAT*gu*subs(s=s*q^khe,POL):
		RAT:=normal(RAT):
		LU:=coeff(gorem,N,khe)*RAT+LU:
		LU:=normal(LU):
	od:
	lu1a:=LU:
	LU1:=numer(LU):
	LU:=1/denom(LU):
	LU:=LU/subs(t=t/q,LU):
	LU:=normal(LU):
	yakhas:=product('ratk(NUM1[y1],q,t)','y1'=1..nops(NUM1))/product('ratk(DEN1[y1],q,t)','y1'=1..nops(DEN1)):
	yakhas:=normal(yakhas):
	kak:=QPOWER-subs(k=k-1,QPOWER):
	kak:=expand(kak):
	kak:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-n*coeff(kak,n,1)):
	yakhas:=ARG1*LU*yakhas*kak:
	yakhas:=normal(yakhas):
	P:=LU1:
	Q:=numer(yakhas):
	R:=denom(yakhas):
	Q:=expand(Q):
	R:=expand(R):
	j1:=0:
	while j1 <= KAMA do
		kapi1:=numer(Q):
		kapi2:=numer(subs(t=t*q^j1,R)):
		kapi1:=expand(kapi1):
		kapi2:=expand(kapi2):
		for j2 from 1 to nops(para) do
			kapi1:=subs(q^op(j2,para)=misht[j2],kapi1):
			kapi2:=subs(q^op(j2,para)=misht[j2],kapi2):
		od:
		g:=gcd(kapi1,kapi2):
		if g <> 1 then
			for j2 from 1 to nops(para) do
				g:=subs(misht[j2]=q^op(j2,para),g):
			od:
			Q:=normal(Q/g):
			R:=normal(R/subs(t=t/q^j1,g)):
			P:=P*product(subs(t=t/q^ja1,g) , ja1=0..j1-1):
			j1:=-1:
		fi:
		j1:=j1+1:
	od:
	P:=expand(P):
	R:=expand(R):
	Q:=expand(Q):
	if(coeff(Q,t,0)*coeff(R,t,0)=0) then
		L1:=0:
	else
		L:=ln(normal(coeff(Q,t,0)/coeff(R,t,0)))/ln(q):
		if type(L,integer) then
			L1:=max(0,L):
		else
			L1:=0:
		fi:
	fi:
	P:=expand(t^L1*P):
	R:=expand(q^L1*R):
	l1:=degree(R,t):
	l2:=degree(Q,t):
	meka1:=coeff(R,t,l1):
	mekb1:=coeff(Q,t,l2):
	if l1 <>l2 then
		k1:=degree(P,t)-max(l1,l2):
	else
		mumu:= meka1/(mekb1*q^l1):
		mumu:=ln(mumu)/ln(q):
		if type(mumu,integer) then
			k1:=max(mumu, degree(P,t)-l1):
		else
			k1:=degree(P,t)-l1:
		fi:
	fi:
	fu:=0:
	if k1 < 0 then
		RETURN(0):
	fi:
	for ia1 from 0 to k1 do
		fu:=fu+a[ia1]*t^ia1:
	od:
	gugu:=subs(t=t*q,Q)*fu-R*subs(t=t/q,fu)-P:
	gugu:=expand(gugu):
	degg:=degree(gugu,t):
	for ia1 from 0 to degg do
		eq[ia1+1]:=coeff(gugu,t,ia1)=0:
	od:
	for ia1 from 0 to k1 do
		va1[ia1+1]:=a[ia1]:
	od:
	for ia1 from 0 to ORDER do
		va2[ia1+1]:=b[ia1]:
	od:
	eq:=convert(eq,set):
	va1:=convert(va1,set):
	va2:=convert(va2,set):
	va1:=va1 union va2 :

	for j2 from 1 to nops(para) do
		eq:=subs(misht[j2]=j2^2+1,eq):
	od:
	va1q:=solve(subs(q=2,eq),va1):
	goremq:=subs(va1q,gorem):
	if goremq=0 then
		RETURN(0):
	fi:
	va1:=solve(eq,va1):
	fu:=subs(va1,fu):
	gorem:=subs(va1,gorem):
	if gorem=0 then
		RETURN(0):
	else
		RETURN(1):
	fi:
end:
###############################################################################################################################################
qctp:= proc(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para)
local  ARG1,  DEN1,  KAMA, 
      KHE, L, L1, LU, LU1, NUM1, ORDER, P, Q, QPOWER, R, RAT, 
      RATB,  S, degg, eq, fu, g, gorem, gu, 
      gugu, i, ia1, j1, ja1, kak, khe, l1, l2, lu1a, meka1, mekb1, 
      mumu, va1, va2, y1, yakhas,k1,j2,kapi1,kapi2,misht,va11,VAK:
	KAMA:=5:
	NUM1:=NUM:DEN1:=DEN:ARG1:=ARG:QPOWER:=QPOWER1:
	ORDER:=ORDER1:
	gorem:=0:
	for i from 0 to ORDER do
		gorem:=gorem+b[i]*N^i:
	od:
	LU:=0:
	KHE:=degree(gorem,N):
	for khe from 0 to KHE do
		RAT:=product('ratn(NUM1[y1],khe,q,s)','y1'=1..nops(NUM1)):
		RATB:=product('ratn(DEN1[y1],khe,q,s)','y1'=1..nops(DEN1)):
		RAT:=RAT/RATB:
		kak:=subs(n=n+khe,QPOWER)-QPOWER:
		kak:=expand(kak):
		gu:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-n*coeff(kak,n,1)):
		RAT:=RAT*gu*subs(s=s*q^khe,POL):
		RAT:=normal(RAT):
		LU:=coeff(gorem,N,khe)*RAT+LU:
		LU:=normal(LU):
	od:
	lu1a:=LU:
	LU1:=numer(LU):
	LU:=1/denom(LU):
	LU:=LU/subs(t=t/q,LU):
	LU:=normal(LU):
	yakhas:=product('ratk(NUM1[y1],q,t)','y1'=1..nops(NUM1))/product('ratk(DEN1[y1],q,t)','y1'=1..nops(DEN1)):
	yakhas:=normal(yakhas):
	kak:=QPOWER-subs(k=k-1,QPOWER):
	kak:=expand(kak):
	kak:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-n*coeff(kak,n,1)):
	yakhas:=ARG1*LU*yakhas*kak:
	yakhas:=normal(yakhas):
	P:=LU1:
	Q:=numer(yakhas):
	R:=denom(yakhas):
	Q:=expand(Q):
	R:=expand(R):
	j1:=0:
	while j1 <= KAMA do
		kapi1:=numer(Q):
		kapi2:=numer(subs(t=t*q^j1,R)):
		kapi1:=expand(kapi1):
		kapi2:=expand(kapi2):
		for j2 from 1 to nops(para) do
			kapi1:=subs(q^op(j2,para)=misht[j2],kapi1):
			kapi2:=subs(q^op(j2,para)=misht[j2],kapi2):
		od:
		g:=gcd(kapi1,kapi2):
		if g <> 1 then
			for j2 from 1 to nops(para) do
				g:=subs(misht[j2]=q^op(j2,para),g):
			od:
			Q:=normal(Q/g):
			R:=normal(R/subs(t=t/q^j1,g)):
			P:=P*product(subs(t=t/q^ja1,g) , ja1=0..j1-1):
			j1:=-1:
		fi:
		j1:=j1+1:
	od:
	P:=expand(P):
	R:=expand(R):
	Q:=expand(Q):
	if(coeff(Q,t,0)*coeff(R,t,0)=0) then
		L1:=0:
	else
		L:=ln(normal(coeff(Q,t,0)/coeff(R,t,0)))/ln(q):
		if type(L,integer) then
			L1:=max(0,L):
		else
			L1:=0:
		fi:
	fi:
	P:=expand(t^L1*P):
	R:=expand(q^L1*R):
	l1:=degree(R,t):
	l2:=degree(Q,t):
	meka1:=coeff(R,t,l1):
	mekb1:=coeff(Q,t,l2):
	if l1 <>l2 then
		k1:=degree(P,t)-max(l1,l2):
	else
		mumu:= meka1/(mekb1*q^l1):
		mumu:=ln(mumu)/ln(q):
		if type(mumu,integer) then
			k1:=max(mumu, degree(P,t)-l1):
		else
			k1:=degree(P,t)-l1:
		fi:
	fi:
	fu:=0:
	if k1 < 0 then
		RETURN(0):
	fi:
	for ia1 from 0 to k1 do
		fu:=fu+a[ia1]*t^ia1:
	od:
	gugu:=subs(t=t*q,Q)*fu-R*subs(t=t/q,fu)-P:
	gugu:=expand(gugu):
	degg:=degree(gugu,t):
	for ia1 from 0 to degg do
		eq[ia1+1]:=coeff(gugu,t,ia1)=0:
	od:

	for ia1 from 0 to k1 do
		va1[ia1+1]:=a[ia1]:
	od:
	for ia1 from 0 to ORDER do
		va2[ia1+1]:=b[ia1]:
	od:
	eq:=convert(eq,set):
	va1:=convert(va1,set):
	va2:=convert(va2,set):
	va1:=va1 union va2 :
	va1:=solve(eq,va1):
	VAK:={}:
	for va11 in va1 do
		if op(1,va11)=op(2,va11) then
			VAK:=VAK union {op(1,va11)}:
		fi:
	od:
	if VAK={} then
		RETURN(FAIL):
	fi:

	fu:=subs(va1,fu):
	fu:=subs({seq(va11=1,va11 in VAK)},fu):
	gorem:=subs(va1,gorem):
	gorem:=subs({seq(va11=1,va11 in VAK)},gorem):
	fu:=normal(fu):
	S:=subs(t=t*q,Q)*fu*lu1a/P :
	S:=normal(S):
	S:=factor(S):
	gorem,S:
end:
###############################################################################################################################################
ratn:=proc(alpha,khe,q,s) local i1,lu,dg:
	dg:=degree(alpha,s)*khe:
	if dg=0 then
		RETURN(1):
	fi:
	if dg>0 then
		lu:=product(1-alpha*q^i1,i1=1..dg):
	fi:
	if dg <0 then
		lu:=1/product(1-alpha/q^i1,i1=0..(-dg)-1):
	fi:
	lu:=normal(lu):
	RETURN(lu):
end:
###############################################################################################################################################
ratk:=proc(alpha,q,t) local lu,dg,i1:
	dg:=degree(alpha,t):
	if dg=0 then
		RETURN(1):
	elif
		dg > 0 then
		lu:=product(1-alpha/q^i1,i1=0..dg-1):
	else
		lu:=1/product(1-alpha*q^i1,i1=1..(-dg)):
	fi:
	lu:=expand(lu):
	lu:=normal(lu):
	RETURN(lu):
end: