###################################################################### ## HalfLine Save this file as HalfLine to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read HalfLine # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg at math dot rutgers dot edu. # ###################################################################### print(`First Written: Oct. 2007: tested for Maple 10 and Maple 11 `): print(`Version of Oct. 2007: `): print(): print(`This is HalfLine, A Maple package`): print(`To study walks on the half-line`): print(`accompanying Manuel Kauers and Doron Zeilberger's article: `): print(`"The Quasi-Holonomic Ansatz and Restricted Lattice Walks" `): print(): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/HalfLine .`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): #print(`For a list of the supporting functions type: ezra1();`): print(): ezra1:=proc() if args=NULL then print(`The supporting procedures are: Findrec`): print(` GenPol, GuessRec11, GuessRec2old,`): print(` GuessRec21old , GuessRecXslow, HafelOper, Mult`): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` HalfLine: A Maple package for `): print(`The MAIN procedures are: Comm`): print(`F, FundRec, GuessRec1, GuessRec2, GuessRecX, Pnx , SeqP `): print(` `): elif nargs=1 and args[1]=Comm then print(`Comm(ope1,ope2,n,a,N,A): The commutator of operators`): print(`ope1(n,a,N,A) and ope2(n,a,N,A) `): print(` [ i.e. ope1*ope2-ope2*ope1 ] `): print(` in the algebra of`): print(`linear recurrence operators in (n,a) where`): print(`N and A are the fundamental shift operators in n`): print(`and a respectively.`): print(`For example, try:`): print(`Comm(a*A,n*N,n,a,N,A); `): elif nargs=1 and args[1]=F then print(`F(n,a,Steps): Inputs: n and a: non-negative integers ; `): print(`Steps: set of integers. Output: the number of walks`): print(`using n steps from the the set Steps that wind-up at a`): print(`and never venture into the negative half-line. For example,`): print(`try: F(5,0,{1,-1});`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N,MaxC): Given a list f tries to find a `): print(`linear recurrence equation with`): print(`poly coffs.`): print(`of maximum DEGREE+ORDER<=MaxC`): print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`): elif nargs=1 and args[1]=FundRec then print(`FundRec(Steps,A,N): The fundumental recurrence operator annihilating`): print(`F(n,a,Steps) (see entry)`): print(`togetger with an integer a0 indicating `): print(`that it is valid for n>=0 and a>=a0`): print(`For example, try:`): print(`FundRec({1,-1},A,N);`): elif nargs=1 and args[1]=GenPol then print(`GenPol(Lx,d,c): a generic polynomial of degree d`): print(`in the set of variables in Lx, indexed by a`): print(`followed by the set of coefficients. For example,`): print(`try: GenPol([n,x,1/N],3,c); `): elif nargs=1 and args[1]=GuessRec1 then print(`GuessRec1(Steps,n,N,a): Given a set of steps Steps,`): print(`and variables n and a, and a symbol N`): print(` tries to guess a linear recurrence operator`): print(`with polynomial coefficients`): print(`P(n,a,N) that annihilates F(n,a,Steps) (see entry)`): print(`It also returns a0, indicating that it is valid for`): print(`a >= a0 `): print(` For example, try:`): print(`GuessRec1({1,-1},n,N,a);`): elif nargs=1 and args[1]=GuessRec11 then print(`GuessRec11(Steps,n,N,a,deg): Given a set of steps Steps,`): print(`and variables n and a, and a symbol N, and a positive`): print(`integer deg, tries to guess a linear recurrence operator`): print(`with polynomial coefficients of total degree deg,`): print(`P(n,a,1/N) that annihilates F(n,a,Steps) (see entry)`): print(`for n>=deg. For example, try:`): print(`GuessRec11({1,-1},n,N,a,4);`): elif nargs=1 and args[1]=GuessRec2 then print(`GuessRec2(Steps,n,N,a,MaxC,OrderA): Given a set of steps Steps,`): print(`and variables n and a, and a symbol N`): print(`tries to guess a linear recurrence operator,`): print(`P(n,1/A,1/N) that annihilates F(n,a,Steps) (see entry)`): print(`where A and N are the shift operators in a and n resp.`): print(`It is such that degree(P,1/N)+degree(P,n)<=MaxC `): print(`and degree(P,A)<=OrderA. For example, try:`): print(`GuessRec2({1,-1},n,N,A,4,2);`): elif nargs=1 and args[1]=GuessRec2old then print(`GuessRec2old(Steps,n,N,A): Given a set of steps Steps,`): print(`and a variable n, and symbols N and A`): print(` tries to guess a linear recurrence operator`): print(`with polynomial coefficients`): print(`P(n,1/A,1/N) that annihilates F(n,a,Steps) (see entry)`): print(` For example, try:`): print(`GuessRec2old({1,-1},n,N,A);`): elif nargs=1 and args[1]=GuessRec21old then print(`GuessRec21old(Steps,n,N,A,deg): Given a set of steps Steps,`): print(`and variable, and symbols N and A, and a positive`): print(`integer deg, tries to guess a linear recurrence operator`): print(`with polynomial coefficients of total degree deg,`): print(`P(n,1/A,1/N) that annihilates F(n,a,Steps) (see entry)`): print(`for n>=deg. For example, try:`): print(`GuessRec21old({1,-1},n,N,A,4);`): elif nargs=1 and args[1]=GuessRecX then print(`GuessRecX(Steps,x,n,N,MaxC,DegX): guesses a `): print(`linear recurrence equation with polynomial coefficients`): print(`in n and x for SeqP(Steps,x,n0) (see entry).`): print(`Such that degree(ope,n)+ORDER(ope,N)<=MaxC and`): print(`degree in x is <=DegX`): print(`For example, try:`): print(`GuessRecX({1,-1},x,n,N,4,2);`): elif nargs=1 and args[1]=GuessRecXslow then print(`GuessRecXslow(Steps,x,n,N,MaxC): guesses a `): print(`linear recurrence equation with polynomial coefficients`): print(`in n and x for SeqP(Steps,x,n0) (see entry)`): print(`For example, try:`): print(`GuessRecX({1,-1},x,n,N,4);`): elif nargs=1 and args[1]=HafelOper then print(`HafelOper(ope,Steps,n,a,N,A,n0,a0): applies the operator`): print(`ope(n,a,N,A) to F(n,a,Steps) (see entry) at n0, a0`): print(`For example, try: `): print(`HafelOper(FundRec({-1,1},A,N)[1],{-1,1},n,a,N,A,5,4);`): elif nargs=1 and args[1]=Mult then print(`Mult(ope1,ope2,n,a,N,A): The product of operators`): print(`ope1(n,a,N,A) and ope2(n,a,N,A) in the algebra of`): print(`linear recurrence operators in (n,a) where`): print(`N and A are the fundamental shift operators in n`): print(`and a respectively`): print(`For example, try:`): print(`Mult(a*A,n*N,n,a,N,A); `): elif nargs=1 and args[1]=Pnx then print(`Pnx(n,x,Steps): The weight-enumerator for steps`): print(`of length n, using steps from the set of steps Steps`): print(`staying in the non-negative half-line, with weight`): print(`x^(FinalDestination). For example, try:`): print(`Pnx(6,x,{1,-1}); `): elif nargs=1 and args[1]=SeqP then print(`SeqP(Steps,x,n0): The sequence of the first n0 terms (n=1..n0)`): print(`of Pnx(n,x,Steps). For example, try:`): print(`SeqP({1,-1},x,20);`): else print(`There is no such thing as`, args): fi: end: #F(n,a,Steps): Inputs: n and a: non-negative integers ; #Steps: set of integers. Output: the number of walks #using n steps from the the set Steps that wind-up at a #and never venture into the negative half-line. For example, #try: F(5,0,{1,-1}); F:=proc(n,a,Steps) local s: option remember: if n=0 then if a=0 then RETURN(1): else RETURN(0): fi: fi: if a<0 then RETURN(0): fi: add(F(n-1,a-s,Steps), s in Steps): end: #GenPol(Lx,d,c): a generic polynomial of degree d #in the set of variables in Lx, indexed by a #followed by the set of coefficients. For example, #try: GenPol([n,x,1/N],3,c): GenPol:=proc(Lx,d,c) local gu,kv,x,i,Lx1,ku: x:=Lx[nops(Lx)]: if nops(Lx)=1 then RETURN(add(c[i]*x^i,i=0..d),{seq(c[i],i=0..d)}): fi: Lx1:=[op(1..nops(Lx)-1,Lx)]: gu:=0: kv:={}: for i from 0 to d do ku:=GenPol(Lx1,d-i,c[i]): kv:=kv union ku[2]: gu:=expand(gu+x^i*ku[1]): od: gu,kv: end: #GuessRec11(Steps,n,N,a,deg): Given a set of steps Steps, #and variables n and a, and a symbol N, and a positive #integer deg, tries to guess a linear recurrence operator #with polynomial coefficients of total degree deg, #P(n,a,1/N) that annihilates F(n,a,Steps) (see entry) #for n>=deg. For example, try: #GuessRec11({1,-1},n,N,a,4); GuessRec11:=proc(Steps,n,N,a,deg) local ope,c,eq,var,pt,i,K,i1,eq1,var1, ind,ind1,v: ope:=GenPol([n,a,1/N],deg,c) : var:=ope[2]: ope:=ope[1]: eq:={}: for K from deg while nops(eq)0 then eq:=eq union {eq1}: fi: od: od: var1:=solve(eq,var): ope:=subs(var1,ope): if ope=0 then RETURN(0): fi: ind:={}: for v in var1 do if op(1,v)=op(2,v) then ind:=ind union {op(1,v)}: fi: od: if nops(ind)>1 then print(`Too much slack`): RETURN({seq(coeff(ope,ind1,1),ind1 in ind)}): fi: subs(ind[1]=1,ope): end: #GuessRec1(Steps,n,N,a): Given a set of steps Steps, #and variables n and a, and a symbol N #tries to guess a linear recurrence operator #P(n,a,N) that annihilates F(n,a,Steps) (see entry) #For example, try: #GuessRec1({1,-1},n,N,a); GuessRec1:=proc(Steps,n,N,a) local deg,ope,lod,n0,a0,maspik: for deg from 1 do ope:=GuessRec11(Steps,n,N,a,deg): if ope<>0 then lod:=-ldegree(ope,N): ope:=Yafe(subs(n=n+lod,ope)*N^(lod),N): maspik:=min(op(Steps))-2: for a0 from max(op(Steps)) to min(op(Steps))-1 by -1 do if {seq(HafelOper(ope,Steps,n,a,N,A,n0,a0),n0=0..10)}<>{0} then maspik:=a0: break: fi: od: RETURN(ope,maspik): fi: od: end: #GuessRec21old(Steps,n,N,A,deg): Given a set of steps Steps, #and a variable n , and symbols N and A for #the shift opertors in n and a resp. and a positive #integer deg, tries to guess a linear recurrence operator #with polynomial coefficients of total degree deg, #P(n,a,1/N) that annihilates F(n,a,Steps) (see entry) #for n>=deg. For example, try: #GuessRec21old({1,-1},n,N,A,4); GuessRec21old:=proc(Steps,n,N,A,deg) local ope,c,eq,var,pt,i,K,i1,eq1,var1, ind,ind1,v,i2: ope:=GenPol([n,1/A,1/N],deg,c) : var:=ope[2]: ope:=ope[1]: eq:={}: for K from 2*deg while nops(eq)0 then eq:=eq union {eq1}: fi: od: od: var1:=solve(eq,var): ope:=subs(var1,ope): if ope=0 then RETURN(0): fi: ind:={}: for v in var1 do if op(1,v)=op(2,v) then ind:=ind union {op(1,v)}: fi: od: if nops(ind)>1 then print(`Too much slack`): RETURN({seq(coeff(ope,ind1,1),ind1 in ind)}): fi: subs(ind[1]=1,ope): end: #GuessRec2old(Steps,n,N,a): Given a set of steps Steps, #and variables n and a, and a symbol N #tries to guess a linear recurrence operator #P(n,1/A,1/N) that annihilates F(n,a,Steps) (see entry) #where A and N are the shift operators in a and n resp. #For example, try: #GuessRec2old({1,-1},n,N,A); GuessRec2old:=proc(Steps,n,N,A) local deg,ope: for deg from 1 do ope:= GuessRec21old(Steps,n,N,A,deg): if ope<>0 then RETURN(ope): fi: od: ope: end: #Pnx(n,x,Steps): The weight-enumerator for steps #of length n, using steps from the set of steps Steps #staying in the non-negative half-line, with weight #x^(FinalDestination). For example, try: #Pnx(6,x,{1,-1}): Pnx:=proc(n,x,Steps) local s,gu,i: option remember: if n=0 then RETURN(1): else gu:=Pnx(n-1,x,Steps): gu:= expand(add(x^s *(gu-add(coeff(gu,x,i)*x^i,i=0..-s-1)), s in Steps)): fi: gu: end: #SeqP(Steps,x,n0): The sequence of the first n0 terms (n=1..n0) #of Pnx(n,x,Steps). For example, try: #SeqP({1,-1},x,20); SeqP:=proc(Steps,x,n0) local i: [seq(Pnx(i,x,Steps),i=1..n0)]: end: ################Start Findrec and its offsprings with(linalg): #with(SolveTools); #findrecEx(f,DEGREE,ORDER,m1): Explores whether thre #is a good chance that there is a recurrence of degree DEGREE #and order ORDER, using the prime m1 #For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003); findrecEx:=proc(f,DEGREE,ORDER,m1) local ope,var,eq,i,j,n0,eq1,a,A1, D1,E1,Eq,Var,f1,n,N: option remember: f1:=f mod m1: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1: od: eq:= eq union {eq1}: od: Eq:= convert(eq,list): Var:= convert(var,list): D1:=nops(Var): E1:=nops(Eq): if E10 then RETURN(false): fi: if E1-D1>=1 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+1],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: if E1-D1>=2 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+2],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: true: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`More than one solution`): RETURN(Yafe(ope[1],N)): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)): else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i: add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)): end: #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 1 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: ################End Findrec and its offsprings #GuessRecXslow(Steps,x,n,N,MaxC): guesses a #linear recurrence equation with polynomial coefficients #in n and x for SeqP(Steps,x,n0) (see entry) #For example, try: #GuessRecX({1,-1},x,n,N,4); GuessRecXslow:=proc(Steps,x,n,N,MaxC) local gu,gu3,gu5,ope3,ope5,n0,DEG,ORD: n0:=(2+MaxC/2)*(1+MaxC/2)+10: gu:=SeqP(Steps,x,n0): gu3:=subs(x=4/3,gu): ope3:=Findrec(gu3,n,N,MaxC): if ope3=FAIL then RETURN(FAIL): fi: gu5:=subs(x=11/7,gu): ope5:=Findrec(gu5,n,N,MaxC): if ope3=FAIL then RETURN(FAIL): fi: if degree(numer(ope3),n)<>degree(numer(ope5),n) then print(`Can't agree on degree `): RETURN(FAIL): fi: if degree(ope3,N)<>degree(ope5,N) then print(`Can't agree on degree ORDER`): RETURN(FAIL): fi: DEG:=degree(numer(ope5),n): ORD:=degree(ope5,N): findrec(gu,DEG,ORD,n,N): end: ##################Begin GR1 with(SolveTools): #Bdok1(DS1,P,x): checks that P in the single variable x #indeed fits the data set DS1 Bdok1:=proc(DaS,P,x) local i: for i from 1 to nops(DaS) do if subs(x=DaS[i][1],denom(P))<>0 and normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then RETURN(false): fi: od: true: end: #GenRat1(x,dxt,dxb,a): a generic (monic top) #rational function of degree of top dxt #and degree of bottom dxb #with coeffs.parametrized by a[i] followed by the set #of coeffs. GenRat1:=proc(x,dxt,dxb,a) local i: (add(a[i]*x^i,i=0..dxt))/ add(a[i+dxt+1]*x^i,i=0..dxb) ,{seq(a[i],i=0..dxt+dxb+1)}: end: #GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2 #that fits the data in DS1 GR1d1d2:=proc(S1,x,d1,d2) local eq,var,P,a,var1,i: P:=GenRat1(x,d1,d2,a): var:=P[2]: P:=P[1]: if nops(S1)<=d1+d2+3 then # print(`Insufficient data`): RETURN(FAIL): fi: eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}: var1:=Linear(eq,var): if expand(subs(var1,numer(P)))=0 then RETURN(FAIL): fi: P:=subs(var1,P): if Bdok1(S1,P,x) then RETURN(factor(normal(P))): else RETURN(FAIL): fi: end: #GR1(S1,x): guesses a rational function from the data-set S1 GR1:=proc(S1,x) local d1,d2,gu,L,i1: if {seq(S1[i1][2],i1=1..nops(S1))}={0} then RETURN(0): fi: for L from 0 to nops(S1)-4 do for d1 from 0 to L do d2:=L-d1: gu:=GR1d1d2(S1,x,d1,d2): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: ##################End GR1 #FindrecX(f,MaxC,x,n,N,DegX): #guesses a recurrence operator of complexity <=MaxC, #in n, and DegX in x #annihilating the sequence f of polynomials in x . #For example, try: FindrecX([seq(x^i,i=1..20)],2,x,n,N,2); FindrecX:=proc(f,MaxC,x,n,N,DegX) local gu,ope0,S, Deg,Ord,i,S0,j,ope,vu: gu:=subs(x=3/4,f): ope0:=Findrec(gu,n,N,MaxC): if ope0=FAIL then RETURN(FAIL): fi: S:={[3/4,ope0]}: Deg:=degree(numer(normal(ope0)),n): Ord:=degree(ope0,N): for i from 1 to DegX+6 do ope0:=findrec(subs(x=i+1/2,f),Deg,Ord,n,N): if ope0<>FAIL then S:=S union {[i+1/2,ope0]}: fi: od: ope:=0: for i from 0 to Ord do S0:=[seq([S[j][1],coeff(S[j][2],N,i)],j=1..nops(S))]: vu:=GR1(S0,x); if vu=FAIL then RETURN(FAIL): else ope:=ope+vu*N^i: fi: od: ope: end: #GuessRecX(Steps,x,n,N,MaxC,DegX): guesses a #linear recurrence equation with polynomial coefficients #in n of degree in and N <=MaxC amd and degree DetX #in x for SeqP(Steps,x,n0) (see entry) #For example, try: #GuessRecX({1,-1},x,n,N,4,2); GuessRecX:=proc(Steps,x,n,N,MaxC,DegX) local gu,n0: n0:=(2+MaxC/2)*(1+MaxC/2)+10: gu:=SeqP(Steps,x,n0): gu:=FindrecX(gu,MaxC,x,n,N,DegX,DegX): gu:=numer(gu): collect(gu,N): end: #GuessRec2(Steps,n,N,a,MaxC,OrderA): Given a set of steps Steps, #and variables n and a, and a symbol N #tries to guess a linear recurrence operator, #P(n,1/A,1/N) that annihilates F(n,a,Steps) (see entry) #where A and N are the shift operators in a and n resp. #It is such that degree(P,1/N)+degree(P,n)<=MaxC and degree(P,A)<=OrderA #For example, try: #GuessRec2({1,-1},n,N,A,4,2); GuessRec2:=proc(Steps,n,N,A,MaxC,OrderA) local x,ope,a,a0,n0,maspik: ope:=GuessRecX(Steps,x,n,N,MaxC,OrderA): ope:=numer(subs(x=1/A,ope)): maspik:=min(op(Steps))-6: for a0 from max(op(Steps)) to min(op(Steps))-5 by -1 do if {seq(HafelOper(ope,Steps,n,a,N,A,n0,a0),n0=0..10)}<>{0} then maspik:=a0: break: fi: od: RETURN(collect(ope,A),maspik): end: #FundRec(Steps,A,N): The fundumental recurrence operator annihilating #F(n,a,Steps) (see entry) #togetger with an integer a0 indicating where it #that it is valid for n>=0 and a>=a0 #For example, try: #FundRec({1,-1},A,N); FundRec:=proc(Steps,A,N) local s,mi: mi:=min(op(Steps)): numer(expand(1-1/N*add(A^(-s),s in Steps))),mi: end: #HafelOper(ope,Steps,n,a,N,A,n0,a0): applies the operator #ope(n,a,N,A) to F(n,a,Steps) (see entry) at n0, a0 #For example, try: #HafelOper(FundRec({-1,1},A,N)[1],{-1,-1},n,a,N,A,5,4); HafelOper:=proc(ope,Steps,n,a,N,A,n0,a0) local i,j: add(add(subs({n=n0,a=a0},coeff(coeff(ope,N,i),A,j) )*F(n0+i,a0+j,Steps), i=ldegree(ope,N)..degree(ope,N)),j=ldegree(ope,A)..degree(ope,A)): end: #Mult(ope1,ope2,n,a,N,A): The product of operators #ope1(n,a,N,A) and ope2(n,a,N,A) in the algebra of #linear recurrence operators in (n,a) where #N and A are the fundamental shift operators in n #and a respectively #For example, try: #Mult(a*A,n*N,n,a,N,A): Mult:=proc(ope1,ope2,n,a,N,A) local i,j,Ope1,Ope2,Ope,mu: Ope1:=expand(ope1): Ope2:=expand(ope2): Ope:=0: for i from ldegree(Ope1,A) to degree(Ope1,A) do mu:=expand(coeff(Ope1,A,i)): for j from ldegree(mu,N) to degree(mu,N) do Ope:=expand(Ope+coeff(mu,N,j)*subs({a=a+i,n=n+j},Ope2)*A^i*N^j): od: od: Ope: end: #Comm(ope1,ope2,n,a,N,A): The commutator of operators #ope1(n,a,N,A) and ope2(n,a,N,A) (ope1*ope2-ope2*ope1) #in the algebra of #linear recurrence operators in (n,a) where #N and A are the fundamental shift operators in n #and a respectively #For example, try: #Comm(a*A,n*N,n,a,N,A): Comm:=proc(ope1,ope2,n,a,N,A) : expand(Mult(ope1,ope2,n,a,N,A)-Mult(ope2,ope1,n,a,N,A)): end: