print(`This is a copy of the package qEKHAD accompanying the book A=B`):
print(`by M. Petkovsek, H. Wilf an D. Zeilberger`):
print(`Because it was mentioned in Zeilberger's paper published in v.2`):
print(`of the New York J. of Math, it is reproduced here for the`):
print(`convenience of future readers, and to guarantee its permanence.`):
print(`Zeilberger would like to thank Mark Steinberger for his suggestion`):
print(`to install a permanent copy of the package in the NYJM's archives.`):
print(`This version works with Maple V vers. 3 and 2,and hopefully above`):
lprint(``):
lprint(``):



print(`This is qEKHAD, one of the Maple packages`):
print(`accompanying the forthcoming book "A=B"`):
print(`(soon to be published by A.K.Peters, Ltd.)`):
print( `by Marko Petkovsek, Herb Wilf, and Doron Zeilberger. `):
print(`Also try out the other packages: EKHAD and multiEKHAD`):
print(``):
print(`The latest versions are always available by anon. ftp to`):
print(`ftp.math.temple.edu`):
print(` in the directory: /pub/zeilberg/programs `):
print(`or, on the Web, at http://www.math.temple.edu/~zeilberg`):
print(`For general help, and a list of the available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);"`):
lprint(``):
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 1995`):
 
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])=`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(`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,REF)`):
print(`Just like qzeil(SUMMAND,k,n,para) but writes a paper with the proof`):
print(`NAME and REF are optional name and reference`):
fi:
 
 
end:
 
 
 
qpaper:=proc(SUMMAND,k,n,N,ope1,SDZ1,NAME,REF)
local SHEM,IDENTITY:
 
if degree(ope1,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:
lprint(`A PROOF OF THE`,NAME,SHEM):
lprint(``):
lprint(`by Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`):
lprint(``):
lprint(`I will give a short proof of the following result(`,REF,`).`):
lprint(``):
if degree(ope,N)=1 then
lprint(`(Note that since the recurrence below is first order, this`):
lprint(`means that the sum`, SUM(n), `has closed form,and it is`):
lprint(`easily seen to be equivalent.)`):
lprint(``):
fi:
lprint(`Theorem:Let`, F(n,k), `be`):
lprint(``):
print(SUMMAND):
lprint(``):
lprint(`and let`, SUM(n),`be the sum of`, F(n,k),` with`):
lprint(`respect to`, k,`.`):
lprint(``):
lprint(SUM(n),` satisfies the following linear recurrence equation`):
lprint(``):
print(yafe(ope1,N,n,SUM(n))):
print(`=0.`):
lprint(``):
lprint(`PROOF: We cleverly construct`, G(n,k),`:=`):
lprint(``):
print(SDZ1*SUMMAND):
lprint(`with the motive that`):
lprint(``):
print(yafe(ope1,N,n,F(n,k))):
lprint(`=`,G(n,k)-G(n,k-1), `(check!)`):
lprint(``):
lprint(`and the theorem follows upon summing with respect to`, k,`.`):
end:
 
 
qzeilpap:=proc()
local SUMMAND,k,para,N:
 
if nargs<4 or nargs>6 then
 ERROR(` qzeilpap(SUMMAND,k,n,para) or qzeilpap(SUMMAND,k,n,para,NAME,REF)`):
 fi:
 
if nargs>=4 then
SUMMAND:=args[1]:
k:=args[2]:
n:=args[3]:
para:=args[4]:
fi:
 
if nargs>=5  then
 NAME:=args[5]:
fi:
 
 
if nargs=6 then
 REF:=args[6]:
fi:
 
gu:=qzeil(SUMMAND,N,k,n,para):
ope1:=gu[1]:
SDZ1:=gu[2]:
qpaper(SUMMAND,k,n,N,ope1,SDZ1,NAME,REF):
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
:
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
#HERE
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):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
fu:=normal(fu):
 
#gorem:=pashet(gorem,N):
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
:
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
#HERE
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:
 
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:
 
#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
 
hafokhpol:=proc(POL,n,k,q,t,s):
subs({q^k=t,q^n=s},expand(POL)):
end:
 
 
#hafokh  converts a q-expression in qfac(a*n+b*k+c) and rf(a,k)
#to the format required in procedure qctp
 
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:
 
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]):
pip[1],factor(expand(normal(simplify(pip[2])))):
end:
 
qzeil:= proc(SUMMAND,N,k, n,para)
local ORDER1,ope:
 
ope:=0:
 
for ORDER1 from 1 while ope=0 do
ope:=qzeilsp(SUMMAND,ORDER1,N, k, n,para):
od:
 
pashet(ope[1],N),factor(ope[2]):
 
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:
 
bdokcert:=proc(SUMMAND,N,k, n,ope,cert,n1,k1)
local gu,G,SU,i,gu1:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*subs(n=n+i,SUMMAND):
od:
 
G:=SUMMAND*cert:
gu:=gu/(G-subs(k=k-1,G)):
 
  gu1:=subs({qfac=qfac1},gu):
  gu1:=subs({k=k1,n=n1},gu1):
normal(eval(gu1)):
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:

###end of qEKHAD