print(`Rogers-Ramanujan Tools Maple package...version 1.112, November 1, 2005`);
print(`[annihiloper, basic_calcs, bosedual, duallist, factorizeop, fermipoly, fermipolydual, fermipolyevendual, fermipolyodddual, fermipolylong, funcrecur, gp, initialconds, polylist, multops, polyrecur, printpolyseq, qfac, qfact, rightdiv, search, spqfac, t0, t1, T1, tau0, Trb, twovargen, U, V]`);
 
qfact := (a,q,n) -> product( (1-a*q^m),m=0..n-1):
PNS:= sort(add( (-1)^n * q^(n*(3*n-1)/2), n=-10..10),q):
SS:=  sort(add( (-1)^n * q^(n^2), n=-10..10),q):
TNS:= sort(add(    q^(n*(n+1)/2),n=0..10),q):
ATNS:=sort(add( (-q)^(n*(n+1)/2),n=0..10),q):
JTP:= (b,c,d) -> sort(sum( (-1)^(b*n) * q^(c*n^2 +d*n),n=-10..10),q):

gp := proc(A,B,k) local r;
   if type(A,integer) and type(B,integer) 
   then if 0 <= B and B<=A 
        then RETURN(subs(q=q^k, collect(simplify(
          product(1-q^r, r=A-B+1..A)/product(1-q^r, r=1..B)),q))) 
        else RETURN(0)
        fi
   else RETURN('gp(A,B,k)')
   fi
end:

#############################################################
# We define 7 q-analogs of the trinomial coefficient        # 
# plus some related functions.                              #
#############################################################
T0 := proc(m,A,k) 
  if type(m,integer) and type(A,integer) then
  RETURN(subs(q=q^k,sort(simplify(
    add( (-1)^j*gp(m,j,2)*gp(2*m-2*j,m-A-j,1),j=0..m)))))
  else 'T0(m,A,k)'
  fi
end:

T1:= proc(m,A,k) 
  if type(m,integer) and type(A,integer) then
  RETURN(subs(q=q^k,sort(simplify(
    add( (-q)^j*gp(m,j,2)*gp(2*m-2*j,m-A-j,1),j=0..m)))))
  else 'T1(m,A,k)'
  fi
end:


t0 := proc(m,A,k)   if type(m,integer) and type(A,integer) then
  RETURN(subs(q=q^k,sort(simplify(
    add( (-1)^j*q^(j^2)*gp(m,j,2)*gp(2*m-2*j,m-A-j,1),j=0..m)))))
  else 't0(m,A,k)'
  fi
end:

t1 :=proc(m,A,k) 
  if type(m,integer) and type(A,integer) then
  RETURN(subs(q=q^k,sort(simplify(
    add( (-1)^j*q^(j^2-j)*gp(m,j,2)*gp(2*m-2*j,m-A-j,1),j=0..m)))))
  else('t1(m,A,k)')
  fi
end:

tau0 :=proc(m,A,k) 
  if type(m,integer) and type(A,integer) then
  RETURN(subs(q=q^k,sort(simplify(
     add( (-1)^j*q^(m*j - (j^2-j)/2)*gp(m,j,1)*gp(2*m-2*j,m-A-j,1),j=0..m)))))
  else('tau0(m,A,k)')
  fi
end:


V:=proc(m,A,k)
   if type(m,integer) and type(A,integer) then
    RETURN(subs(q=q^k,sort(simplify( T1(m-1,A,1) + q^(m-A) *T0(m-1,A-1,1) ))))
 else 'V(m,A,k)'
 fi
end:

U:=proc(m,A,k)
   if type(m,integer) and type(A,integer) then
    RETURN(subs(q=q^k,sort(simplify( T0(m,A,1) + T0(m,A+1,1) ))))
 else 'U(m,A,k)'
 fi
end:

Trb :=proc(m,B,A,k)
   if type(m,integer) and type(A,integer) and type(B,integer) then
     RETURN(subs(q=q^k,sort(simplify(
   add( q^(j^2+j*B)*gp(m,j,1)*gp(m-j,j+A,1),j=0..m)))))
  else 'Trb(m,B,A,k)'
  fi
end:


##### User friendly front-end

encode_summand:=proc(summand) 
local a,b,c,dd,den, den_list,m, nd,num, num_list, PDF, PNF;
num:= numer(summand);
den:= denom(summand);
num_list:= isolate_factors(num);
den_list:= isolate_factors(den);
PNF:= [parse_numer_facs(num_list)];
PDF:= [parse_denom_facs(den_list)];
a :=op(1,PNF);
nd:=op(2,PNF);
dd:=op(1,PDF);
m :=op(2,PDF);
b :=op(3,PNF);
c :=op(4,PNF);
RETURN(a,nd,dd,m,b,c);
end:

encode_long_summand:=proc(summand) 
local a,b,c,d,dd,den, den_list,ed, nd,num, num_list, PDF, PNF;
num:= numer(summand);
den:= denom(summand);
num_list:= isolate_factors(num);
den_list:= isolate_factors(den);
PNF:= [parse_long_numer_facs(num_list)];
PDF:= [parse_long_denom_facs(den_list)];
a :=op(1,PNF);
nd:=op(2,PNF);
dd:=op(1,PDF);
ed:=op(2,PDF);
d :=op(3,PNF);
b :=op(4,PNF);
c :=op(5,PNF);
RETURN(a,nd,dd,ed,d,b,c);
end:


isolate_factors:=proc(f)
local ctr, fac, num_of_facs;
if type(f,`*`) then
  num_of_facs:= nops(f);
  for ctr from 1 to num_of_facs do
    fac[ctr]:=op(ctr,f)
  od
elif type(f,function) then fac[1]:=f; num_of_facs:=1
elif type(f,`^`) then fac[1]:=f; num_of_facs:=1
else ERROR(`Illegal entry`)
fi;
[seq(fac[ctr],ctr=1..num_of_facs)]
end:

parse_numer_facs:= proc(list_of_numer_factors)
local a, altsign, b, c, ctr, length, nd, qflag, qpwr, stopflag;
length:= nops(list_of_numer_factors);
### scan through for functions other than qfac
for ctr from 1 to length do
  if type(op(ctr,list_of_numer_factors),function) and
     op(0,op(ctr,list_of_numer_factors))<> qfac 
  then ERROR(`The only function type allowed in the numerator is qfac.`)
  fi
od;
altsign:=false;
### scan through for factor of (-1)^j
for ctr from 1 to length do
  if op(ctr,list_of_numer_factors) = (-1)^j then altsign:=true fi;
od;
### scan through for power of q:
ctr:=1;
stopflag:=false;
while ctr<=length and stopflag=false do
  if type(op(ctr,list_of_numer_factors),`^`) and 
     depends(op(ctr,list_of_numer_factors),q)
  then qflag:=ctr; stopflag:=true fi;
ctr:=ctr+1;
od;
### Error if more than one q^pwr is present:
while ctr<=length do
  if type(op(ctr,list_of_numer_factors),`^`) and 
     depends(op(ctr,list_of_numer_factors),q)
  then ERROR(`Only one power of q allowed.`) fi;
ctr:=ctr+1
od;
if altsign=true then a:=1 else a:=0 fi;
qpwr:= simplify(subs(q=1/2,log[q](op(qflag,list_of_numer_factors))));
#qpwr:= op(2,op(qflag,list_of_numer_factors));
if degree(qpwr,j)<>2 then ERROR(`The power of q must be quadratic in j.`) fi;
if coeff(qpwr,j,0)<>0 then ERROR(`The power of q must not contain a constant.`) fi;
b:= coeff(qpwr,j,2);
c:= coeff(qpwr,j,1);
#### Build up encoded q-factorials:
nd:=[];
for ctr from 1 to length do
  if type(op(ctr,list_of_numer_factors),function) 
  then nd:=[op(nd),encode_qfac(op(ctr,list_of_numer_factors))]
  fi
od;
RETURN(a,nd, b, c);
end:

parse_long_numer_facs:= proc(list_of_numer_factors)
local a, altsign, b, c, ctr, d, length, nd, qflag, qpwr, stopflag, tflag, tpwr;
length:= nops(list_of_numer_factors);
### scan through for functions other than qfac
for ctr from 1 to length do
  if type(op(ctr,list_of_numer_factors),function) and
     op(0,op(ctr,list_of_numer_factors))<> qfac 
  then ERROR(`The only function type allowed in the numerator is qfac.`)
  fi
od;
altsign:=false;
### scan through for factor of (-1)^j
for ctr from 1 to length do
  if op(ctr,list_of_numer_factors) = (-1)^j then altsign:=true fi;
od;
### scan through for power of q:
ctr:=1;
stopflag:=false;
while ctr<=length and stopflag=false do
  if type(op(ctr,list_of_numer_factors),`^`) and 
     depends(op(ctr,list_of_numer_factors),q)
  then qflag:=ctr; stopflag:=true fi;
ctr:=ctr+1;
od;
### Error if more than one q^pwr is present:
while ctr<=length do
  if type(op(ctr,list_of_numer_factors),`^`) and 
     depends(op(ctr,list_of_numer_factors),q)
  then ERROR(`Only one power of q allowed.`) fi;
ctr:=ctr+1
od;
### scan through for power of t:
ctr:=1;
stopflag:=false;
while ctr<=length and stopflag=false do
   if type(op(ctr,list_of_numer_factors),`^`) and 
      depends(op(ctr,list_of_numer_factors),t)
  then tflag:=ctr; stopflag:=true fi;
ctr:=ctr+1;
od;
### Error if more than one t^pwr is present:
while ctr<=length do
  if type(op(ctr,list_of_numer_factors),`^`) and 
     depends(op(ctr,list_of_numer_factors),t)
  then ERROR(`Only one power of t allowed.`) fi;
ctr:=ctr+1
od;

if altsign=true then a:=1 else a:=0 fi;
qpwr:= simplify(subs(q=1/2,log[q](op(qflag,list_of_numer_factors))));
#qpwr:= op(2,op(qflag,list_of_numer_factors));
if degree(qpwr,j)<>2 then ERROR(`The power of q must be quadratic in j.`) fi;
if coeff(qpwr,j,0)<>0 then ERROR(`The power of q must not contain a constant.`) fi;
tpwr:= expand(log[t](op(tflag,list_of_numer_factors)));
if degree(tpwr,j)<>1 then ERROR(`The power of t must be linear in j.`) fi;
if coeff(qpwr,j,0)<>0 then ERROR(`The power of t must not contain a constant.`) fi;
b:= coeff(qpwr,j,2);
c:= coeff(qpwr,j,1);
d:= coeff(tpwr,j,1);
#### Build up encoded q-factorials:
nd:=[];
for ctr from 1 to length do
  if type(op(ctr,list_of_numer_factors),function) 
  then nd:=[op(nd),encode_long_qfac(op(ctr,list_of_numer_factors))]
  fi
od;
RETURN(a,nd, d, b, c);
end:


encode_qfac:=proc(A)
local a,a1,a2,a3,a4,base,ind;
#Input : qfac(a,base,ind)
#Output: encoded version.
if nops(A)<> 3 then ERROR(`qfac requires three arguments.`) fi;
a   :=op(1,A);
base:=op(2,A);
ind :=op(3,A);
if type(a,`*`) then a1:= a/op(2,a) 
  elif type(a,integer) then a1:=a
  else a1:=1 fi;
a2:= simplify(subs(q=1/2,log[q](a/a1)));
a3:= simplify(subs(q=1/2,log[q](base)));
if degree(ind,j)<> 1 then ERROR(`Index must be of the form j + integer`) fi;
if coeff(ind,j,1)<>1 then ERROR(`Index must be of the form j + integer`) fi;
a4:= ind - j;
if not type(a4,integer) then ERROR(`Index must be of the form j + integer`) fi;
RETURN(a1,a2,a3,a4);
end:

encode_long_qfac:=proc(A)
local a,a1,a2,a3,a4,a5,a6,base,ctr,ind,qflag,tflag;
#Input : qfac(a,base,ind)
#Output: encoded version.
if nops(A)<> 3 then ERROR(`qfac requires three arguments.`) fi;
a   :=op(1,A);
base:=op(2,A);
ind :=op(3,A);
if op(1,a)=-1 then a1:=-1 else a1:=1 fi;
for ctr from 1 to nops(a) do
  if depends(op(ctr,a),t) then tflag:=ctr fi;
  if depends(op(ctr,a),q) then qflag:=ctr fi;
od;
if depends(a,t) and depends(a,q) then
   a2:= simplify(subs(t=1/2,log[t](a/a1/op(qflag,a))))
   elif depends(a,t) then
   a2:= simplify(subs(t=1/2,log[t](a/a1)))
   else a2:=0 fi;
if depends(a,q) and depends(a,t) then
  a3:= simplify(subs(q=1/2,log[q](a/a1/op(tflag,a))))
  elif depends(a,q) then
  a3:= simplify(subs(q=1/2,log[q](a/a1)))
  else a3:=0 fi;
a4:= simplify(subs(q=1/2,log[q](base)));
a5:=coeff(ind,j,1);
a6:= ind - a5*j;
if not type(a5,integer) then 
   ERROR(`Index must be of the form integer*j + integer`) fi;
RETURN(a1,a2,a3,a4,a5,a6);
end:



parse_denom_facs:= proc(list_of_denom_factors)
local arg_spqfac, ctr,dd,length,m,scount,SQF;
length:= nops(list_of_denom_factors);
### scan through for functions other than qfac/spqfac
for ctr from 1 to length do
  if not type(op(ctr,list_of_denom_factors),function) then
    ERROR(`The only function types allowed in the denominator are qfac and spqfac.`)
  fi;
  if op(0,op(ctr,list_of_denom_factors))<> qfac and
     op(0,op(ctr,list_of_denom_factors))<> spqfac then 
     ERROR(`The only function types allowed in the denominator are qfac and spqfac.`)
  fi
od;
### scan thru to make sure exactly one spqfac is present:
scount:=0;
for ctr from 1 to length do
   if op(0,op(ctr,list_of_denom_factors))=spqfac then scount:=scount+1; arg_spqfac:=ctr fi
od;
if scount<>1 then ERROR(`Exactly one spqfac must be present in the denominator.`) fi;
### Check to see if spqfac is in legal form:
SQF:= op(arg_spqfac,list_of_denom_factors);
if op(3,SQF) <> j then ERROR(`The third argument of spqfac must be j.`) fi;
if op(1,SQF) <> op(2,SQF) then 
    ERROR(`The first and second arguments of spqfac must be identical.`) fi;
if not type(op(1,SQF),`^`) and not op(1,SQF)=q 
then ERROR(`The first and second arguments of spqfac must be a power of q.`) fi;
if not depends(op(1,SQF),q) then 
   ERROR(`The first and second arguents of sqpfac must be a power of q.`) fi;
#### Build up encoded q-factorials:
dd:=[];
for ctr from 1 to arg_spqfac-1 do
  dd:=[op(dd),encode_qfac(op(ctr,list_of_denom_factors))]
od;
for ctr from arg_spqfac+1 to length do
  dd:=[op(dd),encode_qfac(op(ctr,list_of_denom_factors))]
od;
m:=simplify(subs(q=1/2,log[q](op(1,SQF))));
#if op(1,SQF) = q then m:=1 else m:=op(2,op(1,SQF)) fi;
RETURN(dd,m);
end:

parse_long_denom_facs:= proc(list_of_denom_factors)
local arg_spqfac, ctr,dd,ed,length,m,scount,SQF;
length:= nops(list_of_denom_factors);
### scan through for functions other than qfac/spqfac
for ctr from 1 to length do
  if not type(op(ctr,list_of_denom_factors),function) then
    ERROR(`The only function types allowed in the denominator are qfac and spqfac.`)
  fi;
  if op(0,op(ctr,list_of_denom_factors))<> qfac and
     op(0,op(ctr,list_of_denom_factors))<> spqfac then 
     ERROR(`The only function types allowed in the denominator are qfac and spqfac.`)
  fi
od;
### scan thru to make sure exactly one spqfac is present:
scount:=0;
for ctr from 1 to length do
   if op(0,op(ctr,list_of_denom_factors))=spqfac then scount:=scount+1; 
   arg_spqfac:=ctr fi
od;
if scount<>1 then ERROR(`Exactly one spqfac must be present in the denominator.`) 
fi;
SQF:= op(arg_spqfac,list_of_denom_factors);
#### Build up encoded q-factorials:
dd:=[];
for ctr from 1 to arg_spqfac-1 do
  dd:=[op(dd),encode_long_qfac(op(ctr,list_of_denom_factors))]
od;
for ctr from arg_spqfac+1 to length do
  dd:=[op(dd),encode_long_qfac(op(ctr,list_of_denom_factors))]
od;
ed:= [encode_long_qfac(SQF)];
RETURN(dd,ed);
end:


twovargen:= proc( summand,min_ind);
if min_ind=0 then
   tify (encode_summand(summand))
elif min_ind=1 then
   tify0(encode_summand(summand))
else ERROR(`Lower limit of summation must be either 0 or 1.`)
fi
end:

funcrecur:= proc( summand,min_ind, sn);
if min_ind=0 then 
  fundrec (encode_summand(summand),sn)
elif min_ind=1 then
  fundrec0(encode_summand(summand),sn)
else ERROR(`Lower limit of summation must be either 0 or 1.`)
fi
end:

polyrecur:= proc( summand,min_ind, sn);
if min_ind=0 then 
  RETURN(P[n] = polyrec (encode_summand(summand)))
elif min_ind=1 then
  RETURN(P[n] = polyrec0(encode_summand(summand)))
else ERROR(`Lower limit of summation must be either 0 or 1.`)
fi
end:
initialconds:= proc( summand,min_ind, sn);
if min_ind=0 then 
  RETURN(initconds (encode_summand(summand)))
elif min_ind=1 then
  RETURN(initconds0(encode_summand(summand)))
else ERROR(`Lower limit of summation must be either 0 or 1.`)
fi
end:

annihiloper:= proc( summand, min_ind, for_back);
if   min_ind=0 and for_back=back then 
  RETURN(annops (encode_summand(summand))[1])
elif min_ind=0 and for_back=forward then
  RETURN(annops (encode_summand(summand))[2])
elif min_ind=1 and for_back=back then
  RETURN(annops0(encode_summand(summand))[1])
elif min_ind=1 and for_back=forward then
  RETURN(annops0(encode_summand(summand))[2])
else ERROR(`Input error.`)
fi
end:

fermipoly:= proc(summand)
LHS(encode_summand(summand))
end:

fermipolylong := proc(summand)
LHS2(encode_long_summand(summand))
end:

fermipolydual:= proc(summand,deg)
LHSdual(encode_summand(summand),deg)
end:

fermipolyevendual:=proc(summand,deg)
LHSevendual(encode_summand(summand),deg)
end:

fermipolyodddual:=proc(summand,deg)
LHSodddual(encode_summand(summand),deg)
end:

basic_calcs:= proc(summand, min_ind, sn)
local e;
e:= encode_summand(summand);
if   min_ind=0 then theworks (e,sn)
elif min_ind=1 then theworks0(e,sn)
else ERROR(`Lower index of summation must be 0 or 1.`)
fi
end:

polylist:= proc(summand, min_ind, length)
if   min_ind=0 then makepolylist (encode_summand(summand), length)
elif min_ind=1 then makepolylist0(encode_summand(summand), length)
fi
end:

printpolyseq:= proc(summand, min_ind, length)
if   min_ind=0 then printpolylist (encode_summand(summand), length)
elif min_ind=1 then printpolylist0(encode_summand(summand), length)
fi
end:
 

###The following tools aid in making the conjectures of the
###bosonic representations:
gpsearch:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly; 
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   elif p=threezero then a:=3: b:=0
   elif p=threeone  then a:=3: b:=1
   elif p=threetwo  then a:=3: b:=2
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-add(gp(subs(n=i,A),
        subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong GP`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(gp(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * gp(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len); 
s;
end:

T0search:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly;
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
      add(T0(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong T0`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(T0(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * T0(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

t1search:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly;
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
       add(t1(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong t1`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(t1(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * t1(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

tau0search:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly;
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
       add(tau0(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong tau0`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(tau0(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * tau0(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

t0search:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly; 
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
      add(t0(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong t0`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(t0(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * t0(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

Usearch:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly;
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
        add(U(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong U`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(U(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * U(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

Vsearch:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly;
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
        add(V(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong V`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(V(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * V(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

T1search:=proc(PP,A,B,l,len,p)
local a,b,s,i,j,k,new,poly;
if p=all then a:=1: b:=0
   elif p=even then a:=2: b:=0
   elif p=odd then a:=2: b:=1
   else ERROR(`invalid 4th parameter`)
fi;
s:=[op(b+1,PP)]:
 while nops(s)< len do 
   k:=nops(s):
   for i from 1 to k-1 do
     if (collect(op(a*i+b+1,PP)-
        add(T1(subs(n=i,A),subs(n=i,B)+j-1,l)*s[j], j=1..nops(s)),q))<>0 then
      ERROR(`Wrong T1`);
     fi;
   od;
   new:= collect(op(a*k+b+1,PP)-
     add(T1(subs(n=k,A),subs(n=k,B)+j-1,l) * s[j], j=1..nops(s)),q):
   s:=[op(s),new]:
 od:
poly := 0;
for i from 0 to nops(s)-1 do
 poly := poly + op(i+1,s) * T1(A,B+i,l)
od;
print(P[a*n+b]  = poly,`for n <`, len);
s;
end:

####The following procedures insert the `t` into the q-sum and calculate the 
#recurrence from which the conjecture is made.
firstterm := proc(a, nd, dd, m, b, c)
 local tmp,i, g;
 g:=tpgcd(nd,dd,m);
 tmp:=1;
 for i from 0 to (-1+nops(nd)/4) do
  tmp := tmp * qfact(nd[4*i+1] * t^(nd[4*i+3]) * q^nd[4*i+2], 
    q^nd[4*i+3], nd[4*i+4])
 od;
 for i from 0 to (-1+nops(dd)/4) do
  tmp:=tmp/ qfact(dd[4*i+1] * t^(dd[4*i+3]) * q^dd[4*i+2], q^dd[4*i+3],     
    dd[4*i+4])
 od;
 tmp:=tmp/(1-t^(m));
 subs(t=t^(1/g),tmp)
end:

fctr := proc(a, nd, dd, m, b, c)
 local tmp,i, g;
 g:=tpgcd(nd,dd,m);
 tmp:=(-1)^a * t^(2*b) * q^(b+c)/ (1-t^(m));
 for i from 0 to (-1+nops(nd)/4) do
  tmp:=tmp*(1-nd[4*i+1]*t^(nd[4*i+3])*q^nd[4*i+2])
 od;
 for i from 0 to (-1+nops(dd)/4) do
  tmp:=tmp/(1-dd[4*i+1]*t^(dd[4*i+3])*q^dd[4*i+2])
 od;
 subs(t=t^(1/g),tmp);
end:

rewrite:= proc(a, nd, dd, m, b, c)
local tmp,i;
tmp:=1;
for i from 0 to (-1+nops(nd)/4) do
 tmp:=tmp*``(nd[4*i+1]*q^nd[4*i+2],q^nd[4*i+3])[j+nd[4*i+4]]
 od;
for i from 0 to (-1+nops(dd)/4) do   
   tmp:=tmp/``(dd[4*i+1]*q^dd[4*i+2],q^dd[4*i+3])[j+dd[4*i+4]]
 od;
tmp:=tmp*(-1)^(a*j)*q^(b*j^2+c*j)/``(q^m,q^m)[j];
Sum( tmp, j=0..infinity);
end:

tify:= proc(a, nd, dd, m, b, c)
local tmp,i, g;
tmp:=1;
g:=tpgcd(nd,dd,m);
for i from 0 to (-1+nops(nd)/4) do
tmp:=tmp*``(nd[4*i+1]*t^(nd[4*i+3]/g)*q^nd[4*i+2],q^nd[4*i+3])[j+nd[4*i+4]]
 od;
for i from 0 to (-1+nops(dd)/4) do   
  tmp:=tmp/``(dd[4*i+1]*t^(dd[4*i+3]/g)*q^dd[4*i+2],q^dd[4*i+3])[j+dd[4*i+4]]
 od;
tmp:=tmp*(-1)^(a*j)*t^(2*b*j/g)*q^(b*j^2+c*j)/``(t^(m/g),q^m)[j+1];
Sum( tmp, j=0..infinity);
end:

tpgcd := proc(nd,dd,m)
local i,ne,de,exs;
# Given the numerator data(nd), denominator data(dd) and base q-factorial
# (m), returns the GCD of the powers of t.
# This is designed to help eliminate P_odd = 0 syndrome.
ne := {seq( nd[4*i+3], i=0..-1+nops(nd)/4)};de := {seq( dd[4*i+3], i=0..-1+nops(dd)/4)};
exs := ne union de union {m};
igcd(op(exs))
end:

fundrec2:=proc(a, nd, dd, m, b, c,sn)
  local L,ft1,ft2,fctr1,fctr2,g;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm(a,nd,dd,m,b,c);
  ft2:=fctr(a,nd,dd,m,b,c);
  print(f[sn](t) = ft1 + ``(ft2) * subs(q=q^g,f[sn](t*q)));
end:

fundrec3:=proc(a, nd, dd, m, b, c, sn)
  local L,ft1,ft2,fctr1,fctr2,g;
  ft1:=firstterm(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  g:=tpgcd(nd,dd,m);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1);
  fctr2:=simplify(L*ft2);
  print(``(L) * f[sn](t) = ``(fctr1) + ``(fctr2) * subs(q=q^g,f[sn](t*q)));
end:

fundrec:=proc(a,nd,dd,m,b,c,sn)  
  fundrec2(a,nd,dd,m,b,c,sn);
  fundrec3(a,nd,dd,m,b,c,sn);
end:

seppoly:=proc(p,maxt)
#Input: p, a polynomial in t; maxt: a positive integer.
#Output: an expression sequence of length maxt+1, 
#   where the (i+1)st entry is the coefficient of t^i
local i;
seq(coeff(p,t,i),i=0..maxt)
end:  

polyrec:=proc(a, nd, dd, m, b,c)
  local L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,i;
  g := tpgcd(nd,dd,m);
  ft1:=firstterm(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  # solve recurrance found in fundrec3 for f(t)
  fctr1:=simplify(L*ft1); 
  fctr1:=collect(expand(fctr1),t);
    #fctr1 is the terms not involving f(t) or f(tq)
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
    #fctr2 is the coeff of f(tq)
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t); 
    #fctr3 is the coeff of f(t) on RHS after 1*f(t) is isolated on LHS
maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));  
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  add(P[n-i]*(s3[i+1]+combine(s2[i+1]*q^(g*(n-i)),power)),i=1..maxt);   
end:

Ipolyrec:=proc(a, nd, dd, m, b,c)  
local L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1);
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  seq((s3[i+1]+combine(s2[i+1]*q^(g*(n-i)),power)),i=1..maxt);  
end:

initconds:=proc(a, nd, dd, m, b, c)
local D,P,L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,j,i;
  g:=tpgcd(nd,dd,m); 
  ft1:=firstterm(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1);
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  for i from -maxt to -1 do D[i]:=0 od;
  D[0]:= s1[1]:
  print(P[0] = s1[1]);
  for j from 1 to maxt-1 do
  D[j]:=s1[j+1] + add(s3[i+1]*D[j-i],i=1..maxt) + 
      add(s2[i+1]*D[j-i]*q^(g*(j-i)),i=1..maxt);
    print(P[j] = sort(simplify(D[j])))
  od;
  end:

Iinitconds:=proc(a, nd, dd, m, b, c)
local D,P,L,ft1,ft2,fctr1,fctr2,fctr3,g,s1,s2,s3,maxt,j,i;
  g:=tpgcd(nd,dd,m); 
  ft1:=firstterm(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1);
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  for i from -maxt to -1 do D[i]:=0 od;
  D[0]:= s1[1]:
  P:=[s1[1]]:
  for j from 1 to maxt-1 do
    D[j]:=s1[j+1] + add(s3[i+1]*D[j-i],i=1..maxt) +  
      add(s2[i+1]*D[j-i]*q^(g*(j-i)),i=1..maxt);
    P:=[op(P),sort(simplify(D[j]))]
  od;
  P;
  end:

theworks:= proc(a,nd,dd,m,b,c,sn)  
  print(phi[sn](q) = rewrite(a,nd,dd,m,b,c));
  print(f[sn](t) = tify(a,nd,dd,m,b,c));
  print(fundrec(a,nd,dd,m,b,c,sn));
  print(P[n] = polyrec(a,nd,dd,m,b,c));
  print(initconds(a,nd,dd,m,b,c));
  print(annops(a,nd,dd,m,b,c)[1]);print(annops(a,nd,dd,m,b,c)[2]);
end:


###The following procedures make it easy to view any number of polynomials 
###in the sequence generated by the recurrance
makepolylist:=proc(a,nd,dd,m,b,c, len)
  local j,i,Q,new,k,R;
  Q:=Iinitconds(a,nd,dd,m,b,c);     
  R:=[Ipolyrec(a,nd,dd,m,b,c)];   
  k:=nops(R);
  for j from nops(Q) to len-1 do   
    new:=sort(collect(expand(add(subs(n=j, op(j-i+1,Q)*op(i,R)),i=1..k)),q));
    Q:=[op(Q),new]    
  od;
  Q;
end:

evenlist := proc(Q)
local j,R;  R:=[]; for j from 1 to nops(Q) by 2 do R:= [op(R), op(j,Q)] od;
R
end:

oddlist:=proc(Q) local j,R; R:=[];
for j from 2 to nops(Q) by 2 do R:=[op(R), op(j,Q)] od;
R
end:

printpolylist:=proc(a,nd,dd,m,b,c,len)
  local S,i;
  S:=makepolylist(a,nd,dd,m,b,c,len);
  for i from 0 to len-1 do print(P[i] = op(i+1,S)) od
 end:

###q-Product procedures

makejtp:=proc(q1,t1) local ns, qpwr,tpwr, t;
 qpwr:= simplify(subs(q=1/2,log[q](q1)));
 t:= t1;
 if op(1,t)=-1 then t:=-t; ns:=-1 else ns:=1 fi;
 tpwr:= simplify(subs(q=1/2,log[q](t)));
   Sum( (-1*ns)^j * q^(qpwr*j^2 + tpwr*j), j=-infinity..infinity) 
    = ``(t1*q1,q1^2)[infinity] *``(q1/t1,q1^2)[infinity]*``(q1^2,q1^2)[infinity]
 end:

makeqp := proc(q,z);
   print( ``(z^3*q,q^3)[infinity]*``(q^2/z^3,q^3)[infinity]*
    ``(q^3,q^3)[infinity] + 
    z*``(q/z^3,q^3)[infinity] * ``(z^3*q^2,q^3)[infinity]*``(q^3,q^3)[infinity]
    = ``(-q/z,q)[infinity]* ``(-z,q)[infinity] * ``(q/z^2,q^2)[infinity] * 
      ``(z^2*q,q^2)[infinity] * ``(q,q)[infinity])
end:

maketp := proc(q,z)
   print(``
 (-q*z^2,q^4)[infinity] *``(-q^3/z^2,q^4)[infinity] * ``(q^4,q^4)[infinity] +
 z*``(-z^2*q^3,q^4)[infinity] * ``(-q/z^2,q^4)[infinity]* ``(q^4,q^4)[infinity] =
 ``(-z,q)[infinity] * ``(-q/z,q)[infinity] * ``(q,q)[infinity])
end:


####Reciprocal Polynomials
parse_bosonic_facs:= proc(list_of_factors)
local altsign, ctr, func, length, qflag, qpwr, s, stopflag;
length:= nops(list_of_factors);
### scan through for functions
for ctr from 1 to length do
  if type(op(ctr,list_of_factors),function) 
    then func:=op(ctr,list_of_factors)
  fi
od;
altsign:=false;
### scan through for factor of (-1)^j
for ctr from 1 to length do
  if op(ctr,list_of_factors) = (-1)^j then altsign:=true fi;
od;
### scan through for power of q:
ctr:=1;
stopflag:=false;
while ctr<=length and stopflag=false do
  if type(op(ctr,list_of_factors),`^`) and 
     depends(op(ctr,list_of_factors),q)
  then qflag:=ctr; stopflag:=true fi;
ctr:=ctr+1;
od;
### Error if more than one q^pwr is present:
while ctr<=length do
  if type(op(ctr,list_of_factors),`^`) and 
     depends(op(ctr,list_of_factors),q)
  then ERROR(`Only one power of q allowed.`) fi;
ctr:=ctr+1
od;
if altsign=true then s:=j else s:=0 fi;
qpwr:= simplify(subs(q=1/2,log[q](op(qflag,list_of_factors))));
#qpwr:= op(2,op(qflag,list_of_factors));
if degree(qpwr,j)<>2 then ERROR(`The power of q must be quadratic in j.`) fi;
if coeff(qpwr,j,0)<>0 then ERROR(`The power of q must not contain a constant.`) fi;
RETURN(s,qpwr,func);
end:

bosedual :=proc(summand, deg)
local list_of_factors;
list_of_factors:= isolate_factors(summand);
recippoly(parse_bosonic_facs(list_of_factors),deg)
end:

recippoly := proc( s, qexp, func, deg)
  local A,m,k,f,fr,qcoefexp,tmp;
  if op(0,func)=CT 
  then 
    f:=T1; 
    m:=op(1,func); 
    A:=op(2,func); 
    k=1/2 
   else
    f :=op(0,func);
    m :=op(1,func);
    A :=op(2,func);
    k :=op(3,func)
  fi;
  if f=T0 then fr:=t0; qcoefexp:= (A^2-m^2)*k 
   elif f=T1 then fr:=t1; qcoefexp:= (A^2-m^2)*k
   elif f=gp then fr:=gp; qcoefexp:= -A*(m-A)*k  
  else ERROR(`Only gp, T0, and T1 accepted`)  
  fi;
  #tmp:=(-1)^s * q^(deg - qexp) * f(m,A,-k) ;
  tmp:=Sum((-1)^s * q^(expand(deg - qexp + qcoefexp)) * 
    fr(m,A,k),j=-infinity..infinity);  
end:

duallist := proc(P)
local i,recip,Q;
Q:=[];
for i from 1 to nops(P) do
  recip:=sort(expand(q^(degree(op(i,P),q))*subs(q=1/q,op(i,P))));
  Q:=[op(Q),recip]
 od;
 Q
end:
#####################

septerms := proc(L)
local i,j,M;
M:=[];
for i from 1 to nops(L) 
 do 
  if not type(op(i,L),`+`) then
   if op(i,L) <> 0 then M:=[op(M), op(i,L)] fi
  else 
   for j from 1 to nops(op(i,L)) 
      do
       M:=[op(M), op(j, op(i,L))]
      od;
   fi;
od;
M
end:

getsigns := proc(L)
# input list should have no zeros or multiterm entries
  local i,S;
  S :=[];
for i from 1 to nops(L) do
  if   op(i,L) = 1 then S:=[op(S),1]
   elif op(1,op(i,L)) = q then S:=[op(S),1]
   elif op(1,op(i,L)) = -1 then S:=[op(S),-1]
   else ERROR(`Sorry, I can't handle this list.`)
  fi
od;
S;
end:

getexps := proc(L)
# Input list should have no zeros or multiterms sums
local NZ,S,EL,i;
 NZ:=septerms(L);
 S:=getsigns(NZ);
 EL :=[];
 for i from 1 to nops(NZ) do
     EL:=[op(EL),simplify(subs(q=1/2, log[q](op(i,S) * op(i,L))))]
    #EL:=[op(EL),     expand  (       log[q](op(i,S) * op(i,L))]
 od;
sort(EL);
end:

conjexp := proc(Q)
# Input is a list of exponents (numbers)
local EL,expguess,i,r,soln;
#EL:=get;
EL:=Q;
if nops(EL)<3 then ERROR(`Not enough terms in input`) 
  else
   soln:=solve( { a + b + c = op(3,EL), a - b + c = op(2,EL), c= op(1,EL)},{a,b,c})
fi;
expguess := subs(soln, a*j^2 + b*j + c);
if nops(EL) =3 then RETURN(expguess,`Not enough terms to check conjecture.`)
else 
  for i from 4 to nops(EL)
    do
      r:= (-1)^(i+1) * floor(i/2);
      if subs(j=r, expguess)<> op(i,EL) then RETURN(`FAIL`) fi;
    od
fi;
RETURN(expguess,`Checked thru`,r)
end:

conjexps:=proc(Q)
local R;
R := getexps(septerms(Q));
conjexp(R);
end:

conjexpeven:= proc(Q)
local EL,R;
EL := getexps(septerms(Q));
R:=pickoffevenk(EL);
conjexp(R);
end:

conjexpodd:= proc(Q)
local EL,R;
EL := getexps(septerms(Q));
R:=pickoffoddk(EL);
conjexp(R);
end:

pickoffevenk := proc(EL)  
#Takes a list of exponents and picks off j=0,+/- 2, +/- 4,...
local i,R;
if nops(EL)=0 then ERROR(`Empty input`)
else 
  R:=[op(1,EL)];
  i:=1;
  while i<= nops(EL)
    do
      i:=i+3; if i>nops(EL) then RETURN(R) fi;
      R:=[op(R),op(i,EL)];
      i:=i+1;if i>nops(EL) then RETURN(R) fi;
      R:=[op(R),op(i,EL)];if i>nops(EL) then RETURN(R) fi;
    od;
fi;
R
end:   

pickoffoddk := proc(EL)  
#Takes a list of exponents and picks off j=+/- 1, +/- 3,...
local i,R;
if nops(EL)< 3 then ERROR(`Insufficient input`)
else 
  R:=[op(3,EL)];
  i:=3;
  while i<= nops(EL)
    do
      i:=i-1;
      R:=[op(R),op(i,EL)];
      i:=i+5;if i>nops(EL) then RETURN(R) fi;
      R:=[op(R),op(i,EL)];if i>nops(EL) then RETURN(R) fi;
    od;
fi;
R
end:   

search := proc(P,fn,A,B,l,len,p)
if   fn=gp   then   gpsearch(P,A,B,l,len,p)
elif fn=T0   then   T0search(P,A,B,l,len,p)
elif fn=T1   then   T1search(P,A,B,l,len,p)
elif fn=tau0 then tau0search(P,A,B,l,len,p)
elif fn=U    then    Usearch(P,A,B,l,len,p)
elif fn=V    then    Vsearch(P,A,B,l,len,p)
else RETURN(`Illegal function`)
fi
end:

searchconj := proc(fn,P,A,B,l,len,p)
local pref,prefodd,prefeven,Q1,S;
Q1:=search(fn,Q,2*n+1,n+1,1,10,all);
if conjexp(Q1)<>FAIL then
  S:=getsigns(septerms(Q1)):
  #pref:=1:
  if   (op(2,S)=-1 and op(3,S)=-1) then pref:=(-1)^j 
   elif (op(2,S)=1  and op(3,S)= 1) then pref:=1
   else RETURN(`FAIL`)
  fi;
  Sum(pref*conjexp(Q1)*fn(A,B,l),j=-infinity..infinity)
elif (conjexpeven(Q1)<>FAIL and conjexpodd(Q1)<>FAIL) then
  S:=getsigns(septerms(Q1)):
  if (op(2,S)=-1 and op(3,S)=-1 and op(4,S)=1 and op(5,S)=1) then 
     prefodd:=-1:
     prefeven:=1
  else RETURN(`FAIL`)
  fi;
  Sum(prefeven*q^(conjexpeven(Q1))*fn(A,B,l) + 
    prefodd*q^(conjexpodd(Q1))*fn(A,B,l), 
    j=-infinity..infinity)
else RETURN(`FAIL`)
fi
end:

rewrite1 := proc(a,nd,dd,m,b,c)
local tmp,i;
tmp:=1;
for i from 0 to (-1 + nops(nd)/4) do
  tmp:=tmp*``(nd[4*i+1]*q^nd[4*i+2],q^nd[4*i+3])[j+nd[4*i+4]]
od;
for i from 0 to (-1 + nops(dd)/4) do
  tmp:=tmp/``(dd[4*i+1]*q^dd[4*i+2],q^dd[4*i+3])[j+dd[4*i+4]]
od;
tmp := tmp*(-1)^(a*j)*q^(b*j^2+c*j)/``(q^m,q^m)[j];
1+Sum(tmp, j=1..infinity);
end:

tify1:= proc(a, nd, dd, m, b, c)
local tmp,i, g;
tmp:=1;
g:=tpgcd(nd,dd,m);
for i from 0 to (-1+nops(nd)/4) do
tmp:=tmp*``(nd[4*i+1]*t^(nd[4*i+3]/g)*q^nd[4*i+2],q^nd[4*i+3])[j+nd[4*i+4]]
 od;
for i from 0 to (-1+nops(dd)/4) do
  tmp:=tmp/``(dd[4*i+1]*t^(dd[4*i+3]/g)*q^dd[4*i+2],q^dd[4*i+3])[j+dd[4*i+4]]
 od;
tmp:=tmp*(-1)^(a*j)*t^(2*b*j/g)*q^(b*j^2+c*j)/``(t^(m/g),q^m)[j+1];
1/(1-t^g) + Sum( tmp, j=1..infinity);
end:

tify0:= proc(a, nd, dd, m, b, c)
local tmp,i, g;
tmp:=1;
g:=tpgcd(nd,dd,m);
for i from 0 to (-1+nops(nd)/4) do
tmp:=tmp*``(nd[4*i+1]*t^(nd[4*i+3]/g)*q^nd[4*i+2],q^nd[4*i+3])[j+nd[4*i+4]]
 od;
for i from 0 to (-1+nops(dd)/4) do
  tmp:=tmp/``(dd[4*i+1]*t^(dd[4*i+3]/g)*q^dd[4*i+2],q^dd[4*i+3])[j+dd[4*i+4]]
 od;
tmp:=tmp*(-1)^(a*j)*t^(2*b*j/g)*q^(b*j^2+c*j)/``(t^(m/g),q^m)[j+1];
1/(1-t) + Sum( tmp, j=1..infinity);
end:

firstterm1 := proc(a, nd, dd, m, b, c)
 local tmp,i, g;
 g:=tpgcd(nd,dd,m);
 tmp:=(-1)^a;
 for i from 0 to (-1+nops(nd)/4) do
  tmp := tmp * qfact(nd[4*i+1] * t^(nd[4*i+3]) * q^nd[4*i+2],
    q^nd[4*i+3], nd[4*i+4]+1)
 od;
 for i from 0 to (-1+nops(dd)/4) do
  tmp:=tmp/ qfact(dd[4*i+1] * t^(dd[4*i+3]) * q^dd[4*i+2], q^dd[4*i+3],
    dd[4*i+4]+1)
 od;
 tmp:=tmp*t^(2*b)*q^(b+c)/(1-t^m)/(1-t^m*q^m);
 subs(t=t^(1/g),tmp)
end:

firstterm0:=proc(a,nd,dd,m,b,c) firstterm1(a,nd,dd,m,b,c) end:

fundrec1:=proc(a,nd,dd,m,b,c,sn)
  fundrec12(a,nd,dd,m,b,c,sn);
  fundrec13(a,nd,dd,m,b,c,sn);
end:

fundrec0:=proc(a,nd,dd,m,b,c,sn)
  fundrec02(a,nd,dd,m,b,c,sn);
  fundrec03(a,nd,dd,m,b,c,sn);
end:

fundrec12:=proc(a, nd, dd, m, b, c,sn)
  local L,ft1,ft2,fctr1,fctr2,g;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm1(a,nd,dd,m,b,c);
  ft2:=fctr(a,nd,dd,m,b,c);
  print(f[sn](t) = 1/(1-t^g) + ft1 + ``(ft2) * subs(q=q^g,f[sn](t*q) - 1/(1-(t*q)^g)));
end:

fundrec02:=proc(a, nd, dd, m, b, c,sn)
  local L,ft1,ft2,fctr1,fctr2,g;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm0(a,nd,dd,m,b,c);
  ft2:=fctr(a,nd,dd,m,b,c);
  print(f[sn](t) = 1/(1-t) 
     + ft1 + ``(ft2) * subs(q=q^g,f[sn](t*q) - 1/(1-t*q) ));
end:

fundrec13:=proc(a, nd, dd, m, b, c, sn)
  local L,ft1,ft2,fctr1,fctr1A,fctr1B,fctr2,fctr5,g;
  ft1:=firstterm1(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  g:=tpgcd(nd,dd,m);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1);
  fctr2:=simplify(L*ft2);
  fctr1A:=simplify(L/(1-t^g));
  fctr1B:=simplify(L*ft2/subs(q=q^g,(1-(t*q)^g)));
  fctr5:=simplify(fctr1A + fctr1 -fctr1B);
  print(``(L) * f[sn](t) = ``(fctr5) + ``(fctr2) * subs(q=q^g,f[sn](t*q)));
end:

fundrec03:=proc(a, nd, dd, m, b, c, sn)
  local L,ft1,ft2,fctr1,fctr1A,fctr1B,fctr2,fctr5,g;
  ft1:=firstterm0(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  g:=tpgcd(nd,dd,m);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1);
  fctr2:=simplify(L*ft2);
  fctr1A:=simplify(L/(1-t));
  fctr1B:=simplify(L*ft2/subs(q=q^g,(1-t*q)));
  fctr5:=simplify(fctr1A + fctr1 -fctr1B);
  print(``(L) * f[sn](t) = ``(fctr5) + ``(fctr2) * subs(q=q^g,f[sn](t*q)));
end:

polyrec1:=proc(a, nd, dd, m, b,c)
  local L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,i;
  g := tpgcd(nd,dd,m);
  ft1:=firstterm1(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  # solve recurrence found in fundrec3 for f(t)
  fctr1:=simplify(L*ft1 + L/(1-t^g) - L*ft2/subs(q=q^g,(1-(t*q)^g)) );
  fctr1:=collect(expand(fctr1),t);
    #fctr1 is the terms not involving f(t) or f(tq)
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
    #fctr2 is the coeff of f(tq)
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
    #fctr3 is the coeff of f(t) on RHS after 1*f(t) is isolated on LHS
maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  add(P[n-i]*(s3[i+1]+combine(s2[i+1]*q^(g*(n-i)),power)),i=1..maxt);
end:

polyrec0:=proc(a, nd, dd, m, b,c)
  local L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,i;
  g := tpgcd(nd,dd,m);
  ft1:=firstterm0(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  # solve recurrence found in fundrec3 for f(t)
  fctr1:=simplify(L*ft1 + L/(1-t) - L*ft2/subs(q=q^g,(1-t*q) ) );
  fctr1:=collect(expand(fctr1),t);
    #fctr1 is the terms not involving f(t) or f(tq)
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
    #fctr2 is the coeff of f(tq)
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
    #fctr3 is the coeff of f(t) on RHS after 1*f(t) is isolated on LHS
maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  add(P[n-i]*(s3[i+1]+combine(s2[i+1]*q^(g*(n-i)),power)),i=1..maxt);
end:

initconds1:=proc(a, nd, dd, m, b, c)
local D,P,L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,j,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm1(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1 + L/(1-t^g) - L*ft2/subs(q=q^g,(1-(t*q)^g)) );
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  for i from -maxt to -1 do D[i]:=0 od;
  D[0]:= s1[1]:
  print(P[0] = s1[1]);
  for j from 1 to maxt-1 do
  D[j]:=s1[j+1] + add(s3[i+1]*D[j-i],i=1..maxt) +
      add(s2[i+1]*D[j-i]*q^(g*(j-i)),i=1..maxt);
    print(P[j] = sort(simplify(D[j])))
  od;
  end:

initconds0:=proc(a, nd, dd, m, b, c)
local D,P,L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,j,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm0(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1 + L/(1-t) - L*ft2/subs(q=q^g,(1-t*q)) );
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  for i from -maxt to -1 do D[i]:=0 od;
  D[0]:= s1[1]:
  print(P[0] = s1[1]);
  for j from 1 to maxt-1 do
  D[j]:=s1[j+1] + add(s3[i+1]*D[j-i],i=1..maxt) +
      add(s2[i+1]*D[j-i]*q^(g*(j-i)),i=1..maxt);
    print(P[j] = sort(simplify(D[j])))
  od;
  end:

Iinitconds1:=proc(a, nd, dd, m, b, c)
local D,P,L,ft1,ft2,fctr1,fctr2,fctr3,g,s1,s2,s3,maxt,j,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm1(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1 + L/(1-t^g) - L*ft2/subs(q=q^g,(1-(t*q)^g)) );
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  for i from -maxt to -1 do D[i]:=0 od;
  D[0]:= s1[1]:
  P:=[s1[1]]:
  for j from 1 to maxt-1 do
    D[j]:=s1[j+1] + add(s3[i+1]*D[j-i],i=1..maxt) +
      add(s2[i+1]*D[j-i]*q^(g*(j-i)),i=1..maxt);
    P:=[op(P),sort(simplify(D[j]))]
  od;
  P;
  end:

Iinitconds0:=proc(a, nd, dd, m, b, c)
local D,P,L,ft1,ft2,fctr1,fctr2,fctr3,g,s1,s2,s3,maxt,j,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm1(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1 + L/(1-t) - L*ft2/subs(q=q^g,(1-t*q)) );
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  for i from -maxt to -1 do D[i]:=0 od;
  D[0]:= s1[1]:
  P:=[s1[1]]:
  for j from 1 to maxt-1 do
    D[j]:=s1[j+1] + add(s3[i+1]*D[j-i],i=1..maxt) +
      add(s2[i+1]*D[j-i]*q^(g*(j-i)),i=1..maxt);
    P:=[op(P),sort(simplify(D[j]))]
  od;
  P;
  end:

makepolylist1:=proc(a,nd,dd,m,b,c, len)
  local j,i,Q,new,k,R;
  Q:=Iinitconds1(a,nd,dd,m,b,c);
  R:=[Ipolyrec1(a,nd,dd,m,b,c)];
  k:=nops(R);
  for j from nops(Q) to len-1 do
    new:=sort(collect(expand(add(subs(n=j, op(j-i+1,Q)*op(i,R)),i=1..k)),q));
    Q:=[op(Q),new]
  od;
  Q;
end:

makepolylist0:=proc(a,nd,dd,m,b,c, len)
  local j,i,Q,new,k,R;
  Q:=Iinitconds0(a,nd,dd,m,b,c);
  R:=[Ipolyrec0(a,nd,dd,m,b,c)];
  k:=nops(R);
  for j from nops(Q) to len-1 do
    new:=sort(collect(expand(add(subs(n=j, op(j-i+1,Q)*op(i,R)),i=1..k)),q));
    Q:=[op(Q),new]
  od;
  Q;
end:

Ipolyrec1:=proc(a, nd, dd, m, b,c)
local L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm1(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1 + L/(1-t^g) - L*ft2/subs(q=q^g,(1-(t*q)^g)) );
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  seq((s3[i+1]+combine(s2[i+1]*q^(g*(n-i)),power)),i=1..maxt);
end:

Ipolyrec0:=proc(a, nd, dd, m, b,c)
local L,ft1,ft2,fctr1,fctr2,fctr3,g, s1,s2,s3,maxt,i;
  g:=tpgcd(nd,dd,m);
  ft1:=firstterm0(a, nd, dd, m, b,c);
  ft2:=fctr(a,nd,dd, m, b,c);
  L:=lcm(denom(ft1),denom(ft2));
  if coeff(L,t,0)=-1 then L:=-L fi;
  fctr1:=simplify(L*ft1 + L/(1-t) - L*ft2/subs(q=q^g,(1-t*q)) );
  fctr1:=collect(expand(fctr1),t);
  fctr2:=simplify(L*ft2);
  fctr2:=collect(expand(fctr2),t);
  fctr3:=1+simplify(-L);
  fctr3:=collect(expand(fctr3),t);
  maxt:=max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
  s1:=seppoly(fctr1,maxt);
  s2:=seppoly(fctr2,maxt);
  s3:=seppoly(fctr3,maxt);
  seq((s3[i+1]+combine(s2[i+1]*q^(g*(n-i)),power)),i=1..maxt);
end:

theworks1:= proc(a,nd,dd,m,b,c,sn)
  print(phi[sn](q) = rewrite1(a,nd,dd,m,b,c));
  print(f[sn](t) = tify1(a,nd,dd,m,b,c));
  print(fundrec1(a,nd,dd,m,b,c,sn));
  print(P[n] = polyrec1(a,nd,dd,m,b,c));
  print(initconds1(a,nd,dd,m,b,c));
  print(annops1(a,nd,dd,m,b,c)[1]);
  print(annops1(a,nd,dd,m,b,c)[2]);  
end:

theworks0:= proc(a,nd,dd,m,b,c,sn)
  print(phi[sn](q) = rewrite1(a,nd,dd,m,b,c));
  print(f[sn](t) = tify0(a,nd,dd,m,b,c));
  print(fundrec0(a,nd,dd,m,b,c,sn));
  print(P[n] = polyrec0(a,nd,dd,m,b,c));
  print(initconds0(a,nd,dd,m,b,c));
  print(annops0(a,nd,dd,m,b,c)[1]);
  print(annops0(a,nd,dd,m,b,c)[2]);
end:

printpolylist1:=proc(a,nd,dd,m,b,c,len)
  local S,i;
  S:=makepolylist1(a,nd,dd,m,b,c,len);
  for i from 0 to len-1 do print(P[i] = op(i+1,S)) od
 end:

printpolylist0:=proc(a,nd,dd,m,b,c,len)
  local S,i;
  S:=makepolylist0(a,nd,dd,m,b,c,len);
  for i from 0 to len-1 do print(P[i] = op(i+1,S)) od
 end:

###Automated expansion of f(t) by q-binomial theorem:

LHS:=proc(a,nd,dd,m,b,c)
local asp,newasp, ctr, elimvar, foldsum, g,gpp, nnf, ndf,newgpp, newqpwr, varlist, varset, qpwr, tpwr;
nnf :=nops(nd)/4;
ndf :=1 + nops(dd)/4;
foldsum := nnf + ndf;
g:=tpgcd(nd,dd,m);
varset := {j};
varlist :=[j];
# Build up the power on t:
tpwr := (2*b/g)*j;
for ctr from 0 to (-1+nops(nd)/4) do
  tpwr := tpwr + (nd[4*ctr+3]/g) * i[ctr+1];
  varset := varset union {i[ctr+1]};
  varlist :=[op(varlist),i[ctr+1]]
od;
varset := varset union{h};
varlist := [op(varlist),h];
tpwr:=tpwr + (m/g)*h;
for ctr from 0 to (-1 + nops(dd)/4) do
  tpwr := tpwr + (dd[4*ctr+3]/g) * k[ctr+1];
  varset := varset union {k[ctr+1]};
  varlist := [op(varlist),k[ctr+1]]
od; 
# Build up the power on q:
qpwr := b*j^2 + c*j;
for ctr from 0 to (-1+nops(nd)/4) do
  qpwr := qpwr 
    + (nd[4*ctr+3]) * (i[ctr+1]^2 - i[ctr+1]) /2 
    + (nd[4*ctr+2]) * i[ctr+1];
od;
for ctr from 0 to (-1+nops(dd)/4) do
  qpwr := qpwr + (dd[4*ctr+2]) * k[ctr+1];
od;
# Build up the power on (-1):
asp:= (-1)^(a*j);
for ctr from 0 to (-1+nops(nd)/4) do
  asp := asp * ((-1)*(nd[4*ctr+1]))^i[ctr+1]
od;
for ctr from 0 to (-1+nops(dd)/4) do
  asp := asp * dd[4*ctr+1]^k[ctr+1]
od;
asp := combine(asp,power);
# Build up Gaussian polynomials:
gpp :=1;
for ctr from 0 to (-1+nops(nd)/4) do
  gpp := gpp * gp(j + nd[4*ctr + 4], i[ctr+1], nd[4*ctr + 3])
od;
for ctr from 0 to (-1+nops(dd)/4) do
  gpp := gpp * gp(j+ dd[4*ctr + 4] + k[ctr+1]-1, k[ctr+1], dd[4*ctr+3])
od;
gpp := gpp * gp(j+h,j,m);
# Next we choose which of the variables to be eliminated:
ctr := 0;
elimvar := FAIL;
#########
if coeff(tpwr, h) = 1 then elimvar :=h fi; 
while (elimvar = FAIL and ctr <= nnf) do
  if ((coeff(qpwr, i[ctr]^2) = 0) and (coeff(qpwr,i[ctr]) = 0) and (coeff(tpwr,i[ctr]) = 1)) then elimvar := i[ctr] fi;
   ctr := ctr + 1;
od;
ctr := 0;
while (elimvar = FAIL and ctr <= ndf) do
  if((coeff(qpwr, k[ctr]^2) = 0) and (coeff(qpwr,k[ctr]) = 0) and (coeff(tpwr,k[ctr]) = 1)) then elimvar := k[ctr] fi;
   ctr := ctr + 1; 
od;
ctr := 0;
while (elimvar = FAIL and ctr <= nnf) do
  if coeff(tpwr,i[ctr]) = 1 then elimvar := i[ctr] fi;
  ctr := ctr + 1; 
od;
ctr := 0;
while (elimvar = FAIL and ctr <= ndf) do
  if coeff(tpwr,k[ctr]) = 1 then elimvar := k[ctr] fi:
  ctr := ctr + 1;
od;
if (nnf=0 and ndf=1 and coeff(tpwr,h)=1) then elimvar := h fi;
print(f(t)=Sum(asp * t^tpwr * q^qpwr * gpp, op(varset)=0..infinity));
print(`Variable to be eliminated`, elimvar);
newasp :=         subs(elimvar = n-tpwr+elimvar, asp);
newqpwr := expand(subs(elimvar = n-tpwr+elimvar, qpwr));
newgpp :=         subs(elimvar = n-tpwr+elimvar, gpp);
print(P[n]=Sum( newasp * q^newqpwr * newgpp,
op(varset minus {elimvar})=0..infinity)); 
end:

LHS2:=proc(a,nd,dd,ed,d,b,c)
local asp,newasp, ctr, elimvar, foldsum, g,gpp, nnf, ndf,newgpp, 
newqpwr, varlist, varset, qpwr, tpwr;
nnf :=nops(nd)/6;
ndf :=1 + nops(dd)/6;
foldsum := nnf + ndf;
varset := {j};
varlist :=[j];
# Build up the power on t:
tpwr := d*j;
for ctr from 0 to (-1+nops(nd)/6) do
  tpwr := tpwr + nd[6*ctr + 2] * i[ctr+1];
  varset := varset union {i[ctr+1]};
  varlist :=[op(varlist),i[ctr+1]]
od;
varset := varset union{h};
varlist := [op(varlist),h];
tpwr:=tpwr + (ed[2])*h;
for ctr from 0 to (-1 + nops(dd)/6) do
  tpwr := tpwr + dd[6*ctr + 2] * k[ctr+1];
  varset := varset union {k[ctr+1]};
  varlist := [op(varlist),k[ctr+1]]
od; 
# Build up the power on q:
qpwr := b*j^2 + c*j;
for ctr from 0 to (-1+nops(nd)/6) do
  qpwr := qpwr 
    + nd[6*ctr + 4] * (i[ctr+1]^2 - i[ctr+1]) /2 
    + nd[6*ctr + 3] *  i[ctr+1];
od;
for ctr from 0 to (-1+nops(dd)/6) do
  qpwr := qpwr + (dd[6*ctr+3]) * k[ctr+1];
od;
# Build up the power on (-1):
asp:= (-1)^(a*j);
for ctr from 0 to (-1+nops(nd)/6) do
  asp := asp * ((-1)*(nd[6*ctr+1]))^i[ctr+1]
od;
asp := asp * ed[1]^h;
for ctr from 0 to (-1+nops(dd)/6) do
  asp := asp * dd[6*ctr+1]^k[ctr+1]
od;
asp := combine(asp,power);
# Build up Gaussian polynomials:
gpp :=1;
for ctr from 0 to (-1+nops(nd)/6) do
  gpp := gpp * gp(nd[6*ctr + 5]*j + nd[6*ctr + 6], i[ctr+1], nd[6*ctr + 4])
od; 
for ctr from 0 to (-1+nops(dd)/6) do
  gpp := gpp * gp(dd[6*ctr + 5]*j + dd[6*ctr + 6] + k[ctr+1]-1, 
  k[ctr+1], dd[6*ctr+4])
od;
gpp := gpp * gp(ed[5]*j + ed[6] + h -1, ed[5]*j + ed[6]-1 , ed[4]);
elimvar :=h;
print(f(t)=Sum(asp * t^tpwr * q^qpwr * gpp, op(varset)=0..infinity));
print(`Variable to be eliminated`,elimvar);
newasp :=         subs(elimvar = n-tpwr+elimvar, asp);
newqpwr := expand(subs(elimvar = n-tpwr+elimvar, qpwr));
newgpp :=         subs(elimvar = n-tpwr+elimvar, gpp);
print(P[n]=Sum(newasp * q^newqpwr * newgpp, op(varset minus {elimvar})=0..infinity)); 
end:

duallist := proc(P)
 local i,recip,Q;
 Q:=[];
 for i from 1 to nops(P) do
   recip := sort(expand(q^(degree(op(i,P),q))*subs(q=1/q,op(i,P))));
   Q:=[op(Q),recip]
 od;
 Q
 end:

LHSdual:=proc(a,nd,dd,m,b,c, p)
 local asp,newasp, ctr, elimvar, foldsum, g,gpp, nnf, ndf,newgpp, newqpwr, varlist, varset, qpwr, tpwr;
 nnf :=nops(nd)/4;
 ndf :=1 + nops(dd)/4;
 foldsum := nnf + ndf;
 g:=tpgcd(nd,dd,m);
 varset := {j};
 varlist :=[j];
 # Build up the power on t:
 tpwr := (2*b/g)*j;
 for ctr from 0 to (-1+nops(nd)/4) do
   tpwr := tpwr + (nd[4*ctr+3]/g) * i[ctr+1];
   varset := varset union {i[ctr+1]};
   varlist :=[op(varlist),i[ctr+1]]
 od;
 varset := varset union{h};
 varlist := [op(varlist),h];
 tpwr:=tpwr + (m/g)*h;
 for ctr from 0 to (-1 + nops(dd)/4) do
   tpwr := tpwr + (dd[4*ctr+3]/g) * k[ctr+1];
   varset := varset union {k[ctr+1]};
   varlist := [op(varlist),k[ctr+1]]
 od; 
 # Build up the power on q
qpwr := p - b*j^2 - c*j -m*h*j;
 for ctr from 0 to (-1+nops(nd)/4) do
   qpwr := qpwr 
     - (nd[4*ctr+3]) * (i[ctr+1]^2 - i[ctr+1]) /2 
     - (nd[4*ctr+2]) * i[ctr+1]
     + (nd[4*ctr+3]) * i[ctr+1]*(i[ctr+1] - (j+nd[4*ctr + 4]));
 od;
 for ctr from 0 to (-1+nops(dd)/4) do
   qpwr := qpwr - (dd[4*ctr+2]) * k[ctr+1]
    + dd[4*ctr+3] * k[ctr+1] * (k[ctr+1]-(j + dd[4*ctr + 4] + k[ctr+1]-1));
 od;
 qpwr :=expand(qpwr);
 qpwr :=factor(qpwr);
# Build up the power on (-1):
 asp:= (-1)^(a*j);
 for ctr from 0 to (-1+nops(nd)/4) do
   asp := asp * ((-1)*(nd[4*ctr+1]))^i[ctr+1]
 od;
 for ctr from 0 to (-1+nops(dd)/4) do
   asp := asp * dd[4*ctr+1]^k[ctr+1]
 od;
 asp := combine(asp,power);
 # Build up Gaussian polynomials:
 gpp :=1;
 for ctr from 0 to (-1+nops(nd)/4) do
   gpp := gpp * gp(j + nd[4*ctr + 4], i[ctr+1], nd[4*ctr + 3])
 od;
 for ctr from 0 to (-1+nops(dd)/4) do
   gpp := gpp * gp(j+ dd[4*ctr + 4] + k[ctr+1]-1, k[ctr+1], dd[4*ctr+3])
 od;
 gpp := gpp * gp(j+h,h,m);
elimvar := j;
print(f(t)=Sum(asp * t^tpwr * q^qpwr * gpp, op(varset)=0..infinity));
print(`Variable to be eliminated`, elimvar);
newasp :=         subs(elimvar = n-tpwr+elimvar, asp);
newqpwr := expand(subs(elimvar = n-tpwr+elimvar, qpwr));
newgpp :=         subs(elimvar = n-tpwr+elimvar, gpp);
print(P[n]=Sum( newasp * q^newqpwr * newgpp,
op(varset minus {elimvar})=0..infinity)); 
end:

LHS2dual:=proc(a,nd,dd,ed,d,b,c, p)
local asp,newasp, ctr, elimvar, foldsum, g,gpp, nnf, ndf,newgpp,
newqpwr, varlist, varset, qpwr, tpwr;
nnf :=nops(nd)/6;
ndf :=1 + nops(dd)/6;
foldsum := nnf + ndf;
varset := {j};
varlist :=[j];
# Build up the power on t:
tpwr := d*j;
for ctr from 0 to (-1+nops(nd)/6) do
  tpwr := tpwr + nd[6*ctr + 2] * i[ctr+1];
  varset := varset union {i[ctr+1]};
  varlist :=[op(varlist),i[ctr+1]]
od;
varset := varset union{h};
varlist := [op(varlist),h];
tpwr:=tpwr + (ed[2])*h;
for ctr from 0 to (-1 + nops(dd)/6) do
  tpwr := tpwr + dd[6*ctr + 2] * k[ctr+1];
  varset := varset union {k[ctr+1]};
  varlist := [op(varlist),k[ctr+1]]
od;
# Build up the power on q:
qpwr := p - b*j^2 - c*j - h*(ed[5]*j+ed[6]-1);
for ctr from 0 to (-1+nops(nd)/6) do
  qpwr := qpwr 
    - nd[6*ctr + 4] * (i[ctr+1]^2 - i[ctr+1]) /2
    - nd[6*ctr + 3] *  i[ctr+1]
    - nd[6*ctr + 4] *  i[ctr+1]* ( nd[6*ctr+5]*j - i[ctr+1]);
od;
for ctr from 0 to (-1+nops(dd)/6) do
  qpwr := qpwr -  dd[6*ctr+3] * k[ctr+1]
               -  dd[6*ctr+4] * k[ctr+1] * 
           (dd[6*ctr+5]*j+dd[6*ctr+6]+ k[ctr+1]-1 - k[ctr+1]);
od;
qpwr :=expand(qpwr);
# Build up the power on (-1):
asp:= (-1)^(a*j);
for ctr from 0 to (-1+nops(nd)/6) do
  asp := asp * ((-1)*(nd[6*ctr+1]))^i[ctr+1]
od;
asp := asp * ed[1]^h;
for ctr from 0 to (-1+nops(dd)/6) do
  asp := asp * dd[6*ctr+1]^k[ctr+1]
od;
asp := combine(asp,power);
# Build up Gaussian polynomials:
gpp :=1;
for ctr from 0 to (-1+nops(nd)/6) do
  gpp := gpp * gp(nd[6*ctr + 5]*j + nd[6*ctr + 6], i[ctr+1], nd[6*ctr + 4])
od;
for ctr from 0 to (-1+nops(dd)/6) do
  gpp := gpp * gp(dd[6*ctr + 5]*j + dd[6*ctr + 6] + k[ctr+1]-1,
  k[ctr+1], dd[6*ctr+4])
od;
gpp := gpp * gp(ed[5]*j + ed[6] + h -1, h , ed[4]);
elimvar :=j;
#print(foldsum+1,`fold sum over`,varset,`>=0:`,asp * t^tpwr * q^qpwr * gpp,
print(f(t)=Sum(asp * t^tpwr * q^qpwr * gpp, op(varset)=0..infinity));
print(`Variable to be eliminated`, elimvar);
newasp :=         subs(elimvar = n-tpwr+elimvar, asp);
newqpwr := expand(subs(elimvar = n-tpwr+elimvar, qpwr));
newgpp :=         subs(elimvar = n-tpwr+elimvar, gpp);
print(P[n]=Sum( newasp * q^newqpwr * newgpp,
op(varset minus {elimvar})=0..infinity)); 
end:

LHSevendual :=proc(a, nd, dd, m, b, c, p)
local asp, newasp, newerasp, ctr, elimvar, foldsum, g, gpp, nnf, ndf,
newgpp, newergpp, newqpwr, newerqpwr, varlist, varset, qpwr, tpwr;
    subs(n=2*n,p);
    nnf := 1/4*nops(nd);
    ndf := 1 + 1/4*nops(dd);
    foldsum := nnf + ndf;
    g := tpgcd(nd, dd, m);
    varset := {j};
    varlist := [j];
    tpwr := 2*b*j/g;
    for ctr from 0 to -1 + 1/4*nops(nd) do
        tpwr := tpwr + nd[4*ctr + 3]*i[ctr + 1]/g;
        varset := varset union {i[ctr + 1]};
        varlist := [op(varlist), i[ctr + 1]]
    od;
    varset := varset union {h};
    varlist := [op(varlist), h];
    tpwr := tpwr + m*h/g;
    for ctr from 0 to -1 + 1/4*nops(dd) do
        tpwr := tpwr + dd[4*ctr + 3]*k[ctr + 1]/g;
        varset := varset union {k[ctr + 1]};
        varlist := [op(varlist), k[ctr + 1]]
    od;
   qpwr := p - b*j^2 - c*j - m*h*j;
    for ctr from 0 to -1 + 1/4*nops(nd) do qpwr := qpwr
         - 1/2*nd[4*ctr + 3]*(i[ctr + 1]^2 - i[ctr + 1])
         - nd[4*ctr + 2]*i[ctr + 1] +
        nd[4*ctr + 3]*i[ctr + 1]*(i[ctr + 1] - j - nd[4*ctr + 4])
    od;
    for ctr from 0 to -1 + 1/4*nops(dd) do qpwr := qpwr
         - dd[4*ctr + 2]*k[ctr + 1]
         + dd[4*ctr + 3]*k[ctr + 1]*(-j - dd[4*ctr + 4] + 1)
    od;
    qpwr := expand(qpwr);
    asp := (-1)^(a*j);
    for ctr from 0 to -1 + 1/4*nops(nd) do
        asp := asp*(-nd[4*ctr + 1])^i[ctr + 1]
    od;
    for ctr from 0 to -1 + 1/4*nops(dd) do
        asp := asp*dd[4*ctr + 1]^k[ctr + 1]
    od;
    asp := combine(asp, power);
    gpp := 1;
    for ctr from 0 to -1 + 1/4*nops(nd) do gpp :=
        gpp*gp(j + nd[4*ctr + 4], i[ctr + 1], nd[4*ctr + 3])
    od;
    for ctr from 0 to -1 + 1/4*nops(dd) do gpp := gpp*gp(
        j + dd[4*ctr + 4] + k[ctr + 1] - 1, k[ctr + 1],
        dd[4*ctr + 3])
    od;
    gpp := gpp*gp(j + h, j, m);
    elimvar := h;
#       print(foldsum + 1, `fold sum over`, varset, `>=0`,
#        asp * t^tpwr *q^qpwr*gpp, `Var to be eliminated:`, elimvar);
#    newasp := subs(elimvar = 2*n - tpwr + elimvar, asp);
#    newqpwr := expand(subs(elimvar = 2*n - tpwr + elimvar, qpwr));
#    newgpp := subs(elimvar = 2*n - tpwr + elimvar, gpp);
#    print(foldsum, `fold sum over`, varset minus {elimvar},
#        `>=0:`, newasp*q^newqpwr*newgpp);
#    newerasp  := expand(subs(j=n-J,newasp));
#    newerqpwr := expand(subs(j=n-J,newqpwr));
#    newergpp  := expand(subs(j=n-J,newgpp));
#    print(foldsum, `fold sum over`, varset minus {elimvar} minus {j} union {J},
#        `>=0:`, newerasp*q^newerqpwr*newergpp);
print(f(t)=Sum(asp * t^tpwr * q^qpwr * gpp, op(varset)=0..infinity));
print(`Variable to be eliminated`, elimvar);
newasp :=         subs(elimvar = 2*n-tpwr+elimvar, asp);
newqpwr := expand(subs(elimvar = 2*n-tpwr+elimvar, qpwr));
newgpp :=         subs(elimvar = 2*n-tpwr+elimvar, gpp);
print(P[n]=Sum( newasp * q^newqpwr * newgpp,
op(varset minus {elimvar})=0..infinity));
    newerasp  := expand(subs(j=n-J,newasp));
    newerqpwr := expand(subs(j=n-J,newqpwr));
    newergpp  := expand(subs(j=n-J,newgpp));
print(P[n] = Sum( newerasp * q^newerqpwr * newergpp, op(varset minus {elimvar}
  minus {j} union {J})=0..infinity)); 
end:

LHSodddual :=proc(a, nd, dd, m, b, c, p)
local asp, newasp, newerasp, ctr, elimvar, foldsum, g, gpp, nnf, ndf,
newgpp, newergpp, newqpwr, newerqpwr, varlist, varset, qpwr, tpwr;
    subs(n=2*n+1,p);
    nnf := 1/4*nops(nd);
    ndf := 1 + 1/4*nops(dd);
    foldsum := nnf + ndf;
    g := tpgcd(nd, dd, m);
    varset := {j};
    varlist := [j];
    tpwr := 2*b*j/g;
    for ctr from 0 to -1 + 1/4*nops(nd) do
        tpwr := tpwr + nd[4*ctr + 3]*i[ctr + 1]/g;
        varset := varset union {i[ctr + 1]};
        varlist := [op(varlist), i[ctr + 1]]
    od;
    varset := varset union {h};
    varlist := [op(varlist), h];
    tpwr := tpwr + m*h/g;
    for ctr from 0 to -1 + 1/4*nops(dd) do
        tpwr := tpwr + dd[4*ctr + 3]*k[ctr + 1]/g;
        varset := varset union {k[ctr + 1]};
        varlist := [op(varlist), k[ctr + 1]]
    od;
   qpwr := p - b*j^2 - c*j - m*h*j;
    for ctr from 0 to -1 + 1/4*nops(nd) do qpwr := qpwr
         - 1/2*nd[4*ctr + 3]*(i[ctr + 1]^2 - i[ctr + 1])
         - nd[4*ctr + 2]*i[ctr + 1] +
        nd[4*ctr + 3]*i[ctr + 1]*(i[ctr + 1] - j - nd[4*ctr + 4])
    od;
    for ctr from 0 to -1 + 1/4*nops(dd) do qpwr := qpwr
         - dd[4*ctr + 2]*k[ctr + 1]
         + dd[4*ctr + 3]*k[ctr + 1]*(-j - dd[4*ctr + 4] + 1)
    od;
    qpwr := expand(qpwr);
    asp := (-1)^(a*j);
    for ctr from 0 to -1 + 1/4*nops(nd) do
        asp := asp*(-nd[4*ctr + 1])^i[ctr + 1]
    od;
    for ctr from 0 to -1 + 1/4*nops(dd) do
        asp := asp*dd[4*ctr + 1]^k[ctr + 1]
    od;
    asp := combine(asp, power);
    gpp := 1;
    for ctr from 0 to -1 + 1/4*nops(nd) do gpp :=
        gpp*gp(j + nd[4*ctr + 4], i[ctr + 1], nd[4*ctr + 3])
    od;
    for ctr from 0 to -1 + 1/4*nops(dd) do gpp := gpp*gp(
        j + dd[4*ctr + 4] + k[ctr + 1] - 1, k[ctr + 1],
        dd[4*ctr + 3])
    od;
    gpp := gpp*gp(j + h, j, m);
    elimvar := h;
#     print(foldsum + 1, `fold sum over`, varset, `>=0`,
#        asp * t^tpwr *q^qpwr*gpp, `Var to be eliminated:`, elimvar);
#    newasp := subs(elimvar = 2*n+1 - tpwr + elimvar, asp);
#    newqpwr := expand(subs(elimvar = 2*n+1 - tpwr + elimvar, qpwr));
#    newgpp := subs(elimvar = 2*n+1 - tpwr + elimvar, gpp);
#    print(foldsum, `fold sum over`, varset minus {elimvar},
#        `>=0:`, newasp*q^newqpwr*newgpp);
#    newerasp  := expand(subs(j=n-J,newasp));
#    newerqpwr := expand(subs(j=n-J,newqpwr));
#    newergpp  := expand(subs(j=n-J,newgpp));
#    print(foldsum, `fold sum over`, varset minus {elimvar} minus {j} union {J},
#        `>=0:`, newerasp*q^newerqpwr*newergpp);
print(f(t)=Sum(asp * t^tpwr * q^qpwr * gpp, op(varset)=0..infinity));
print(`Variable to be eliminated`, elimvar);
newasp :=         subs(elimvar = 2*n+1 -tpwr+elimvar, asp);
newqpwr := expand(subs(elimvar = 2*n+1 -tpwr+elimvar, qpwr));
newgpp :=         subs(elimvar = 2*n+1 -tpwr+elimvar, gpp);
print(P[n]=Sum( newasp * q^newqpwr * newgpp,
op(varset minus {elimvar})=0..infinity));
    newerasp  := expand(subs(j=n-J,newasp));
    newerqpwr := expand(subs(j=n-J,newqpwr));
    newergpp  := expand(subs(j=n-J,newgpp));
print(P[n] = Sum( newerasp * q^newerqpwr * newergpp, op(varset minus {elimvar}
  minus {j} union {J})=0..infinity)); 

end:

### Tools for annihilating operators:

multops :=proc(A, B)
local antidegA, antidegB, degA, degB, j, k, L, M, prod;
    degA := degree(A, N);
    degB := degree(B, N);
    antidegA := degree(subs(N = 1/N, A), N);
    antidegB := degree(subs(N = 1/N, B), N);
    L := max(antidegA,antidegB);
    M := max(degA, degB);
    prod := add(add(
       coeff(A, N, j)
      *subs({n= n+j, s=s*q^j}, coeff(B, N, k - j)
      *N^k),
        j = -L .. M), k = -2*L .. 2*M);
    sort(collect(prod,N))
end:

factorizeop := proc(R)
 local antidegR, tmp, ctr, degR;
 degR := degree(R,N);
 antidegR := degree( subs(N=1/N, R),N); 
 tmp:=0:
 for ctr from antidegR to degR do
   tmp:= tmp + factor( coeff(R,N,ctr)) * N^ctr;
 od;
 sort(tmp,N)
end: 

rightdiv := proc(AB,B)
local a, antidegAB, antidegB, b, ctr, degAB, degB, eqns, J, K, L, M, neweqn, neweqnLHS, soln, varnames, vars;
# input: AB, B; output A.
    degAB := degree(AB, N);
    degB := degree(B, N);
    antidegAB := degree(subs(N = 1/N, AB), N);
    antidegB := degree(subs(N = 1/N, B), N);
    L := max(antidegAB,antidegB);
    M := max(degAB, degB);
varnames:={};
for ctr from -L to M do varnames := varnames union {a[ctr]} od;
for ctr from -2*L-M to 2*M+L do b[ctr]:= coeff(B,N,ctr) od;
eqns :={};
for K from -2*L to 2*M do
   neweqnLHS := 0;
   for J from -L to M do
     neweqnLHS := neweqnLHS + a[J] * subs({n=n+J, s=s*q^J},b[K-J])
   od;
   neweqn := neweqnLHS = coeff(AB,N,K);
   eqns := eqns union {neweqn}
od;
soln :=solve(eqns,varnames);
RETURN(add(subs(soln, a[J]*N^J), J=-L..M))
end:

annops := proc(a, nd, dd, m, b, c)
   local A, backward, ctr, fctr1, fctr2, fctr3, forward, ft1, ft2, g,L,  maxt;
   A := polyrec(a,nd,dd,m,b,c);
   g := tpgcd(nd,dd,m);
   ft1 := firstterm(a,nd,dd,m,b,c);
   ft2 := fctr(a,nd,dd,m,b,c);
   L := lcm(denom(ft1),denom(ft2));
   if coeff(L,t,0) =-1 then L := -L fi;
   fctr1 := simplify(L*ft1);
   fctr1 := collect(expand(fctr1),t);
   fctr2 := simplify(L*ft2);
   fctr2 := collect(expand(fctr2),t);
   fctr3 := 1+simplify(-L);
   fctr3 := collect(expand(fctr3),t);
   maxt := max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
   forward := 1 - add(coeff(A,P[n-ctr])*N^(-ctr), ctr=1..maxt);
   backward := collect(combine(expand(subs(n=n+maxt,N^maxt *forward)),power),N);
   RETURN(forward,backward)
 end:

annops1 := proc(a, nd, dd, m, b, c)
   local A, backward, ctr, fctr1, fctr2, fctr3, forward, ft1, ft2, g,L,  maxt;
   A := polyrec1(a,nd,dd,m,b,c);
   g := tpgcd(nd,dd,m);
   ft1 := firstterm(a,nd,dd,m,b,c);
   ft2 := fctr(a,nd,dd,m,b,c);
   L := lcm(denom(ft1),denom(ft2));
   if coeff(L,t,0) =-1 then L := -L fi;
   fctr1 := simplify(L*ft1 + L/(1 - t^g) - L*ft2/subs(q=q^g, 1-(t*q)^g));
   fctr1 := collect(expand(fctr1),t);
   fctr2 := simplify(L*ft2);
   fctr2 := collect(expand(fctr2),t);
   fctr3 := 1+simplify(-L);
   fctr3 := collect(expand(fctr3),t);
   maxt := max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
   forward := 1 - add(coeff(A,P[n-ctr])*N^(-ctr), ctr=1..maxt);
   backward := collect(combine(expand(subs(n=n+maxt,N^maxt *forward)),power),N);
   RETURN(forward,backward)
 end:

annops0 := proc(a, nd, dd, m, b, c)
   local A, backward, ctr, fctr1, fctr2, fctr3, forward, ft1, ft2, g,L,  maxt;
   A := polyrec0(a,nd,dd,m,b,c);
   g := tpgcd(nd,dd,m);
   ft1 := firstterm0(a,nd,dd,m,b,c);
   ft2 := fctr(a,nd,dd,m,b,c);
   L := lcm(denom(ft1),denom(ft2));
   if coeff(L,t,0) =-1 then L := -L fi;
   fctr1 := simplify(L*ft1 + L/(1 - t) - L*ft2/subs(q=q^g, 1-t*q));
   fctr1 := collect(expand(fctr1),t);
   fctr2 := simplify(L*ft2);
   fctr2 := collect(expand(fctr2),t);
   fctr3 := 1+simplify(-L);
   fctr3 := collect(expand(fctr3),t);
   maxt := max(degree(fctr3,t),degree(fctr1,t),degree(fctr2,t));
   forward := 1 - add(coeff(A,P[n-ctr])*N^(-ctr), ctr=1..maxt);
   backward := collect(combine(expand(subs(n=n+maxt,N^maxt *forward)),power),N);
   RETURN(forward,backward)
 end: