##################################################################### ## This program is written to study the Generalized Rondom Walk # ## also known as Correlated Rondom walk. The first version of # ## the program was written as a class excercise in Dr. Z's # ## Bioinformatics class of 2006. This version contains more # ## functions and prepares the ground for further investigation. # ## tHIS VERSION OF mY 2006. # ## # ## CRW: Save this file as CRW To use it, stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: CRW: # ## Then follow the instructions given there # ## # ## Written by Moa Apagodu # ## mohamudm@math.rutgers.edu # ##################################################################### with(SolveTools): _EnvAllSolutions:=true: print(`Tested for Maple 8, 9, and 10`): print(): print(`This is CRW, A Maple package`): print(`written by Moa Apagodu`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(): ezra:=proc() if args=NULL then print(`This program is written to study the Generalized Rondom`): print(`Walk also known as Correlated Rondom walk. The first version`): print(`of the program was written as a class excercise in Dr. Z`): print(`Bioinformatics class of 2006. This is improved`): print(`version of May 2006.`): print(`The main procedures are `): print(` GFOne, GFOnep, GFOnepRA, GFTwo, GFTwop, GFTwopRA, Asy, MAZ`): print(): elif nargs=1 and args[1]=Asy then print(`Asy(ope,n,N,K): the asymptotic expansion of solutions `): print(`to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator`): print(`up to the K's term`): elif nargs=1 and args[1]=MAZ then print(`MAZ(POL,H,rat,x,n,N,pars): Given a polynomial, POL, an hyper-exponential function, H, and`): print(`a rational function, rat, of the continuous variables in the list x, and a discrete variable`): print(`n, and the shift-operator in n, N, and a set of auxiliary parameters, par`): print(`ouputs a recurrence operator, let's call it ope(n,n), and a multi-certificate, let's call it R`): print(`such that the function F(n;x):=POL*H*rat^n satisfies`): print(` ope(n,N)F=sum(D_{x[i]} (R[i]H*rat^n),i=1..nops(R)`): print(`In particular, try:`): print(`MAZ(1,((1-1/x)*(1-y))^b*((1-1/x/y)*(1-1/y))^c/x/y,(1-x)*(1-x*y),[x,y],a,A,{b,c});`): print(`MAZ(1,exp(-x^2/2-y^2/2),(x-y)^2,[x,y],n,N,{});`): elif nargs=1 and args[1]=GFOne then print(`GFOne(x,t): the weight-enumerator of the set`): print(`of all words in {[1,0],[-1,0]} where the weight of`): print(`a word is t^lenght(w)*x^(sum of the entry)`): print(`# This outputs the generating function SUM(SUM(a_mx^mt^n,n=0..infinity),m=-infinity..infinity).`): print(`For example try the Polya case: GFOne(x,t)`): elif nargs=1 and args[1]=GFOnep then print(`GFOnep(x,t,a,p,q): the probability weight-enumerator of the set`): print(`Given that a=[a1,a2] initial prob., p=[p1,p2],q=[q1,q2]`): print(` are the adv. and retreat probabilty in `): print(`the +ve x and -ve x-direc respectively,`): print(`of all words in {[1,0],[-1,0]} where the weight of`): print(`a word is t^lenght(w)*x^(sum of the entry)`): print(`# This outputs the generating function SUM(SUM(P(m,k)x^my^kt^n,n=0..infinity),m,k=-infinity..infinity).`): print(`For example try the Polya case: GFOnep(x,t,[1/2,1/2],[1/2,1/2],[1/2,1/2])`): elif nargs=1 and args[1]=GFOnepRA then print(`GFOnepRA(px,a,p,q,K): the probability weight-enumerator of the set`): print(`of all words in {[1,0],[-1,0]} where the weight of`): print(`Given that a=[a1,a2] initial prob., p=[p1,p2],q=[q1,q2]`): print(` are the adv. and retreat probabilty in `): print(`the +ve x and -ve x-directions respectively. The weight of`): print(`a word is t^lenght(w)*x^(sum of the entry)`): print(`This outputs the recurrence satisfied by P(2n,px), and the asymptotic of order K`): print(`For example try the Polya's case: GFOnepRA(1,[1/2,1/2],[1/2,1/2],[1/2,1/2],2)`): elif nargs=1 and args[1]=GFTwo then print(`GFTwo(x,y,t): the weight-enumerator of the set`): print(`of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of`): print(`a word is t^lenght(w)*x^(current nubmber of #D's)*y^(#of E's)`): print(`For example try the Polya's case: GFTwo(x,y,t)`): elif nargs=1 and args[1]=GFTwop then print(`GFTwop(x,y,t,p): the probability weight-enumerator of the set`): print(`of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of`): print(`a word is t^lenght(w)*x^(current nubmber of #D's)*y^(#of E's)`): print(`For example try the Polya's case: GFTwopRA(x,y,t,[1/4,1/4,1/4,1/4]$5)`): elif nargs=1 and args[1]=GFTwopRA then print(`GFTwopRA(px,py,a,p): the probability weight-enumerator of the set`): print(`of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of`): print(`a word is t^lenght(w)*x^(current nubmber of #D's)*y^(#of E's)`): print(`The out put is the recurrence satisfied by P(2n,(px,py)), Asumptotic of order K`): print(`For example try the Polya's case: GFTwopRA(1,1,[1/4,1/4,1/4,1/4]$5,2)`): elif nargs=1 and args[1]=GFThree then print(`GFThree(x,y,z,t): the weight-enumerator of the set`): print(`of all words in {[1,0,0],[-1,0,0],[0,-1,0],[0,1,0],[0,0,1],[0,0,-1]}`): print(`where the weight of a word is t^lenght(w)*x^(current nubmber of #D)*y^(#of E's)*z^(# of O's)`): elif nargs=1 and args[1]=GFThreep then print(`GFThree(x,y,z,t): the weight-enumerator of the set`): print(`of all words in {[1,0,0],[-1,0,0],[0,-1,0],[0,1,0],[0,0,1],[0,0,-1]}`): print(`where the weight of a word is t^lenght(w)*x^(current nubmber of #D)*y^(#of E's)*z^(# of O's)`): elif nargs=1 and args[1]=GFThreepRA then print(`GFThreepRA(px,py,pzt,a,p,K): the probability weight-enumerator of the set`): print(`of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of`): print(`a word is t^lenght(w)*x^(current nubmber of #D's)*y^(#of E's)*z^(# of B's)`): print(`The out put is the recurrence satisfied by P(2n,(px,py,pz)), Asumptotic of degree K`): print(`For example try the Polya's case: GFThreepRA(1,1,1,[1/6,1/6,1/6,1/6,1/6,1/6],[[1/6,1/6,1/6]$6],2)`): fi: end: ##############Begining of the mian program ######################## #GFOne(x,t,p,q): the weight-enumerator of the set #of all words in {[1,0],[-1,0]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #) GFOne:=proc(x,t,p,q) local fD,fmD,sol: sol:= solve({ fD=t*x+ fD*t*x+ fmD*t*x, fmD=t/x+ fD*t/x+ fmD*t/x }, {fD,fmD}); factor(normal(1+subs(sol,fD)+subs(sol,fmD) )); end: #GFOnep(x,t,a,p,q): the probabilty weight-enumerator of the set #of all words in {[1,0],[-1,0]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #) #Moa's hopeful colleague's way GFOnep:=proc(x,t,a,p,q) local fD,fmD,sol: sol:= solve({ fD=t*x*a[1]+ fD*t*x*p[1]+ fmD*t*x*q[2], fmD=t/x*a[2]+ fD*t/x*p[2]+ fmD*t/x*q[1] }, {fD,fmD}); factor(normal(1+subs(sol,fD)+subs(sol,fmD) )); end: #GFOnepRA(px,a,p,q,K): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) #a,p,q are a list of lists of four entry each. GFOnepRA:=proc(px,a,p,q,K) local ff,ope,sol,x,t: sol:= solve({ fD=t*x*a[1]+ fD*t*x*p[1]+ fmD*t*x*q[2], fmD=t/x*a[2]+ fD*t/x*p[2]+ fmD*t/x*q[1] }, {fD,fmD}); ff:=factor(normal(1+subs(sol,fD)+subs(sol,fmD) )); ope:=MAZ(1,ff/x^(px+1)/t,1/t^2,[x,t],n,N,{})[1]: ope,Asy(ope,n,N,K): end: ###################################### END of 1-D ############## #GFTwo(x,y,t): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) GFTwo:=proc(x,y,t) local fD,fmD,fE,fmE,sol: sol:= solve({ fD=t*x+ fD*t*x+ fmD*t*x+fE*t*x+fmE*t*x, fmD=t/x+ fD*t/x+ fmD*t/x+fE*t/x+fmE*t/x, fE=t*y+ fD*t*y+ fmD*t*y+fE*t*y+fmE*t*y, fmE=t/y+ fD*t/y+ fmD*t/y+fE*t/y+fmE*t/y }, {fD,fmD,fE,fmE}); factor(normal(1+subs(sol,fD)+subs(sol,fmD)+ subs(sol,fE)+subs(sol, fmE) )); end: #GFTwop(x,y,t,a,p,q,r,s): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) #a,p,q,r,s are a list of lists of four entry each. GFTwop:=proc(x,y,t,a,p,q,r,s) local fD,fmD,fE,fmE,sol: sol:= solve({ fD=t*x*a[1]+ fD*t*x*p[1]+ fmD*t*x*p[2]+fE*t*x*p[3]+fmE*t*x*p[4], fmD=t/x*a[2]+ fD*t/x*q[1]+ fmD*t/x*q[2]+fE*t/x*q[3]+fmE*t/x*q[4], fE=t*y*a[3]+ fD*t*y*r[1]+ fmD*t*y*r[2]+fE*t*y*r[3]+fmE*t*y*r[4], fmE=t/y*a[4]+ fD*t/y*s[1]+ fmD*t/y*s[2]+fE*t/y*s[3]+fmE*t/y*s[4] }, {fD,fmD,fE,fmE}); factor(normal(1+subs(sol,fD)+subs(sol,fmD)+ subs(sol,fE)+subs(sol, fmE) )); end: #GFTwopRA(px,py,a,p,q,r,s,K): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) #a,p,q,r,s are a list of lists of four entry each. GFTwopRA:=proc(px,py,a,p,q,r,s,K) local fD,fmD,fE,fmE,sol,ff,ope,x,y,t: sol:= solve({ fD=t*x*a[1]+ fD*t*x*p[1]+ fmD*t*x*p[2]+fE*t*x*p[3]+fmE*t*x*p[4], fmD=t/x*a[2]+ fD*t/x*q[1]+ fmD*t/x*q[2]+fE*t/x*q[3]+fmE*t/x*q[4], fE=t*y*a[3]+ fD*t*y*r[1]+ fmD*t*y*r[2]+fE*t*y*r[3]+fmE*t*y*r[4], fmE=t/y*a[4]+ fD*t/y*s[1]+ fmD*t/y*s[2]+fE*t/y*s[3]+fmE*t/y*s[4] }, {fD,fmD,fE,fmE}); ff:=factor(normal(1+subs(sol,fD)+subs(sol,fmD)+ subs(sol,fE)+subs(sol, fmE) )); ope:=MAZ(1,ff/x^(px+1)/t/y^(py+1),1/t^2,[x,y,t],n,N,{})[1]: ope,Asy(ope,n,N,K): end: ###################################### END of 2-D ################## #GFThree(x,y,z,t): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) GFThree:=proc(x,y,z,t) local fD,fmD,fE,fmE,fO,fmO,sol: sol:= solve({ fD=t*x+ fD*t*x+ fmD*t*x+fE*t*x+fmE*t*x+fO*t*x+fmO*x*t, fmD=t/x+ fD*t/x+ fmD*t/x+fE*t/x+fmE*t/x+fO*t/x+fmO*t/x, fE=t*y+ fD*t*y+ fmD*t*y+fE*t*y+fmE*t*y+fO*t*y+fmO*y*t, fmE=t/y+ fD*t/y+ fmD*t/y+fE*t/y+fmE*t/y+fO*t/y+fmO*t/y, fO=t*z+fD*t*z+ fmD*t*z+fE*t*z+fmE*t*z+fO*t*z+fmO*z*t, fmO=t/z+fD*t/z+ fmD*t/z+fE*t/z+fmE*t/z+fO*t/z+fmO*t/z }, {fO,fmO,fD,fmD,fE,fmE}); factor(normal(1+subs(sol,fD)+subs(sol,fmD)+ subs(sol,fE)+subs(sol, fmE)+subs(sol,fO)+subs(sol, fmO) )); end: #GFThreep(x,y,z,t,a,pro): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) #Moa's hopeful colleague's way GFThreep:=proc(x,y,z,t,a,pro) local fD,fmD,fE,fmE,sol,p1,p2,q1,q2,r1,r2: p1:=pro[1]:p2:=pro[2]: #x-direction q1:=pro[3]:q2:=pro[4]: #y-direction r1:=pro[5]:r2:=pro[6]: #z-direction sol:= solve({ fD=t*x*a[1]+ fD*t*x*p1[1]+ fmD*t*x*p2[2]+fE*t*x*q1[3]+fmE*t*x*q2[3]+fO*t*x*r1[3]+fmO*x*t*r2[3], fmD=t/x*a[2]+ fD*t/x*p1[2]+ fmD*t/x*p2[1]+fE*t/x*q1[3]+fmE*t/x*q2[3]+fO*t/x*r1[3]+fmO*t/x*r2[3], fE=t*y*a[3]+ fD*t*y*p1[3]+ fmD*t*y*p2[3]+fE*t*y*q1[1]+fmE*t*y*q2[2]+fO*t*y*r1[3]+fmO*y*t*r2[3], fmE=t/y*a[4]+ fD*t/y*p1[3]+ fmD*t/y*p2[3]+fE*t/y*q2[2]+fmE*t/y*q2[1]+fO*t/y*r1[3]+fmO*t/y*r2[3], fO=t*z*a[5]+fD*t*z*p1[3]+ fmD*t*z*p2[3]+fE*t*z*q1[3]+fmE*t*z*q2[3]+fO*t*z*r1[1]+fmO*z*t*r2[2], fmO=t/z*a[6]+fD*t/z*p1[3]+ fmD*t/z*p2[3]+fE*t/z*q1[3]+fmE*t/z*q2[3]+fO*t/z*r1[2]+fmO*t/z*r2[1] }, {fO,fmO,fD,fmD,fE,fmE}); factor(normal(1+subs(sol,fD)+subs(sol,fmD)+ subs(sol,fE)+subs(sol, fmE)+subs(sol,fO)+subs(sol, fmO) )); end: #GFThreepRA(px,py,pz,a,pro,K): the weight-enumerator of the set #of all words in {[1,0],[-1,0],[0,-1],[0,1]} where the weight of #a word is t^lenght(w)*x^(current nubmber of #)*y^(#of E's) #Moa's hopeful colleague's way GFThreepRA:=proc(px,py,pz,a,pro,K) local fD,fmD,fE,fmE,sol,x,y,z,t, p1,p2,q1,q2,r1,r2,ff,ope: p1:=pro[1]:p2:=pro[2]: q1:=pro[3]:q2:=pro[4]: r1:=pro[5]:r2:=pro[6]: sol:= solve({ fD=t*x*a[1]+ fD*t*x*p1[1]+ fmD*t*x*p2[2]+fE*t*x*q1[3]+fmE*t*x*q2[3]+fO*t*x*r1[3]+fmO*x*t*r2[3], fmD=t/x*a[2]+ fD*t/x*p1[2]+ fmD*t/x*p2[1]+fE*t/x*q1[3]+fmE*t/x*q2[3]+fO*t/x*r1[3]+fmO*t/x*r2[3], fE=t*y*a[3]+ fD*t*y*p1[3]+ fmD*t*y*p2[3]+fE*t*y*q1[1]+fmE*t*y*q2[2]+fO*t*y*r1[3]+fmO*y*t*r2[3], fmE=t/y*a[4]+ fD*t/y*p1[3]+ fmD*t/y*p2[3]+fE*t/y*q2[2]+fmE*t/y*q2[1]+fO*t/y*r1[3]+fmO*t/y*r2[3], fO=t*z*a[5]+fD*t*z*p1[3]+ fmD*t*z*p2[3]+fE*t*z*q1[3]+fmE*t*z*q2[3]+fO*t*z*r1[1]+fmO*z*t*r2[2], fmO=t/z*a[6]+fD*t/z*p1[3]+ fmD*t/z*p2[3]+fE*t/z*q1[3]+fmE*t/z*q2[3]+fO*t/z*r1[2]+fmO*t/z*r2[1] }, {fO,fmO,fD,fmD,fE,fmE}); ff:=factor(normal(1+subs(sol,fD)+subs(sol,fmD)+ subs(sol,fE)+subs(sol, fmE)+subs(sol,fO)+subs(sol, fmO) )); ope:=MAZ(1,ff/x^(px+1)/t/y^(py+1)/z^(pz+1),1/t^2,[x,y,z,t],n,N,{})[1]: ope,Asy(ope,n,N,K): end: ######################### Asy from GuessHolo ################# #Aope1(ope,n,N): the asymptotics of the difference operator #with poly. coeffs. ope(n,N) Aope1:=proc(ope,n,N) local gu: gu:=expand(numer(normal(ope))): gu:=coeff(gu,n,degree(gu,n)): gu:=expand(gu/coeff(gu,N,degree(gu,N))): factor(gu): end: C3:=proc(a,b,c) option remember: if [a,b,c]=[0,0,0] then RETURN(1): fi: if a<0 or b<0 or c<0 then RETURN(0): fi: C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1): end: #Aluf1(pit): is it positive dominant? Aluf1:=proc(pit) local aluf,shi,i,makom: shi:=evalf(abs(pit[1])): aluf:={evalf(evalc(pit[1]))}: makom:={1}: for i from 2 to nops(pit) do if evalf(abs(pit[i]))>evalf(shi) then shi:=abs(evalf(evalc(pit[i]))): aluf:={pit[i]}: makom:={i}: elif evalf(abs(pit[i]))=shi then aluf:=aluf union {evalf(evalc(pit[i]))}: makom:=makom union {i}: fi: od: aluf,shi,makom: end: #OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepAS1:=proc(ope1,n,N,alpha,f,S1) local x,f1,L,F,A,mu,i,A1: f1:=subs(n=1/x,f): L:=degree(f1,x): F:=f1+A*x^(L+S1): mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F) ,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): mu:=taylor(mu,x=0,L+11): mu:=simplify(mu): for i from L+1 to L+9 while coeff(mu,x,i)=0 do od: if i=L+10 then RETURN(FAIL): fi: mu:=coeff(mu,x,i): A1:=subs(Linear({mu},{A}),A): if A1=0 then RETURN(FAIL): fi: subs({x=1/n,A=A1},F): end: #OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepA:=proc(ope1,n,N,alpha,f) local S1: for S1 from 1 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od: OneStepAS1(ope1,n,N,alpha,f,S1): end: #Asy(ope,n,N,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the K's term Asy:=proc(ope,n,N,K) local gu,pit,lu,makom,alpha,mu,ope1,ku,i,f,x,ka: gu:=Aope1(ope,n,N): if degree(gu,N)=0 then print(`Not of exponential growth, case not implemented`): RETURN(FAIL): fi: if not type(expand(normal(ope/gu)),`+`) then pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: RETURN(pit[makom[1]]^n*expand((nops(makom)-1)!*binomial(nops(makom)-1+n,nops(makom)-1))): fi: pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Not only that but Dominant roots are complex`): RETURN(FAIL): fi: fi: lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Dominant root is complex`): RETURN(FAIL): fi: lu:=pit[makom[1]]: ope1:=Yafeasy(subs(N=lu*N,ope),N)[2]: mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,nops(makom)+1),x,nops(makom))): ka:=simplify([solve(ku,alpha)]): alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN(lu^n*n^alpha): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: lu^n*n^alpha*f: end: Yafeasy:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: ########################### End of Asy ################################# ############################### AlmkvistZeilberger algorithm ############## #IV1(d,n): all the vectors of non-negative integres of length d whose #sum is n IV1:=proc(d,n) local gu,i,gu1,i1: if d=0 then if n=0 then RETURN({[]}): else RETURN({}): fi: fi: gu:={}: for i from 0 to n do gu1:=IV1(d-1,n-i): gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}: od: gu: end: yafe:=proc(ope,N) local i: add(factor(coeff(ope,N,i))*N^i,i=ldegree(ope,N)..degree(ope,N)):end: #IV(d,n): all the vectors of non-negative integres of length d whose #sum is <=n IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=0..n)}: end: #GenPol(kList,a,deg): The generic polynomial in kList of #degree deg, using the indexed variable a, followed by the set #of coeffs. GenPol:=proc(kList,a,deg) local gu,i,i1: gu:=IV(nops(kList),deg): convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`), i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}: end: #Extract1(POL,kList): extracts all the coeffs. of a POL in kList Extract1:=proc(POL,kList) local k1,kList1,i: k1:=kList[nops(kList)]: kList1:=[op(1..nops(kList)-1,kList)]: if nops(kList)=1 then RETURN({seq(coeff(POL,k1,i),i=0..degree(POL,k1))}): else RETURN({seq( op(Extract1(coeff(POL,k1,i),kList1)), i=0..degree(POL,k1))}): fi: end: FindDeg:=proc(POL,H,rat,x,xSet,n,L) local s,t,h,e,i,Hbar,q,r,gu: s:=numer(rat): t:=denom(rat): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=normal(simplify(diff(Hbar,x)/Hbar)): q:=numer(gu): r:=denom(gu): max(degree(h,xSet)-degree( diff(r,x)+q,xSet)+degree(r,xSet) ,degree(h,xSet)-degree(r,xSet)+1+degree(r,xSet)): end: #MAZ1(POL,H,rat,x,n,N,L,pars): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1:=proc(POL,H,rat,x,n,N,L,pars) local deg,deg1,m,gu,i,i1: deg:=[]: for i from 1 to nops(x) do deg:=[op(deg),FindDeg(POL,H,rat,x[i],{seq(x[i1],i1=1..nops(x))},n,L)]: od: m:=min(op(deg)): for i from 0 to m do deg1:=[seq(deg[i1]-m+i,i1=1..nops(deg))]: gu:=MAZ1deg(POL,H,rat,x,n,N,L,pars,deg1): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: MAZ:=proc(POL,H,rat,x,n,N,pars) local gu,L: for L from 0 do print(`trying L=`,L): gu:=MAZ1(POL,H,rat,x,n,N,L,pars): if gu<>FAIL then RETURN(gu): fi: od: end: CheckMAZ:=proc(POL,H,rat,x,n,N,pars) local F,gu,luL,luR,R,ope,i: F:=H*rat^n: gu:=MAZ(POL,H,rat,x,n,N,pars): if gu=FAIL then RETURN(FAIL): fi: ope:=gu[1]: R:=gu[2]: luL:=normal(add(subs(n=n+i,POL)*coeff(ope,N,i)*normal(simplify(subs(n=n+i,F)/F)),i=0..degree(ope,N))): luR:=normal(add( normal(diff(R[i]*F,x[i])/F ), i=1..nops(R))): evalb(normal(luL/luR)=1): end: #MAZ1deg(POL,H,rat,x,n,N,L,pars,deg): Given a polynomial POL of the discrete variable n #and the continuous variables in x #and a pure-hyperexponential function H of the variables in x, and a rational function #Rat of x, and a shift operator N, and a non-neg. integer L #returns FAIL or, if there is one, an operator ope(N,n) of order L and #a list of rational functions R, such that F:=POL*H*rat^n satisfies #ope(N,n)F=sum(diff(F*R[i],x[i]),i=1..nops(R)) . #For example, try MAZ1(1,exp(-x),x,[x],n,N,1,{}) MAZ1deg:=proc(POL,H,rat,x,n,N,L,pars,deg) local s,t,h,e,i,Hbar,gu,X,a,var,q,r,eq,ope,var1,ku,i1,eqN,var1N,opeN,bilti,meka: s:=numer(rat): t:=denom(rat): ope:=add(e[i]*N^i,i=0..L): h:=convert([seq(e[i]*subs(n=n+i,POL)*s^i*t^(L-i),i=0..L)],`+`): Hbar:=H*s^n/t^(n+L): gu:=[]: for i from 1 to nops(x) do gu:=[op(gu),normal(simplify(diff(Hbar,x[i])/Hbar))]: od: q:=[]: r:=[]: for i from 1 to nops(gu) do q:=[op(q),numer(gu[i])]: r:=[op(r),denom(gu[i])]: od: var:={}: X:=[]: for i from 1 to nops(x) do ku:=GenPol(x,a[i],deg[i]): X:=[op(X),ku[1]]: var:=var union ku[2]: od: var:=var union {seq(e[i],i=0..L)}: gu:=h: for i from 1 to nops(x) do gu:=expand(normal(gu-(diff(r[i],x[i])+q[i])*X[i]/r[i]-r[i]*diff(X[i]/r[i],x[i]))): od: eq:=Extract1(numer(gu),x): eqN:=subs(n=9/17,eq): eqN:=subs({seq(pars[i]=1/(i^2+3),i=1..nops(pars))},eqN): var1N:=solve(eqN,var): opeN:=subs(var1N,ope): if opeN=0 then RETURN(FAIL): else print(`There is hope for a recurrence of order`,L): fi: var1:=solve(eq,var): ope:=subs(var1,ope): X:=subs(var1,X): if ope=0 then RETURN(FAIL): fi: bilti:={}: for i from 1 to nops(var1) do if op(1,op(i,var1))=op(2,op(i,var1)) then bilti:=bilti union {op(1,op(i,var1))}: fi: od: if nops(bilti)>1 then print(`There is some slack`): fi: for i from 1 to nops(bilti) do if coeff(ope,bilti[i],1)<>0 then ope:=coeff(ope,bilti[i],1): X:=[seq(coeff(X[i1],bilti[i],1),i1=1..nops(X))]: meka:=coeff(ope,N,degree(ope,N)): ope:=expand(normal(ope/meka)): ope:=yafe(ope,N): X:=[seq(normal(X[i1]/meka/t^L),i1=1..nops(X))]: RETURN(ope,X): fi: od: end: ########################### End of AlmkvistZeilberger ########## #Khalek(ope1,ope2,n,N): given recurrence #operators ope1 and ope2, both monic in #finds the operators ope3 such that ope1=ope3*ope2 #or returns FAIL Khalek:=proc(ope1,ope2,n,N) local ope3,d,lu,mu: if degree(ope1,N)=degree(ope2,N) then if ope1=ope2 then RETURN(1): else RETURN(FAIL): fi: fi: d:=degree(ope1,N)-degree(ope2,N): ope3:=expand(ope1-N^d*subs(n=n+d,ope2)): lu:=yafe11(ope3,N): mu:=Khalek(lu[2],ope2,n,N): if mu=FAIL then RETURN(FAIL): fi: N^d+lu[1]*mu: end: yafe11:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: ################################## END of Khalek ####################