module zeilberg; % An implementation of the Gosper and Zeilberger
                 % algorithms.
% zeilberg.red
% authors: Gregor Stoelting & Wolfram Koepf
% version 1.2, April 1995
% Reduce version 3.6
%
%
% References:
%
% R. W. Gosper, Jr.:
% Decision procedure for indefinite hypergeometric summation,
% Proc. Nat. Acad. Sci. USA 75 (1978), 40-42.
%
% Koornwinder, T. H.:
% On Zeilberger's algorithm and
% its q-analogue: a rigorous description.
% J. of Comput. and Appl. Math. 48 (1993), 91-111.
%
% Zeilberger, D.:
% A fast algorithm for proving terminating hypergeometric identities,
% Discrete Math. 80 (1990), 207-211.
%
% Koepf, W.:
% Algorithms for the indefinite and definite summation.
% Konrad-Zuse-Zentrum Berlin (ZIB), Preprint SC 94-33, 1994.
%
%
fluid '(zb_version);
zb_version:="package zeilberg, version 1.1, Feb. 15, 1995"$
global '(inconsistent!*);

create!-package('(zeilberg),'(contrib misc));

algebraic;
share zb_order$
zb_order:=5 $
gosper_representation:=nil;
zeilberger_representation:=nil;

operator gamma,pochhammer,binomial;
operator hypergeom;
operator summ,zb_f,zb_sigma;
operator local_gamma,local_prod;
gamma1!*rules:={
gamma(~k)=> factorial(k-1)
 when fixp(k) and k>0
};
let gamma1!*rules;

pochhammer!*rules:={
pochhammer(~z,~k)  =>  ( for i:=0:(k-1) product(z + i))
   when fixp k and k < 20 and k > 0,
pochhammer(~z,~k)  => factorial(z+k-1)/factorial(z-1)
   when fixp z and z > 0
};

let pochhammer!*rules;

onerules:=
{gamma(~zb_x)=>1,
 binomial(~zb_x,~zb_y)=> 1,
 factorial(~zb_x)=> 1,
 pochhammer(~zb_x,~zb_y)=> 1
};

onerules2:={summ(~zb_x)=>1,hypergeom(~zb_x1,~zb_x2,~zb_x3)=>1};

gammatofactorial:={gamma(~a) => factorial(a-1)};

zb_binomialrules:={binomial(~n,0)=>1,
binomial(~n,~n)=>1,
binomial(~n,~k)=>0 when (fixp k and k<0),
binomial(~n,~k)=>0 when (fixp (n-k) and (n-k)<0),
binomial(~n,~k)=>factorial(n)/(factorial(k)*factorial(n-k)) when
(fixp k and fixp n)
};

let zb_binomialrules;




switch zb_factor, zb_timer,zb_proof, zb_trace,zb_inhomogeneous;
lisp setq(!*zb_factor,t);
% zb_factor:=1;
zb_direction:=down;

symbolic  procedure gosper!-eval u;
% gosper(f,k) searches for a hpergeometric term = closed form
% antidifference.
% form solution does not exist.
% gosper(f,k,m,n) determines
%
%  __n___
%  \
%   \   f(k)
%   /
%  /
%  -----
%   k=m
%
% using Gosper's algorithm. This is only successful if
% Gosper's algorithm applies.
(<< abc:=
<<
if length u = 2 then gosper0(aeval car u, aeval cadr u)
else
if length u = 3 then gosperborders(aeval car u, aeval cadr u,
0,caddr u)
else if  length u = 4 then gosperborders(aeval car u, aeval cadr u,
  aeval caddr u, aeval cadddr u)
else rederr("illegal number of arguments")>>;
if eqcar (abc,'!*sq) then
      list('!*sq,cadr abc,nil) else abc>>) where abc=nil;
put ('gosper, 'psopfn,'gosper!-eval);




algebraic procedure gosperborders(func,k,k0,k1);
% gosperborders(func,k,k0,k1) = gosper(func,k,k0,k1)
begin
scalar tmp,gosper2,!*factor;
gosper2:=gosper0(func,k);
tmp:=sub(k=k1,gosper2)-sub(k=k0-1,gosper2);
if lisp !*zb_factor then
 <<
  on factor;
  return num(tmp)/den(tmp)
 >>
else
 return tmp;
end;


algebraic procedure gosper0(func,k);
begin
scalar tmp;
tmp:=gosper1(func,k);
if tmp = zb_gancfse then
  rederr("Gosper algorithm: no closed form solution exists")
else
 return (tmp);
end;

algebraic procedure gosper1(func,k);
% gosper1(func,k) = gosper(func,k)
begin
scalar dexp,gexp,d,g,dg,degree, downmax,facj,partj, j,jj, p,r1,r2,
polynomials,sol,!*exp,!*factor,equations,equationlist,f,varlist,s,l;
clear(gosper_representation,
 zeilberger_representation,rational_certificate);
on exp;
%tester:=func;
%tester:=(tester where onerules);
%if not(type_ratpoly(tester,k)) then
%<<
% tester:=tester/sub(k=k-1,tester);
% if not(type_ratpoly(tester,k))  then
%  rederr("tester: algorithm not applicable")
%>>;
if polynomqq(func,k) then
  <<
  if lisp !*zb_trace then write "Gosper algorithm applicable";
  polynomials:={func,1,1}
  >>
else
 <<
 on factor;
 off exp;
  dg:=simplify_combinatorial(func/sub(k=k-1,func));
  %on exp;
  if dg = 0 then return 0;
  d:=num(dg);
  g:=den(dg);
 %Dexp:=D;
 %Gexp:=G;
 if lisp !*zb_trace then
   write "a(",k,")/a(",k,"-1):=",dg;
 %off factor;
  %on exp;
 %if not ratpoly2(D,G,k) then
 %rederr("Gosper algorithm not applicable");
 if not (polynomq4(d,k) and polynomq4(g,k)) then
  rederr("Gosper algorithm not applicable");
 % Gexp:=NIL;
 % Dexp:=NIL;
  if lisp !*zb_trace then write "Gosper algorithm applicable";
 %off exp;
 %on factor;
 % write "D:=",D;
 %write "G:=",G;
  polynomials:=determine!_polynomials2(d,g,k);
  d:=nil; g:=nil;
  %on exp;
 >>;
if lisp !*zb_timer then
 << write "flag: determine polynomials"; showtime>>;
p:=first(polynomials);
r1:=second(polynomials);
r2:=third(polynomials);
if lisp !*zb_trace then
  << write  "p:=",p; write  "q:=",r1; write  "r:=",r2>>;
off factor;
on exp;
% compute the maximum degree for the polynomial f.
% (Lemma 3.6 in Koornwinder)
degree:=maxdegf(r1,r2,p,k);
if lisp !*zb_trace then write  "degreebound:=",degree;
if lisp !*zb_timer then
 << write "flag: maxdegf"; showtime>>;
if (degree< 0) then
 return (zb_gancfse);
f:=(for j:=0 :degree sum (zb_f(j) * k^j));
equations:=(sub(k=k+1,r1) * f - r2 *  sub(k=k-1,f) - p);
on exp;
equationlist:=coeff(equations,k);
varlist:=(for j:=0 : degree collect (zb_f(j)));
l:=arglength(equationlist);
downmax:=max(degree -l +1 ,0);
sol:={};
sol2:={};
for j:=degree step -1 until downmax   do
<<
 off factor;
 partj:=sub(sol,part(equationlist,l -degree  +j));
 on exp;
 jj:=degree;
 while freeof(partj,zb_f(jj)) and jj geq 0 do
  jj:=jj-1;
 facj:=coeff(partj,zb_f(jj));
 off exp;
 on factor;
 if arglength(facj) = 2 then
 <<
  solj:={zb_f(jj)= -part(facj,1)/part(facj,2)};
  off factor;
  sol:=append(solj,sub(solj,sol));
  sol2:=append({zb_f(jj)},sol2)
 >>;
>>;
 if arglength(sol) = degree  then
    <<
     tmp:=t; j:=0;
     while tmp and j leq degree do
      <<
       if freeof (sol2,zb_f(j)) then
         <<
          sol:=append({zb_f(j) = 0},sub(zb_f(j) = 0,sol));
          tmp:=nil
         >>;
      j:=j +1;
     >>
    >>;
if lisp !*zb_timer then
 << write "flag: sol";showtime>>;
if sub(sol,equations) neq 0 then
 <<
  if lisp !*zb_proof then gosper_representation:={p,r1,r2,nil};
  return(zb_gancfse)
 >>;
f:=sub(sol,f);
if lisp !*zb_proof then gosper_representation:={p,r1,r2,f};
if lisp !*zb_proof then
 if zb_direction = down  then
  <<
   rational_certificate:=sub(k=k+1,r1)/p*f;
   s:=rational_certificate*func
  >>
 else
  <<
   rational_certificate:=sub(k=k+1,r2/p)*f;
   s:=rational_certificate*sub(k=k+1,func)
  >>
else
 s:=sub(k=k+1,r1)/p*f*func;
if lisp !*zb_trace then write  "f:=",f;
off factor;
if lisp !*zb_timer then
<< write "flag: simplify comb von sol:=";showtime>>;
if lisp !*zb_factor then
 <<
  on factor;
  s:=num(s)/den(s);
 >>;
if lisp !*zb_timer then
<< write "flag: num(s)/den(s) under factor";showtime>>;
if lisp !*zb_trace then write  "Gosper algorithm successful";
if zb_direction = down  then
 return s
else
 return sub(k=k-1,s);
end; %gosper1


symbolic procedure sumrecursion!-eval u;
% sumrecursion(f,k,n,j) determines a holonomic recurrence equation
% of order less or equal j for
%
%  _____
%  \
%   \   f(n,k)   with respect to n
%   /
%  /
%  -----
%   k
% sumrecursion(f,k,n) = sumrecursion(f,k,n,zb_order), where zb_order
% is a global variable (default = 5)
<<
 if length u = 3 then
  sumrecursion0(aeval car u,aeval cadr u,aeval caddr u,1,zb_order)
 else
  if  length u = 4 then
   sumrecursion0(aeval car u, aeval cadr u, aeval caddr u,
    aeval cadddr u, aeval cadddr u)
  else rederr("illegal number of arguments")
>>;

put ('sumrecursion, 'psopfn,'sumrecursion!-eval);

algebraic procedure sumrecursion0(func,secundus,n,mini, maxi);
begin
scalar !*factor,!*exp;
if lisp !*zb_factor then on factor;
return part(sumrecursion1(func,secundus,n,mini, maxi),1);
end;

algebraic procedure sumrecursion1(func,secundus,n,mini, maxi);
begin
scalar result1,ank,b,c,d,g,bc,dg, j,jj, inhomogeneous,k,aa,bb,
!*factor,!*exp,order1;
clear(gosper_representation,
 zeilberger_representation,rational_certificate);
result1:=-1;
on exp;
inhomogeneous:=nil;
% provisorisch
if (arglength(secundus )= 3) and (part(secundus,0)=list) then
 rederr("not yet implemented.");
if (arglength(secundus )= 3) and (part(secundus,0)=list) then
  <<
   k:=part(secundus,1);
   aa:=part(secundus,2);
   bb:=part(secundus,3);
   inhomogeneous:=t
  >>
 else
  <<
   k:=secundus;
   aa:=0;
   bb:=0
  >>;
ank:=func;
tester:=func;
%Write("tester:=",tester);
tester:=(tester where onerules);
%Write("tester:=",tester);
if tester neq 0  then
<<
 tester2:=tester;
 tester:=tester/sub(k=k-1,tester);
 %Write("tester:=",tester);
 if not(type_ratpoly(tester,k))  then
  rederr("algorithm not applicable");
 %Write("tester2:=",tester);
 tester2:=tester2/sub(n=n-1,tester2);
 %Write("tester2:=",tester);
 if not(type_ratpoly(tester2,n))  then
  rederr("algorithm not applicable");
>>;
bc:=simplify_combinatorial(ank/sub(n=n-1,ank));
b:=num(bc);
c:=den(bc);
if lisp !*zb_trace then
  write "F(",n,",",k,")/F(",n,"-1,",k,"):=",bc;
on exp;
if not type_ratpoly(bc,n) then
 <<
 if lisp !*zb_trace then
   << write "not rational";
      write "Zeilberger algorithm not applicable";
   >>;
 return extended_sumrecursion1(ank,k,n);
 >>;
dg:=simplify_combinatorial(ank/sub(k=k-1,ank));
d:=num(dg);
g:=den(dg);
if lisp !*zb_trace then
  write "F(",n,",",k,")/F(",n,",",k,"-1):=",dg;
on exp ; % achtung
if not type_ratpoly(dg,k) then
 <<
  if lisp !*zb_trace then
   << write "not rational";
      write "Zeilberger algorithm not applicable";
   >>;
return extended_sumrecursion1(ank,k,n);
 >>;
if lisp !*zb_trace then write "Zeilberger algorithm applicable";
on factor;
order1:=0;
for j:=mini : maxi  do
 if result1 = -1 then
   result1:=sumrecursion2(ank,b,c,d,g,k,n,aa,bb,inhomogeneous,j)
 else
   if order1= 0 then order1:=j -1;
   if result1 = -1 then
   rederr("Zeilberger algorithm fails. Enlarge zb_order");
off factor;
if lisp !*zb_factor then
<<
 on factor;
 if lisp minusp (lc numr cadr aeval prepsq cadr result1) then
   result1:=-result1;
>>;
return {result1,order1};
end;
% sumrecursion1

algebraic procedure
   sumrecursion2(ank,b,c,d,g,k,n,aa,bb,inhomogeneous,order1);
% applies Zeilberger algorithm for order order1
begin
scalar j,jj, p,r10,r20,r1,r2,p1,polynomials, gg, recursion,recursion2,
equations,inhomogeneous,f,varlist,r12,summe,s,k,z,!*factor,!*exp;
if lisp !*zb_factor then
on factor;
if lisp !*zb_trace then
  write  "applying Zeilberger algorithm for order:=",order1;
p0:=
 (for j:=0 : (order1-1) product sub(n=n-j,b) )+
 (for j:=1 : order1 sum (zb_sigma(j) *
 (for jj:=0 : (j-1) product sub(n=n-jj,c)) *
 (for jj:=j : (order1-1) product sub(n=n-jj,b)) ));
r12:=
 d * (for j:=0 : (order1-1) product sub({n=n-j,k=k-1},b))/
 (g * (for j:=0 : (order1-1) product sub(n=n-j,b)));
r12:=simplify_combinatorial(r12);
r10:=num(r12);
r20:=den(r12);
off factor;
polynomials:=
  determine!_polynomials(r10,r20,k);
p1:=first(polynomials);
p:=p1 *p0;
r1:=second(polynomials);
r2:=third(polynomials);
if lisp !*zb_trace then
   <<write  "p:=",p; write  "q:=",r1; write  "r:=",r2>>;
% compute the maximum degree for the polynomial f.
% (Lemma 3.6 in Koornwinder)
degree:=maxdegf(r1,r2,p,k);
if lisp !*zb_trace then write  "degreebound:=",degree;
if (degree< 0) then
  return (-1);
f:=(for j:=0 :degree sum (zb_f(j) * k^j));
equations:=(sub(k=k+1,r1) * f - r2 *  sub(k=k-1,f) - p);
on exp;
equationlist:=coeff(equations,k);
va:=(for j:=0 : degree collect zb_f(j));
vb:=(for j:=1 : order1 collect zb_sigma(j));
varlist:=append(va,vb);
sol:=savesolve(equationlist,varlist);
if sub(sol,equations) neq 0 then return -1;
f:=sub(sol,f);
%if arglength(sol)<degree + order1 then
%write "warning: solution not unique";
f:=(f where arbcomplex(~x)=>0);
if lisp !*zb_trace then write  "f:=",f;
p:=sub(sol,p);
p:=(p where arbcomplex(~x)=>0);
if lisp !*zb_trace then write  "p:=",p;
if lisp !*zb_proof then zeilberger_representation:={p,r1,r2,f};
if zb_direction = down then
  <<if lisp !*zb_proof then rational_certificate:=sub(k=k+1,r1)/p *f>>
else
  <<if lisp !*zb_proof then rational_certificate:=sub(k=k+1,r2/p)*f>>;
va:=sub(sol,va);
vb:=sub(sol,vb);
n0:=order1-1;
for j:=1 : degree do
    n0:=max(testnonnegintroots(den(part(va,j)),n),n0);
%write "n0:=",n0;
%write "first va " , testnonnegintroots(den(part(va,1)),n);
%write "last va " , testnonnegintroots(den(part(va,degree)),n);
for j:=1 :order1 do
    n0:=max(testnonnegintroots(den(part(vb,j)),n),n0);
%write "n0:=",n0;
%write "first vb " , testnonnegintroots(den(part(vb,1)),n);
%write "last vb " , testnonnegintroots(den(part(vb,order1)),n);
n0:=max(testnonnegintroots(den(sub(k=n+1,q)),n),n0);
n0:=max(testnonnegintroots(num(sub(k=n+1,p)),n),n0);
%write "n0:=",n0;
if n0>=order1 then write "recursion valid for n>=",n0+1;
%zb_testnonnegintroots:=n0+1;
zb_testnonnegintroots:=n0-order1+1;
recursion:=summ(n) +
(for j:=1 : order1 sum part(vb,j)* summ(n-j));
 recursion:=num(recursion);
recursion:=(recursion where arbcomplex(~local_x) => 1);
if lisp !*zb_proof or inhomogeneous then
 <<
  gg:=
  f*sub(k=k+1,r2*ank/
    (p1*(for j:=0: order1-1 product sub(n=n-j ,b))));
  gg:=sub(sol,gg);
  proof:=gg;
  gg:=sub({k=k-1,n=n+1},gg);
  gg:=simplify_combinatorial(gg);
  if lisp !*zb_trace then write "G:=",gg;
 >>;
if inhomogeneous then
<<
 on factor;
 if lisp !*zb_inhomogeneous then
 <<
  recursion:=
  {recursion ,
  simplify_combinatorial(sub(k=bb+1,gg) - sub(k=aa,gg))};
  tempo:=
   simplify_combinatorial(sub(k=bb,gg) - sub(k=aa-1,gg));
 >>
 else
 <<
  recursion:=
   gg * sub(k=k+1,recursion) -sub(k=k+1,gg)*recursion;
  recursion:=simplify_combinatorial(recursion);
  recursion:=num(recursion)
 >>;
>>;
%if inhomogeneous
if lisp !*zb_trace then write  "Zeilberger algorithm successful";
if zb_direction = down then
 return recursion
else
 return sub(n=n+order1,recursion);
end; %sumrecursion2

algebraic procedure testnonnegintroots(term1,n);
begin
scalar  n0,l,j,n1;
n0:=-1;
term1:=factorize(term1);
l:=arglength(term1);
for j:=1:l do
 <<
 f:=part(term1,j);
 if  deg(f,n )= 1 then
  n1:=-part(coeff(f,n),1)/part(coeff(f,n),2);
   if fixp(n1) then n0:=max(n1,n0)
 >>;
%write "returning",n0;
return n0;
end;




symbolic  procedure hypersum!-eval u;
 <<
  if length u = 4 then
   hypersum1(aeval car u, aeval cadr u, aeval
       caddr u,aeval cadddr u,1,zb_order)
  else
   if  length u = 5 then hypersum1(aeval car u, aeval cadr u,
      aeval caddr u, aeval cadddr u,car cddddr u,car cddddr u)
   else rederr("illegal number of arguments")
 >>;

put ('hypersum, 'psopfn,'hypersum!-eval);

algebraic procedure recursion_to_closed_form(recursion,startl,n,m);
begin
scalar aj,recj,list1,p,q,tmp,j,nonhyp,order1,!*factor,!*exp;
on exp;
list1:={};
order1:= arglength(startl);
p:=part(coeff(recursion,summ(n)),2);
q:=part(coeff(recursion,summ(n-order1)),2);
nonhyp:=0;
if not(freeof(summ,p) and freeof(summ,q)) then
  <<
  nonhyp:=1;
  write "no hypergeometric solution found";
  return  recursion;
   >>;
for j:=1:order1 do
  <<
   aj:=part(startl,j);
   if aj=0 then
    list1:=append(list1,{0})
   else
    <<
     recj:=sub(n= n*order1 +j-1,p) * summ(n) +
        sub(n= n*order1 +j-1,q) *summ(n-1);
      tmp:=rectopoch(recj,n,1,m);
      tmp:={aj * sub(n=(n-j+1)/order1,tmp)};
      list1:=append(list1,tmp);
    >>;
   >>;%for
if order1 = 1 then
 return part(list1,1)
else
 return list1;
end;%recursion_to_closed_form



algebraic procedure hypersum1(upper,lower,z,n,mini, maxi);
begin
scalar tmp1,tmp,j,jj,aj,order1,recursion,term1,!*exp,startl;
off exp;
tmp:=hyperrecursion1(upper,lower,z,n,mini,maxi);
recursion:=part(tmp,1);
order1:=part(tmp,2);
%order1:=recorder(f,n);
if lisp !*zb_trace then
 write "recursion for underlying hypergeometric term:=",recursion;
startl:={1};
if order1 > 1 then
<<
 for j:=1: order1-1 do
 <<
  aj:=sub(n=j,(for jj:=0 :j sum hyperterm(upper,lower,z,jj)));
  aj:=simplify_combinatorial(aj);
  %write "aj:=",aj;
   startl:=append(startl,{aj});
 >>;
>>;
% write "startl in hyp1:=",startl;
return recursion_to_closed_form(recursion, startl,n,0);
end;%hypersum1



% sumtohyper(hyperterm({-a,b},{c},z,k),k);
% sumtohyper(hyperterm({-a,b},{c},z,k),k);
algebraic procedure summation(f,k,n);
begin
scalar l,localhypersum ,upper,lower,z,term,i,tmp,
startl, aj,piecewiseterm,piecewiseseq,f1,partj,
recursion,counter,m,tmpterm,upper,lower,z,prefactor,init,ht,
initial,initialnumber,summand,j,gammasummand,!*exp;
on exp;
ht:=sumtohyper(f,k);
prefactor:=(ht where onerules2);
%write "prefactor:=",prefactor;
ht:=ht/prefactor;
upper:=part(ht,1);
lower:=part(ht,2);
z:=part( ht,3);
f1:=simplify_combinatorial(f);
tmp:=sumrecursion1(f1,k,n,1,zb_order);
%write("zb_testnonnegintroots:=",zb_testnonnegintroots);
recursion:=part(tmp,1);
order1:=part(tmp,2);
if (order1 = 1) and (zb_testnonnegintroots= 0)  then
<<
 return recursion_to_closed_form(recursion,{prefactor},n,0);
>>;
% evaluate upper border
l:=arglength(upper);
initialnumber:=0;
%write "UPPER:=",upper;
for j:=1:l do
 <<
 partj:=part(upper,j);
 %write "partj :=",partj;
 tmp:=coeff(partj,n);
 if arglength(tmp)=2 then
  if fixp(part(tmp,2)) and part(tmp,2)<0 and fixp(part(tmp,1)) then
    initialnumber:=-partj;
 >>;
if initialnumber = 0 then
  rederr("no reccurent evaluation possible");
startl:={};
 for j:=zb_testnonnegintroots: order1-1+ zb_testnonnegintroots do
 <<
  write "prefactor:=",prefactor;
  write "sum(hyperterm(UPPER,LOWER,z,k),k,0,initialnumber):=",
   sum(hyperterm(upper,lower,z,k),k,0,initialnumber);
  aj:=sub(n=j,prefactor*
   sum(hyperterm(upper,lower,z,k),k,0,initialnumber));
  aj:=simplify_combinatorial(aj);
   startl:=append(startl,{aj});
 >>;
write "startl:=",startl;
   term:=recursion_to_closed_form(recursion,startl,n,
            zb_testnonnegintroots);
write "term:=",term;
if freeof(term,summ) then
  return(term)
else
 if freeof(prefactor,n) then
  recursion:=term
 else
  recursion:=sumrecursion(f,k,n);
if lisp !*zb_trace then
 %write "recursion:=",recursion;
counter:=0;
l:=arglength(recursion);
for i:=1 : l do
 <<
  term:=
   part(recursion,i)/(part(recursion,i) where onerules2);
  term:=part(term,1);
   m:=part(term,1);
    counter:=max(counter,m-term);
 >>;
%initial values, depend on testnonnegintroots
if lisp !*zb_trace then
 write "calculating initial values";
 initialnumber:=0;
 l:=arglength(upper);
 for i:=1 : l do
  <<
  tmp:=part(coeff(part(i,upper),n),2);
  if fixp(tmp) and (tmp <0) then
     initialnumber:=part(upper,i);
  % still to implement: rational case
 if initialnumber=0 then
   rederr("no initialization found");
 if zb_testnonnegintroots=0 then
  errorset;
>>;
tmp:=sub(n=0,prefactor);
end; %summation


%in "zeilberger.red"$ hypersum({-n,-n},{1},-1,n);
%hypersum({-n,n+3*a,a},{3*a/2,(3*a+1)/2},3/4,n);
%hyperrecursion({-n,-n},{1},-1,n);
%hypersum({-2 *n ,-2 *n},{1},-1,n);
%sub(n=n/2,hypersum({-2 *n ,-2 *n},{1},-1,n));
%sumtohyper((-1)^k*binomial(n,k)^2,k);
% boolsche Variable
% Polynomgeschichten revisted
% w. schickt
% t durch 1 erstezen.



algebraic procedure recorder(f,n);
begin
pa:=patternarguments(f,summ,{});
pa:=sub(n=0,pa);
return (-min(pa));
end;

algebraic procedure rectopoch(f,n,order1,m);
begin
scalar dennum,denden,cases1,k,nume,deno,!*exp,!*gcd;
on exp; on gcd;
%write "order1:=",order1;
%order1:=recorder(f,n);
%write "f:=",f;
deno:=-part(coeff(f, summ(n)),2);
%if freeof(f, summ(n-order1)) then rederr("not yet implemented");
nume:=part(coeff(f, summ(n-order1)),2);
if order1 >1 then
 <<
 for j:=1 : order1 -1 do
  if not freeof(f,summ(n-j)) then rederr("no hypergeometric solution");
 cases1:={};
 for j:=0 :order1-1 do
 cases1:=append({sub(n=(n-j)/order1,
  rectopoch(summ(n) * sub(n=order1*n+j,nume) +
  summ(n-1) * sub(n=order1*n+j,deno),n,1,m))},cases1);
 return cases1
 >>;
%if not freeof(deno,summ) or not freeof(nume,summ) then
% rederr("no hypergeometric solution");
lcr2:=first(reverse(coeff(deno,n)));
lcr1:=first(reverse(coeff(nume,n)));
nume:=nume/lcr1;
dennum:=den(nume);
nume:=num(nume);
deno:=deno/lcr2;
denden:=den(deno);
deno:=num(deno);
deno:=factorize(deno);
nume:=factorize(nume);
deno:=(part(deno,1):=part(deno,1)/denden);
nume:=(part(nume,1):=part(nume,1)/dennum);
deno:=refactors(deno,n);
nume:=append({1},refactors(nume,n));
tmp:={};l:=arglength(nume);
for j:=1:l do
 tmp:=append(tmp, {part(nume,j)+m});
nume:=tmp;
tmp:={};l:=arglength(deno);
for j:=1:l do
 tmp:=append(tmp, {part(deno,j)+m});
deno:=tmp;
%write "deno:=",deno;
%write "nume:=",nume;
return hyperterm(nume,deno,lcr1/lcr2,n-m)*factorial(n-m)/
pochhammer(m+1,n-m);
end;
%extended_sumrecursion((pochhammer( - n,k)* pochhammer(b,k)*
%pochhammer(c,k))/(factorial(k)*pochhammer((b - n + 1)/2,k)*
%pochhammer(2*c,k)), k, n);
%hypersum({-n,b,c},{1/2*(1-n+b),2*c},1,n);
% hypersum({a,b},{c},1,a);
%hypersum({-n,b},{c},1,n);
%zeilb([-n,b],[c],1,n,1);




algebraic procedure refactors(term1,n);
begin
scalar a, l,i,c,d,g,pol,degree ,!*exp, !*factor,
denpol,numpol;denpol;
on exp;
g:={};
l:=arglength(term1);
%p1:=part(term1,1);
for i:=1:l do
<<
   pol:=part(term1,i);
   on exp;
  if not freeof(pol,n) then
   <<
   numpol:=num(pol);
   denpol:=den(pol);
   degree:=deg(numpol,n);
   if degree=1 then
     <<
      d:=part(coeff(numpol,n),2)/denpol;
      c:=part(coeff(numpol,n),1)/d/denpol+1;
      <<for j:=1 :degree do g:=append(g,{c})>>;
       >>
    else
      newrederr{pol," does not factorize."}
   >>
>>;
return g;
end;

symbolic  procedure hyperrecursion!-eval u;
% hyperrecursion({a_1,...,a_p},{b_1,...,b_q},x,n,j) determines
% a holonomic recurrence equation (up to order j)
% with respect to n for the peneralized hypergeometric function
%  F (a_1,...,a_p;b_1,...,b_q;x)
% hyperrecursion(upper,lower,x,n) =
% hyperrecursion(upper,lower,x,n,zb_order)
% where zb_order is a global variable (default = 5)
 <<
  if length u = 4 then
   hyperrecursion0(aeval car u, aeval cadr u, aeval
       caddr u,aeval cadddr u,1,zb_order)
  else
   if  length u = 5 then hyperrecursion0(aeval car u, aeval cadr u,
      aeval caddr u, aeval cadddr u,car cddddr u,car cddddr u)
   else rederr("illegal number of arguments")
 >>;

put ('hyperrecursion, 'psopfn,'hyperrecursion!-eval);


algebraic procedure hyperrecursion0(upper,lower,z,n,mini, maxi);
begin
scalar !*factor,!*exp;
if lisp !*zb_factor then on factor;
return part(hyperrecursion1(upper,lower,z,n,mini, maxi),1)
end;


algebraic procedure hyperrecursion1(upper,lower,z,n,mini, maxi);
begin
scalar tester,result1,b,c,d,g,bc,dg, upl,lol,func,j,goon,x,
liste,!*factor,!*exp,order1;
clear(gosper_representation,
zeilberger_representation,rational_certificate);
result1:=-1;
upl:=arglength(upper);
lol:=arglength(lower);
%% test if some upper index is a nonnegative integer
%goon:=t;
%liste:=upper;
%while goon do
% <<
%  if liste = {} then
%   goon:=NIL
%  else
%  <<
%   x:=first(liste);
%   if fixp(x) and x>-1 then
%     goon:=NIL
%   else
%     liste:=rest(liste)
%  >>
% >>;
%if arglength(liste)>0 then
% rederr ("some upper index is a nonnegative integer");
goon:=t;
liste:=lower;
while goon do
 <<
  if liste = {} then
   goon:=nil
  else
  <<
   x:=first(liste);
   if fixp(x) and x<1 then
     goon:=nil
   else
     liste:=rest(liste)
  >>
 >>;
if arglength(liste)>0 then
  rederr ("some lower index is a nonpositive integer");
func:=hyperterm(upper,lower,z,local_k);
tester:=func;
%Write("tester:=",tester);
tester:=(tester where onerules);
%Write("tester:=",tester);
if tester neq 0  then
<<
 tester:=tester/sub(n=n-1,tester);
 %Write("tester:=",tester);
 if not(type_ratpoly(tester,n))  then
  rederr("algorithm not applicable");
>>;
bc:=simplify_combinatorial(func/sub(n=n-1,func));
on factor;
b:=num(bc);
c:=den(bc);
if lisp !*zb_trace then
  write "F(",n,",local_k)/F(",n,"-1,local_k):=",b/c;
%off factor;
on exp;
if not type_ratpoly(bc,n) then
 <<
 if lisp !*zb_trace then
  <<
   write "not rational";
   write "Zeilberger algorithm not applicable"
  >>;
 return extended_hyperrecursion1(upper,lower,z,n);
 >>;
dg:=(for j:=1 : upl product(local_k - 1 + part(upper,j))) * z/
  ((for j:=1 :lol product(local_k - 1 + part(lower,j)))*(local_k));
d:=num(dg);
g:=den(dg);
if lisp !*zb_trace then
 <<
  write "F(",n,",local_k)/F(",n,",local_k-1):=",d/g;
  write "Zeilberger algorithm applicable"
 >>;
order1:=0;
for j:=mini:maxi  do
 if result1 = -1 then
  result1:=sumrecursion2(func,b,c,d,g,local_k,n,0,0,nil,j)
 else
  if order1= 0 then order1:=j -1;
if result1 = -1 then
 rederr("Zeilberger algorithm fails. Enlarge zb_order");
if lisp !*zb_factor then
<<
 on factor;
 if lisp minusp (lc numr cadr aeval prepsq cadr result1) then
   result1:=-result1;
>>;
return {result1,order1};
end;
% hyperrecursion1

algebraic procedure determine!_polynomials(r10,r20,k);
% determines polynomials p(k),r1(k),r2(k) as in Lemma 3.1 in
% Koornwinder,  or p_k,r_k,q_k as in Gosper,
% respectively.
begin
scalar tmp,r1divr2,p,r1,r2,j,jj, gamma1,!*exp,!*factor;
on exp;
%write "enter maxshift with ",{r10,r20};
%globalns:={r10,r20};
maxshift1:=maxshift(r10,r20,k);
%write "maxshift:=",maxshift1;
p:=1;
r1:=r10;
r2:=r20;
for jj:=0: maxshift1 do
<<
 %write "jj:=",jj;
 gamma1:=gcd(r1,sub(k= k+jj,r2));
 %write "jj:=",jj;
 if gamma1 neq 1 then
  <<
    r1:=r1/gamma1;
    r2:=r2/sub(k= k-jj,gamma1);
    p:=p * (for j:=0 : (jj-1) product sub(k=k-j,gamma1))
  >>; % if
>>;
return {p,r1,r2};
end;

algebraic procedure determine!_polynomials2(r10,r20,k);
% determines polynomials r1(k),r2(k),p(k) as in Lemma 3.1 in
% Koornwinder
begin
scalar
!*exp,!*factor,
f1,f2,order1,order2,ma,leadj,leadjj,jj,j,r1,r2,p;
on factor;off exp;
p:=1;
r1:=r10;
r2:=r20;
f1:=factorize(r1);
f2:=factorize(r2);
order1:=arglength(f1);
order2:=arglength(f2);
for j:=1 : order1 do
 for jj:=1 : order2 do
  <<
  complist:=comppol(part(f1,j),part(f2,jj),k);
  comp:=part(complist,1);
  leadj:=part(complist,2);
  leadjj:=part(complist,3);
  if comp> -1 then
  <<
   gamma1:=part(f1,j);
   gamma2:=part(f2,jj);
  %if  gamma1 neq  sub(k=k+j,gamma2) then
    r1:=r1/gamma1;
    r2:=r2/sub(k= k-comp,gamma1);
    p:=p * (for jj:=0 : (comp-1) product sub(k=k-jj,gamma1));
    f1:=(part(f1,j):=1);
    f2:=(part(f2,jj):=1);
 % neu
  >> % if
 >>;
on exp;
return {p,r1,r2};
end; % determine!_polynomials2

algebraic procedure maxshift(p1,p2,k);
% computes the maximal j with
% gcd(p1(k),p2(k+j)) neq 1
begin
scalar f1,f2,order1,order2,ma,j,jj;
ma:=-1;
f1:=factorize(p1);
f2:=factorize(p2);
order1:=arglength(f1);
order2:=arglength(f2);
for j:=1 : order1 do
 for jj:=1 : order2 do
  ma:=max(ma,comppol(part(f1,j),part(f2,jj),k));
return ma;
end;

algebraic procedure maxdegf(r1,r2,p,k);
% evalutes an upper bound for the degree of f
% with respect to variable k
% (Lemma 3.6 in Koornwinder)
begin
scalar l,dp, hold,hold2,!*exp,!*factor;
on exp;
pminus:=sub(k= k+1,r1) - r2;
pplus:=sub(k= k+1,r1) + r2;
lplus:=deg( pplus,k);
lminus:=deg( pminus,k);
if pminus=0 then lminus:=-1;
dp:=deg(p,k);
if (lplus leq lminus) then
 return max(dp - lminus,0)
else
 <<
 el:=
 part(coeff(pplus,k),lplus+1);
>>;
 if arglength(coeff(pminus,k))<lplus then
   dlminus1:=0
 else
   dlminus1:=part(coeff(pminus,k),lplus);
 hold:=-2 * dlminus1/el;
 hold2:=dp -lplus +1;
 if fixp(hold) and (hold geq  0) then
 % case b2 in Koornwinder
   return max(hold,hold2)
 else
 % case b1 in Koornwinder
 return max(hold2,0);
;
end;

algebraic procedure comppol(f,g,k);
% Tests for polynomials f and g in k if f(k)=g(k+j)
% for some nonnegative
% integer j, while f and g are not constant in k
% if this is the case comppol returns that j and
% the leading coefficiants of the polynomials
begin
 scalar nn,a,b,c,d,j, !*exp;
on exp;
nn:=deg(f,k);
if (nn=0) or (nn neq deg(g,k)) then return {-1,1,1};
a:=part(coeff(f,k),nn+1);
b:=part(coeff(f,k),nn);
c:=part(coeff(g,k),nn+1);
d:=part(coeff(g,k),nn);
j:=(b*c-a*d)/nn/a/c;
if not fixp(j) then return {-1,1,1}
else if j<0 then return {-1,1,1};
if (c * f - a * sub(k=k+j,g) = 0) then return {j,a,c}
else return {-1,1,1};
end;

algebraic procedure hyperterm(upper,lower,z,k);
% converts the representation of a hypergeomeric term
begin
scalar lu,ll;
lu:=arglength(upper);
ll:=arglength(lower);
return ((for j:=1:lu product(pochhammer(part(upper,j),k)))*z^k/
((for j:=1:ll product(pochhammer(part(lower,j),k)))*factorial(k)));
end;

algebraic procedure simplify_combinatorial(term1);
% converts binomials, products, pochhammers, and
% factorials in term1 into gammas, and
% applies simplify_gamma to the modified term1
   simplify_gamma(togamma(term1));

% write *mode;

% in "zeilberg.red"$
% gosper(gamma(k+n+2)*n/((k+n+1)*gamma(k+2)*gamma(n+1)),k);

algebraic procedure togamma(term1);
% converts binomials, products, pochhammers, and
% factorials in term1 into gammas
begin
term1:=(term1 where prod(~term,~k,~m1,~m2)=>
   producttopochhammer(term,k,m1,m2));
term1:=(term1 where local_prod(~term,~k,~m1,~m2)=>
   prod(~term,~k,~m1,~m2));
term1:=(term1 where pochhammer(0,~k)=>0);
term1:=(term1 where pochhammer(~n,~k)=> gamma(~n+~k)/gamma(~n));
term1:=(term1 where binomial(~n,~k) =>
        factorial(~n)/(factorial(~n - ~k)*factorial(~k)));
term1:=(term1 where factorial(~k)=> gamma(~k+1));
return term1;
end;

%ratsimplify_gamma(gamma(n) *n);

algebraic procedure ratsimplify_gamma(term1);
begin
scalar !*exp,!*factor,deno,nume, ln,ld,dega,nuga,derest,nurest,
lnurest,lnuga,lderest,ldega,jj,j,sp,term2,tmp;
on factor;
deno:=den(term1);
nume:=num(term1);
nurest:={};nuga:={};
derest:={};dega:={};
% construct two lists
% dega with parts that are gamma terms
% and derest with the others.
if arglength(deno) >0 then
<<
 if not(part(deno,0)= times) then
  if not freeof(deno,gamma) then
   <<
   deno:=strip_power(deno);
   tmp:=part(deno,1);
   if not(part( tmp,0) = gamma) then
    return term1
   else
    dega:=deno
   >>
  else
   derest:=strip_power(deno)
 else
 <<
  ld:=arglength(deno);
  for j:=1: ld do
  <<
   sp:=strip_power(part(deno,j));
   tmp:=part(sp,1);
   if not freeof(tmp,gamma) and part(tmp,0) = gamma then
     dega:=append(dega,sp)
   else
     derest:=append(derest,sp);
  >>; %for
 >>; %else
>> % if
else
 derest:={deno};
%ende
if arglength(nume) >0 then
<<
 if not(part(nume,0)= times) then
  if not freeof(nume,gamma) then
   <<
   nume:=strip_power(nume);
   tmp:=part(nume,1);
   if not(part( tmp,0) = gamma) then
    return term1
   else
    nuga:=nume
   >>
  else
   nurest:=strip_power(nume)
 else
 <<
  ln:=arglength(nume);
  for j:=1: ln do
  <<
   sp:=strip_power(part(nume,j));
   tmp:=part(sp,1);
   if not freeof(tmp,gamma) and part(tmp,0) = gamma then
     nuga:=append(nuga,sp)
   else
     nurest:=append(nurest,sp);
  >>; %for
 >>; %else
>> % if
else
 nurest:={nume};
%ende
% dega with parts that are gamma terms
% and derest with the others.
ldega:=arglength(dega);
lderest:=arglength(derest);
lnuga:=arglength(nuga);
lnurest:=arglength(nurest);
if ldega>0 then
<<
 for j:=1 : ldega  do
 <<
  tmp:=part(dega ,j);
  tmp:=part(tmp,1);
  for jj:=1 : lderest do
   if (part(derest,jj) - tmp) = 0  then
    <<
     derest:=(part(derest,jj):=1);
     tmp:=tmp+1;
     dega:=(part(dega,j):=gamma(tmp));
    >>;
  for jj:=1 : lnurest do
   if (part(nurest,jj) - tmp) = -1  then
    <<
     nurest:=(part(nurest,jj):=1);
     tmp:=tmp-1;
     dega:=(part(dega,j):=gamma(tmp));
    >>;
 >> %for j
>>; %ldega>0
if lnuga>0 then
<<
 for j:=1 : lnuga  do
 <<
  tmp:=part(nuga ,j);
  tmp:=part (tmp,1);
  for jj:=1 : lnurest do
   if (part(nurest,jj) - tmp) = 0  then
    <<
     nurest:=(part(nurest,jj):=1);
     tmp:=tmp+1;
     nuga:=(part(nuga,j):=gamma(tmp));
    >>;
  for jj:=1 : lderest do
   if (part(derest,jj) - tmp) = -1  then
    <<
     derest:=(part(derest,jj):=1);
     tmp:=tmp-1;
     nuga:=(part(nuga,j):=gamma(tmp));
    >>;
 >>% for j;
>>; %lnuga>0
 term2:=1;
%if lnuga>0 then
for j:=1 : lnuga do term2:=term2 *part(nuga,j);
%if lnurest>0 then
for j:=1 : lnurest  do term2:=term2 * part(nurest,j);
%if ldega>0 then
for j:=1 : ldega do term2:=term2 /part(dega,j);
%if lderest>0 then
for j:=1 : lderest  do term2:=term2 / part(derest,j);
if term2 = term1 then
 return term2
else
 return ratsimplify_gamma(term2);
end; %ratsimplify_gamma

algebraic procedure strip_power(term1);
begin
scalar j,!*factor,list1;
on factor;
list1:={};
if (arglength(term1)<2) or
   (part(term1,0) neq expt) or
   not fixp(part(term1,2))  then
  return {term1}
else
  for j:=1: part(term1,2) do
   list1:=append(list1,{part(term1,1)});
return list1;
end;





% ratsimplify_gamma(gamma(n)/n);
% ratsimplify_gamma(gamma(n)/(n-1));
% ratsimplify_gamma(gamma(n)^2/(n-1)^2);
% ratsimplify_gamma(gamma(n)^2*n^2);
% ratsimplify_gamma((n+1) * gamma(n)^2*n^2);


algebraic procedure simplify_gamma(term1);
% converts all subexpressions
% gamma(xi) -> gamma(xi + m)/((xi)*(xi+1)*...* (xi+m-1))
% where m is the largest integer , so that a subexpression
% gamma(xj) of term1 exists with xj = xi + m.
%
begin
scalar !*exp,!*factor,!*gcd,high,highl,highlength,j;
%on gcd;
%on factor;
if freeof(term1,gamma) then return term1;
highl:={};
highl:=highest_gamma_order(term1,highl);
if lisp !*zb_timer then
 << write "flag:highl:=",highl, "at "; showtime>>;
if highl = {} then return term1;
%term1:=gammashift(term1,highl);
term1:=matchgammashift(term1,highl);
%globalterm3:=term1;
%globalterm1:=gammashift(term1,highl);
%highlength:=arglength(highl);
%for j:=1:highlength do
%<<
%high:=part(highl,j);
%term1:=gammashift(term1,high);
%term1:=(term1 where gamma(~local_x)=>shift_gamma(~local_x,high));
%term1:=(term1 where local_gamma(~local_x)=>gamma(~local_x));
%>>;
%globalterm2:=term1;
%on exp;
return term1;
end;

algebraic procedure matchgammashift(term1,highl);
begin
scalar deno,nume,!*factor;
%on factor;
nume:=num(term1);
deno:=den(term1);
nume:=(nume where gamma(~local_x)=>listshift_gamma(~local_x,highl));
nume:=(nume where local_gamma(~local_x)=>gamma(~local_x));
if nume=0 then return 0;
deno:=(deno where gamma(~local_x)=>listshift_gamma(~local_x,highl));
deno:=(deno where local_gamma(~local_x)=>gamma(~local_x));
return nume/deno;
end;


algebraic procedure highest_gamma_order(term1,highl);
% produces a list of maximal xi for which
% exist subexpressions gamma(xi) of term1, and
% xi-xj is no integer iff xi neq xj
begin
scalar jjj,jj,j,max, term1length,localhighl,localhighllength,new;
term1length:=arglength(term1);
if (term1length<1) or freeof(term1,gamma)  then return highl;
new:=1;
highllength:=arglength(highl);
if part(term1,0) = gamma then
 <<
  if term1length neq 1 then
   rederr("gamma has illegal number of arguments");
  for j:=1 : highllength do
   if fixp(part(highl,j) - part(term1,1)) then
     <<
      if (part(highl,j) - part(term1,1)<0) then
        highl:=(part(highl,j):=part(term1,1));
      new:=0;
     >>;
  if new = 1 then highl:=append(highl,{part(term1,1)});
 >>
else
 for j:=1 : term1length do
 <<
  localhighl:=highest_gamma_order(part(term1,j),{});
  localhighllength:=arglength(localhighl);
  for jjj:=1:localhighllength do
   <<
    highllength:=arglength(highl);
    new:=1;
    for jj:=1 :highllength do
      if fixp(part(highl,jj) - part(localhighl,jjj)) then
        <<
        if (part(highl,jj) - part(localhighl,jjj)<0)  then
         highl:=(part(highl,jj):=part(localhighl,jjj));
        new:=0
        >>;
    if new = 1 then highl:=append(highl,{part(localhighl,jjj)})
   >>;
 >>; % for j:=1 : term1length
% if new = 1
return highl;
end;

algebraic procedure gammashift(term1,highl);
begin
scalar lhighl,term2,xx,nminusxx, j,jj,n;
if freeof(term1,gamma) then
  return term1;
if (arglength(term1)>1) then
  return map(gammashift(~zbglobal,highl),term1);
if (part(term1,0) = gamma) then
<<
 lhighl:=arglength(highl);
 term2:=term1;jj:=1;
 while (term1=term2) and (jj leq lhighl) do
 <<
  xx:=part(term1,1);
  n:=part(highl,jj);
  nminusxx:=n-xx;
  if (nminusxx = 0) then term1:=0;
  if fixp(nminusxx) and (nminusxx neq 0) then
   if nminusxx>0 then
    term2:=(gamma(n) / (for j:=1: nminusxx product(n-j)))
   else
    term2:=gamma(n) * (for j:=1: -nminusxx product(xx-j));
  jj:=jj+1;
  >>;
 return term2
>>
else
 return map(gammashift(~zbglobal,highl),term1);
end;




algebraic procedure shift_gamma(xx,n);
% shifts gamma-expression if n - xx is an integer
% warning: returns operator local_gamma instead of gamma
begin
scalar nminusx,j;
nminusx:=n-xx;
if not fixp(nminusx) then return local_gamma(xx);
if nminusx>0 then return local_gamma(n) /
(for j:=1: nminusx product(n-j)) else
return local_gamma(n) * (for j:=1: -nminusx product(xx-j));
end;


algebraic procedure listshift_gamma(xx,highl);
begin
scalar lhighl,nminusx,j,n,ret;
lhighl:=arglength(highl);
ret:=local_gamma(xx);
for j:=1 :lhighl  do
<<
 n:=part(highl,j);
 nminusx:=n-xx;
 if fixp(nminusx) then
 <<
  if nminusx>0 then
   ret:=local_gamma(n) / (for j:=1: nminusx product(n-j))
  else
   ret:=local_gamma(n) * (for j:=1: -nminusx product(xx-j));
  % j:=highl
 >>
>>;
return ret;
end;

algebraic procedure producttopochhammer(term,k,m1,m2);
% converts products into pochhammers
begin
 scalar fehler,ar,co,aa,bb,liste,tlength,j,pa;
 fehler:=nil;
if (den(term) neq 1) then return
  producttopochhammer(num(term),k,m1,m2)/
producttopochhammer(den(term),k,m1,m2);
 liste:=factorize(term);
 %gets initialized with factors of term
 %during the procedure I exchange them with pochhammer terms
 tlength:=arglength(liste);
 for j:=1 : tlength do
 <<
  pa:=part(liste,j);
  co:=coeff(pa,k);
  ar:=arglength(co);
  if ar>2 then fehler:=t
  else if ar<2 then
     liste:=(part(liste,j):= pa^(m2-m1+1))
  else
   <<
    aa:=part(co,2);
    bb:=pa/aa -k;
    if bb = 0 then
     liste:=(part(liste,j):=pochhammer(m1+ part(co,1),m2-m1+1))
    else
     liste:=(part(liste,j):=
      aa^(m2-m1)*pochhammer(bb,m2+1)/pochhammer(bb,m1))
   >>
 >>;
 if fehler then return local_prod(term,k,m1,m2);
 return (for j:=1: tlength  product( part(liste,j)));
end;

% extended
% authors: Gregor Stoelting & Wolfram Koepf

symbolic  procedure extended_gosper!-eval u;
(<< abc:=
<<
if length u = 2 then extended_gosper1(aeval car u, aeval cadr u)
else  if length u = 3 then extended_gosper2(aeval car u, aeval cadr u,
  aeval caddr u)
else  if length u = 4
 then extended_gosperborders(aeval car u, aeval cadr u,
  aeval caddr u,aeval cadddr u)
else rederr("illegal number of arguments")>>;
if eqcar (abc,'!*sq) then
      list('!*sq,cadr abc,nil) else abc>>) where abc=nil;
put ('extended_gosper, 'psopfn,'extended_gosper!-eval);

algebraic procedure extended_gosperborders(term1,k,k0,k1);
begin
scalar tmp,gosper2,!*factor;
gosper2:=extended_gosper1(term1,k);
if zb_direction = up then gosper2:=sub(k=k+1,gosper2);
tmp:=sub(k=k1,gosper2)-sub(k=k0-1,gosper2);
if lisp !*zb_factor then
 <<
  on factor;
  return num(tmp)/den(tmp)
 >>
else
 return tmp;
end;% extended_gosperborders


algebraic procedure extended_gosper2(term1,k,m);
begin
scalar !*exp,!*factor,s,tmp;
tmp:=gosper1(sub(k=k*m,term1),k);
if tmp = zb_gancfse then
  newrederr {"extended Gosper algorithm (Koepf): no ",m,
  "-fold hypergeometric solution"};
s:=sub(k=k/m,tmp);
if lisp !*zb_factor then on factor;
return s;
end; %extended_gosper2

algebraic procedure extended_gosper1(term1,k);
begin
scalar sol,!*factor,j,l,partj,s,m,tmp;
if lisp !*zb_trace then
     write "Koepf extension of Gosper algorithm entered...";
list1:=argumentlist(term1,{});
if list1 = {} then return gosper0(term1,k);
list2:=foreach partj in list1 collect linearfactor(partj,k);
m:=lcml(list2);
if lisp !*zb_trace then
  write "linearizing integer with respect to ",k," is ",m;
s:=extended_gosper2(term1,k,m);
if lisp !*zb_trace then write "s(",k,"):=",s;
sol:=(for j:=0:m-1 sum(sub(k=k-j,s)));
%if m>1 then sol:=simplify_combinatorial(sol);
if zb_direction = up then
 sol:=sub(k=k+1,sol);
if lisp !*zb_factor then on factor;
return sol
end; %extended_gosper1


symbolic  procedure extended_sumrecursion!-eval u;
(<< abc:=
<<
if length u = 3 then
  extended_sumrecursion0(aeval car u, aeval cadr u,aeval caddr u)
else  if length u = 5
 then extended_sumrecursion20(aeval car u, aeval cadr u,
  aeval caddr u,aeval cadddr u,car cddddr u)
else rederr("illegal number of arguments")>>;
if eqcar (abc,'!*sq) then
      list('!*sq,cadr abc,nil) else abc>>) where abc=nil;
put ('extended_sumrecursion, 'psopfn,'extended_sumrecursion!-eval);

algebraic procedure extended_sumrecursion0(term1,k,n);
begin
scalar !*factor,!*exp;
if lisp !*zb_factor then on factor;
return part(extended_sumrecursion1(term1,k,n),1);
end;

%extended_hyperrecursion1({ - n,b,c},{(b - n + 1)/2,2*c},1,n);

algebraic procedure extended_sumrecursion1(term1,k,n);
begin
scalar m,j,l,partj,s,!*exp,dg,bc;
on exp;
if lisp !*zb_trace then
     write "Koepf extension of Zeilberger algorithm entered...";
list1:=argumentlist(term1,{});
if list1 = {} then return sumrecursion1(term1,k,n,1,zb_order);
listk:=foreach partj in list1 collect linearfactor(partj,k);
listn:=foreach partj in list1 collect linearfactor(partj,n);
l:=lcml(listk);
m:=lcml(listn);
if lisp !*zb_trace then
   <<
     write "linearizing integer with respect to ",k," is ",l;
     write "linearizing integer with respect to ",n," is ",m;
   >>;
if m=1 and l=1 then
<<
 bc:=simplify_combinatorial(term1/sub(n=n-1,term1));
 globalbc:=bc;
 if not type_ratpoly(bc,n) then
 <<
  if lisp !*zb_trace then
    write "F(",n,",local_k)/F(",n,"-1,local_k):=",bc;
  rederr("Zeilberger algorithm not applicable")
 >>;
 dg:=simplify_combinatorial(term1/sub(k=k-1,term1));
 on exp;
 if not type_ratpoly(dg,k) then
 <<
  if lisp !*zb_trace then
    write "F(",n,",",k,")/F(",n,",",k,"-1):=",dg;
  rederr("Zeilberger algorithm not applicable")
 >>;
 return(sumrecursion1(term1,k,n,1,zb_order))
>>;
return extended_sumrecursion2(term1,k,n,m,l);
end; %extended_sumrecursion1

algebraic procedure extended_sumrecursion20(term1,k,n,m,l);
begin
scalar !*factor,!*exp;
if lisp !*zb_factor then on factor;
return part(extended_sumrecursion2(term1,k,n,m,l),1);
end;

algebraic procedure extended_sumrecursion2(term1,k,n,m,l);
begin
scalar term2,tmpterm,rule,!*factor,!*exp,order1;
term2:=sub({k=k*l,n=n*m},term1);
if lisp !*zb_trace then
  write  "applying Zeilberger algorithm to F(",n,",",k,"):=",term2;
tmpterm:=sumrecursion1(term2,k,n,1,zb_order);
order1:=m* part(tmpterm,2);
tmpterm:=part(tmpterm,1);
tmpterm:=sub({n=n/m},tmpterm);
rule:={summ(~nn/~mm)=>summ(nn) when mm=m};
tmpterm:=num(tmpterm where rule);
off factor;
tmpterm:=tmpterm;
if lisp !*zb_factor then on factor;
return({tmpterm,order1})
end;

symbolic  procedure extended_hyperrecursion!-eval u;
(<< abc:=
<<
if length u = 4 then
  extended_hyperrecursion0(aeval car u, aeval cadr u,aeval caddr u,
   aeval cadddr u)
%else  if length u = 5
%        then extended_hyperrecursion2(aeval car u, aeval cadr u,
%  aeval caddr u,aeval cadddr u,car cddddr u)
else rederr("illegal number of arguments")>>;
if eqcar (abc,'!*sq) then
      list('!*sq,cadr abc,nil) else abc>>) where abc=nil;
put ('extended_hyperrecursion, 'psopfn,'extended_hyperrecursion!-eval);


algebraic procedure extended_hyperrecursion0(upper,lower,x,n);
   part(extended_hyperrecursion1(upper,lower,x,n),1);

algebraic procedure extended_hyperrecursion1(upper,lower,x,n);
  extended_sumrecursion1(hyperterm(upper,lower,x,local_k),local_k,n);


algebraic procedure linearfactor(term1,n);
begin
scalar p,co;
co:=coeff(term1,n);
if arglength(co) = 1 then return 1;
p:=den(part(co,2));
if arglength(co) > 2 or (not fixp(p)) then
   rederr("Extended Gosper algorithm not applicable");
return p;
end;

algebraic procedure lcml(list1);
begin
% finds least common multiple of a list of integers
scalar p1,l;
l:=arglength(list1);
p1:=part(list1,1);
if l = 1 then return p1;
if l = 2 then return lcm( p1, part(list1,2));
return lcm(p1,lcml(rest(list1)));
end;

algebraic procedure argumentlist(term1, list1);
begin
scalar head1,j,l;
l:=arglength(term1);
if l<1 then return list1;
head1:=part(term1,0);
if head1 = gamma or
   %head1 = expt or
   head1 = factorial then
   list1:=append(list1,{part(term1,1)})
else if
   head1 = pochhammer or
   head1 = binomial then
   list1:=append(list1,{part(term1,1),part(term1,2)})
else
 for j:=1:l do
  list1:=argumentlist(part(term1,j),list1);
return list1;
end;

operator hypergeometric;

%let {gamma(~n)=>factorial(n-1) when (fixp(n) and n>0)};


%sumtohyper(hyperterm({a,a,a,a,a,b},{c},x,k),k);


algebraic procedure negintoccurs(list1);
begin
scalar l,tmp,tmp2,j;
 tmp2:=nil;
 l:=arglength(list1);
 if l = 0 then return nil;
 for j:=1 : l do
 <<
  tmp:=part(list1,j)   ;
  if fixp(tmp) and tmp<0 then tmp2:=t
 >>;
 return tmp2 ;
end; % negintoccurs

algebraic procedure sumtohyper(ank,k);
begin
scalar de,rat,numerator,denominator,numfactors,denfactors,lc,l,numlist,
oldnumlist, olddenlist,tmp,tmp2,
numdegree,denfactors,denlist, dendegree,i,j,lcden,
lcnum,!*exp,!*factor,!*gcd, gcdterm;
on exp;on gcd;
ank:=simplify_combinatorial(ank);
de:=simplify_combinatorial(sub(k=k+1,ank)/ank);
 if lisp !*zb_trace then
   write "a(",k,"+1)/a(",k,"):=",de;
numerator:=num(de);
denominator:=den(de);
if not polynomq4(numerator,k) then
 rederr("cannot be converted into hypergeometric form");
if not polynomq4(denominator,k) then
 rederr("cannot be converted into hypergeometric form");
numerator:=numerator;
denominator:=denominator;
numfactors:=factorize(numerator);
denfactors:=factorize(denominator);
lcnum:=lcof(numerator,k);
lcden:=lcof(denominator,k);
if lcnum = 0 then lcnum:=numerator;
if lcden = 0 then lcden:=denominator;
lc:=lcnum/lcden;
if freeof(first(numfactors),k) then numfactors:=rest(numfactors);
numlist:={};
len:=length(numfactors);
for i:=1:len do
<<
 fir:=first(numfactors);
 if not freeof(fir,k) then
 <<
  new:=-part(first(solve(fir,k)),2);
  numlist:=append(numlist,{new});
 >>;
 numfactors:=rest(numfactors);
>>;
maxint:=maxposint(numlist);
len:=length(denfactors);
denlist:={};
for j:=1:len do
<<
 fir:=first(denfactors);
 if not freeof(fir,k) then
 <<
  if not polynomq4(fir,k) or deg(fir,k)>2  then
    rederr("not yet implemented")
  else
   tmp:=solve( fir,k);
  for jj:=1: arglength(tmp) do
   denlist:=append(denlist,{-part(part(tmp,jj),2)});
 >>;
 denfactors:=rest(denfactors);
>>;
minint:=minnegint(denlist);
if minint leq 0 then
<<
 if lisp !*zb_trace then
    write "shifting by ",1-minint;
 numlist:=sub(k= k+1-minint,numlist);
 if numberofzeros(numlist)>0 then
  rederr("not yet implemented")
>>
else
<<
 if maxint geq 0 then
 <<
 if lisp !*zb_trace then
    write "shifting by ",1-maxint;
 denlist:=sub(k= k+1-maxint,denlist);
 if numberofzeros(denlist)>0 then
  rederr("not yet implemented");
  minint:=maxint;
 >>
>>;
shiftnumber:=1-minint;
if lisp !*zb_trace then
    write "calculating initial value";
olddenlist:=denlist;
denlist:={};
for j:=1: arglength(olddenlist) do
 denlist:=append({part(olddenlist,j ) + 1-minint},denlist);
oldnumlist:=numlist;
numlist:={};
for j:=1: arglength(oldnumlist) do
 numlist:=append({part(oldnumlist,j ) + 1-minint},numlist);
if sub(k=1-minint,den(ank)) = 0 or sub(pochhammer= poch, den(ank)) = 0
  then tmp:=limit(ank,k,1-minint)
else
 tmp:=sub(k=1-minint,ank);
if  member(1,denlist) then
<<
 tmplist:={};
 done:=0;
 for i:=1 : arglength(denlist) do
  if not(part(denlist,i)=1) or done then
   tmplist:=append(tmplist,{ part(denlist,i)})
  else
   done:=1;
 denlist:=tmplist;
>>
else
 numlist:=append(numlist,{1});
tmp:=simplify_combinatorial(tmp)*hypergeom(numlist,denlist,lc);
if lisp !*zb_trace then
  <<
    write "finished conversion in hypergeometric notation";
    write tmp;
  >>;
return tmp;
end$  % sumtohyper


%remove_reduntant_elements({1,3,6},{1,1,1});
%remove_reduntant_elements({1,3,6},{1,1,3});


algebraic procedure remove_reduntant_elements(denlist,numlist);
begin
scalar j,jj,jjj,ln,ld,tmp;
ln:=arglength(numlist);
ld:=arglength(denlist);
if (ln>0) and (ld>0) then
 <<
  for j:=1:arglength(numlist) do
   for jj:=1 : arglength(denlist) do
     if part(numlist,j) =  part(denlist,jj) then
       <<
        tmp:=denlist; denlist:={};
        for jjj:=1 : jj-1 do denlist:=append(denlist,{part(tmp,jjj)});
        for jjj:=jj+1 :arglength(tmp) do
         denlist:=append(denlist,{part(tmp,jjj)});
         tmp:=numlist; numlist:={};
        for jjj:=1 : j-1 do numlist:=append(numlist,{part(tmp,jjj)});
        for jjj:=j+1 :arglength(tmp) do
         numlist:=append(numlist,{part(tmp,jjj)});
       jj:=arglength(denlist)
        >>
 >>;
return {denlist,numlist};
end;



algebraic procedure trim (u);
  if u = {} then {} else
      if member(first u,rest u) then trim rest u
         else first u . trim rest u;



algebraic procedure maxposint(list1);
begin
scalar partj, l,j,tmp;
tmp:=-1;
l:=arglength(list1);
for j:=1 : l do
<<
 partj:=part(list1,j);
 if fixp(partj) and (partj geq 0)  then
   tmp:=max(tmp,partj);
>>;
return tmp;
end;


algebraic procedure minnegint(list1);
begin
scalar partj, l,j,tmp;
tmp:=1;
l:=arglength(list1);
for j:=1 : l do
<<
 partj:=part(list1,j);
 if fixp(partj) and (partj leq 0)  then
   tmp:=min(tmp,partj);
>>;
return tmp;
end;





algebraic procedure binom(n,k);
begin
scalar i;
 if fixp(n) then
  if n>0  then
   return factorial(n)/(factorial(k)*factorial(n-k))
  else
   if n<0  then
    rederr("negative integer argument")
   else return delta(0,k)
 else
  if fixp(k) then
   return (for i:=0:k-1 product(n-i))/factorial(k)
  else
   return binomial(n,k);
end;



algebraic procedure numberofzeros(list1);
begin scalar c,l,j;
c:=0;
l:=arglength(list1);
for j:=1 :l do
 if part(list1,j) = 0 then c:=c+1;
return c;
end;


algebraic procedure patternarguments(term1,pattern,list1);
begin
scalar j,l;
if freeof(term1,pattern) then return list1;
l:=arglength(term1);
if part(term1,0) = pattern then
 return append(list1 ,{part(term1,1)})
else
 for j:=1:l do
  list1:=patternarguments(part(term1,j),pattern,list1);
return list1;
end;

algebraic procedure remove_part(list1,j);
begin
scalar jj,l,list2;
list2:={};
l:=arglength(list1);
for jj:=1 :j-1 do
list2:=append(part(list1,jj),list2);
for jj:=j+1 : l do
 list2:=append(part(list1,jj),list2);
return list2;
end;

algebraic procedure remove_nonlinear_parts(list1,k);
begin
scalar j,list2,!*exp;
on exp;
list2:=list1;
while list1 neq {} do
<<
 if deg(first(list1),k) > 1 then
  rederr("nonlinear argument in gamma")
 else
  if deg(first(list1),k) = 0 then
   list2:=rest(list2);
 list1:=rest(list1)
>>;
return list2;
end;



algebraic procedure closedform_initialization(f,k,n);
begin
scalar co,j,l,ga,mini,maxi,!*exp,ba,b,a,tmpmax,tmpmin;
on exp;
f:=den(simplify_combinatorial(f));
mini:=nil;
maxi:=nil;
ga:=patternarguments(f,gamma,{});
ga:=remove_nonlinear_parts(ga,k);
l:=arglength(ga);
for j:=1 :l do
<<
 co:=coeff(part(ga,j),k);
 a:=part(co,2);
 b:=part(co,1);
 ba:=-b/a;
 if numberp(a) and fixp(ba) then
  if a >0 then
   if maxi = nil then
    maxi:=ba
   else
    maxi:=max(maxi,ba)
  else % a <0
   if mini = nil then
    mini:=ba
   else
    mini:=min(mini,ba)
 else
  if not freeof(ba,n) then
   <<
   if a >0 then
    tmpmax:=ba
   else
    tmpmin:=ba;
   >>;
>>;
if maxi = nil then maxi:=tmpmax;
if mini = nil then mini:=tmpmin;
return {maxi,mini};
end; %closedform_initialization

%f:=(2*pi)^(-1/2)* 2^(2*w2-1/2)*gamma(w2) * gamma(w2+1/2)-gamma(w2*2);
%simplify_gamma2(f);
%simplify_gamman(f,2);


algebraic procedure simplify_gamma2(term1);
begin
scalar p,l,j,jj,jjj,w1,w2,list1,changed;
list1:=patternarguments(term1,gamma,{});
l:=arglength(list1);
changed:={};
for j:=1:l do
<<
 changed:=nil;
 w1:=part(list1,j);
 jj:=0;
 while not changed and jj < l do
  <<
  jj:=jj+1;
  w2:=part(list1,jj);
  p:=w1 - 2* w2;
  if fixp(p) then
   <<
    if p = 0 then
     term1:=sub(gamma(w1)= (2*pi)^(-1/2)* 2^(2*w2-1/2)*
             gamma(w2) * gamma(w2+1/2),term1)
    else
     if p>0 then
      term1:=sub(gamma(w1)= (for jjj:=1 : p product w1-jjj) *
             (2*pi)^(-1/2)* 2^(2*w2-1/2)*
             gamma(w2) * gamma(w2+1/2),term1)
     else
      term1:=sub(gamma(w1)= 1/(for jjj:=0 : (-p-1) product w1 -jjj)*
             (2*pi)^(-1/2)* 2^(2*w2-1/2)*
              gamma(w2) * gamma(w2+1/2),term1);
    changed:=1
   >> % if
  >> %while
>>; % for
return simplify_combinatorial(term1);
end; %simplify_gamma2

%f:=(2*pi)^(-1/2)* 2^(2*w2-1/2)*gamma(w2) * gamma(w2+1/2)-gamma(w2*2);
%simplify_gamma2(f);
%simplify_gamman(f,2);
%simplify_gamman(gamma(3*w2) -subst,3);
%ff:=( - 2*sqrt(3)*gamma(3*w2)*pi + 3**(3*w2)*gamma((3*w2 + 2)/3)*
%    gamma((3*w2 + 1)/3)*gamma(w2))/(2*sqrt(3)*pi)$
%sub(w2=w2+1,ff);
%simplify_gamman(ws,3);
%simplify_gamman(ff,3);



algebraic procedure simplify_gamman(term1,n);
% applies rule 6.1.20 p 77 in Abramowitz
begin
scalar subst,p,l,j,jj,jjj,jjjj,w1,w2,list1,changed;
list1:=patternarguments(term1,gamma,{});
l:=arglength(list1);
changed:={};
for j:=1:l do
<<
 changed:=nil;
 w1:=part(list1,j);
 jj:=0;
 while not changed and jj < l do
  <<
  jj:=jj+1;
  w2:=part(list1,jj);
  p:=w1 - n* w2;
  if fixp(p) then
   <<
    subst:=(2*pi)^(1/2*(1-n))*n^(n*w2-1/2)*(for jjjj:=0:(n-1)
        product (gamma(w2+ jjjj/n)));
    if p = 0 then
     term1:=sub(gamma(w1)=subst ,term1)
    else
     if p>0 then
      term1:=sub(gamma(w1)= (for jjj:=1 : p product w1-jjj) *
             subst,term1)
     else
      term1:=sub(gamma(w1)= 1/(for jjj:=0 : (-p-1) product w1 +jjj)*
             subst,term1);
    changed:=1
   >> % if
  >> %while
>>; % for
return simplify_combinatorial(term1);
end; %simplify_gamman

operator zb_subst;

% simplify_gamma3(f);
algebraic procedure simplify_gamma3(term1);
begin
scalar subst,p,l,j,jj,jjj,w1,w2,list1,changed;
list1:=patternarguments(term1,gamma,{});
l:=arglength(list1);
changed:={};
for j:=1:l do
<<
 changed:=nil;
 w1:=part(list1,j);
 jj:=0;
 while not changed and jj < l do
  <<
  jj:=jj+1;
  w2:=part(list1,jj);
  p:=w1 - 3* w2;
  if fixp(p) then
   <<
   subst:=
  (2*pi)^(-1) * 3^(3*w2-1/2)* gamma(w2) * gamma(w2+1/3)* gamma(w2+2/3);
    if p = 0 then
     term1:=sub(gamma(w1)= zb_subst(j) ,term1)
    else
     if p>0 then
      term1:=sub(gamma(w1)= (for jjj:=1 : p product w1-jjj) *
              zb_subst(j)  ,term1)
     else
      term1:=sub(gamma(w1)= 1/(for jjj:=0 : (-p-1) product w1 +jjj)*
               zb_subst(j)  ,term1);
    term1:=sub(zb_subst(j)=subst,term1);
    changed:=1
   >> % if
  >> %while
>>; % for
return simplify_combinatorial(term1);
end; %simplify_gamma3



% auxiliary functions

symbolic procedure newrederr(u);
   <<terpri!* t;
     prin2!* "***** ";
     if eqcar(u,'list) then foreach xx in cdr u do newrederr1(xx)
       else  newrederr1 u;
     terpri!* nil; erfg!*:=t; error1()>>;

symbolic procedure newrederr1(u);
     if not atom u and atom car u and cdr u and atom cadr u
        and null cddr u
       then <<prin2!* car u; prin2!* " "; prin2!* cadr u>>
      else maprin u;

flag('(newrederr),'opfn);



% some compatibility functions for Maple sources.
% by Winfried Neun

put('polynomqq,'psopfn,'polynomqqq);

algebraic procedure polynomq4(expr1,k);
begin scalar !*exp;
on exp;
return polynomqq(expr1,k);
end;


% checks if expr is rational in var
algebraic procedure type_ratpoly(expr1,var);
begin
scalar deno, nume;
deno:=den expr1;
nume:=num expr1;
  if (polynomqq (deno,var) and polynomqq (nume,var))
    then return t else return nil;
end;
flag ('(type_ratpoly),'boolean);

symbolic procedure tttype_ratpoly(u,xx);
  ( if fixp xx then t else
        if not eqcar (xx , '!*sq) then  nil
          else and(polynomqqq(list(mk!*sq (numr cadr xx ./ 1),
                                  reval cadr u))
                 ,polynomqqq(list(mk!*sq (denr cadr xx ./ 1),
                                  reval cadr u)))
 ) where xx = aeval(car u);

flag ('(tttype_ratpoly),'boolean);





symbolic flag('(savesolve ),'opfn);

symbolic procedure savesolve (x,y);
<< switch solveinconsistent;
   on solveinconsistent;
   inconsistent!*:=nil;
   if pairp (x:=errorset!*(list ('solveeval,mkquote list(x,y)),nil))
        and not inconsistent!* and not (x = '((list)))
        then << x:=car x;
                if eqcar(cadr x,'equal) then % one element solution
                        list('list,x) else x>>
        else <<erfg!*:=nil;list('list)>> >>;


%checks if x is polynomial in var
symbolic procedure polynomq (x,var);

 if not fixp denr simp x then nil else
 begin scalar kerns,kern,aa;

 kerns:=kernels !*q2f simp x;

 aa: if null kerns then return t;
     kern:=first kerns;
     kerns:=cdr kerns;
     if not(eq (kern, var)) and depends(kern,var)
                then return nil else go aa;
end;

flag('(polynomq),'opfn);

flag ('(polynomq type_ratpoly),'boolean);



symbolic procedure polynomqqq (x);

(
 if not fixp denr (xx:=cadr aeval car x) then nil else
 begin scalar kerns,kern,aa,var;

 var:=reval cadr x;
 kerns:=kernels !*q2f xx;

 aa: if null kerns then return t;
     kern:=first kerns;
     kerns:=cdr kerns;
     if not(eq (kern, var)) and depends(kern,var)
                then return nil else go aa;
end) where xx = x;


put('polynomqq,'psopfn,'polynomqqq);

symbolic procedure polynomqqq (x);

(if fixp xx then t else
 if not onep denr (xx:=cadr xx) then nil
 else begin scalar kerns,kern,aa,var,fform,mvv,degg;

 fform:=sfp  mvar  numr xx;
 var:=reval cadr x;
 if fform then << xx:=numr xx;
    while (xx neq 1) do
     << mvv:=mvar  xx;
        degg:=ldeg  xx;
        xx:=lc  xx;
        if domainp mvv then <<if not freeof(mvv,var) then
                << xx:=1 ; kerns:=list list('sin,var) >> >> else
        kerns:=append ( append (kernels mvv,kernels degg),kerns) >> >>
   else kerns:=kernels !*q2f xx;

 aa: if null kerns then return t;
     kern:=first kerns;
     kerns:=cdr kerns;
     if not(eq (kern, var)) and depends(kern,var)
                then return nil else go aa;
end) where xx = aeval(car x);

put('polynomqq,'psopfn,'polynomqqq);

symbolic procedure ttttype_ratpoly(u);
  ( if fixp xx then t else
        if not eqcar (xx , '!*sq) then nil
          else and(polynomqqq(list(mk!*sq (numr cadr xx ./ 1),
                              reval cadr u)),
                   polynomqqq(list(mk!*sq (denr cadr xx ./ 1),
                              reval cadr u)))
 ) where xx = aeval(car u);

flag ('(type_ratpoly),'boolean);

put('type_ratpoly,'psopfn,'ttttype_ratpoly);

endmodule;


end;
