#    file qzeilb for summation of terminating q-hypergeometric series
#    by the q-version of Zeilberger's algorithm.
#
#    Copyright 1992, 1996, 1998, 2000 by Tom H. Koornwinder (thk@science.uva.nl)
#
#    The 1992 version was a competely rewritten version of a Maple procedure
#    q_ident_prover.maple by D. Zeilberger.
#
#    The 1996 version was adapted to Maple V, Release 3 and 4.
#     
#    The 1998 version, July 15, 1998, was adapted to Maple V, Release 4 and 5.
#
#    Thanks to Peter Paule, Axel Riese and Wolfram Koepf for helpful remarks.
#
#    The present version, 11 December 2000, works in Maple 6, as well as in
#    Maple V, Release 4 and 5.
#
#    Thanks to Harald Boeing & Wolfram Koepf for the procedure q_exponent
#    from their package qsum.


`type/homogeneous`:=proc(f,z)
local c,d:
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
d:=degree(f,z):
c:=coeff(f,z,d):
f=c*z^d and not has(c,z):
end:

`type/intqpower`:=proc(x,q)
option `Copyright 1996 by Tom H. Koornwinder`:
x<>0 and type(simplify(ln(x)/ln(q)),integer)
end:

`type/nonnegintqpower`:=proc(x,q)
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
x<>0 and type(simplify(ln(x)/ln(q)),nonnegint)
end:

`type/nonposintqpower`:=proc(x,q)
local b:
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
if x=0 then
    RETURN(false):
fi:
b:=simplify(ln(x)/ln(q)):
b=0 or type(b,negint)
end:

qfac:=proc(a,q,n) local k:
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
if n=0 then
    1
elif type(n,posint) then
    product('1-a*q^k','k'=0..n-1)
elif type(n,negint) then
    product('(1-a*q^(-k))^(-1)','k'=1..-n)
else
    RETURN('procname(args)')
fi
end:

qcomppol:=proc(f1,g1,q,x)
local f,g,n,m,a,b,c,d,j:
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
# Tests for polynomials f1 and g1 in x if f1(x)=f2(q^j*x) for some
# nonnegative integer j, while f1 and f2 are not monomial in x.
# Returns j if test succeeds, returns -1 if test fails.
f:=collect(f1,x):
g:=collect(g1,x):
n:=degree(f,x):
if n=0 or n<>degree(g,x) then
    RETURN(-1):
fi:
m:=ldegree(f,x):
if m=n or m<>ldegree(g,x) then
    RETURN(-1):
fi:
a:=coeff(f,x,n):
b:=coeff(f,x,m):
c:=coeff(g,x,n):
d:=coeff(g,x,m):
j:=collect(a*d/c/b,q):
if not type(j,homogeneous(q)) then
    RETURN(-1):
fi:
j:=degree(j,q)/(n-m):
if not type(j,nonnegint) then
    RETURN(-1):
fi:
g:=collect(subs(x=x*q^j,g),x):
c:=coeff(g,x,n):
if collect(c*f-a*g,x)=0 then
    RETURN(j):
else
    RETURN(-1):
fi:
end:

qmaxshift:=proc(r1,r2,q,x)
local fr1,fr2,i,j,m,ma,n:
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
ma:=-1:
fr1:=factors(r1):
fr2:=factors(r2):
m:=nops(op(2,fr1)):
n:=nops(op(2,fr2)):
for i from 1 to m do
    for j from 1 to n do
        ma:=max(ma,qcomppol(op(1,op(i,op(2,fr1))),op(1,op(j,op(2,fr2))),q,x)):
    od:
od:
RETURN(ma):
end:

testnonnegintqpowerroots:=proc(p,q,x) local f,fp,i,a,b,c,k0:
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
k0:=-1:
fp:=factors(p):
for i from 1 to nops(op(2,fp)) do
    f:=collect(op(1,op(i,op(2,fp))),x):
    if degree(f,x)=1 then
        a:=coeff(f,x,1):
        b:=coeff(f,x,0):
        c:=collect(-b/a,q):
        if type(c,nonnegintqpower(q)) then
            k0:=max(k0,simplify(ln(c)/ln(q))):
        fi:
    fi:
od:
RETURN(k0):
end:

qrefactors:=proc(f,x,q,n)
local A,i,c,d,g,qc,pol,deg;
option `Copyright 1992, 1996 by Tom H. Koornwinder`:
A:=op(1,f);
g:=1;
qc:=0;
for i from 1 to nops(op(2,f)) do
    pol:=op(1,op(i,op(2,f)));
    deg:=op(2,op(i,op(2,f)));
    if degree(pol,x)<2 then
        c:=coeff(pol,x,0);
        d:=q*coeff(pol,x,1);
        if c=0 then
            A:=A*d^deg;
            qc:=qc+deg
        elif d=0 then
            A:=A*c^deg
        else
            A:=A*c^deg;
            d:=-d/c;
            g:=g*qfac(d,q,n)^deg
        fi
    else
        RETURN(FAIL)
    fi
od;
A^n*q^(n*(n-1)/2*qc)*g
end:

q_exponent:= proc(ff, q)
local f;
    option `Copyright (c) 1997 by Harald Boeing & Wolfram Koepf.`;
    f:= combine(factor(ff), power);
    if type(f, identical(q)^anything) then
        op(2,f);
    elif (f = q) then
        1;
    elif (f = 1) then
        0;
    else 
        FAIL;
    fi;
end:

qzeilb:=proc()
local A, B, BO, BOdeg, C, D, E, F, F1, F2, INQH, K, N, P, P0, P1, ORDER, R1,
R2, S, TO, TOdeg, a, b, co, co1, co2, deg, deg1, deg2, eq,
g, i, inqh, j, k, m, ma, n, n0, nonunique, nvb, p, q, r, ra1, ra2,
s, sol, t, talklevel, time0, va, vb, var, z;
option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`;
time0:=time():
if nargs<6 then
    ERROR(`too few parameters in function call`):
fi:
q:=args[3]:
if not type(q,name) then
    ERROR(`third argument is not a name`):
fi:
if not type(args[5],function) or nops(args[5])<>1 then
    ERROR(`fifth argument is not a function of one variable`)
fi:
n:=op(1,args[5]):
if not type(n,name) then
    ERROR(`fifth argument is not of form f(n) with n being a name`):
fi:
ORDER:=args[6]:
if not type(ORDER,nonnegint) then
    ERROR(`sixth argument is not a nonnegative integer`):
fi:
n0:=ORDER-1:
if nargs>6 then
    talklevel:=args[7]:
else
    talklevel:=max(1,printlevel):
fi:
if not type(talklevel,posint) then
    ERROR(`seventh argument is not a positive integer`):
fi:
TO:=simplify(subs(n=ln(N)/ln(q),args[1])):
if not type(TO,list(ratpoly)) then
    ERROR(`first argument is not a list of rational functions`):
fi:
BO:=simplify(subs(n=ln(N)/ln(q),args[2])):
if not type(BO,list(ratpoly)) then
    ERROR(`second argument is not a list of rational functions`):
fi:
z:=simplify(subs(n=ln(N)/ln(q),args[4])):
if not type(z,ratpoly) then
    ERROR(`fourth argument is not a rational function`):
fi:
TO:=map(normal,TO,expanded):
BO:=map(normal,BO,expanded):
z:=normal(z):
if z=0 then
    RETURN(1):
fi:
for i while i <= nops(TO) do
    if member(TO[i], BO, 'j') then
        TO := subsop(i = NULL, TO):
        BO := subsop(j = NULL, BO):
        i := i-1:
    fi:
od:
r:=nops(TO):
s:=nops(BO):
if has(map(type,TO,intqpower(q)),true) then
    ERROR(cat(`some upper index is integer  power of `,q)):
fi:
if has(map(type,BO,nonposintqpower(q)),true) then
    ERROR(cat(`some lower index is nonpositive integer of `,q)):
fi:
if ORDER=0 then
    if has(z,N) or has(TO,N) or has(BO,N) then
        ERROR(cat(`order=0 while argument or some index depends on `,n)):
    fi:
fi:
if ORDER>0 then
    if not member(1/N,TO) then
        ERROR(cat(`order>0 and no upper index is equal to `,
          convert(q^(-n),string))):
    fi:
    TO:=map(collect,TO,N):
    BO:=map(collect,BO,N):
    z:=collect(z,N):
    if not type(z,homogeneous(N)) then
        ERROR(cat(`fourth argument is not homogeneous in `,
          convert(q^(-n),string))):
    fi:
    if has(map(type,TO,homogeneous(N)),false) then
        ERROR(cat(`some upper index is not homogeneous in `,
          convert(q^(-n),string))):
    fi:
    if has(map(type,BO,homogeneous(N)),false) then
        ERROR(cat(`some lower index is not homogeneous in `,
          convert(q^(-n),string))):
    fi:
    TOdeg:=map(degree,TO,N):
    BOdeg:=map(degree,BO,N):
    deg:=degree(z,N):
    for i from 1 to s do
        if BOdeg[i]<0 and type(coeff(BO[i],N,BOdeg[i]),intqpower(q)) then
            ERROR(cat(`some lower index equals integer power of `,q,` times \
              negative integer power of`, convert(q^(-n),string))):
        fi:
    od:
fi:
if talklevel>2 then
    print(`qzeilb: computing P, R1 and R2`):
fi:
if talklevel>3 then
    print(cat(`time=`,convert(time()-time0,string))):
fi:
B:=product('qfac(TO[t]/q^TOdeg[t]*K,q,TOdeg[t])/
        qfac(TO[t]/q^TOdeg[t],q,TOdeg[t])','t'=1..r)*
        product('qfac(BO[t]/q^BOdeg[t],q,BOdeg[t])/
        qfac(BO[t]/q^BOdeg[t]*K,q,BOdeg[t])','t'=1..s)*K^deg:
B:=normal(B):
C:=denom(B):
B:=numer(B):
D:=product('1-TO[t]*K/q','t'=1..r)/
    product('1-BO[t]*K/q','t'=1..s)/(1-K)*(-K/q)^(s-r+1)*z:
D:=normal(D):
E:=denom(D):
D:=numer(D):
P0:=product('subs(N=N/q^i,B)','i'=0..ORDER-1):
for j from 1 to ORDER do
    P0:=P0+b[j]*product('subs(N=N/q^i,C)','i'=0..j-1)*
        product('subs(N=N/q^i,B)','i'=j..ORDER-1):
od:
R1:=D*product('subs(N=N/q^i,K=K/q,B)','i'=0..ORDER-1)/
E/product('subs(N=N/q^i,B)','i'=0..ORDER-1):
R1:=normal(R1):
R2:=denom(R1):
R1:=numer(R1):
ma:=qmaxshift(R1,R2,q,K):
j:=0:
P1:=1:
while j<= ma do
    g:=gcd(R1,subs(K=K*q^j,R2)):
    if g <> 1 then
        R1:=normal(R1/g):
        R2:=normal(R2/subs(K=K/q^j,g)):
        P1:=P1*product('subs(K=K/q^p,g)' , 'p'=0..j-1):
    fi:
    j:=j+1:
od:
P:=expand(P0*P1,K):
R1:=expand(R1,K):
deg1:=ldegree(R1,K):
co1:=coeff(R1,K,deg1):
R2:=expand(R2,K):
deg2:=ldegree(R2,K):
co2:=coeff(R2,K,deg2):
if deg1 <> deg2 then
    deg:=ldegree(P,K)-min(deg1,deg2):
else
    co:= co2/co1:
#    co:=simplify(ln(co)/ln(q)):   (causes wrong results in Maple V4)
#    co:=expand(ln(co)/ln(q)):     (causes wrong results in Maple V5)
    co:=q_exponent(co2/co1,q):
    if type(co,integer) then
        deg:=min(co,ldegree(P,K))-deg1:
    else
        deg:=ldegree(P,K)-deg1:
    fi:
fi:
P:=expand(K^(-deg)*P,K):
if talklevel>2 then
    print(`qzeilb: computing maximal degree of P`):
fi:
if talklevel>3 then
    print(cat(`time=`,convert(time()-time0,string))):
fi:
if ldegree(P,K)<0 then
    print(`Algorithm fails`):
    if talklevel>1 then
        print(`ldegree(P)<0`):
    fi:
    RETURN():
fi:
P1:=K^(-deg)*P1:
deg1:=degree(R1,K):
co1:=coeff(R1,K,deg1):
R2:=expand(q^(-deg)*R2,K):
deg2:=degree(R2,K):
co2:=coeff(R2,K,deg2):
if deg1 <>deg2 then
    deg:=degree(P,K)-max(deg1,deg2):
else
    co:= co2/co1:
#    co:=simplify(ln(co)/ln(q)): (causes wrong results in Maple V4)
#    co:=expand(ln(co)/ln(q)):   (causes wrong results in Maple V5)
    co:=q_exponent(co2/co1,q):
    if type(co,integer) then
#        deg:=max(co,degree(P,K)-deg1) (was bug in earlier version) 
        deg:=max(co,degree(P,K))-deg1
    else
        deg:=degree(P,K)-deg1:
    fi:
fi:
if deg < 0 then
    print(`Algorithm fails`):
    if talklevel>1 then
        print(`degree(F)<0`):
    fi:
    RETURN():
fi:
if talklevel>2 then
    print(`qzeilb: solving equations`):
fi:
if talklevel>3 then
    print(cat(`time=`,convert(time()-time0,string))):
fi:
F:=sum('a[i]*K^i','i'=0..deg):
F1:=expand(subs(K=K*q,R1)*F-R2*subs(K=K/q,F)-P):
deg1:=degree(F1,K):
eq:=seq(coeff(F1,K,i),i=0..deg1):
va := seq(a[i],i=0..deg):
vb := seq(b[i],i=1..ORDER):
var:=vb,va:
sol:=solve({eq},{var}):
#     eq := Groebner['gbasis']([eq], 'plex'(var) );
#     eq := convert(eq, 'set'): 
#     # select equations involving b[i]'s
#     eqb := select(has,eq,{vb}):
#     eqa := eq minus eqb: # equation with a[i]'s only
#     sola := solve(eqa,{va}):
#     solb := solve(eqb,{vb}):
#     sol:=sola union solb:
if sol=NULL then
    print(`Algorithm fails`):
    if talklevel>1 then
        print(`System of equations has no solutions`):
    fi:
    RETURN():
fi:
if ORDER=0 then
    F:=subs(sol,F):
    if testnonnegintqpowerroots(numer(normal(F)),q,K)>-1 then
        print(`Algorithm fails`):
        if talklevel>1 then
            print(`F(k)=0 at some nonnegative integer`):
        fi:
        RETURN():
    fi:
    A:=product('qfac(TO[j],q,k)','j'=1..r)/
        product('qfac(BO[j],q,k)','j'=1..s)/qfac(q,q,k)*
            ((-1)^k*q^(k*(k-1)/2))^(s-r+1)*z^k:
    F2:=subs(K=q^(n+1),R1)*subs(k=n,A)*subs(K=q^n,F)/subs(K=q^n,P):
    RETURN(factor(F2)):
fi:
nonunique:=has({seq(op(2,sol[i]),i=1..deg+ORDER+1)},{var});
va:=subs(sol,[va]):
vb:=subs(sol,[vb]):
if normal(1+sum('vb[i]','i'=1..ORDER))=0 then
    print(`Algorithm fails`):
    if talklevel>1 then
        print(`a(0) is equal to 0`):
    fi:
    RETURN():
fi:
P:=subs(sol,P):
if testnonnegintqpowerroots(numer(normal(P)),q,k)>-1 then
    print(`Algorithm fails`):
    if talklevel>1 then
        print(`P(k)=0 at some nonnegative integer`):
    fi:
    RETURN():
fi:
F:=subs(sol,F):
if testnonnegintqpowerroots(numer(normal(F)),q,k)>-1 then
    print(`Algorithm fails`):
    if talklevel>1 then
        print(`F(k)=0 at some nonnegative integer`):
    fi:
    RETURN():
fi:
for j from 1 to nops(va) do
    n0:=max(testnonnegintqpowerroots(denom(normal(va[j])),q,N),n0):
od:
for j from 1 to nops(vb) do
    n0:=max(testnonnegintqpowerroots(denom(normal(vb[j])),q,N),n0):
od:
n0:=max(testnonnegintqpowerroots(denom(normal(subs(K=q*N,R2))),q,N),n0):
n0:=max(testnonnegintqpowerroots(numer(normal(subs(K=q*N,P))),q,N),n0):
if ORDER =1 and not nonunique then
    nvb:=-normal(vb[1]):
    ra1:=factors(numer(nvb)):
    ra2:=factors(denom(nvb)):
    ra1:=qrefactors(ra1,N,q,n):
    ra2:=qrefactors(ra2,N,q,n)
fi:
INQH:=op(0,args[5]):
inqh:='INQH':
if talklevel>2 then
    S:=Sum(inqh(n,k)+sum('subs(N=q^n,vb[j])*inqh(n-j,k)','j'=1..ORDER),
        'k'=0..m):            
    print(`Write the input series as`):
    print(INQH(n)=Sum(inqh(n,k),'k'=0..n)):
    print(`Then `):
    print(S=factor(subs(N=q^n,K=q^m*q,R2)*subs(N=q^n,K=q^m,F)/
        subs(N=q^n,K=q^m*q,P1)/
        product('subs(N=q^n/q^i,K=q^m*q,B)','i'=0..ORDER-1))*
        inqh(n,m+1),n>n0):
    if talklevel>3 then
    print(cat(`time=`,convert(time()-time0,string))):
    fi:
    print(cat(`Hence `,INQH,`(`,n,`) equals`)):
fi:
if talklevel=2 then
    print(cat(`The input series `,INQH,`(`,n,`) equals`)):
fi:
if ORDER>1 or nonunique or n0>=ORDER or ra1=FAIL or ra2=FAIL then
    RETURN(-sum('factor(subs(N=q^n,vb[i]))*INQH(n-i)','i'=1..ORDER),n>n0):
fi:
if talklevel>1 then
    print(-sum('factor(subs(N=q^n,vb[i]))*INQH(n-i)','i'=1..ORDER),n>n0):
    print(cat(`Hence `,INQH,`(`,n,`) equals`)):
fi:
RETURN(simplify(ra1/ra2)):
end: