#    file zeilb for summation of terminating hypergeometric series
#    by Zeilberger's algorithm.
#
#    Copyright 1992, 1996, 1998, 2000 by Tom H. Koornwinder (thk@science.uva.nl)
#
#    The 1992 version was a completely rewritten version of a Maple procedure
#    ident_prover.maple by D. Zeilberger.
#
#    The 1996 version was adapted to Maple V, Release 3 and 4.
#     
#    The 1998 version was adapted to Maple V, Release 4 and 5.
#
#    The present version, 9 December 2000, works in Maple 6, as well as in
#    Maple V, Release 4 and 5.

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

comppol:=proc(f1,g1,x)
local f,g,n,a,b,c,d,j:
option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`:
# Tests for polynomials f1 and g1 in x if f1(x)=g1(x+j) for some nonnegative
# integer j, while f1 and g1 are not constant in x
f:=collect(f1,x):
g:=collect(g1,x):
n:=degree(f,x):
if n=0 or n<>degree(g,x) then
    RETURN(-1):
fi:
a:=coeff(f,x,n):
b:=coeff(f,x,n-1):
c:=coeff(g,x,n):
d:=coeff(g,x,n-1):
j:=normal((b*c-a*d)/n/a/c):
if not type(j,nonnegint) then
    RETURN(-1):
fi:
if collect(c*f-a*subs(x=x+j,g),x)=0 then
    RETURN(j):
else
    RETURN(-1):
fi:
end:

maxshift:=proc(r1,r2,x)
local fr1,fr2,i,j,m,ma,n:
option `Copyright 1992, 1996, 1998 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,comppol(op(1,op(i,op(2,fr1))),op(1,op(j,op(2,fr2))),x)):
    od:
od:
RETURN(ma):
end:

testnonnegintroots:=proc(p,k) local f,fp,i,k0:
option `Copyright 1992, 1996, 1998 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))),k):
    if degree(f,k)=1 and coeff(f,k)=1 and type(-subs(k=0,f),nonnegint) then
        k0:=max(k0,-subs(k=0,f))
    fi:
od:
RETURN(k0):
end:

refactors:=proc(f,n)
local A,i,c,d,g,pol,deg:
option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`:
A:=op(1,f):
g:=1:
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,n)<2 then
        d:=coeff(pol,n,1):
        c:=coeff(pol,n,0)+d:
        if degree(pol,n)=1 then
            A:=A*d^deg:
            c:=c/d:
            g:=g*fac(c,n)^deg:
        else
            A:=A*c^deg:
        fi:
    else
        RETURN(FAIL):
    fi:
od:
A^n*g:
end:

zeilb:=proc()
local A, B, BO, BOdeg, C, D, EE, F, F1, INH, ORDER, P, P0, P1, R1, R2, R3, R4, S,
TO, TOdeg, a, b, co, co1, co2, deg, deg1, deg2, eq, 
g, i, inh, j, k, m, ma, n, n0, nonunique, p, r, ra1, ra2,
s, sol, t, talklevel, time0, va, vb, vbn, var, z;
options `Copyright 1992, 1996, 1998, 2000 by Tom H. Koornwinder`;
time0:=time():
if nargs<5 then
    ERROR(`too few parameters in function call`):
fi:
if not type(args[4],function) or nops(args[4])<>1 then
    ERROR(`fourth argument is not a function of one variable`)
fi:
n:=op(1,args[4]):
if not type(n,name) then
    ERROR(`fourth argument is not of form f(n) with n being a name`):
fi:
ORDER:=args[5]:
if not type(ORDER,nonnegint) then
    ERROR(`fifth argument is not a nonnegative integer`):
fi:
n0:=ORDER-1:
if nargs>5 then
    talklevel:=args[6]:
else
    talklevel:=max(1,printlevel):
fi:
if not type(talklevel,posint) then
    ERROR(`sixth argument is not a positive integer`):
fi:
TO:=args[1]:
if not type(TO,list(ratpoly)) then
    ERROR(`first argument is not a list of rational functions`):
fi:
BO:=args[2]:
if not type(BO,list(ratpoly)) then
    ERROR(`second argument is not a list of rational functions`):
fi:
z:=args[3]:
if not type(z,ratpoly) then
    ERROR(`third 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:
if has(z,n) then
    ERROR(cat(`third argument depends on `,n)):
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,integer),true) then
    ERROR(`some upper index is integer`):
fi:
if has(BO,0) or has(map(type,BO,negint),true) then
    ERROR(`some lower index is nonpositive integer`):
fi:
if ORDER=0 then
    if has(TO,n) or has(BO,n) then
        ERROR(cat(`order=0 while some index depends on `,n)):
    fi:
fi:
if ORDER>0 then
    if not member(-n,TO) then
        ERROR(cat(`order>0 and no upper index is equal to `,`-n`)):
    fi:
    for i from 1 to r do
        if has(TO[i],n) and not type(TO[i],linear(n)) then
            ERROR(cat(`some upper index is not polynomial of degree <=1 in `,n)):
        fi:
        if not type(coeff(TO[i],n),integer) then
            ERROR(cat(`some upper index has non-integer coefficient of `,n)):
        fi:
    od:
    for i from 1 to s do
        if has(BO[i],n) and not type(BO[i],linear(n)) then
            ERROR(cat(`some lower index is not polynomial of degree <=1 in `,n)):
        fi:
        if not type(coeff(BO[i],n),integer) then
            ERROR(cat(`some lower index has non-integer coefficient of `,n)):
        fi:
        if type(coeff(BO[i],n),negint) and type(subs(n=0,BO[i]),integer) then
            ERROR(cat(`some lower index equals integer minus negative integer \
times `,n)):
        fi:
    od:
fi:
if talklevel>2 then
    print(`zeilb: computing P, R1 and R2`):
fi:
if talklevel>3 then
    print(cat(`time=`,convert(time()-time0,string))):
fi:
TOdeg:=map(coeff,TO,n):
BOdeg:=map(coeff,BO,n):
B:=product('fac(TO[t]-TOdeg[t]+k,TOdeg[t])/
        fac(TO[t]-TOdeg[t],TOdeg[t])','t'=1..r)*
        product('fac(BO[t]-BOdeg[t],BOdeg[t])/
        fac(BO[t]-BOdeg[t]+k,BOdeg[t])','t'=1..s):
B:=normal(B):
C:=denom(B):
B:=numer(B):
D:=product('TO[t]+k-1','t'=1..r)/
    product('BO[t]+k-1','t'=1..s)/k*z:
D:=normal(D):
EE:=denom(D):
D:=numer(D):
P0:=product('subs(n=n-i,B)','i'=0..ORDER-1):
for j from 1 to ORDER do
    P0:=P0+b[j]*product('subs(n=n-i,C)','i'=0..j-1)*
        product('subs(n=n-i,B)','i'=j..ORDER-1):
od:
R1:=D*product('subs(n=n-i,k=k-1,B)','i'=0..ORDER-1)/
    EE/product('subs(n=n-i,B)','i'=0..ORDER-1):
R1:=normal(R1):
R2:=denom(R1):
R1:=numer(R1):
ma:=maxshift(R1,R2,k):
j:=0:
P1:=1:
while j <= ma do
    g:=gcd(R1,subs(k=k+j,R2)):
    if g <> 1 then
        R1:=normal(R1/g):
        R2:=normal(R2/subs(k=k-j,g)):
        P1:=P1*product('subs(k=k-p,g)' , 'p'=0..j-1):
    fi:
    j:=j+1:
od:
P:=expand(P0*P1,k):
if talklevel>2 then
    print(`zeilb: computing maximal degree of P`):
fi:
if talklevel>3 then
    print(cat(`time=`,convert(time()-time0,string))):
fi:
R3:=expand(subs(k=k+1,R1)+R2,k):
deg1:=degree(R3,k):
co1:=coeff(R3,k,deg1):
R4:=expand(subs(k=k+1,R1)-R2,k):
if R4=0 then
    deg2:=-1:
else
    deg2:=degree(R4,k):
fi:
co2:=coeff(R4,k,deg1-1):
if deg1 <= deg2 then
    deg:=degree(P,k)-deg2
else
    co:= -2*co2/co1:
    co:=simplify(co):
    if type(co,nonnegint) then
        deg:=max(co,degree(P,k)-deg1+1)
    else
        deg:=degree(P,k)-deg1+1:
    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(`zeilb: 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+1,R1)*F-R2*subs(k=k-1,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 testnonnegintroots(numer(normal(F)),k)>-1 then
        print(`Algorithm fails`):
        if talklevel>1 then
            print(`F(k)=0 at some nonnegative integer`):
        fi:
        RETURN():
    fi:
    A:=product('fac(TO[j],k)','j'=1..r)/
        product('fac(BO[j],k)','j'=1..s)/fac(1,k)*z^k:
    RETURN(factor(subs(k=n+1,R1)*subs(k=n,A)*subs(k=n,F)/subs(k=n,P))):
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 testnonnegintroots(numer(normal(P)),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 testnonnegintroots(numer(normal(F)),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(testnonnegintroots(denom(normal(va[j])),n),n0):
od:
for j from 1 to nops(vb) do
    n0:=max(testnonnegintroots(denom(normal(vb[j])),n),n0):
od:
n0:=max(testnonnegintroots(denom(normal(subs(k=n+1,R2))),n),n0):
n0:=max(testnonnegintroots(numer(normal(subs(k=n+1,P))),n),n0):
if ORDER =1 and not nonunique then
    vbn:=-normal(vb[1]):
    ra1:=factors(numer(vbn)):
    ra2:=factors(denom(vbn)):
    ra1:=refactors(ra1,n):
    ra2:=refactors(ra2,n)
fi:
INH:=op(0,args[4]):
inh:='INH':
if talklevel>2 then
    S:=Sum(inh(n,k)+sum('vb[j]*inh(n-j,k)','j'=1..ORDER),'k'=0..m):
    print(`Write the input series as`):
    print(INH(n)=Sum(inh(n,k),'k'=0..n)):
    print(`Then `):
    print(S=subs(k=m+1,R2)*subs(k=m,F)/subs(k=m+1,P1)/
        product('subs(n=n-i,k=m+1,B)','i'=0..ORDER-1)*inh(n,m+1),n>n0):
   if talklevel>3 then
        print(cat(`time =`,convert(time()-time0,string))):
    fi:
    print(cat(`Hence `,INH,`(`,n,`) equals`)):
fi:
if talklevel=2 then
    print(cat(`The input series `,INH,`(`,n,`) equals`)):
fi:
if ORDER>1 or nonunique or n0>=ORDER or ra1=FAIL or ra2=FAIL then
    RETURN(-sum('factor(vb[i])*INH(n-i)','i'=1..ORDER),n>n0):
fi:
if talklevel>1 then
    print(-sum('factor(vb[i])*INH(n-i)','i'=1..ORDER),n>n0):
    print(cat(`Hence `,INH,`(`,n,`) equals`)):
fi:
RETURN(simplify(ra1/ra2)):
end: