Artifact cda853c7f93611dbc592ed0c6e616589fc0a7cc973f63ef58a8948f1ea355cd0:
- Executable file
r36/src/zeilberg.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 58077) [annotate] [blame] [check-ins using] [more...]
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;