######################################################################
##qEKHAD: Save this file as  qEKHAD                                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read qEKHAD<Enter>                                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################

print(`This is qEKHAD, one of the Maple packages`):
print(`accompanying the classic book "A=B"`):
print(`(published by A.K.Peters Ltd. (that now became part of CRC press) )`):
print( `by Marko Petkovsek, Herb Wilf, and Doron Zeilberger. `):
print(``):
print(`The book is available freely from `):
print(`https://www.math.upenn.edu/~wilf/AeqB.html`):
print(``):
print(`The latest version of this package is available from`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/qEKHAD`):
print(``):

print(`---------------------------------------------`):
print(`For  help with supporting procedures, and a list of the available functions,`):
print(` type "ezra1();". For specific help type "ezra(procedure_name);"`):
lprint(``):

print(`---------------------------------------------`):
print(`For general help, and a list of the available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);"`):
lprint(``):

print(`---------------------------------------------`):
print(`Warning: q is a global variable`):
print(`qfac(k) means (1-q)*(1-q^2)*...*(1-q^k), and `):
print(`NOT (1)(1+q)*...*(1+q+...+q^(k-1))`):
print(` `):
print(`Version of  July 17, 2015.`):

print(`Thanks to Brunu Gauthier for spotting a bug in an earlier version`):

ezra1:=proc()
if args=NULL then
print(`The supporting procedures are: `):
print(` hafokh `):
else
ezra(args):
fi:

end:
 
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(`qzeil,qzeilpap`):
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`):
fi:



if 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);`):
fi:
  
if 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-Saalschutz do:`):
print(` qzeil( q^(n^2-k^2)*qfac(a+b+c-n)/qfac(a-n)/qfac(b-n)/qfac(c-n)/qfac(n+k)/qfac(n-k)/qbin(a+b,a+k)/qbin(a+c,c+k)/qbin(b+c,b+k),K,n,k,[a,b,c]);`):

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); `):
fi:
 
if nops([args])=1 and op(1,[args])=`qzeilpap` then
print(` qzeilpap(SUMMAND,k,n,para) or qzeilpap(SUMMAND,k,n,para,NAME)`):
print(`Just like qzeil(SUMMAND,N,k,n,para) but writes a paper with the proof`):
print(`NAME is optional`):
print(`Examples: The left side of the Finite-Form Rogers-Ramanujan is given`):
print(`by qzeilpap(q^(k^2)/qfac(k)/qfac(n-k),k,n,[]);`):
print(`( while the right side (after Peter Paule's symmetrization) is:`):
print(` qzeilpap(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,[]);`):
print(`Since the second-order recurrences outputted are identical (try it!)`):
print(`we have a proof of the First Rogers-Ramanujan identity.`):
print(`For the two above examples, with the nake indicated, try:`):
print(` qzeilpap(q^(k^2)/qfac(k)/qfac(n-k),k,n,[],"Left Side of Finite Rogers-Ramanujan"); `);
print(``):
print(` qzeilpap(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,[], "Right Side of Finite Rogers-Ramanujan" );`):


print(`Another example: To prove q-Chu do:`):
print(``):
print(` qzeilpap(q^(k^2)*qbin(n,k)*qbin(m,k)/qbin(n+m,m),k,n,[m]);`):
print(``):
print(`and with the name`):
print(``):
print(` qzeilpap(q^(k^2)*qbin(n,k)*qbin(m,k)/qbin(n+m,m),k,n,[m], "Chu-Vandermonde" );`):
print(``):
print(`Another example: for q-Dixon (after the Paule-symmetrization)`):
print(``):
print(` qzeilpap((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(2*n,n+k)^3/qfac(3*n)*qfac(n)^3,k,n,[]);`):


print(`Yet another example: for q-Dixon with three parameters(after the Paule-symmetrization)`):
print(``):
print(`qzeilpap((-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),k,a,[b,c]);`):
fi:




 
end:
 
 
 



qpaper4:=proc(SUMMAND,k,n,para) local SHEM,IDENTITY,RECURRENCE, gu, N,ope,SDZ1:


gu:=qzeil(SUMMAND,N,k, n,para):

if gu=FAIL then 
  RETURN(FAIL):
fi:

ope:=gu[1]:
SDZ1:=gu[2]:

if degree(ope,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:

print(``):
print(`A PROOF OF a q-Series `, SHEM):
lprint(``):
print(`by Shalosh B. Ekhad, Rutgers University`):
lprint(``):
print(`I will give a short proof of the following result.`):
print(``):
if degree(ope,N)=1 then
print(`(Note that since the recurrence below is first order, this`):
print(`means that the sum`, SUM(n), `has closed form,and it is`):
print(`easily seen to be equivalent.)`):
print(``):
fi:
print(`Theorem: Let `, F(n,k), ` be given by `):
lprint(``):
print(SUMMAND):
print(``):
print(`and let`, SUM(n),` be the sum of `, F(n,k),` with respect to `, k):
print(``):
print(SUM(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`PROOF: We cleverly construct`, G(n,k),`:=`):
print(``):
print(SDZ1*SUMMAND):
print(``):
print(`and in Maple format `):
print(``):
lprint(SDZ1*SUMMAND):
print(``):
lprint(`with the motive that`):
print(``):
print(yafe(ope,N,n,F(n,k))=G(n,k)-G(n,k-1) ):
print(`(check!)`):
print(``):
print(`and the theorem follows upon summing with respect to `, k, ` . `):
end:
 


qpaper5:=proc(SUMMAND,k,n,para,NAME) local SHEM,IDENTITY,RECURRENCE, gu, N,ope,SDZ1:


gu:=qzeil(SUMMAND,N,k, n,para):

if gu=FAIL then 
  RETURN(FAIL):
fi:

ope:=gu[1]:
SDZ1:=gu[2]:

if degree(ope,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:

print(``):
print(`A PROOF OF The`, NAME , ` q-Series  `, SHEM):
lprint(``):
print(`by Shalosh B. Ekhad, Rutgers University`):
lprint(``):
print(`I will give a short proof of the following result.`):
print(``):
if degree(ope,N)=1 then
print(`(Note that since the recurrence below is first order, this`):
print(`means that the sum`, SUM(n), `has closed form,and it is`):
print(`easily seen to be equivalent.)`):
print(``):
fi:
print(`Theorem: Let `, F(n,k), ` be given by `):
lprint(``):
print(SUMMAND):
print(``):
print(`and let`, SUM(n),` be the sum of `, F(n,k),` with respect to `, k):
print(``):
print(SUM(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`PROOF: We cleverly construct`, G(n,k),`:=`):
print(``):
print(SDZ1*SUMMAND):
print(``):
print(`and in Maple format `):
print(``):
lprint(SDZ1*SUMMAND):
print(``):
lprint(`with the motive that`):
print(``):
print(yafe(ope,N,n,F(n,k))=G(n,k)-G(n,k-1) ):
print(`(check!)`):
print(``):
print(`and the theorem follows upon summing with respect to `, k, ` . `):
end:
 


 




 
qzeilpap:=proc(SUMMAND,k,n,para,NAME):
 
if not (nargs=4 or nargs=5) then
 ERROR(` qzeilpap(SUMMAND,k,n,para) or qzeilpap(SUMMAND,k,n,para,NAME)`):
 fi:
 
if nargs=5 then
qpaper5(SUMMAND,k,n,para,NAME):
else
qpaper4(SUMMAND,k,n,para):
fi:
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:
 
pashet:=proc(p,N)
local i,gu,p1:
p1:=expand(p):
 
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu):
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:
 
 
 
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:
 
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:
 
#qbin converts the q-binomial coefficient to q-factorials
 
qbin:=proc(a,b):
qfac(a)/qfac(b)/qfac(a-b):
end:
 

hafokhpol:=proc(POL,n,k,q,t,s):
subs({q^k=t,q^n=s},expand(POL)):
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:
 
 
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:
 
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:
 
 
 
qfac1:=proc(r)
option remember:
 
if r=0 then
 1:
elif r<0 then
 0:
else
 expand(qfac1(r-1)*normal((1-q^r)/(1-q))):
fi:
end:
 
 
 
 
 
yafe:=proc(ope,N,n,SUM)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,N) do
 gu:=gu+coeff(ope,N,i)*subs(n=n+i,SUM):
od:
 
gu:
end: