print(`This is the GUESS WZ CERTIFICATE PACKAGE for Maple`);
print(`Written by Drew Sills and Yuriy Choliy`);
print(`Version of July 14, 2004`);
print(`The main procedure is GuessCert`);

with(LinearAlgebra):
Verbose:=false:

## Rising factorial function
rf:=proc(a,n) local i;
   (a+n-1)!/(a-1)!
end:

GuessCert:=proc() local Args,F, n,k;
Args:=[args];
if type(Args[1],list) then
  if nops(Args)=3 then GuessCertSemiNumHG(Args[1],Args[2],Args[3],10)
  else GuessCertSemiNumHG(Args[1],Args[2],Args[3],Args[4])
  fi
else
  if nops(Args)=3 then GuessCertSemiNum(Args[1],Args[2],Args[3],10)
  else GuessCertSemiNum(Args[1],Args[2],Args[3],Args[4])
  fi
fi
end:

GenRatForm:=proc(x,c,d,NumDeg,DenDeg) local i;
add(c[i]*x^i, i=0..NumDeg)/add(d[i]*x^i, i=0..DenDeg)
end:

GenRatForm2Var:=proc(x,y,c,d,NumDegx,NumDegy,DenDegx,DenDegy) local i,j:
add(add( c[i,j]*x^i*y^j, i=0..NumDegx),j=0..NumDegy)/
add(add( d[i,j]*x^i*y^j, i=0..DenDegx),j=0..DenDegy)
end:

GenSpecPol:=proc(k,x,degk) local i;
  add(k^i*x[i],i=0..degk)
end:

GenSpecPol2Var:=proc(n,k,x,degn,degk) local i,j:
   add(add( x[i,j]*n^i *k^j,i=0..degn),j=0..degk)
end:

GuessCertSemiNumOld:=proc(F,n,k,min,MaxDeg) local i,j, v; global r_R;
### find r_R[min] here
r_R[min]:=0;
for i from 0 to MaxDeg do
  ###****find r_R[min+i+1]
   r_R[min+i+1]:=
     normal(
       expand(
         F(n,min+i)*r_R[min+i]/F(n,min+i+1) + F(n+1,min+i)/F(n,min+i+1) -
         F(n,min+i)/F(n,min+i+1)
              )
            );
   if Verbose then print(R(n,min+i+1) = r_R[min+i+1]) fi;
   for j from 0 to i do
     v:=GuessCertSemiNum1(F,n,k,min,i-j,j);
     if v<>0 then RETURN(v) fi
   od
od;
RETURN(Fail)
end:

GuessCertSemiNum:=proc(F,n,k,MaxDeg) local kval,i, v; global r_R;
r_R[0]:=0;
for kval from 0 to 1 do
   r_R[kval+1]:=
     normal(
       expand(
         F(n,kval)*r_R[kval]/F(n,kval+1) + F(n+1,kval)/F(n,kval+1) -
         F(n,kval)/F(n,kval+1)
              )
            );
   if Verbose then print(R(n,kval+1) = r_R[kval+1]) fi;
od;
for i from 1 to MaxDeg do
  for kval from 2*i to 2*i+1 do
    r_R[kval+1]:=
     normal(
       expand(
         F(n,kval)*r_R[kval]/F(n,kval+1) + F(n+1,kval)/F(n,kval+1) -
         F(n,kval)/F(n,kval+1)
              )
            );
    if Verbose then print(R(n,kval+1) = r_R[kval+1]) fi;
  od;
v:=GuessCertSemiNum1(0,i,i);
if v<>0 then RETURN(v) fi
od;
RETURN(Fail)
end:

GuessCertSemiNumHG:=proc(F,RHS1,n,MaxDeg)
local a,b,z,r,s,k,kval,i,RHS, v,MainRatio,NRatio; global r_R; #, r_Rnum,
r_Rden;
RHS:=RHS1;
a:=op(1,F); b:=op(2,F); z:=op(3,F);
r:=nops(a); s:=nops(b);
MainRatio:= (k)->
   convert([seq(b[i]+k,i=1..s)],`*`)
  /convert([seq(a[i]+k,i=1..r)],`*`)*(k+1)/z;
if RHS = 0 then RHS:=1 fi;
NRatio:= (k)->
normal(expand( RHS/subs(n=n+1,RHS))*
expand(subs(n=n+1, z^k*
  convert([seq(rf(a[i],k),i=1..r)],`*`)
/convert([seq(rf(b[i],k),i=1..s)],`*`))
/z^k/
convert([seq(rf(a[i],k),i=1..r)],`*`)
*convert([seq(rf(b[i],k),i=1..s)],`*`)));
if Verbose then print(`Main ratio`,MainRatio(k));
                print(`N ratio`,NRatio(k))
fi;
r_R[0]:=0;
#for kval from 0 to 1 do
for kval from 0 to r+s+1 do
   r_R[kval+1]:=
     normal(expand( MainRatio(kval)* (NRatio(kval)-1+r_R[kval])));
#  r_Rnum[kval+1]:= numer(r_R[kval+1]);
#  r_Rden[kval+1]:= denom(r_R[kval+1]);
   if Verbose then print(R(n,kval+1) = r_R[kval+1]) fi;
od;
for i from r to MaxDeg do
  for kval from 2*i+1 to 2*i+2 do
     r_R[kval+1]:=
     normal(expand( MainRatio(kval)* (NRatio(kval)-1+r_R[kval])));
  #  r_Rnum[kval+1]:= numer(r_R[kval+1]);
  #  r_Rden[kval+1]:= denom(r_R[kval+1]);
    if Verbose then print(R(n,kval+1) = r_R[kval+1]) fi;
  od;
v:=GuessCertSemiNum1(0,i,i-1);
if v<>0 then RETURN(v) fi
od;
RETURN(Fail)
end:

GuessCertFullNum:=proc(F,n,k,MaxDeg) local i,j,kval,v,min; global r_R;
### find r_R[min] here
r_R[0]:=0;
#for i from 0 to 2+2*MaxDeg do
for kval from 0 to 1 do
    r_R[kval+1]:=
      normal(
    expand(F(n,kval)/F(n,kval+1) *
    (F(n+1,kval)/F(n,kval)-1+r_R[kval])
            ));
  if Verbose then print(R(n,kval+1)=r_R[kval+1]) fi
od;
#   r_R[kval+1]:=
#     normal(
#       expand(
#         F(n,kval)*r_R[kval]/F(n,kval+1) + F(n+1,kval)/F(n,kval+1) -
#         F(n,kval)/F(n,kval+1)
#              )
#            );
for i from 1 to MaxDeg do
   for kval from 2*i to 2*i+1 do
     r_R[kval+1]:= normal(expand(
       F(n,kval)/F(n,kval+1) *
       (F(n+1,kval)/F(n,kval)-1+r_R[kval])));
      if Verbose then print(R(n,kval+1)=r_R[kval+1]) fi
   od;
  v:=GuessCertFullNum1(F,n,k,i,i,i,i);
  if v<>0 then RETURN(v) fi
od;
RETURN(Fail)
end:

GuessCertSemiNum1:=proc(min,NumDeg,DenDeg)
local a,b,i,max,vars,eqns,soln,test;
vars:={seq(a[i],i=0..NumDeg)} union {seq(b[i],i=0..DenDeg)};
if Verbose then print(`==========================================`) fi;
if Verbose then print(`Trying with numerator degree =`,NumDeg);
                print(`denominator degree = `,DenDeg)
fi;
if Verbose then print(vars) fi;
max:=min+NumDeg+DenDeg+1;
eqns:={seq(r_R[i]*GenSpecPol(i,b,DenDeg) =
GenSpecPol(i,a,NumDeg),i=min..max)};
#eqns:={seq(r_Rnum[i]*GenSpecPol(i,b,DenDeg) -
#r_Rden[i]*GenSpecPol(i,a,NumDeg), i=min..max)};
if Verbose then print(eqns) fi;
soln:=SolveTools[Linear](eqns,vars);
if Verbose then print(soln) fi;
if Verbose then print(vars) fi;
test:= op(subs(soln,vars));
if   test=0 then 0
else factor(normal(subs(soln, GenRatForm(k,a,b,NumDeg,DenDeg))))fi

end:

GuessCertFullNum1:=proc(F,n,k,NumDegn,NumDegk,DenDegn,DenDegk)
local
a,b,i,max,vars,eqns,soln,test,NumbUnknowns,nval,kval,mink,maxk,minn,maxn,tmp,ctr,A,B,Bnew,Rnumer,Rdenom,ZeroVec,VectorSoln,ListSoln,j;
#global r_R;
vars:={seq(seq(a[i,j],i=0..NumDegn),j=0..NumDegk)} union
{seq(seq(b[i,j],i=0..DenDegn),j=0..DenDegk)};
if Verbose then print(`==========================================`) fi;
#if Verbose then print(`Trying with numerator degree for `,n=NumDegn);
#                print(`numerator degree for`,k=NumDegk);
#                print(`denominator degree for`,n=DenDegn);
#                print(`denominator degree for`,k=DenDegk)
#fi;
if Verbose then print(`Attempting with max degrees`, NumDegn) fi;

#if Verbose then print(`Set of Variables:`,vars) fi;
NumbUnknowns:= (NumDegn+1)*(NumDegk+1) + (DenDegn+1)*(DenDegk+1);
#if Verbose then print(`number of unknowns,`=NumbUnknowns) fi;
mink:=0: maxk:=nval:
minn:=0: maxn:= ceil(-3/2 + sqrt(1+8*NumbUnknowns)/2 );
#print(mink, maxk, minn, maxn);
#eqns:={seq(seq(
#    subs(n=nval,r_R[kval])*GenSpecPol2Var(nval,kval,b,DenDegn,DenDegk)
#      =    GenSpecPol2Var(nval,kval,a,NumDegn,NumDegk),kval=mink..nval),
#nval=minn..maxn)};
#if Verbose then print(`Original system of`,nops(eqns), `equations`,eqns)

A:=
[seq(
   seq(
        [seq(
            seq( subs(n=nval,r_R[kval]) * nval^i * kval^j,
                j=0..NumDegk),
         i=0..NumDegn),
         seq(
            seq(
                 - (nval^i * kval^j),
            j=0..DenDegk),
         i=0..DenDegn)],
      kval=mink..nval),
nval=minn..maxn)];
#B:= [seq(convert(op(i,A),Vector),i=1..nops(A))];
B:=  Basis([seq( convert(op(i,A),Vector[row]), i=1..nops(A))]);
if Verbose then print(`A matrix`,A) fi;
if Verbose then print(`B matrix`,B) fi;
ZeroVec:=Vector(nops(B));
Bnew:=  Matrix([seq(convert(op(i,B),list),i=1..nops(B))]);
if Verbose then print(`B matrix`,Bnew) fi;
VectorSoln:=LinearSolve(Bnew,ZeroVec);
ListSoln:=convert(VectorSoln,list);
if Verbose then
   print(`Solution to system:`,ListSoln);fi;

if {op(ListSoln)}={0} then RETURN(0) fi;
tmp:=0;
ctr:=1;
for i from 0 to DenDegn do
  for j from 0 to DenDegk do
     tmp:= tmp+ op(ctr,ListSoln)*n^i* k^j;
     ctr:=ctr+1
  od
od;
Rdenom:=tmp;
tmp:=0;
for i from 0 to NumDegn do
  for j from 0 to NumDegk do
     tmp:= tmp+ op(ctr,ListSoln)*n^i*k^j;
     ctr:=ctr+1
  od
od;
Rnumer:=tmp;
factor(normal(Rnumer/Rdenom))
## Remove extra equations:
#if nops(eqns)<NumbUnknowns then RETURN(0) fi;
#eqns:= {seq( op(i,eqns), i=1..NumbUnknowns)};
#if Verbose then print(`Systemof`,nops(eqns),`equationsto besolved`,eqns) 
#fi;
#soln:=SolveTools[Linear](eqns,vars);
#if Verbose then print(`Solution to the system`,soln) fi;
#if Verbose then print(`Set of variables`,vars) fi;
#test:= op(subs(soln,vars));
#if   test=0 then 0
#else factor(normal(subs(soln,
#   GenRatForm2Var(n,k,a,b,NumDegn,NumDegk,DenDegn,DenDegk))))
#fi
end:

################################################

###This is a procedure for comparison of speeds
###between GuessWZCert method and zeil algorithm.
###z_input should be a list of parameters in form
###suitable for input to zeil procedure.
timing:=proc(code,hl,hr,z_input,n) local ctr,
ourtime,zeiltime,ourAv,zeilAv,ourMaxDev,zeilMaxDev:
print(code);
ourtime:=[seq(time(GuessCert(hl,hr,n) ),ctr=1..10)]:
zeiltime:=[seq(time(zeil(z_input,k,n,N)),ctr=1..10)]:
#ourAv:=((add(op(ctr,ourtime)),ctr=1..10))/10:
ourAv:= (add(ourtime[ctr],ctr=1..10))/10:
#zeilAv:=((add(op(ctr,zeiltime)),ctr=1..10))/10:
zeilAv:= (add(zeiltime[ctr],ctr=1..10))/10:
#ourMaxDev:=max(abs(op(ctr,ourtime)-ourAv)):
#zeilMaxDev:=max(abs(op(ctr,zeiltime)-zeilAv)):
print(`Our times are:`);
print(ourtime);
#print [op(ctr,ourtime),ctr=1..10];
print (`Zeilberger's times are:`);
print(zeiltime);
#print [op(ctr,zeiltime),ctr=1..10];
print (`Our average time is:`, ourAv);
print (`Zeilberger's average time is:`, zeilAv);
print (`Ratio (EKHAD's average time)/(our average time) is:`);
print ( zeilAv/ourAv);
end:




################################################


Print (`Written by Andrew Sills and Yuriy Choliy`);
print (`version of July 13 2004`);

print(`This file contains examples of hypergeometric identities`);
print ( `It is primarily used to test new methods for finding`);
print (`and/or verifying these identities`);
print (` For every identity there are three variables:`);
print (`left-hand side(hl stands for`);
print (`hypergeometric left) right-hand side(rh) and input form for`);
print( `Zeilberger's algorithm `);
print (` from EKHAD8 packege(z). Each variable's name begins with`);
print(`abbreviated name of `);
print (`corresponding identity using capital letters. For example,`);
print(`right-hand side of Pfaff-Saalschutz's `);
print (`identity has a name PShr. Below is list of all identities and
their`);
print( `abbreviation:`);
print (`Chu-Vandermonde's: CV`);
print (`Pfaff-Saalschultz's: PS`);
print (`Gauss's Theorem(2F1(: GT`);
print (`Kummer's Theorem(2F1): K`);
print (`Gauss's Second Theorem(2F1): GT2`);
print (`Bailey's Theorem: B`);
print (`Dixon's Theorem I (3F2) : D1;`);
print (`Dixon's Theorem IA (3F2): D1A;`);
print (`No name 4F3 identity: NN1;`);
print (`MRR 5F4 identity: MRR;`);
print (`No name 4F3 (A) (if c=-n): NN1A;`);
print (`No name 5F4: NN2;`);
print (`No name 5F4 (A): NN2A;`);
print (`Dougall's Theorem (7F6): DOUG;`);
print (`Nearly-poised summation t-m; NP1;`);


###  Chu-Vandermonde
CVhl:=[ [a,-n],[c],1]:
CVhr:= rf(c-a,n)/rf(c,n):
CVz:= (rf(a,k)*rf(-n,k)/(k!*rf(c,k)))/CVhr:


###  Pfaff-Saalschutz
PShl:=[ [a,b,-n], [c,a+b-c-n+1], 1]:
PShr:= rf(c-a,n)*rf(c-b,n)/ (rf(c,n)*rf(c-a-b,n)):
PSz:= (rf(a,k) * rf(b,k) * rf(-n,k)/ (k!*rf(c,k)*rf(a+b-c-n+1,k))) /PShr:

###  Gauss's Theorem
GThl:=[ [n,b],[c],1]:
GThr:= rf(1,c-n-b-1)*rf(1,c-1)/(rf(1,c-n-1)*rf(1,c-b-1)):
GTz:= (rf(n,k)*rf(b,k)/(k!*rf(c,k)))/GThr:

###  Kummer's Theorem
Khl:=[ [a,n],[1+a-n],-1]:
Khr:=rf(1,a-n)*rf(1,a/2)/(rf(1,a)*rf(1,a/2-n)):
Kz:= (rf(a,k)*rf(n,k)*(-1)^k/(k!*rf(1+a-n,k)))/Khr:

###  Gauss's Second Theorem
GT2hl:=subs(a=2*n,[ [a,b],[1/2+(a+b)/2],1/2]):
GT2hr:=subs(a=2*n,
Gamma(1/2)*rf(1,a/2+b/2-1/2)/(rf(1,a/2-1/2)*rf(1,b/2-1/2))):
GT2z:=subs(a=2*n,(rf(a,k)*rf(b,k)*(1/2)^k/(k!*rf(1/2+a/2+b/2,k)))/GT2hr):


### Bailey's Theorem
Bhl:=subs(a=2*n,[ [a,1-a],[c],1/2]):
Bhr:=subs(a=2*n,rf(1,c/2-1)*rf(1,c/2-1/2)/(rf(1,c/2+a/2-1)*rf(1,c/2-a/2-1/2))):
Bz:=subs(a=2*n,(rf(a,k)*rf(1-a,k)*(1/2)^k/(k!*rf(c,k)))/Bhr):


###  Dixon's Theorem I
D1hl:=subs(c=n,[[a,b,c],[1+a-b,1+a-c],1]):
D1hr:=subs(c=n,rf(1,a/2)*rf(1,a-b)*rf(1,a/2-b-c)
  /(rf(1,a)*rf(1,a/2-b)*rf(1,a/2-c)*rf(1,a-b-c))):
D1z:=subs(c=n,(rf(a,k)*rf(b,k)*rf(c,k)/(k!*rf(1+a-b,k)*rf(1+a-c,k)))/D1hr):


###Dixon's Theorem IA (We assume c=-n)
D1Ahl:=[[a,b,-n],[1+a-b,1+a+n],1]:
D1Ahr:=rf(1+a,n)*rf(1+a/2-b,n)/(rf(1+a/2,n)*rf(1+a-b,n)):
D1Az:=(rf(a,k)*rf(b,k)*rf(-n,k)/(k!*rf(1+a-b,k)*rf(1+a+n,k)))/D1Ahr:


###No name 4F3
NN1hl:=subs(c=n,[[a,1+a/2,b,c],[a/2,1+a-b,1+a-c],-1]):
NN1hr:=subs(c=n,rf(1,a-b)*rf(1,a-c)/(rf(1,a)*rf(1,a-b-c))):
NN1z:=subs(c=n,
  (rf(a,k)*rf(1+a/2,k)*rf(b,k)*rf(c,k)*(-1)^k/(k!*rf(a/2,k)*rf(1+a-b,k)
  *rf(1+a-c,k)))/NN1hr):



### MRR
MRRhl:= [[-2*n-1,x+2*n+2,x-z+1/2,x+n+1,z+n+1],[x/2+1/2, x/2+1, 2*z+2*n+2,
2*x-2*z+1],1]:
MRRhr:= 0:


###No name 4F3 (A) (if c=-n)
NN1Ahl:=[[a,1+a/2,b,-n],[a/2,1+a-b,1+a+n],-1]:
NN1Ahr:=rf(1+a,n)/rf(1+a-b,n):
NN1Az:=
(rf(a,k)*rf(1+a/2,k)*rf(b,k)*rf(-n,k)*(-1)^k/(k!*rf(a/2,k)*rf(1+a-b,k)*rf(1+a+n,k)))/NN1Ahr:

###############################################fixed above.

###No name 5F4
NN2hl:=subs(d=n,[[a,a/2,b,c,d],[a/2,1+a-b,1+a-c,1+a-d],1]):
NN2hr:=subs(d=n,rf(1,a-b)*rf(1,a-c)*rf(1,a-d)*rf(1,a-b-c-d)/(rf(1,a)*rf(1,a-b-c)*rf(1,a-b-d)*rf(1,a-c-d))):
NN2z:=subs(d=n,(rf(a,k)*rf(1+a/2,k)*rf(b,k)*rf(c,k)*rf(d,k)/(k!*rf(a/2,k)*rf(1+a-b,k)*rf(1+a-c,k)*rf(1+a-d,k)))/NN2hr):


###No name 5F4 (A) (similar to NN2 where d=-n)
NN2Ahl:= [[a,1+a/2,b,c,-n],[a/2,1+a-b,1+a-c,1+a+n],1]:
NN2Ahr:=rf(1+a,n)*rf(1+a-b-c,n)/(rf(1+a-b,n)*rf(1+a-c,n)):
NN2Az:=(rf(a,k)*rf(1+a/2,k)*rf(b,k)*rf(c,k)*rf(-n,k)/(k!*rf(a/2,k)*rf(1+a-b,k)*rf(1+a-c,k)*rf(1+a+n,k)))/NN2Ahr:


###Dougall's Theorem (7F6)
DOUGhl:=[[a,a/2,b,c,d,1+2*a-b-c-d+n,-n],[a/2,1+a-b,1+a-c,1+a-d,-a+b+c+d-n,1+a+n],1]:
DOUGhr:=rf(1+a,n)*rf(1+a-b-c,n)*rf(1+a-b-d,n)*rf(1+a-c-d,n)/(rf(1+a-b,n)*rf(1+a-c,n)*rf(1+a-d,n)*rf(1+a-b-c-d,n)):
DOUGz:=(rf(a,k)*rf(a,2,k)*rf(b,k)*rf(c,k)*rf(d,k)*rf(1+2*a-b-c-d+n,k)*rf(-n,k)/(k!*rf(a/2,k)*rf(1+a-b,k)*rf(1+a-c,k)*
rf(1+a-d,k)*rf(-a+b+c+d-n,k)*rf(1+a+n,k)))/Dhr:


###Nearly-poised summation t-m
NP1hl:=[[a,1+a/2,-n],[a/2,b],1]:
NP1hr:=(b-a-1-n)*rf(b-a,n-1)/rf(b,n):
NP1z:=(rf(a,k)*rf(1+a/2,k)*rf(-n,k)/(k!*rf(a/2,k)*rf(b,k)))/NP1hr: