# multsum11.mpl # (C) Torsten Sprenger, University of Kassel, 2003-2007 # # Torsten Sprenger # sprenger@mathematik.uni-kassel.de # http://www.mathematik.uni-kassel.de/~sprenger ##PACKAGE multsum ## ##HALFLINE Overview of the multsum package ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- multsum['command']('arguments') ##- 'command'('arguments') ## ##DESCRIPTION ##- The `multsum` package provides the following functions ##-- "multsum[multsumrecursion]" : determine a linear and homogeneous ## recurrence equation with polynomial coefficients for the multiple (or single) sum over a given ## hypergeometric term ##-- "multsum[findrec]" : find a linear and homogeneous ## recurrence equation with polynomial coefficients for a given ## hypergeometric term ## ##EXAMPLES ##> with(multsum); ## ##SEEALSO ##- "multsum[multsumrecursion]" ##- "multsum[findrec]" ##- "sumtools" ##- "SumTools" ## ##ENDPACKAGE unprotect('multsum'): multsum:=module() option package; export MULTSUMSmax, MULTSUMSolve, MULTSUMFilter, properterm, multsumrecursion, multsumdiffeq, findrec, checkhyper, checkproper, checkrec, sumrec, orderrec, shiftrec, solverec, rect, opt, plotopt, rectocer, certorec, _pexports; local _properterm,_checkhyper,_checkproper,convertgam,convertfacex,facdiff,structurefct,_boundptsSstd,orderrec_, boundptsS,inequat,isfinite,Scov,selectSmax,_krec,_rec,_findrec,maxrec,_sumrec,_sumcer,simplifysol,selectSstd, decomposeS,replacevar,_orderrec,_ordercer,minrec,Jmax,_rectocer,_certorec,_shiftcer,_shiftrec,_checkcer, addrec,prinpart,maxorder,filter,updateF,rat_linear,replacevar_sum,boundptsSstd,_findkfreerec,_checkrec, contentrec_,contentrec,_contentrec,updateF_back,filter_back,simpcomb,ratio,hyperterm,Sstd,Smax,_orderrec_, _ordercer_; #_pexports:=()->[op({exports(thismodule)} minus {':-_pexports'})]; _pexports:=()->[':-properterm',':-multsumrecursion',':-multsumdiffeq',':-findrec',':-checkhyper',':-checkproper',':-checkrec',':-sumrec', ':-orderrec',':-shiftrec',':-solverec',':-rect',':-opt',':-plotopt',':-rectocer',':-certorec']; MULTSUMSmax:=10: MULTSUMSolve:=SolveTools[Linear]: MULTSUMFilter:=0: simpcomb:=sumtools[simpcomb]: hyperterm:=sumtools[hyperterm]; ratio:=(f,n)->simpcomb(subs(n=n+1,f)/f): ########################################################### # # # '_checkhyper' überprüft, ob ein eingegebener Term t # # hypergeometrisch ist. Ist das nicht der Fall, so wird # # ein Fehler ausgegeben, andernfalls passiert nichts. # # # ########################################################### _checkhyper:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local i,var,hyper; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; var:=[n,op(k)]; hyper:=andmap(evalb, {seq(type(ratio(t,var[i]),ratpoly(anything,var[i])), i=1..nops(var)) }); if not(hyper) then error "Term is not hypergeometric" end if; end proc: ##PROCEDURE(help) checkhyper ## ##HALFLINE determine whether a term is hypergeometric or not ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- checkhyper('F', 'k', 'n') ## ##PARAMETERS ##- 'F' : algebraic expression ##- 'k' : summation variable or list of summation variables ##- 'n' : recurrence variable ## ##DESCRIPTION ##- This command decides whether the input term is a hypergeometric term in 'k1',...,'kr' (k=[k1,...,kr]) and 'n', i.e. ## whether the ratios _F(n+1,k)/F(n,k)_ and _F(n,k1,...,kj+1,...,kr)/F(n,k1,...,kj,...,kr)_ are rational functions in 'k1',...,'kr' and 'n'. ## ##EXAMPLES ##> with(multsum): ##> checkhyper(binomial(n,k), k, n); ##< true ##> checkhyper(pochhammer(n,k), k, n); ##< true ##> checkhyper(k/(n^2+k^2), k, n); ##< true ##> checkhyper(n^k, k, n); ##< false ## ##SEEALSO ##- "multsum" ##- "multsum[checkproper]" ##- "multsum[properterm]" ##ENDPROCEDURE checkhyper:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local check; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; check:=_checkhyper(t,k,n); if check=NULL then return true else return false end if; end proc: ########################################################### # # # '_checkproper' überprüft, ob ein eingegebener Term t # # zulässig (proper) ist. Ist das nicht der Fall, so # # wird ein Fehler ausgegeben, andernfalls wird das # # Polynom (pol) und der GAMMA-Ausdruck (gam) des PHT's # # (proper hypergeometric term) zurückgegeben. # # # ########################################################### _checkproper:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local i,x,term,kk,pol_den,pol_den_,den_set,var,rules,gam_den,gam_num,gam,pol; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; _checkhyper(t,k,n); kk:=[op(k)]; var:=[n,op(k)]; rules:=[seq((x::freeof(var))^var[i]=1,i=1..nops(kk)+1)]; term:=applyrule(rules,simpcomb(t)); gam_num:=numer(term); pol_den:=factor(applyrule(GAMMA(x::anything)=1,denom(term))); pol_den_:=1; if type(pol_den,`*`) then for i from 1 to nops(pol_den) do if type(op(i,pol_den),freeof(var)) then gam_num:=gam_num/op(i,pol_den); pol_den_:=pol_den_*op(i,pol_den); end if; end do; end if; pol_den:=pol_den/pol_den_; gam_den:=1; if not(type(pol_den,freeof(var))) then den_set:=convert(factor(pol_den),`multiset`); for i from 1 to nops(den_set) do if degree(den_set[i,1],{op(var)})=1 then gam_den:=gam_den*(GAMMA(den_set[i,1]+1)/GAMMA(den_set[i,1]))^(den_set[i,2]) else error "Hypergeometric term is not proper" end if end do else pol_den:=1; end if; gam_den:=gam_den*(denom(term)/pol_den); gam:=gam_num/gam_den; pol:=applyrule(GAMMA(x::anything)=1,gam); gam:=normal(gam/pol); return [pol,gam]; end proc: ##PROCEDURE(help) checkproper ## ##HALFLINE determine whether a term is a proper hypergeometric term or not ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- checkproper('F', 'k', 'n') ## ##PARAMETERS ##- 'F' : algebraic expression ##- 'k' : summation variable or list of summation variables ##- 'n' : recurrence variable ## ##DESCRIPTION ##- This command decides whether the input term 'F' is a proper hypergeometric term in 'k1',...,'kr' (k=[k1,...,kr]) and 'n', i.e. ## whether the ratios _F(n+1,k)/F(n,k)_ and _F(n,k1,...,kj+1,...,kr)/F(n,k1,...,kj,...,kr)_ are rational functions in 'k1',...,'kr' and 'n' and ## 'F' has a representation as ## _P(n,k)*product((a[p]*n+add(b[l,p]*k[l],l=1..r)+c[p])!,p=1..P)/product((d[q]*n+add(e[l,q]*k[l],l=1..r)+f[q])!,q=1..Q)*product(x[l]^k[l],l=1..r)*y^n_, ## where 'P(n,k)' is a polynomial, 'x[l]','y' are complex numbers and 'a[p]','b[l,p]','c[p]','d[q]','e[r,q]','f[q]','P','Q' are integers. ## ##EXAMPLES ##> with(multsum): ##> checkproper(binomial(n,k), k, n); ##< true ##> checkproper(pochhammer(n,k), k, n); ##< true ##> checkproper(n^k, k, n); ##< false ##> checkproper(k/(n^2+k^2), k, n); ##< false ## ##SEEALSO ##- "multsum" ##- "multsum[checkhyper]" ##- "multsum[properterm]" ##ENDPROCEDURE checkproper:=proc(t::algebraic,k::Or(symbol,list),n::symbol) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; try _checkproper(t,k,n); catch: return false; end try; return true; end proc: ########################################################### # # # 'convertgam' extrahiert die Argumente der GAMMA- # # Funktionen des GAMMA-Ausdrucks und überführt sie in # # eine Liste für die Zähler-Ausdrücke (num_fac_ex) und # # in eine Liste für die Nenner-Ausdrücke (den_fac_ex). # # # ########################################################### convertgam:=proc(gam::algebraic) local i,j,num,den,upper,lower,num_fac_ex,den_fac_ex; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; num:=convert(numer(gam),`multiset`); den:=convert(denom(gam),`multiset`); upper:=applyrule(GAMMA(x::anything)=x,num); lower:=applyrule(GAMMA(x::anything)=x,den); num_fac_ex:=[]; den_fac_ex:=[]; if not(type(num,freeof(GAMMA))) then for i from 1 to nops(upper) do num_fac_ex:=[op(num_fac_ex),seq(upper[i,1],j=1..upper[i,2])] end do; end if; if not(type(den,freeof(GAMMA))) then for i from 1 to nops(lower) do den_fac_ex:=[op(den_fac_ex),seq(lower[i,1],j=1..lower[i,2])] end do; end if; return [num_fac_ex,den_fac_ex]; end proc: ########################################################### # # # 'convertfacex' bestimmt # # die Koeffizienten der Zähler-Ausdrücke (a,b,c) und # # die Koeffizienten der Nenner-Ausdrücke (u,v,w) und # # gibt sie jeweils als Listen aus. # # # ########################################################### convertfacex:=proc(num_fac_ex::list,den_fac_ex::list,k::Or(symbol,list),n::symbol,{facoutput::boolean:=false}) local i,j,kk,pp,qq,a,b,c,u,v,w,r,x,var; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; kk:=[op(k)]; var:=[n,op(k)]; r:=nops(k); pp:=nops(num_fac_ex); qq:=nops(den_fac_ex); a:=[seq(coeff(num_fac_ex[i],n,1),i=1..pp)]; b:=[seq([seq(coeff(num_fac_ex[i],var[j+1],1),i=1..pp)],j=1..r)]; c:=[]; for i from 1 to pp do if ldegree(num_fac_ex[i],var)=1 then c:=[op(c),0] else c:=[op(c),tcoeff(num_fac_ex[i],var)] end if end do; u:=[seq(coeff(den_fac_ex[i],n,1),i=1..qq)]; v:=[seq([seq(coeff(den_fac_ex[i],var[j+1],1),i=1..qq)],j=1..r)]; w:=[]; for i from 1 to qq do if ldegree(den_fac_ex[i],var)=1 then w:=[op(w),0] else w:=[op(w),tcoeff(den_fac_ex[i],var)] end if end do; if facoutput then c:=map(x->x-1,c); w:=map(x->x-1,w) end if; return [a,b,c,u,v,w] end proc: ########################################################### # # # 'properterm' bestimmt die Bestandteile eines einge- # # gebenen Terms t. Ausgegeben wird eine Liste mit drei # # Elementen: # # (1): Polynom (pol) # # (2): die Koeffizienten der Zähler- und Nenner- # # Ausdrücke (a,b,c,u,v,w) # # (3): restliche Potenzen. # # Die lokale Prozedur gibt die Koeffizienten bzgl. GAMMA # # aus, die exportierte Prozedur bzgl. Fakultäten. # # # ########################################################### _properterm:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local pol_gam,abc_uvw,var; option remember, `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; var:=[n,op(k)]; pol_gam:=_checkproper(t,k,n); abc_uvw:=convertgam(pol_gam[2]); return [pol_gam[1],convertfacex(abc_uvw[1],abc_uvw[2],k,n),simpcomb(t/(pol_gam[1]*pol_gam[2]))]; end proc: ##PROCEDURE(help) properterm ## ##HALFLINE determine the parameters of a proper hypergeometric term ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- properterm('F', 'k', 'n') ## ##PARAMETERS ##- 'F' : algebraic expression ##- 'k' : summation variable or list of summation variables ##- 'n' : recurrence variable ## ##DESCRIPTION ##- This command determines a list of upper and lower parameters for a proper hypergeometric term 'F' in 'k1',...,'kr' (k=[k1,...,kr]) and 'n'. ## If 'F' is proper then 'F' has the representation ## _P(n,k)*product((a[p]*n+add(b[l,p]*k[l],l=1..r)+c[p])!,p=1..P)/product((d[q]*n+add(e[l,q]*k[l],l=1..r)+f[q])!,q=1..Q)*product(x[l]^k[l],l=1..r)*y^n_ of 'F', ## where 'P(n,k)' is a polynomial, 'x[l]','y' are complex numbers and 'a[p]','b[l,p]','c[p]','d[q]','e[r,q]','f[q]','P','Q' are integers. ## The output of 'properterm' is the following list ## _[P(n,k),[[a[1],...,a[r]],[[b[1,1],...,b[1,P]],...,[b[r,1],...,b[r,P]]],[c[1],...,c[r]], ## [d[1],...,d[r]],[[e[1,1],...,e[1,Q]],...,[e[r,1],...,e[r,Q]]],[f[1],...,f[r]]],product(x[l]^k[l],l=1..r)*y^n]_. Note: This representation is not unique ## ##EXAMPLES ##> with(multsum): ##> properterm(binomial(n,k), k, n); ##< [1, [[1], [[0]], [0], [1, 0], [[-1, 1]], [0, 0]], 1] ##> properterm(pochhammer(n,k), k, n); ##< [1, [[1], [[1]], [-1], [1], [[0]], [-1]], 1] ## ##SEEALSO ##- "multsum" ##- "multsum[checkhyper]" ##- "multsum[checkproper]" ##ENDPROCEDURE properterm:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local pol_gam,abc_uvw,var; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; var:=[n,op(k)]; pol_gam:=_checkproper(t,k,n); abc_uvw:=convertgam(pol_gam[2]); return [pol_gam[1],convertfacex(abc_uvw[1],abc_uvw[2],k,n,facoutput),simpcomb(t/(pol_gam[1]*pol_gam[2]))]; end proc: ########################################################### # # # 'facdiff' bestimmt die Fakultäts-Differenz # # und gibt sie als [A_t,B_t] aus. # # # ########################################################### facdiff:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local i,j,r,propterm,abc_uvw,p,q,pp,qq,A,B; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; propterm:=_properterm(t,k,n); r:=nops(k); pp:=nops(propterm[2,1]); qq:=nops(propterm[2,4]); abc_uvw:=propterm[2]; A:=0; for p from 1 to pp do if {seq(abc_uvw[2,j,p],j=1..r)}<>{0} then A:=A+abc_uvw[1,p] end if end do; for q from 1 to qq do if {seq(abc_uvw[5,j,q],j=1..r)}<>{0} then A:=A-abc_uvw[4,q] end if end do; for i from 1 to r do B[i]:=0; for p from 1 to pp do B[i]:=B[i]+abc_uvw[2,i,p] end do; for q from 1 to qq do B[i]:=B[i]-abc_uvw[5,i,q] end do; end do; return [A,seq(B[i],i=1..r)]; end proc: ########################################################### # # # 'structurefct' bestimmt die Parameter der Struktur- # # funktionen eines PHT's t. # # # ########################################################### structurefct:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local i,j,r,propterm,abc_uvw,AB,p,q,pp,qq,gh_abc,gh_uvw,gh_AB; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; propterm:=_properterm(t,k,n); r:=nops(k); pp:=nops(propterm[2,1]); qq:=nops(propterm[2,4]); abc_uvw:=propterm[2]; for p from 1 to pp do gh_abc[p]:=[]; if {seq(abc_uvw[2,j,p],j=1..r)}<>{0} then gh_abc[p]:=[op(gh_abc[p]),abc_uvw[1,p]/igcd(abc_uvw[1,p],seq(abc_uvw[2,j,p],j=1..r))]; for i from 1 to r do gh_abc[p]:=[op(gh_abc[p]),abc_uvw[2,i,p]/igcd(abc_uvw[1,p],seq(abc_uvw[2,j,p],j=1..r))]; end do end if end do; for q from 1 to qq do gh_uvw[q]:=[]; if {seq(abc_uvw[5,j,q],j=1..r)}<>{0} then gh_uvw[q]:=[op(gh_uvw[q]),-abc_uvw[4,q]/igcd(abc_uvw[4,q],seq(abc_uvw[5,j,q],j=1..r))]; for i from 1 to r do gh_uvw[q]:=[op(gh_uvw[q]),-abc_uvw[5,i,q]/igcd(abc_uvw[4,q],seq(abc_uvw[5,j,q],j=1..r))]; end do end if end do; gh_AB:=NULL; AB:=facdiff(t,k,n); if {op(AB)}<>{0} then gh_AB:=[seq(-AB[j]/igcd(seq(AB[i],i=1..r+1)),j=1..r+1)] end if; return remove(x->(x=[]),[seq(gh_abc[p],p=1..pp), seq(gh_uvw[q],q=1..qq), gh_AB]); end proc: ########################################################### # # # '_boundptsSstd' bestimmt Randpunkte der # # Standard-Strukturmenge Sstd (bp) und gibt zusätzlich # # die zugrunde liegenden Parameter der Struktur- # # funktionen (sf) aus. # # # ########################################################### _boundptsSstd:=proc(t::algebraic,k::Or(symbol,list),n::symbol,IJ::list) local i,j,num_IJ,sf,num_sf,bp; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; sf:=structurefct(t,k,n); num_IJ:=nops(IJ); num_sf:=nops(sf); for i from 1 to num_sf do for j from 1 to num_IJ do if type(sf[i,j],nonnegative) then bp[i,j]:=IJ[j] else bp[i,j]:=0 end if; end do; end do; return [sf,[seq([seq(bp[i,j],j=1..num_IJ)],i=1..num_sf)]]; end proc: boundptsSstd:=proc(t::algebraic,k::Or(symbol,list),n::symbol,IJ::list) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; return _boundptsSstd(t,k,n,IJ)[2]; end proc: ########################################################### # # # 'boundptsS' bestimmt Randpunkte einer # # Strukturmenge S (bp) und gibt zusätzlich # # die zugrunde liegenden Parameter der Struktur- # # funktionen (sf) aus. # # # ########################################################### boundptsS:=proc(t::algebraic,k::Or(symbol,list),n::symbol,S::set) local i,j,num_S,sf,num_sf,bp; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; sf:=structurefct(t,k,n); num_S:=nops(S); num_sf:=nops(sf); for i from 1 to num_sf do bp[i]:=[S[1,1],S[1,2]]; for j from 2 to num_S do if bp[i][1]*sf[i,1]+bp[i][2]*sf[i,2] with(multsum): ##> rect([1,2]); ##< {[1, 1], [0, 0], [1, 0], [0, 1], [2, 0], [2, 1]} ##> rect([[0,2],1]); ##< {[0, 0, 0], [0, 0, 1], [0, 0, 2], [1, 0, 0], [1, 0, 1], [1, 0, 2]} ## ##SEEALSO ##- "multsum" ##- "multsum[opt]" ##- "multsum[plotopt]" ##ENDPROCEDURE rect:=proc(IJ::list) local newIJ; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; newIJ:=ListTools[Flatten](IJ); Sstd([newIJ[-1],op(1..-2,newIJ)]) end proc: ########################################################### # # # 'Scov' erzeugt aus einer vorgegebenen Struktur- # # menge Sstd eine neue Strukturmenge Scov, in der # # die P-maximale Strukturmenge Smax bezüglich Sstd # # liegt. In Scov werden nur die Grenzen gespeichert. # # # ########################################################### Scov:=proc(t::algebraic,k::Or(symbol,list),n::symbol,IJ::list) local i,num_IJ,Scov_,min_max_IJ,extend; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; num_IJ:=nops(IJ); min_max_IJ:=isfinite(t,k,n,IJ); if type(min_max_IJ,`list`) then Scov_:=min_max_IJ; else if 4 with(multsum): ##> opt(binomial(n,k), k, n, [1,2]); ##< {[2, 2], [1, 1], [0, 0], [0, 1], [1, 0], [1, 2], [2, 0], [2, 1], [2, 3], [-1, 0]} ##> opt(binomial(n,i)*binomial(i,j), [i,j], n, [[0,1],1]); ##< {[0, 0, 1], [-1, -1, 0], [1, 1, 0], [1, 0, 0], [0, -1, 0], [1, 1, 2], [0, 0, 0], [1, -1, 0], [1, 0, 1], [1, 1, 1]} ## ##SEEALSO ##- "multsum" ##- "multsum[plotopt]" ##- "multsum[rect]" ##ENDPROCEDURE opt:=proc(t::algebraic,k::Or(symbol,list),n::symbol,IJ::list) local newIJ; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; newIJ:=ListTools[Flatten](IJ); Smax(t,k,n,[newIJ[-1],op(1..-2,newIJ)]) end proc: ########################################################### # # # 'plotopt' zeichnet die Punkte der P-maximalen # # Strukturmenge eines PHT's in ein zwei- bzw. drei- # # dimensionales Koordinatensystem. # # # ########################################################### ##PROCEDURE(help) plotopt ## ##HALFLINE plot an optimal structure set ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- plotopt('F', 'k', 'n', 'IJ') ## ##PARAMETERS ##- 'F' : algebraic expression ##- 'k' : summation variable or list of summation variables ##- 'n' : recurrence variable ##- 'IJ' : list of upper bounds for each component of the rectangular structure set. The first entry must be a list of nonnegative integers ## according to the upperbounds w.r.t. the summation variables. The second entry represents the upper bound w.r.t. the recurrence vairable. If there's only one summation variable, ## then the list parentheses of the first list can be omitted. 'IJ' represents the underlying rectangular structure set of the optimal structure set. ## ##DESCRIPTION ##- This command plots an optimal structure set for the rectangular structure set represented as IJ=[[x1],y] or IJ=[[x1,x2],y]. The underlying rectangular structure set and the boundary points are also plotted. ## If `plotopt` creates a three-dimensional plot, it will be useful to rotate the whole plot by the mouse to get an overview of the optimal structure set. ## ##OPTIONS ##- `border`: nonnegative integer (default: `3`). This option defines the size of the border, i.e. the distance between the plotted area and the end of the plot. ##- `hyperplane`: boolean (default: `true`). If set to `true` the corresponding structure hyperplanes (in the case of one summation variable: structure lines) ## are plotted. ## ##EXAMPLES ##> with(multsum): ##> plotopt(binomial(n,k), k, n, [2,3]); ##> plotopt(binomial(n,i)*binomial(i,j), [i,j], n, [[1,2],2]); ##> plotopt(binomial(n,i)*binomial(i,j), [i,j], n, [[1,2],2], border=0, hyperplane=false); ## ##SEEALSO ##- "multsum" ##- "multsum[opt]" ##- "multsum[rect]" ##ENDPROCEDURE plotopt:=proc(t::algebraic,k::Or(symbol,list),n::symbol,IJ::list,{border::nonnegint:=3,hyperplane::boolean:=true}) local i,j,Smax_pts,Sstd_pts,extend,extend_,equat,lines,Smax_,Sstd_,num_IJ,Scov_,bdpts,bdpts_, plotlist,planes,newIJ; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; extend:=border; newIJ:=ListTools[Flatten](IJ); newIJ:=[newIJ[-1],op(1..-2,newIJ)]; Sstd_pts:=Sstd(newIJ); Smax_pts:=Smax(t,k,n,newIJ,extend) minus Sstd_pts; bdpts:=_boundptsSstd(t,k,n,newIJ)[2]; num_IJ:=nops(newIJ); Scov_:=Scov(t,k,n,newIJ,0); equat:=convert(inequat(t,k,n,newIJ,x),equality); if num_IJ=2 then lines:=plots[implicitplot](equat, x[1]=Scov_[1,1]-extend..Scov_[1,2]+extend, x[2]=Scov_[2,1]-extend..Scov_[2,2]+extend,color=black,linestyle=DOT); if not hyperplane then lines:=NULL end if; Smax_:=seq(plottools[point](Smax_pts[i], symbol=circle,color=grey,symbolsize=15), i=1..nops(Smax_pts)); Sstd_:=seq(seq(plottools[point](Sstd_pts[i], symbol=circle,color=black,symbolsize=15-j), i=1..nops(Sstd_pts)),j=0..14); bdpts_:=seq(plottools[point](bdpts[i], symbol=circle,color=black,symbolsize=22), i=1..nops(bdpts)); return plots[display](Sstd_,Smax_,bdpts_,lines, axes=normal,scaling=constrained, view=[Scov_[1,1]-extend..Scov_[1,2]+extend, Scov_[2,1]-extend..Scov_[2,2]+extend]); end if; if num_IJ=3 then for j from 1 to 3 do plotlist[j]:=[] end do; for i from 1 to nops(equat) do for j from 1 to 3 do if has(equat[i],x[j]) then plotlist[j]:=[op(plotlist[j]),solve(equat[i],x[j])]; break; end if; end do; end do; if Scov_[1][1]=0 and Scov_[1][2]=0 then extend_:=extend; else extend_:=0; end if; planes[1]:=plot3d({seq([plotlist[1][j],x[2],x[3]],j=1..nops(plotlist[1]))}, x[2]=Scov_[2][1]-extend_..Scov_[2][2]+extend_, x[3]=Scov_[3][1]-extend_..Scov_[3][2]+extend_,style=wireframe,color=grey,grid=[2,2]): planes[2]:=plot3d({seq([x[1],plotlist[2][j],x[3]],j=1..nops(plotlist[2]))}, x[1]=Scov_[1][1]-extend_..Scov_[1][2]+extend_, x[3]=Scov_[3][1]-extend_..Scov_[3][2]+extend_,style=wireframe,color=grey,grid=[2,2]): planes[3]:=plot3d({seq([x[1],x[2],plotlist[3][j]],j=1..nops(plotlist[3]))}, x[1]=Scov_[1][1]-extend_..Scov_[1][2]+extend_, x[2]=Scov_[2][1]-extend_..Scov_[2][2]+extend_,style=wireframe,color=grey,grid=[2,2]): if not hyperplane then planes[1]:=NULL;planes[2]:=NULL;planes[3]:=NULL; end if; Smax_:=seq(plottools[point](Smax_pts[i], symbol=circle,color=grey,symbolsize=10), i=1..nops(Smax_pts)); Sstd_:=seq(seq(plottools[point](Sstd_pts[i], symbol=circle,color=black,symbolsize=10-j), i=1..nops(Sstd_pts)),j=0..9); bdpts_:=seq(plottools[point](bdpts[i], symbol=box,color=black,symbolsize=15), i=1..nops(bdpts)); return plots[display3d](planes[1],planes[2],planes[3],Sstd_,Smax_,bdpts_, axes=normal,scaling=constrained, view=[Scov_[1,1]-extend..Scov_[1,2]+extend, Scov_[2,1]-extend..Scov_[2,2]+extend, Scov_[3,1]-extend..Scov_[3,2]+extend]); end if; end proc: ########################################################### # # # 'selectSmax' sortiert die P-maximalen Struktur- # # Mengen aufsteigend nach der Grösse. # # Als obere Schranke dient K. # # # ########################################################### selectSmax:=proc(t::algebraic,k::Or(symbol,list),n::symbol,K::posint) local i,r,x,y,IJ,Smax_,Smax_sorted,num_IJ; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r:=nops(k); IJ:=Sstd([seq(K,i=1..r+1)]) minus {[seq(0,i=1..r+1)]}; num_IJ:=nops(IJ); Smax_:=[seq([IJ[i],Smax(t,k,n,IJ[i])],i=1..num_IJ)]; Smax_:=sort(Smax_,(x,y)->evalb(nops(x[2])<=nops(y[2]))); Smax_sorted:=[Smax_[1]]; for i from 2 to num_IJ do if evalb(nops(Smax_[i,2])<>nops(Smax_[i-1,2])) then Smax_sorted:=[op(Smax_sorted),Smax_[i]] end if; end do; return Smax_sorted; end proc: ########################################################### # # # 'selectSsstd' sortiert die rechteckigen Struktur- # # Mengen aufsteigend nach der Grösse. # # Als obere Schranke dient upperbound. # # # ########################################################### selectSstd:=proc(k::Or(symbol,list),K::posint,{upperbound::integer:=-1}) local i,r,x,y,IJ,Sstd_,num_IJ; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r:=nops(k); IJ:=Sstd([seq(K,i=1..r+1)]) minus {[seq(0,i=1..r+1)]}; num_IJ:=nops(IJ); if upperbound<=0 then Sstd_:=[seq([IJ[i],Sstd(IJ[i])],i=1..num_IJ)]; else Sstd_:=[]; for i from 1 to num_IJ do if convert(map(x->x+1,IJ[i]),'`*`')<=upperbound and convert(IJ[i],'`*`')<>0 then Sstd_:=[op(Sstd_),[IJ[i],Sstd(IJ[i])]] end if; end do; end if; Sstd_:=sort(Sstd_,(x,y)->evalb(x[1][1]<=y[1][1])); Sstd_:=sort(Sstd_,(x,y)->evalb(nops(x[2])<=nops(y[2]))); return Sstd_; end proc: ########################################################### # # # 'replacevar' ersetzt die auftretenden Variablen in # # einer Rekursionsgleichung durch möglichst einfache # # Zahlen (0 oder 1). # # # ########################################################### replacevar:=proc(recurrence,F::symbol,n::symbol,a::symbol) local i,j,var,num_var,subs_S,rec_wvar,rules,recurrences; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; var:=select(has,indets(recurrence),a); num_var:=nops(var); recurrences:={}; if num_var>0 then for i from 1 to num_var do rec_wvar:=subs(var[i]=1,recurrence); rec_wvar:=subs({seq(var[j]=0,j=1..num_var)},rec_wvar); if orderrec_(rec_wvar=0)>=1 then rec_wvar:=numer(rec_wvar); rec_wvar:=contentrec_(rec_wvar,F,n,a); rec_wvar:=collect(rec_wvar,F,factor); recurrences:=recurrences union {rec_wvar=0}; end if; end do; end if; return recurrences; end proc: ########################################################### # # # 'replacevar_sum' ersetzt die auftretenden Variablen # # in einer Rekursionsgleichung so, daß möglichst viele # # Summanden wegfallen. # # # ########################################################### replacevar_sum:=proc(recurrence,S::symbol,n::symbol,a::symbol) local i,rec,max_r,min_r,eq,temp_rec; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; rec:=recurrence; min_r:=minrec(rec)[1]; max_r:=maxrec(rec)[1]; for i from min_r to max_r-2 do try eq:=coeff(collect(lhs(rec)-rhs(rec),S),S(n+i)); temp_rec:=subs(solve(eq,select(has,indets(rec),a)),rec); temp_rec:=numer(lhs(temp_rec)-rhs(temp_rec))=0; if orderrec_(temp_rec)>=1 then rec:=temp_rec; end if; catch: end try; end do; temp_rec:=replacevar(lhs(rec),S,n,a); if nops(temp_rec)>0 then rec:=temp_rec[1]; end if; return op(2,contentrec(rec)); end proc: ########################################################### # # # 'contentrec' bestimmt den Inhalt einer Rekursion # # und teilt die Rekursionsgleichung durch diesen. # # # ########################################################### _contentrec:=proc(recurrence) local i,content,rec,F_list,F,recurrence_; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if type(recurrence,equation) then recurrence_:=lhs(recurrence)-rhs(recurrence); else recurrence_:=recurrence; end if; F_list:=select(type,indets(recurrence_),function); F:=op(0,F_list[1]); rec:=factor(recurrence_); content:=1; if type(rec,`*`) then for i from 1 to nops(rec) do if type(op(i,rec),freeof(F)) then content:=content*op(i,rec) end if; end do; end if; return [content,collect(normal(rec/content),F,factor)=0]; end proc: contentrec:=proc(recurrence) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if type(recurrence,list) or type(recurrence,set) then return map(x->_contentrec(x),[op(recurrence)]); else return _contentrec(recurrence); end if; end proc: ########################################################### # # # 'contentrec_' bestimmt den Inhalt einer Rekursion # # und teilt die Rekursionsgleichung durch Faktoren von # # diesen, dessen Nullstellen nichtnegativ sind. # # # ########################################################### contentrec_:=proc(recurrence,F::symbol,n::name,a::name) local i,content,rec,order; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; rec:=recurrence=0; order:=orderrec_(rec)-maxrec(rec)[1]; rec:=factor(lhs(rec)); content:=1; if type(rec,`*`) then for i from 1 to nops(rec) do if type(op(i,rec),freeof(F)) then if has(op(i,rec),a) then content:=content*op(i,rec); else if nops(select(type,[solve(subs(n=n+order,op(i,rec)),n)],nonnegint))=0 then content:=content*op(i,rec); end if; end if; end if; end do; end if; return normal(rec/content); end proc: ########################################################### # # # 'rat_linear' gibt die Zahlen [m,l1,...,lr] bei Eingabe # # eines (m,l1,...lr)-fachen hypergeometr. Terms aus. # # # ########################################################### rat_linear:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local term,i,m,p,q,pp,qq,r,denom_list; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r:=nops(k); term:=convertfacex(op(convertgam(select(has,simpcomb(t),GAMMA))),k,n); pp:=nops(term[1]); qq:=nops(term[4]); for i from 1 to r+1 do denom_list[i]:=[]; end do; for p from 1 to pp do denom_list[1]:=[op(denom_list[1]),denom(term[1][p])]; for i from 1 to r do denom_list[i+1]:=[op(denom_list[i+1]),denom(term[2][i][p])]; end do; end do; for q from 1 to qq do denom_list[1]:=[op(denom_list[1]),denom(term[4][q])]; for i from 1 to r do denom_list[i+1]:=[op(denom_list[i+1]),denom(term[5][i][q])]; end do; end do; return [seq(lcm(op(denom_list[m])),m=1..r+1)]; end proc: ########################################################### # # # 'updateF' wandelt einen (m,l1,...lr)-fachen hyper- # # geometr. Terms in einen hypergeometr. Term um. # # # ########################################################### updateF:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local rat,i,var,r; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r:=nops(k); var:=[n,op(k)]; rat:=rat_linear(t,k,n); return subs(seq(var[i]=rat[i]*var[i],i=1..r+1),t); end proc: ########################################################### # # # 'updateF_back' wandelt eine Rekursion in eine # # Rekursion bzgl. des Ausgangsterms um. # # # ########################################################### updateF_back:=proc(recurrence,t::algebraic,k::Or(symbol,list),n::symbol,F::symbol) local i,m,new_rec,rat_,tt,pol,var,r; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r:=nops(k); var:=[n,op(k)]; rat_:=rat_linear(t,k,n); new_rec:=0; for i from 1 to nops(recurrence) do tt:=select(has,op(i,recurrence),F); pol:=subs({seq(var[m]=var[m]/rat_[m],m=1..r+1)},normal(op(i,recurrence)/tt)); new_rec:=new_rec+pol*subs({seq(op(m,tt)=var[m]+expand(rat_[m]*(op(m,tt)-var[m])),m=1..r+1)},tt); end do; return new_rec; end proc: ########################################################### # # # '_krec' gibt eine k-freie Rekursionsgleichung aus, # # sofern eine existiert. # # # ########################################################### _krec:=proc(t_::algebraic,k::Or(symbol,list),n::symbol,IJ::list) local t,tt,i,j,m,num_IJ,num_Smax,F,a,Smax_,subs_S,var,variables,eq,IJ_, sol,recurrence,print,rat_term,arg,newIJ; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if 41 then if not(type(IJ[2],`list`)) then Smax_:=IJ[2]; else Smax_:=IJ; print:=0; end if; else Smax_:=IJ; print:=0; end if; else Smax_:=Smax(t,k,n,IJ); IJ_:=IJ; num_IJ:=nops(IJ); end if; num_Smax:=nops(Smax_); for i from 1 to num_Smax do rat_term[i]:=simpcomb(subs({seq(var[m]=var[m]-Smax_[i,m],m=1..num_IJ)},t)/t); end do; eq:=add(a[op(Smax_[i])]*rat_term[i],i=1..num_Smax); variables:=select(has,indets(eq),a); eq:={coeffs(expand(numer(eq)),k)}; if print=1 then newIJ:=ListTools[Flatten](IJ_); userinfo(3,multsum,`structure set:`,[[op(2..-1,newIJ)],newIJ[1]],`\t`,`number of equations:`,nops(eq), `\t`, `number of variables:`,nops(variables)); else userinfo(3,multsum,`number of equations:`,nops(eq), `\t`, `number of variables:`,nops(variables)); end if; sol:=MULTSUMSolve(eq,variables); if {seq(op(2,sol[i]),i=1..nops(sol))}={0} then error "No k-free recurrence equation exists" end if; for i from 1 to num_Smax do arg[i]:=seq(var[m]-Smax_[i,m],m=1..num_IJ); end do; recurrence:=add(a[op(Smax_[i])]*F(arg[i]),i=1..num_Smax); if MULTSUMFilter=1 then recurrence:=filter_back(recurrence,tt,k,n,F) end if; recurrence:=subs(sol,recurrence); recurrence:=updateF_back(recurrence,t_,k,n,F); if nargs<=5 then recurrence:=replacevar(recurrence,F,n,a); else recurrence:=numer(normal(recurrence)); recurrence:=contentrec_(recurrence,F,n,a); recurrence:=collect(recurrence,F,factor); recurrence:={recurrence=0}; end if; return recurrence; end proc: ########################################################### # # # '_findkfreerec' versucht eine k-freie Rekursion # # zu finden. Es werden keine P-maximalen Strukturmengen # # verwendet. Die Strukturmengen werden dabei solange # # vergrössert, bis eine Rekursionsgleichung gefunden # # wird. # # # ########################################################### _findkfreerec:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local i,j,r,IJ,recurrence,F,a,S; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; F:=NULL; a:=NULL; S:=NULL; if 41 then i:=args[2]; for j from 1 to nops(F_list[1]) do if op(indets(op(j,F_list[1])))=i then nvar:=j; end if; end do; v_list:=map(x->(op(nvar,x)-i),F_list); else v_list:=map(x->(op(1,x)-op(indets(op(1,F_list[1])))),F_list); end if; return max(op(v_list))-min(op(v_list)); end proc: _orderrec_:=proc(recurrence) local i,var,r,F_list,v_list; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; F_list:=[op(select(type,indets(lhs(recurrence)-rhs(recurrence)),function))]; r:=nops([op(F_list[1])]); var:=map(op@indets,F_list[1]); for i from 1 to r do v_list[i]:=map(x->(op(i,x)-op(i,var)),F_list); end do; if r<>1 then return [[seq(max(op(v_list[i]))-min(op(v_list[i])),i=2..r)],max(op(v_list[1]))-min(op(v_list[1]))] else return max(op(v_list[1]))-min(op(v_list[1])) end if; #if r<>1 then # return [seq(max(op(v_list[i]))-min(op(v_list[i])),i=1..r)] #else # return max(op(v_list[1]))-min(op(v_list[1])) #end if; end proc: ########################################################### # # # '_ordercer' gibt die Ordnungen von Zertifikats- # # Rekursionen in der ersten Variablen an (optional: # # in einer bestimmten Variablen). # # # ########################################################### _ordercer:=proc(certificate) local var; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if nargs>1 then var:=args[2] else var:=NULL end if; _orderrec(certorec(certificate),var) end proc: _ordercer_:=proc(certificate) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; _orderrec_(certorec(certificate)) end proc: ########################################################### # # # 'orderrec' kombiniert _orderrec und _ordercer. # # # ########################################################### orderrec_:=proc(recurrence::Or(equation,set,list),{variable:=default}) local var; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if variable<>default then var:=variable else var:=NULL end if; if map2(op,0,select(type,map2(op,0,indets(recurrence,function)),indexed))={} then if type(recurrence,list) or type(recurrence,set) then return map(_orderrec,[op(recurrence)],var); else return _orderrec(recurrence,var); end if; else if type(recurrence,list) or type(recurrence,set) then return map(_ordercer,[op(recurrence)],var); else return _ordercer(recurrence,var); end if; end if; end proc: ##PROCEDURE(help) orderrec ## ##HALFLINE determine the order of a recurrence equation (certificate or normal) ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- orderrec('rec') ## ##PARAMETERS ##- 'rec' : recurrence equation or certificate recurrence equation or a list of them ## ##DESCRIPTION ##- Given a recurrence equation or recurrence equations in 'F(n,k1,...kr)' ('r' can be zero), this command determines the order of the recurrence equation(s) (either certificate or normal) in all occuring variables 'n', 'k1', ..., 'kr'. ## The output is a single nonnegative integer for every recurrence equation in the case of one variable and otherwise a list [[x1,...,xr],y], where ## the 'xi''s according to 'ki''s and the 'y' to 'n'. ## ##EXAMPLES ##> with(multsum): ##> orderrec(S(n+1) - (n+1)*S(n) = 0); ##< 1 ##> orderrec(k*F(n,k-1) + (n-1)*F(n+2,k) - (n+k)*F(n-1,k-2) = 0); ##< [[2],3] ## ##SEEALSO ##- "multsum" ##- "multsum[solverec]" ##- "multsum[shiftrec]" ##- "multsum[checkrec]" ##- "multsum[certorec]" ##- "multsum[rectocer]" ## ##ENDPROCEDURE orderrec:=proc(recurrence::Or(equation,set,list)) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if map2(op,0,select(type,map2(op,0,indets(recurrence,function)),indexed))={} then if type(recurrence,list) or type(recurrence,set) then return map(_orderrec_,[op(recurrence)]); else return _orderrec_(recurrence); end if; else if type(recurrence,list) or type(recurrence,set) then return map(_ordercer_,[op(recurrence)]); else return _ordercer_(recurrence); end if; end if; end proc: ########################################################### # # # '_checkrec' überprüft, ob t die eingegebene # # Rekursionsgleichung erfüllt. # # # ########################################################### _checkrec:=proc(recurrence,t::Or(algebraic,equation)) local i,j,var,l_rec,F_list,fulfilled,F,rules,subs_S,rec1,rec2,err; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_rec:=lhs(recurrence)-rhs(recurrence); if type(t,equation) then if orderrec_(recurrence)>orderrec_(t) then rec1:=t; rec2:=recurrence; else rec1:=recurrence; rec2:=t; end if; err:=0; while orderrec_(rec2)>orderrec_(rec1) do rec2:=addrec(rec1,rec2); err:=err+1; if err>200 then error "Recurrences must contain the same variables" end if; end do; rec2:=addrec(rec1,rec2); if lhs(rec2)=0 and rhs(rec2)=0 then return true else return false end if; end if; if type(l_rec,`+`) then F_list:=[op(select(type,indets(l_rec),function))]; var:=[]; for i from 1 to nops(F_list[1]) do var:=[op(var),op(indets(op(i,F_list[1])))]; end do; rules:={}; for i from 1 to nops(F_list) do subs_S:={}; for j from 1 to nops(F_list[i]) do subs_S:=subs_S union {var[j]=op(j,F_list[i])}; end do; rules:=rules union {F_list[i]=subs(subs_S,t)}; end do; fulfilled:=simpcomb(subs(rules,l_rec))=0; if fulfilled then return true else return false end if; end if; end proc: ########################################################### # # # '_checkcer' überprüft, ob t die eingegebene # # Rekursionsgleichung erfüllt. # # # ########################################################### _checkcer:=proc(recurrence,t::Or(algebraic,equation)) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; return _checkrec(certorec(recurrence),t) end proc: ########################################################### # # # 'checkrec' kombiniert _checkrec und _checkcer. # # # ########################################################### ##PROCEDURE(help) checkrec ## ##HALFLINE determine whether a term fulfills a recurrence equation or not or whether two recurrence equations have the same solutions or not ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- checkrec('rec', 'F') ##- checkrec('rec1', 'rec2') ## ##PARAMETERS ##- 'rec' : recurrence equation or certificate recurrence equation ##- 'F' : a hypergeometric term ##- 'rec1' : recurrence equation or certificate recurrence equation ##- 'rec2' : recurrence equation or certificate recurrence equation ## ##DESCRIPTION ##- Given a recurrence equation in 'F(n,k1,...kr)' ('r' can be zero) and a hypergeometric term 'G(n,k1,...,kr)', this command checks ## whether 'G(n,k1,...,kr)' fulfills the (certificate) recurrence equation (output: 'true') or not (output: 'false'). ##- Given two recurrence equations, this command checks whether one of the recurrence equation is a multiple of the other or not. ## ##EXAMPLES ##> with(multsum): ##> checkrec(S(n+1) - (n+1)*S(n) = 0, n!); ##< true ##> checkrec(S(n+1) - (n+1)*S(n) = 0, -(n+1)*S(n) - (n+1)*S(n+1) + S(n+2) = 0); ##< true ## ##SEEALSO ##- "multsum" ##- "multsum[solverec]" ##- "multsum[shiftrec]" ##- "multsum[orderrec]" ##- "multsum[certorec]" ##- "multsum[rectocer]" ## ##ENDPROCEDURE checkrec:=proc(recurrence::Or(equation,set,list),t::Or(algebraic,equation)) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if map2(op,0,select(type,map2(op,0,indets(recurrence,function)),indexed))={} then if type(recurrence,list) or type(recurrence,set) then return map(_checkrec,[op(recurrence)],t); else return _checkrec(recurrence,t); end if; else if type(recurrence,list) or type(recurrence,set) then return map(_checkcer,[op(recurrence)],t); else return _checkcer(recurrence,t); end if; end if; end proc: ########################################################### # # # 'prinpart' gibt den Hauptteil eines Zertifikats aus. # # # ########################################################### prinpart:=proc(certificate) local F_list,delta,l_cer,i,k; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_cer:=lhs(certificate)-rhs(certificate); F_list:=select(type,indets(l_cer),function); delta:={}; for i from 1 to nops(F_list) do if type(op(0,F_list[i]),indexed) then delta:=delta union {F_list[i]} end if; end do; return l_cer-add(delta[k],k=1..nops(delta)); end proc: ########################################################### # # # 'minrec' bestimmt die Minima der Indizes der # # zulässigen hypergeometrischen Terme. # # # ########################################################### minrec:=proc(recurrence) local i,j,mini,r,l_rec,var,F_list; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_rec:=lhs(recurrence)-rhs(recurrence); F_list:=select(type,indets(prinpart(l_rec=0)),function); var:=[]; for i from 1 to nops(F_list[1]) do var:=[op(var),op(indets(op(i,F_list[1])))]; end do; r:=nops(var)-1; for j from 1 to r+1 do mini[j]:=infinity end do; for i from 1 to nops(F_list) do for j from 1 to nops(op(i,F_list)) do if evalb(op(j,op(i,F_list))-var[j]maxi[j]) then maxi[j]:=op(j,op(i,F_list))-var[j] end if; end do; end do; return [seq(maxi[i],i=1..r+1)]; end proc: ########################################################### # # # '_shiftrec' führt eine Indexverschiebung aus. # # # ########################################################### _shiftrec:=proc(recurrence,{direction::identical(up,down):=up}) local i,j,max,min,r,l_rec,var,F_list,a; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_rec:=lhs(recurrence)-rhs(recurrence); F_list:=select(type,indets(l_rec),function); var:=[]; for i from 1 to nops(F_list[1]) do var:=[op(var),op(indets(op(i,F_list[1])))]; end do; r:=nops(var)-1; if direction=up then min:=minrec(recurrence); l_rec:=subs(seq(var[j]=var[j]-min[j],j=1..r+1),l_rec); else max:=maxrec(recurrence); l_rec:=subs(seq(var[j]=var[j]-max[j],j=1..r+1),l_rec); end if; l_rec:=contentrec_(l_rec,op(0,F_list[1]),op(1,var),a); l_rec:=collect(l_rec,op(0,F_list[1]),factor); return l_rec=0; end proc: ########################################################### # # # '_shiftcer' führt eine Indexverschiebung aus. # # # ########################################################### _shiftcer:=proc(recurrence,{direction::identical(up,down):=up}) local i,j,min,r,l_rec,var,F_list; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_rec:=lhs(recurrence)-rhs(recurrence); if direction=up then min:=maxrec(recurrence); min[1]:=minrec(recurrence)[1]; else min:=minrec(recurrence); min[1]:=maxrec(recurrence)[1]; end if; F_list:=select(type,indets(l_rec),function); var:=[]; j:=1; while j with(multsum): ##> shiftrec((n-1)*S(n+1) + S(n) - n*S(n-1) = 0); ##< (-n-1)*S(n)+n*S(n+2)+S(n+1) = 0 ##> shiftrec((n-1)*S(n+1) + S(n) - n*S(n-1) = 0, direction=down); ##< (n-2)*S(n)+S(n-1)+(1-n)*S(n-2) = 0 ## ##SEEALSO ##- "multsum" ##- "multsum[solverec]" ##- "multsum[checkrec]" ##- "multsum[orderrec]" ##- "multsum[certorec]" ##- "multsum[rectocer]" ## ##ENDPROCEDURE shiftrec:=proc(recurrence::Or(equation,list,set),{direction::identical(up,down):=up}) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if map2(op,0,select(type,map2(op,0,indets(recurrence,function)),indexed))={} then if type(recurrence,list) or type(recurrence,set) then return map(_shiftrec,{op(recurrence)},_options); else return _shiftrec(recurrence,_options); end if; else if type(recurrence,list) or type(recurrence,set) then return map(_shiftcer,{op(recurrence)},_options); else return _shiftcer(recurrence,_options); end if; end if; end proc: ########################################################### # # # 'rectocer' wandelt eine Rekursionsgleichung # # (nicht notwendig k-frei) in eine Rekursionsgleichung # # mit Delta-Ausdrücken (Zertifikats-Rekursion) um. # # # ########################################################### _rectocer:=proc(recurrence) local i,j,r,kk,l_rec,max,rec_list,F,F_list,cer_recurrence, delta,Delta,exp,p_part,p_part_test,coeff_,m; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_rec:=lhs(recurrence)-rhs(recurrence); if not(type(l_rec,`+`)) then error "No non-trivial certificate recurrence exists" end if; F_list:=[]; for i from 1 to nops(l_rec) do F_list:=[op(F_list),op(select(type,indets(op(i,l_rec)),function))]; end do; F:=op(0,F_list[1]); kk:=[]; for i from 2 to nops(F_list[1]) do kk:=[op(kk),op(indets(op(i,F_list[1])))]; end do; r:=nops(kk); rec_list:=[]; for i from 1 to nops(l_rec) do rec_list:=[op(rec_list),[op(i,l_rec)/F_list[i], [op(F_list[i])]]]; end do; max:=-minrec(recurrence); cer_recurrence:=0; for i from 1 to r do delta:=0; for j from 1 to nops(rec_list) do exp:=[seq(-(rec_list[j,2,l+1]-kk[l]),l=1..r)]; delta:=delta+add(subs({seq(kk[m]=kk[m]-(max[m+1]-exp[m]),m=1..i-1),kk[i]=kk[i]-l}, rec_list[j,1]*F(op(rec_list[j,2]))),l=0..max[i+1]-exp[i]-1); end do; delta:=subs(kk[i]=kk[i]-1,delta); delta:=collect(delta,F,factor); if delta<>0 then cer_recurrence:=cer_recurrence+Delta[kk[i]](delta); end if; end do; p_part:=0; p_part_test:=0; for j from 1 to nops(rec_list) do exp:=[seq(-(rec_list[j,2,l+1]-kk[l]),l=1..r)]; coeff_:=subs({seq(kk[m]=kk[m]-(max[m+1]-exp[m]),m=1..r)},rec_list[j,1]); p_part:=p_part+coeff_*F(rec_list[j,2,1],seq(kk[m]-max[m+1],m=1..r)); p_part_test:=p_part_test+coeff_*F(rec_list[j,2,1]); end do; p_part:=collect(p_part,F,factor); p_part_test:=collect(p_part_test,F,factor); if p_part=0 or not(type(p_part_test,freeof(kk))) then error "No non-trivial certificate recurrence exists" end if; cer_recurrence:=p_part+cer_recurrence; if nargs>1 then return p_part else return cer_recurrence=0 end if; end proc: rectocer:=proc(recurrence::Or(equation,list,set)) local opt; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if nargs>1 then opt:=1 else opt:=NULL end if; if type(recurrence,list) or type(recurrence,set) then return map(x->_rectocer(x,opt),{op(recurrence)}) else return _rectocer(recurrence,opt) end if; end proc: ########################################################### # # # 'certorec' wandelt eine Rekursionsgleichung mit Delta- # # Ausdrücken (Zertifikats-Rekursion) in eine # # Rekursionsgleichung (nicht notwendig k-frei) um. # # # ########################################################### _certorec:=proc(certificate) local i,m,l_cer,F_list,delta_list,F,recurrence,index; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; l_cer:=lhs(certificate)-rhs(certificate); F_list:=select(type,indets(l_cer),function); delta_list:={}; for i from 1 to nops(F_list) do if type(op(0,F_list[i]),indexed) then delta_list:=delta_list union {F_list[i]} end if; end do; F_list:=F_list minus delta_list; F:=op(0,F_list[1]); recurrence:=subs({seq(delta_list[m]=0,m=1..nops(delta_list))},l_cer); for i from 1 to nops(delta_list) do index:=op(1,op(0,delta_list[i])); recurrence:=recurrence-op(1,delta_list[i]); recurrence:=recurrence+subs(index=index+1,op(1,delta_list[i])); end do; return collect(recurrence,F,factor)=0; end proc: certorec:=proc(recurrence::Or(equation,list,set)) option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if type(recurrence,list) or type(recurrence,set) then return map(x->_certorec(x),{op(recurrence)}) else return _certorec(recurrence) end if; end proc: ########################################################### # # # 'Jmax' bestimmt die Maxima der Komponenten der Punkte # # einer beliebigen Strukturmenge S, wobei die erste # # Komponente nicht betrachtet wird. # # # ########################################################### Jmax:=proc(S) local i,j,m,r,Jmax_; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r:=nops(S[1])-1; Jmax_:=[seq(-infinity,m=1..r)]; for i from 1 to nops(S) do for j from 1 to r do if S[i][j+1]>Jmax_[j] then Jmax_[j]:=S[i][j+1] end if; end do; end do; return Jmax_; end proc: ########################################################### # # # '_rec' versucht eine Rekursionsgleichung zu finden # # (bzgl.der P-maximalen Strukturmenge der Standard- # # Strukturmenge Sstd), bei der der Hauptteil k-frei ist. # # ( 1. Verallgemeinerung ) # # # ########################################################### _rec:=proc(t_::algebraic,k::Or(symbol,list),n::symbol,IJ,M,{structureset::identical(rectangular,optimal):=optimal}) local t,tt,i,j,r,kk,l,m,variables,eq,var,F,a,Smax_,num_Smax,num_IJ,pol,pol_set,Jmax_,newIJ, ansatz,rec_eq,reduced_eq,reduced_variables,recurrence,reduced_sol,sol,IJ_; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if 0<_nrest then F:=_rest[1] end if; if 1<_nrest then a:=_rest[2] end if; tt:=updateF(t_,k,n); if MULTSUMFilter=1 then t:=filter(tt,k,n); else t:=tt; end if; kk:=[op(k)]; var:=[n,op(k)]; IJ_:=IJ; num_IJ:=nops(IJ); if type(op(1,IJ),'list') then Smax_:=op(2,IJ); IJ_:=op(1,IJ); num_IJ:=nops(op(1,IJ)); else if structureset=rectangular then Smax_:=Sstd(IJ) else Smax_:=Smax(t,k,n,IJ) end if; end if; num_Smax:=nops(Smax_); Jmax_:=Jmax(Smax_); ########################################################### for i from 1 to num_Smax do ansatz[i]:=a[op(Smax_[i]),seq(l[m],m=1..num_IJ-1)]* mul(kk[m]^l[m],m=1..num_IJ-1)* simpcomb(subs({seq(var[m]=var[m]-Smax_[i,m],m=1..num_IJ)},t)/t); recurrence[i]:=a[op(Smax_[i]),seq(l[m],m=1..num_IJ-1)]* mul(kk[m]^l[m],m=1..num_IJ-1)* F(op(var-Smax_[i])); for j from 1 to num_IJ-1 do ansatz[i]:=sum(ansatz[i],l[j]=0..M); recurrence[i]:=sum(recurrence[i],l[j]=0..M); end do; end do; ansatz:=add(ansatz[r],r=1..num_Smax); recurrence:=add(recurrence[r],r=1..num_Smax); ########################################################### for i from 1 to num_Smax do pol[op(Smax_[i])]:=a[op(Smax_[i]),seq(l[m],m=1..num_IJ-1)]* mul((kk[m]+Smax_[i,m+1]-Jmax_[m])^l[m],m=1..num_IJ-1); pol[Smax_[i,1]]:=0; for j from 1 to num_IJ-1 do pol[op(Smax_[i])]:=sum(pol[op(Smax_[i])],l[j]=0..M); end do; end do; for i from 1 to num_Smax do pol[Smax_[i,1]]:=expand(pol[Smax_[i,1]]+pol[op(Smax_[i])]) end do; pol_set:=[op({seq(pol[Smax_[i,1]],i=1..num_Smax)})]; for i from 1 to nops(pol_set) do if ldegree(pol_set[i],{op(kk)})=0 then pol_set[i]:=pol_set[i]-tcoeff(pol_set[i],kk) end if; reduced_variables[i]:=select(has,indets(pol_set[i]),a); reduced_eq[i]:={coeffs(pol_set[i],k)}; reduced_sol[i]:=solve(reduced_eq[i],reduced_variables[i]); end do; ########################################################### eq:=subs({seq(op(reduced_sol[i]),i=1..nops(pol_set))},ansatz); variables:=select(has,indets(eq),a); eq:={coeffs(numer(eq),k)}; newIJ:=ListTools[Flatten](IJ_); userinfo(3,multsum,`structure set:`,[[op(2..-1,newIJ)],newIJ[1]],`\t`,`number of equations:`,nops(eq), `\t`, `number of variables:`,nops(variables)); sol:=MULTSUMSolve(eq,variables); if {seq(op(2,sol[i]),i=1..nops(sol))}={0} then error "No non-trivial certificate recurrence exists" end if; recurrence:=subs({seq(op(reduced_sol[i]),i=1..nops(pol_set))},recurrence); if MULTSUMFilter=1 then recurrence:=filter_back(recurrence,tt,k,n,F) end if; recurrence:=subs(sol,recurrence); recurrence:=updateF_back(recurrence,t_,k,n,F); if _nrest<=1 then recurrence:=replacevar(recurrence,F,n,a) else recurrence:=numer(normal(recurrence)); recurrence:=contentrec_(recurrence,F,n,a); recurrence:=collect(recurrence,F,factor); recurrence:={recurrence=0}; end if; return recurrence; end proc: ########################################################### # # # '_findrec' versucht eine Rekursionsgleichung zu # # finden (M=0: k-frei und sonst: nicht k-frei). # # Die Strukturmengen werden dabei solange vergrössert, # # bis eine Rekursionsgleichung gefunden wird. # # # ########################################################### _findrec:=proc(t::algebraic,k::Or(symbol,list),n::symbol,K,M,{structureset::identical(rectangular,optimal):=optimal}) local i,j,r,IJ,recurrence,cer,F,a,Smax_,temp_rec; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; cer:=1; F:=NULL; a:=NULL; if 0<_nrest then if type(_rest[1],name) then F:=_rest[1]; if 1<_nrest then a:=_rest[2]; else a:=NULL; end if; else cer:=_rest[1]; if 1<_nrest then F:=_rest[2]; if 2<_nrest then a:=_rest[3]; else a:=NULL; end if; else F:=NULL; a:=NULL; end if; end if; end if; recurrence:={}; if type(K,'list') then if M=0 then recurrence:=_krec(t,k,n,K,F,a); else recurrence:=_rec(t,k,n,K,M,F,a); end if; if cer=1 then temp_rec:={}; for j from 1 to nops(recurrence) do try temp_rec:=temp_rec union {rectocer(recurrence[j])}; catch "No non-trivial certificate recurrence exists": end try; end do; recurrence:=temp_rec; end if; else if K<>0 then if structureset=optimal then Smax_:=selectSmax(updateF(t,k,n),k,n,K) else Smax_:=selectSstd(k,K,upperbound=K) end if; i:=1; if M=0 then while i<=nops(Smax_) do try recurrence:=_krec(t,k,n,Smax_[i],F,a); if cer=1 then recurrence:=map(x->rectocer(x),recurrence); end if; break; catch "No k-free recurrence equation exists": i:=i+1; end try; end do; else while i<=nops(Smax_) do try recurrence:=_rec(t,k,n,Smax_[i],M,F,a); if cer=1 then temp_rec:={}; for j from 1 to nops(recurrence) do try temp_rec:=temp_rec union {rectocer(recurrence[j])}; catch "No non-trivial certificate recurrence exists": end try; end do; recurrence:=temp_rec; end if; if recurrence={} then i:=i+1; else break; end if; catch "No non-trivial certificate recurrence exists": i:=i+1; end try; end do; end if; else r:=nops(k); IJ:=[1,seq(0,i=1..r)]; i:=1; if M=0 then while i<=r+1 do try recurrence:=_krec(t,k,n,IJ,F,a); if cer=1 then recurrence:=map(x->rectocer(x),recurrence); end if; break; catch "No k-free recurrence equation exists": if i=r+1 then i:=1 else i:=i+1 end if; IJ[i]:=IJ[i]+1; end try; end do; else while i<=r+1 do try recurrence:=_rec(t,k,n,IJ,M,F,a,_options); if cer=1 then temp_rec:={}; for j from 1 to nops(recurrence) do try temp_rec:=temp_rec union {rectocer(recurrence[j])}; catch "No non-trivial certificate recurrence exists": end try; end do; recurrence:=temp_rec; end if; if recurrence={} then if i=r+1 then i:=1 else i:=i+1 end if; IJ[i]:=IJ[i]+1; else break; end if; catch "No non-trivial certificate recurrence exists": if i=r+1 then i:=1 else i:=i+1 end if; IJ[i]:=IJ[i]+1; end try; end do; end if; end if; end if; if recurrence={} then if M=0 then error "No k-free recurrence equation exists" else error "No non-trivial certificate recurrence exists" end if; end if; return recurrence; end proc: findrec:=proc(t::algebraic,k::Or(symbol,list),n::symbol,{structureset::Or(identical(rectangular,optimal),set):=optimal,strategy::nonnegint:=2, certificate::boolean:=true,upperkbound::nonnegint:=1,direction::identical(up,down):=up}) local rec,F,a,K,M,i; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if _nrest>0 then F:=_rest[1] else `tools/genglobal`('F',[],'reset'); F:=`tools/genglobal`('F',indets(t,'name')); end if; if _nrest>1 then a:=_rest[2] else a:=NULL end if; M:=upperkbound; K:=strategy; if type(structureset,'set') then if certificate then rec:=_findrec(t,k,n,[[seq(`-`,i=1..nops(k)+1)],structureset],M,1,F,a) else rec:=_findrec(t,k,n,[[seq(`-`,i=1..nops(k)+1)],structureset],M,0,F,a) end if; elif structureset=optimal then if certificate then rec:=_findrec(t,k,n,K,M,1,F,a) else rec:=_findrec(t,k,n,K,M,0,F,a) end if; elif structureset=rectangular then if M=0 and K=0 then if certificate then rec:=rectocer(_findkfreerec(t,k,n,F,a)) else rec:=_findkfreerec(t,k,n,F,a) end if; else if certificate then rec:=_findrec(t,k,n,K,M,1,F,a,_options[4]) else rec:=_findrec(t,k,n,K,M,0,F,a,_options[4]) end if; end if; end if; shiftrec(rec,_options) end proc: ########################################################### # # # '_sumrec' summiert über alle Summationsvariablen # # der Rekursionsgleichungen in der Menge recurrence. # # # ########################################################### _sumrec:=proc(recurrence) local i,j,m,S,rules,sum_recurrence,recurrences,p_part,F_list, minimal,sum_rec,ord,min_ord,new_rec,a; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; minimal:=true; `tools/genglobal`('S',[],'reset'); S:=`tools/genglobal`('S',indets(recurrence,'name')); if 11 and min_ord>1 then for i from 1 to nops(sum_rec) do for j from i+1 to nops(sum_rec) do new_rec:=addrec(sum_rec[1],sum_rec[2]); if new_rec<>(0=0) then return new_rec end if; end do; end do; end if; return sum_rec[1]; end if; return sum_recurrence; end proc: ########################################################### # # # '_sumcer' summiert über alle Summationsvariablen # # der Zertifikats-Rekursionen in der Menge certificate. # # # ########################################################### _sumcer:=proc(certificate) local i,j,m,S,rules,subs_S,sum_certificate,certificates,p_part,F_list,l_cer,delta_list, minimal,sum_rec,ord,min_ord,new_rec,a; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; minimal:=true; `tools/genglobal`('S',[],'reset'); S:=`tools/genglobal`('S',indets(certificate,'name')); if 11 and min_ord>1 then for i from 1 to nops(sum_rec) do for j from i+1 to nops(sum_rec) do new_rec:=addrec(sum_rec[1],sum_rec[2]); if new_rec<>(0=0) then return new_rec; end if; end do; end do; end if; return sum_rec[1]; end if; return sum_certificate; end proc: ########################################################### # # # 'sumrec' kombiniert _sumrec und _sumcer. # # # ########################################################### sumrec:=proc(recurrence::Or(equation,list,set)) local var; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; var:=NULL; if nargs>1 then if type(args[2],'name') then var:=args[2] end if; end if; if map2(op,0,select(type,map2(op,0,indets(recurrence,function)),indexed))={} then if type(recurrence,list) or type(recurrence,set) then return map(_sumrec,{op(recurrence)},var); else return _sumrec(recurrence,var); end if; else if type(recurrence,list) or type(recurrence,set) then return map(_sumcer,{op(recurrence)},var); else return _sumcer(recurrence,var); end if; end if; end proc: ########################################################### # # # 'solverec' bestimmt die Lösung einer Rekursion vom # # hypergeometrischen Typ. # # # ########################################################### solverec:=proc(recurrence::equation) local A_list,A,k,i,j,mn,mx,start,temp,RE,RE2,arg1,arg2,sol,rat,num,den,upp,low,minupp,minlow,coe0,coe1,x,order,shift,piecewisesol,content,n,m; A_list:=select(type,indets(recurrence),function); if nops(A_list)>2 then error "Input is no recurrence equation of hypergeometric type" end if; k:=op(indets(op(1,A_list[1]))); A:=op(0,A_list[1]); RE:=lhs(recurrence)-rhs(recurrence); # determine the two arguments of A if nops(op(1,RE))=1 then arg1:=op(1,op(1,RE))-k; else arg1:=op(1,select(has,op(1,RE),A))-k; end if; if nops(op(2,RE))=1 then arg2:=op(1,op(2,RE))-k; else arg2:=op(1,select(has,op(2,RE),A))-k; end if; mx:=max(arg1,arg2); mn:=min(arg1,arg2); order:=mx-mn; RE:=subs(k=k-mn,RE); # call solverec recursively if order>1 then RE:=subs(A(k+order)=A(temp+1),A(k)=A(temp),RE); # forward substitution of args (to get a recurrence of order 1) piecewisesol:=[]; for i from 0 to order-1 do RE2:=RE; RE2:=subs(k=order*k+i,RE2); # substitute coefficients of RE RE2:=subs(temp=k,RE2); # backward substitution sol:=solverec(RE2=0); if op(0,sol)=A then start:=op(1,sol); else start:=op(1,select(has,sol,A)); end if; sol:=subs(A(start)=A(order*start+i),sol); sol:=subs(k=(k-i)/order,sol); piecewisesol:=[op(piecewisesol),irem(k,order)=i,sol] end do; if nargs>1 then return subs(args[2],piecewise(op(piecewisesol))) else return piecewise(op(piecewisesol)) end if; end if; # determine the ratio A(k+order)/A(k) sol:=normal(solve(RE,A(k+order))/A(k)); rat:=simplify(expand(sol)); num:=ListTools[Flatten](convert(factor(numer(rat)),multiset)); den:=ListTools[Flatten](convert(factor(denom(rat)),multiset)); x:=1; # determine upper pochhammer symbols upp:=[];minupp:=infinity; #x:=x*op(1,num); for i from 1 by 2 to nops(num) do coe1:=coeff(op(i,num),k,1); coe0:=coeff(op(i,num),k,0); if coe1<>0 then x:=x*coe1^op(i+1,num); if type(coe0/coe1,integer) then minupp:=min(minupp,coe0/coe1) end if; for j from 1 to op(i+1,num) do upp:=[op(upp),coe0/coe1] end do; else x:=x*op(i,num)^op(i+1,num) end if; end do; # determine lower pochhammer symbols low:=[];minlow:=infinity; #x:=x/op(1,den); for i from 1 by 2 to nops(den) do coe1:=coeff(op(i,den),k,1); coe0:=coeff(op(i,den),k,0); if coe1<>0 then x:=x/coe1^op(i+1,den); if type(coe0/coe1,integer) then minlow:=min(minlow,coe0/coe1) end if; for j from 1 to op(i+1,den) do low:=[op(low),coe0/coe1] end do; else x:=x/op(i,den)^op(i+1,den) end if; end do; # determine shift content:=op(1,contentrec(RE)); shift:=min(minlow,minupp,op(map(z->-z,solve(content,[k])))); if type(shift,negint) or shift=0 then shift:=shift-1; upp:=map(y->y-shift,upp); low:=map(y->y-shift,low); else shift:=0 end if; sol:=subs(k=k+shift,A(-shift)*x^k* convert(map(y->pochhammer(y,k),upp),'`*`')/ convert(map(y->pochhammer(y,k),low),'`*`')); if nargs>1 then sol:=subs(args[2],sol) end if; sol:=simplifysol(sol); sol end proc: ########################################################### # # # 'multsumrecursion' bestimmt zu einem zulässigen Term t # # eine Rekursion minimaler Ordnung für die Summe. # # # ########################################################### ##PROCEDURE(help,testlimit=560) multsumrecursion ## ##HALFLINE determine a linear and homogeneous ## recurrence equation with polynomial coefficients for the multiple (or single) sum over a given ## hypergeometric term ## ##AUTHOR Torsten Sprenger ## ##DATE October 2006 ## ##CALLINGSEQUENCE ##- multsumrecursion('F', 'k', 'n', 'opts') ##- multsumrecursion('F', 'k', 'S(n)', 'opts') ## ##PARAMETERS ##- 'F' : algebraic expression (must be a hypergeometric or a m-fold hypergeometric term) ##- 'k' : summation variable or list of summation variables ##- 'n' : recurrence variable ##- 'S(n)' : name of the sum ##- 'opts' : sequence of options of the form 'keyword'='value'; ## possible keywords are `direction`, `hypersol`, `strategy`, ## `structureset` and `upperkbound` ## ##DESCRIPTION ##- This command tries to determine a linear and homogeneous ## recurrence equation in 'n' with polynomial coefficients for the multiple (or particularly single) sum 'S(n)' over a given ## hypergeometric term 'F(n,k)', where 'k' is a vector of 'r' summation variables. ## ##OPTIONS ##- `direction`: either `up` or `down` (default: `up`). This option defines whether the resulting recurrence ## equation contains upward or downward shifts. ##- `hypersol`: boolean (default: `false`). If set to `true` and a recurrence equation is found, then the equation ## is tried to be solved (by direct formula for recurrences of hypergeometric type and Van Hoeij algorithm). ##- `strategy`: one of `coarse`, `medium` or `fine` (default: `medium`). This option specifies, in which way the structure sets are increased. ## If `strategy` is set to `fine` normally more linear equation systems are considered. However, no system is ignored, so that, ## if the procedure proceeds, the solved system is the smallest linear equation system which leads to a solution. `coarse` ## is the counterpart, so the structure sets are increased coarsely and the number of the considered structure sets is smaller. One of the main ## disadvantages of `coarse` is that often a smaller successful equation system is ignored. If that happens, the computation of the recurrence lasts longer. ## `medium` is the default value (similar to `fine`) and should be used, because in most cases it is sufficient to get the recurrence equation. In some cases, when there could be no recurrence computed with `medium`, ## `fine` can be chosen. ##- `structureset`: one of `rectangular` or `optimal` or a specific structure set (default: `optimal`). The option `rectangular` is only for, because it is ## not efficient as `optimal`. It is also possible to enter a specific structure set in form of a set of lists, where each list consists of 'r'+1 nonnegative integers ('r' ## is the number of the summation variables plus the recurrence variable). The first entry of each list is the integer which corresponds to the recurrence variable. ##- `upperkbound`: nonnegative integer (default: `1`). The option `upperkbound` represents the upper bound of the degree of the coefficient polynomials in 'k' in the approach of the underlying algorithm. ## If 'k=0' then only 'k'-free recurrence equations are computed. ## ##SECTION(nohelp) Error conditions ##-(nolead) `multsumrecursion` may generate one of the ## following errors: ##- \"No non-trivial certificate recurrence exists\". This happens if ## the considered structure sets leads to no recurrence equation ## ##REFERENCES ##- Kurt Wegschaider, Computer Generated Proofs of Binomial Multi-Sum Identities. Diploma Thesis, 1997 ##- Torsten Sprenger, Algorithmen für mehrfache Summmen. Diplomarbeit, 2004 ## ##EXAMPLES ##> with(multsum): ##> multsumrecursion(binomial(n,k)^2, k, S(n)); ##< (-1-n)*S(n+1)+(2+4*n)*S(n) = 0 ##> multsumrecursion(binomial(n,k)^2, k, S(n), hypersol); ##< S(0)*binomial(2*n, n) ##> multsumrecursion() ## ##SEEALSO ##- "multsum" ##- "multsum[findrec]" ##- "SumTools" ##- "LREtools[hypergeomsols]" ##ENDPROCEDURE multsumrecursion:=proc(t::algebraic,k::Or(symbol,list),sn::Or(function,symbol),{structureset::Or(identical(rectangular,optimal),set):=optimal,strategy::nonnegint:=2, hypersol::boolean:=false,upperkbound::nonnegint:=1,direction::identical(up,down):=up}) local S,F,n,a,rec,sol; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; if type(sn,function) then S:=op(0,sn); n:=op(1,sn); else `tools/genglobal`('S',[],'reset'); S:=`tools/genglobal`('S',indets(t,'name')); n:=sn; end if; #_checkproper(t,k,n); try #rec:=_sumcer(findrec(t,k,n),S); rec:=_sumcer(findrec(t,k,n,F,a,_options),S); rec:=replacevar_sum(rec,S,n,a); catch: error "No non-trivial certificate recurrence exists" end try; if not hypersol then return op(2,contentrec(shiftrec(rec))) else try sol:=solverec(rec) catch: sol:=LRETools[hypergeomsols](rec,S(n),{}) end try; return sol end if; end proc: ########################################################### # # # 'multsumdiffeq' bestimmt zu einem hyperexponentiellen # # Term F eine Differentialgleichung für die Summe. # # # ########################################################### multsumdiffeq:=proc(F::algebraic,k::Or(symbol,list),sx::Or(symbol,function),{difforder::posint:=5}) local S,x,var,nvar,nstr,coe,i,j,approach,pol,diffeq,sol,polcoe,s,p,str,a; option `Copyright (c) 2006 by Torsten Sprenger, University of Kassel`; description `determines a holonomic differential equation`; if type(sx,function) then S:=op(0,sx); x:=op(1,sx) else `tools/genglobal`('S',[],'reset'); S:=`tools/genglobal`('S',indets(F,'name')); x:=sx end if; if type(k,list) then var:=k else var:=[k] end if; nvar:=nops(var); for j from 1 to nvar+1 do s[j]:=0; end do; p:=1; while p<=difforder do for j from 1 to nvar+1 do s[nvar+2-j]:=s[nvar+2-j]+1; str:=Sstd([seq(s[m],m=1..nvar+1)]); nstr:=nops(str); coe:=map(x->a[op(x)],str); approach:=add(coe[i]*simpcomb(diff(subs([seq(var[j]=var[j]+str[i,j+1],j=1..nvar)],F),[x$str[i,1]])/F),i=1..nstr); pol:=collect(numer(normal(approach)),var,distributed); polcoe:={coeffs(pol,var)}; sol:=MULTSUMSolve(polcoe,coe); if {op(subs(sol,coe))}<>{0} then #diffeq:=subs(sol,add(coe[i]*diff(S(x),[x$str[i,1]]),i=1..nstr)); diffeq:=subs(sol,add(coe[i]*S(x+str[i,1]),i=1..nstr)); diffeq:=collect(numer(normal(diffeq)),S,factor)=0; diffeq:=replacevar_sum(diffeq,S,x,a); diffeq:=subs({seq(S(x+m)=diff(S(x),x$m),m=1..p)},diffeq); return diffeq; end if; if j=nvar+1 then p:=p+1 end if: end do; end do end proc: ########################################################### # # # 'addrec' addiert zwei Rekursionsgleichungen. Dabei # # wird ein Summand eliminiert. # # # ########################################################### addrec:=proc(rec1,rec2) local i,r1,r2,cf1,cf2,F_list,var,num_IJ,F,new_rec,a; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; r1:=shiftrec(rec1,direction=down); r2:=shiftrec(rec2,direction=down); F_list:=select(type,indets(r1),function); num_IJ:=nops(F_list[1]); var:=seq(op(indets(op(i,F_list[1]))),i=1..num_IJ); F:=op(0,F_list[1]); cf1:=coeff(lhs(r1),F(var)); cf2:=coeff(lhs(r2),F(var)); new_rec:=lhs(expand(cf2*r1-cf1*r2)); if new_rec=0 then return 0=0 else return shiftrec(contentrec_(new_rec,F,var[1],a)=0) end if; end proc: ########################################################### # # # 'filter' filtert alle Faktoren von einem zulässigen # # Term aus, die nicht von den Summationsvariablen k # # abhängen. # # # ########################################################### filter:=proc(t::algebraic,k::Or(symbol,list),n::symbol) local kfree,tt; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; tt:=simpcomb(t); kfree:=select(type,tt,freeof(k)); return normal(tt/kfree); end proc: ########################################################### # # # 'filter_back' macht das Filtern in der eingegebenen # # Rekursion wieder rückgängig. # # # ########################################################### filter_back:=proc(recurrence,t::algebraic,k::Or(symbol,list),n::symbol,F::symbol) local i,m,kfree,new_rec; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; kfree:=select(type,simpcomb(t),freeof(k)); new_rec:=0; for i from 1 to nops(recurrence) do m:=op(1,select(has,op(i,recurrence),F)); new_rec:=new_rec+op(i,recurrence)/normal(simpcomb(subs(n=m,kfree)/kfree)); end do; return new_rec; end proc: ########################################################### # # # 'simplifysol' vereinfacht einen hypergeometrischen # # Term. # # # ########################################################### simplifysol:=proc(t_) local t,n,m; option `Copyright 2003-2007 Torsten Sprenger, University of Kassel`; t:=t_; t:=applyrule(pochhammer(1,m::anything)=m!,t); t:=applyrule(pochhammer(2,m::anything)=(m+1)!,t); t:=applyrule(pochhammer(1/2,m::anything)=(2*m)!/4^m/m!,t); t:=applyrule([(-1)^(m::anything)*pochhammer(n::anything,m::anything)/(m::anything)!=binomial(-n,m), (-1)^(m::anything)*pochhammer(n::anything,m::anything)/pochhammer(1,m::anything)=binomial(-n,m)],t); t:=convert(t,binomial); t:=simplify(t); t end proc: end module: #savelib('multsum');