with(SolveTools): ###################################################################### ## GuessHolo: Save this file as GuessHolo To use it, stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read GuessHolo: # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## zeilberg@math.rutgers.edu. # ####################################################################### print(`First Written: Feb. 15, 2005: tested for Maple 8 and 9 `): print(`Version of Feb. 15, 2005: `): print(): print(`This is GuessHolo, A Maple package`): print(`accompanying the article `): print(`It guesses recurrence equations for discrete variables of several variables`): print(): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg .`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): print(): print(`For a list of all functions type: ezra1();`): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(): ezra1:=proc() if args=NULL then print(` GuessHolo: A Maple package for empirically guessing partial recurrence`): print(`equations satisfied by multi-sequences`): print(): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: ApplyOper `): print(`Asy, BaW, BaWSimp,FindGRL1L2, FindGR, FindGRd, findrec, Findrec, FindDegOrd`): print(`GenOper, GenOperU, GuessOperU, GuessRat, RecLaW, RecLaW1, `): print(`RecLaW1D , RecLaWD `): print(`LaW , LaWCube`): print(` SeqFromRec `): print(): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` GuessHolo: A Maple package for empirically guessing partial recurrence`): print(`equations satisfied by multi-sequences`): print(): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: `): print(`Asy, BaW, FindGR, FindGRd, Findrec, GuessRat,`): print(`LaW , LaWCube, RecLaW, RecLaWD, SeqFromRec `): print(): elif nargs=1 and args[1]=ApplyOper then print(`ApplyOper(DataSet,Oper,n,Pt): Given a data-set `): print(`describing a discrete function, let's call it f,`): print(`in terms of pairs {[Pt,value]}, a partial-recurrence `): print(`operator in Ugly notation`): print(`For example, try: `): print(`ApplyOper(LaWCube({[0,1],[1,0]},10),{[1,[0,0]]},[m,n],[1,1]);`): 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]=BaW then print(`BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps ending at the point a`): print(`For example, try: BaW({[1,0],[0,1]},[1,1]);`): elif nargs=1 and args[1]=BaWSimp then print(`BaWSimp(Steps,K): all the values of the Lattice walks with steps Steps`): print(`in the simplex of size K`): print(`For example, try: BaWSimp({[1,0],[0,1]},2);`): elif nargs=1 and args[1]=FindDegOrd then print(`FindDegOrd(DataSet,Ini,Dir1,MaxC): finds the degree and order of a recurrence`): print(`satisfied in the direction Dir1, looking at Ini, and its neigborhood`): print(`For example, try:`): print(`FindDegOrd({seq{seq([[a,b],(a+b)!/a!/b!),a=1..10),b=1..10)},[1,1],[1,0],4);`): elif nargs=1 and args[1]=FindGR then print(`FindGR(DataSet1,L2,Dir1,MaxC,MaxD,n,N): A general recurrence in direction`): print(`Dir1 for a discrete function given by a cubical-data set, with domain`): print(`[0,L2]^k`): print(`and maximum complexity MaxC. `): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`FindGR(LaWCube({[0,1],[1,0]},20),20,2,12,10,n,N);`): elif nargs=1 and args[1]=FindGRd then print(`FindGRd(DataSet1,L2,Dir1,MaxC,MaxD,n,N): A general recurrence in diagonal`): print(` direction, for a discrete function given by a cubical-data set, with domain`): print(`[0,L2]^k`): print(`and maximum complexity MaxC. `): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`FindGRd(LaWCube({[0,1],[1,0]},40),40,4,4,n,N);`): elif nargs=1 and args[1]=FindGRL1L2 then print(`FindGRL1L2(DataSet1,L1,L2,Dir1,MaxC,MaxD,n,N): A general recurrence in direction`): print(`Dir1 for a discrete function given by a cubical-data set, with domain`): print(`[0,L2]^k`): print(`and maximum complexity MaxC. `): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`FindGRL1L2(LaWCube({[0,1],[1,0]},20),10,20,2,12,10,n,N);`): elif nargs=1 and args[1]=findrec then print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`): print(`the sequence f of degree DEGREE and order ORDER.`): print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`): elif nargs=1 and args[1]=Findrec then print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`): print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `): 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]=GenOper then print(`GenOper(n,deg,Shifts,a): a generic linear-partial `): print(`recurrence operator with polynomial coefficients `): print(`(of degree k) using`): print(`the shifts in Shifts, of k variables with`): print(`coefficents indexed a, followed by the set of coeffs.`): print(`For example, try: GenOper([n],1,{1,N},a);`): elif nargs=1 and args[1]=GenOperU then print(`GenOper(n,deg,Shifts,a): a generic linear-partial `): print(`recurrence operator with polynomial coefficients `): print(`, in Ugly notation, (of degree k) using`): print(`the shifts in Shifts, of k variables with`): print(`coefficents indexed a, followed by the set of coeffs.`): print(`For example, try: GenOperU([n],1,{[0,1]},a);`): elif nargs=1 and args[1]=RecLaW then print(`RecLaW(Steps,Dir1,MaxC,n,N): A general recurrence in direction`): print(`Dir1 for LaW(Steps,a), of max degree (in the non-active variables) MaxD`): print(`and maximum complexity MaxC. `): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`RecLaW({[0,1],[1,0]},1,5,n,N);`): elif nargs=1 and args[1]=RecLaWD then print(`RecLaWD(Steps,MaxC,n,N): A general recurrence in direction`): print(`of the diag. for LaW(Steps,a), of max degree `): print(`(in the non-active variables) MaxD`): print(`and maximum complexity MaxC. `): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`RecLaWD({[0,1],[1,0]},5,n,N);`): elif nargs=1 and args[1]=RecLaW1 then print(`RecLaW1(Steps,Dir1,MaxC,MaxD,L1,L2,n,N): A general recurrence in direction`): print(`Dir1 for LaW(Steps,a), using as data-set a cube [10,10+L]^(k-1)`): print(`and maximum complexity MaxC`): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`RecLaW1({[0,1],[1,0]},2,12,10,20,40,n,N);`): elif nargs=1 and args[1]=RecLaW1D then print(`RecLaW1D(Steps,MaxC,MaxD,L1,L2,n,N): A general recurrence in direction`): print(`of the diag. for LaW(Steps,a), using as data-set a cube [10,10+L]^(k-1)`): print(`and maximum complexity MaxC`): print(`n is the (indexed symbol for the variable)`): print(`For example, try: `): print(`RecLaW1({[0,1],[1,0]},12,10,20,40,n,N);`): elif nargs=1 and args[1]=GuessOperU then print(`Do not use yet!, does not seem to be so good!`): print(`GuessOperU(DataSet1,n,deg,Shifts):Given a data-set DataSet1,`): print(`and a list of variables n, and a non-neg. integer deg, and`): print(`a set of Shifts (vectors of the same lengh as the first-componets`): print(`of DataSet1, guesses an operator of degree deg with shifts Shifts`): print(`For example, try:`): print(`GuessOperU(LaWCueb({[1,0],[0,1]},10),[n,m],0,{[1,0],[0,1],[1,1]});`): elif nargs=1 and args[1]=GuessRat then print(`GuessRat(F,Resh,deg): Given a set of pairs [vect,value]`): print(`and a symbolic list, Resh, guesses a rational function`): print(`with total degrees of sums of degrees of top and bottom up to`): print(`deg that seems fits it, or RETURNS FAIL`): print(`For example, try: GuessRat({seq([[i],i]],i=1..10)},[a],1);`); elif nargs=1 and args[1]=LaW then print(`LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps`): print(`in the set of steps Steps ending at the point a`): print(`For example, try: LaW({[1,0],[0,1]},[1,1]);`): elif nargs=1 and args[1]=LaWCube then print(`LaWCube(Steps,K): all the values of the Lattice walks with steps Steps`): print(`in the cube of size K`): print(`For example, try: LaWCube({[1,0],[0,1]},2);`): elif nargs=1 and args[1]=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): print(`For example, try:`): print(`SeqFromRec(N-n-1,n,N,[1],10);`): fi: end: ############From GuessRat #GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a, #and a non-negative integer deg, outputs a generic polynomial #of total degree deg,in the variables Resh, indexed by the letter a of #followed by the set of coeffs. For example, try: #GenPol([x,y],a,1); GenPol:=proc(Resh,a,deg) local var,POL1,i,x,Resh1,POL: if not type(Resh,list) or nops(Resh)=0 then ERROR(`Bad Input`): fi: x:=Resh[1]: if nops(Resh)=1 then RETURN(convert([seq(a[i]*x^i,i=0..deg)],`+`),{seq(a[i],i=0..deg)}): fi: Resh1:=[op(2..nops(Resh),Resh)]: var:={}: POL:=0: for i from 0 to deg do POL1:=GenPol(Resh1,a[i],deg-i): var:=var union POL1[2]: POL:=expand(POL+POL1[1]*x^i): od: POL,var: end: #GenRat(Resh,a,b,degT,degB): The generic form of a rational function #in the list of variables Resh indexed by a on top and b at bottom #of total degree degT on top and degB at bottom #for example, try: GenRat([x,y],a,b,1,2) GenRat:=proc(Resh,a,b,degT,degB) local TOP,BOT: TOP:=GenPol(Resh,a,degT): BOT:=GenPol(Resh,b,degB): TOP[1]/BOT[1],TOP[2] union BOT[2]: end: #GuessRat1(F,Resh,degT,degB): Given a set of pairs [vect,value] #and a symbolic lis Resh, nops(Resh) arguments (that are integers), guesses a rational function #with total degree degT and degB for numerator and denominator respectively #that seems fits it GuessRat1:=proc(F,Resh,degT,degB) local eq,var,Rat,a,b,i1,j1,k: k:=nops(Resh): Rat:=GenRat(Resh,a,b,degT,degB): var:=Rat[2]: Rat:=Rat[1]: eq:= {seq( numer(normal(F[i1][2]-subs({seq(Resh[j1]=F[i1][1][j1],j1=1..k)},Rat))), i1=1..nops(F))}: eq:=eq minus {0}: if nops(eq)FAIL then RETURN(factor(lu)): fi: od: FAIL: end: #####End GuessRat ###Findrec #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); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+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)+2 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:=Linear(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(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: 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: 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:=Linear(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 RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrecK(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); #If there is slack it also returns FAIL findrecK:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+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)+2 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:=Linear(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 RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: Yafe:=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: #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 0 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: CheckRec:=proc(L,ope,n,N) local i,n0,lu: lu:=normal([seq(add(subs(n=n0,coeff(ope,N,i))*L[n0+i],i=0..degree(ope,N)),n0=1..nops(L)-degree(ope,N))]): evalb(lu=[0$nops(lu)]): end: #findrecK(f,DEGREE,ORDER,n,N): Is there a unique recurrence if degree DEGREE and order ORDER #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecK:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu: if (1+DEGREE)*(1+ORDER)+3+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)+2 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:=Linear(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(false): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)=1 then RETURN(true): else RETURN(false): fi: end: #end Findrec #FindDegOrd1(DataSet,Ini,Dir1,MaxC): finds the degree and order of a recurrence #satisfied in the direction Dir1, looking at Ini, and its neigborhood #For example, try: #FindDegOrd1({seq{seq([[a,b],(a+b)!/a!/b!),a=1..10),b=1..10)},[1,1],[1,0],4); FindDegOrd1:=proc(DataSet1,Ini,Dir1,MaxC) local i,Domain1,T,n,N,ope1,j,i1,f1: Domain1:={seq(DataSet1[i][1],i=1..nops(DataSet1))}: for i from 1 to nops(DataSet1) do T[DataSet1[i][1]]:=DataSet1[i][2]: od: for i from 0 while member([seq(Ini[j]+Dir1[j]*i,j=1..nops(Ini))],Domain1) do od: i:=i-1: f1:=[seq(T[[seq(Ini[j]+Dir1[j]*i1,j=1..nops(Ini))]],i1=0..i)]: ope1:=Findrec(f1,n,N,MaxC): if ope1<>FAIL then RETURN([ degree(numer(normal(ope1)),n),degree(ope1,N)]): else RETURN(FAIL): fi: end: ####Begin Asymptotics #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): 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,x,mu,ope1,ku,i,f: gu:=Aope1(ope,n,N): pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then print(`There is no unique dominant root`): RETURN(FAIL): 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:=Yafe(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:=coeff(taylor(mu,x=0,2),x,1): alpha:=subs(Linear({ku},{alpha}),alpha); f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: lu^n*n^alpha*f: end: ###EndAsymptotics #LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps ending at the point a #For example, try: LaW({[1,0],[0,1]},[1,1]); LaW:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: add(LaW(Steps,Prev[i]),i=1..nops(Prev)): end: #BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps #in the set of steps Steps ending at the point a #For example, try: BaW({[1,0],[0,1]},[1,1]); BaW:=proc(Steps,a) local i,j,k,Prev: option remember: k:=nops(a): if a=[0$k] then RETURN(1): fi: if min(op(a))<0 then RETURN(0): fi: if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then RETURN(0): fi: Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}: add(BaW(Steps,Prev[i]),i=1..nops(Prev)): end: #Cube1(k,K): [0,K]^k Cube1:=proc(k,K) local gu,i,j: if k=0 then RETURN({[]}): fi: gu:=Cube1(k-1,K): {seq(seq([op(gu[i]),j],j=0..K),i=1..nops(gu))}: end: #Cube2(k,K1,K2): [K1,K2]^k Cube2:=proc(k,K1,K2) local gu,i,j: if k=0 then RETURN({[]}): fi: gu:=Cube2(k-1,K1,K2): {seq(seq([op(gu[i]),j],j=K1..K2),i=1..nops(gu))}: end: #Simp1(k,K): all sequenes K>=a1..>=ak>=0 Simp1:=proc(k,K) local gu,i,j,j1,lu: if k=0 then RETURN({[]}): fi: gu:={}: for j from 0 to K do lu:=Simp1(k-1,j): gu:=gu union {seq(seq([j1,op(lu[i])],j1=j..K),i=1..nops(lu))}: od: gu: end: #Simp2(k,K1,K2): all sequenes K2>=a1..>=ak>=K1 Simp2:=proc(k,K1,K2) local gu,i,j,j1,lu: if k=0 then RETURN({[]}): fi: gu:={}: for j from K1 to K2 do lu:=Simp2(k-1,K1,j): gu:=gu union {seq(seq([j1,op(lu[i])],j1=j..K2),i=1..nops(lu))}: od: gu: end: #LaWCube(Steps,K): all the values of the Lattice walks with steps Steps #in the cube of size K #For example, try: LaWCube({[1,0],[0,1]},2); LaWCube:=proc(Steps,K) local lu,a,k,gu,i: k:=nops(Steps[1]): lu:=Cube1(k,K): gu:={}: for i from 1 to nops(lu) do a:=lu[i]: gu:=gu union {[a,LaW(Steps,a)]}: od: gu: end: #BaWSimp(Steps,K): all the values of the Ballot Lattice walks with steps Steps #in the simplex of size K #For example, try: BaWSimp({[1,0],[0,1]},2); BaWSimp:=proc(Steps,K) local lu,a,k,gu,i: k:=nops(Steps[1]): lu:=Simp1(k,K): gu:={}: for i from 1 to nops(lu) do a:=lu[i]: gu:=gu union {[a,BaW(Steps,a)]}: od: gu: end: #FindDegOrd(DataSet,Ini,Dir1,MaxC): finds the degree and order of a recurrence #satisfied in the direction Dir1, looking at Ini, and its neigborhood #For example, try: #FindDegOrd({seq{seq([[a,b],(a+b)!/a!/b!),a=1..10),b=1..10)},[1,1],[1,0],4); FindDegOrd:=proc(DataSet1,Ini,Dir1,MaxC) local i,gu1,gu2,gu3: if not Dir1=[1,0$(nops(Dir1)-1)] then gu1:=FindDegOrd1(DataSet1,Ini,Dir1,MaxC): gu2:=FindDegOrd1(DataSet1,[Ini[1]+1,op(2..nops(Ini),Ini)],Dir1,MaxC): gu3:=FindDegOrd1(DataSet1,[Ini[1]+2,op(2..nops(Ini),Ini)],Dir1,MaxC): if member(FAIL,{gu1,gu2,gu3}) then RETURN(FAIL): fi: for i from 3 while nops({gu1,gu2,gu3})<>1 do gu1:=gu2: gu2:=gu3: gu3:=FindDegOrd1(DataSet1,[Ini[1]+i,op(2..nops(Ini),Ini)],Dir1,MaxC): if gu3=FAIL then RETURN(FAIL): fi: od: RETURN(gu1,[Ini[1]+i-2,op(2..nops(Ini),Ini)]): else gu1:=FindDegOrd1(DataSet1,Ini,Dir1,MaxC): gu2:=FindDegOrd1(DataSet1,[op(1..nops(Ini)-1,Ini),Ini[nops(Ini)]+1],Dir1,MaxC): gu3:=FindDegOrd1(DataSet1,[op(1..nops(Ini)-1,Ini),Ini[nops(Ini)]+2],Dir1,MaxC): if member(FAIL,{gu1,gu2,gu3}) then RETURN(FAIL): fi: for i from 3 while nops({gu1,gu2,gu3})<>1 do gu1:=gu2: gu2:=gu3: gu3:=FindDegOrd1(DataSet1,[op(1..nops(Ini)-1,Ini),Ini[nops(Ini)]+i],Dir1,MaxC): if gu3=FAIL then RETURN(FAIL): fi: od: RETURN(gu1,[op(1..nops(Ini)-1,Ini),Ini[nops(Ini)]+i-2]): fi: end: #GenOper(n,deg,Shifts,a): a generic linear-partial #recurrence operator with polynomial coefficients #(of degree k) in the list of variables n using #the shifts in Shifts, of k variables with #coefficents indexed a, followed by the set of coeffs. GenOper:=proc(n,deg,Shifts,a) local gu,kv,i,lu: kv:={}: gu:=0: for i from 1 to nops(Shifts) do lu:=GenPol(n,a[i],deg): gu:=gu+lu[1]*Shifts[i]: kv:=kv union lu[2]: od: gu,kv: end: #GenOperU(n,deg,Shifts,a): a generic linear-partial #recurrence operator with polynomial coefficients #(of degree k) in the list of variables n using #the shifts in Shifts, in ugly notation #coefficents indexed a, followed by the set of coeffs. #For example, try: #GenOperU([n],1,{[0,1]},a); GenOperU:=proc(n,deg,Shifts,a) local gu,kv,i,lu: kv:={}: gu:={}: for i from 1 to nops(Shifts) do lu:=GenPol(n,a[i],deg): gu:=gu union {[lu[1],Shifts[i]]}: kv:=kv union lu[2]: od: gu,kv: end: AddPts:=proc(Pt1,Pt2) local i: [seq(Pt1[i]+Pt2[i],i=1..nops(Pt1))]:end: MonomToVect:=proc(Mono,N) local i: [seq(degree(Mono,N[i]),i=1..nops(N))]: end: Domain:=proc(DataSet) local i: option remember: {seq(DataSet[i][1],i=1..nops(DataSet))}: end: #ApplyOper(DataSet,Oper,n,Pt): Given a data-set #describing a discrete function, let's call it f, #in terms of pairs {[Pt,value]}, a partial-recurrence #operator in ugly notation, where n is a list of discrete variables #point Pt, evaluates Oper(n,N)f(Pt). #For example, try: #ApplyOper(LaWCube({[0,1],[1,0]},10),{[1,[0,0]]},[m,n],[1,1]); ApplyOper:=proc(DataSet,Oper,n,Pt) local i,i1,Domain1,Neighs1,Shifts1,T: Domain1:=Domain(DataSet): Shifts1:={seq(Oper[i][2],i=1..nops(Oper))}: Neighs1:={seq(AddPts(Shifts1[i],Pt),i=1..nops(Shifts1))}: for i from 1 to nops(DataSet) do T[DataSet[i][1]]:=DataSet[i][2]: od: if Neighs1 minus Domain1<>{} then RETURN(FAIL): fi: add(subs({seq(n[i1]=Pt[i1],i1=1..nops(n))},Oper[i][1])* T[AddPts(Oper[i][2],Pt)],i=1..nops(Oper)): end: #GuessOperU(DataSet1,n,deg,Shifts):Given a data-set DataSet1, #and a list of variables n, and a non-neg. integer deg, and #a set of Shifts (vectors of the same lengh as the first-componets #of DataSet1, guesses an operator of degree deg with shifts Shifts #For example, try: #GuessOperU(LaWCueb({[1,0],[0,1]},10),[n,m],0,{[1,0],[0,1],[1,1]}); GuessOperU:=proc(DataSet1,n,deg,Shifts) local gu,var,Oper,i,var1,a, Domain1,Pt,eq,eq1,atsma: gu:=GenOperU(n,deg,Shifts,a): var:=gu[2]: Oper:=gu[1]: eq:={}: Domain1:=Domain(DataSet1): for i from 1 to nops(Domain1) while nops(eq)<=nops(var)+30 do Pt:=Domain1[i]: eq1:=ApplyOper(DataSet1,Oper,n,Pt): if eq1<>FAIL then eq:= eq union {eq1}: fi: od: if nops(eq)1 then print(`There is too much slack`): print(`Try lower degree`): RETURN(FAIL): fi: Oper:=subs(var1,Oper): if {seq(Oper[i][1],i=1..nops(Oper))}={0} then RETURN(FAIL): else Oper:=subs(atsma[1]=1,Oper): RETURN(Oper): fi: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #GuessRatP(F,Resh,Deg,N): If the data are polynomials in N of degree degN #does GuesRat separately for each coeffs. GuessRatP:=proc(F,Resh,Deg,N,degN) local F1,i,gu,gu1,i1: gu:=0: for i from 0 to degN do F1:={seq([F[i1][1],coeff(F[i1][2],N,i)],i1=1..nops(F))}: if nops({seq(F1[i1][2],i1=1..nops(F1))})=1 then gu1:=F1[1][2]: else gu1:=GuessRat(F1,Resh,Deg): if gu1=FAIL then RETURN(FAIL): fi: fi: gu:=gu+gu1*N^i: od: gu: end: ################Parallel to the axes ###Unrestriced Walks #RecLaW1(Steps,Dir1,MaxC,MaxD,L1,L2,n,N): A general recurrence in direction #Dir1 for LaW(Steps,a), using as data-set a cube [10,10+L]^(k-1) #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #RecLaW1({[0,1],[1,0]},2,12,10,530n,N); RecLaW1:=proc(Steps,Dir1,MaxC,MaxD,L1,L2,n,N) local gu,ope,lu,k,Pt1,Pt2,Pt3,ope1,ope2,ope3,i,gu1,gu2,gu3, Deg,Ord,Pt,i1,ku,Resh: k:=nops(Steps[1]): lu:=Cube2(k-1,L1,L2): gu:={}: Pt1:=lu[1]: gu1:=[seq(LaW(Steps,[op(1..Dir1-1,Pt1),i ,op(Dir1..k-1,Pt1)]),i=1..MaxC^2+10)]: ope1:=Findrec(gu1,n[Dir1],N[Dir1],MaxC): if ope1=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt2:=lu[2]: gu2:=[seq(LaW(Steps,[op(1..Dir1-1,Pt2),i,op(Dir1..k-1,Pt2)]), i=1..MaxC^2+10)]: ope2:=Findrec(gu2,n[Dir1],N[Dir1],MaxC): if ope2=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt3:=lu[3]: gu3:=[seq(LaW(Steps,[op(1..Dir1-1,Pt3),i,op(Dir1..k-1,Pt3)]), i=1..MaxC^2+10)]: ope3:=Findrec(gu3,n[Dir1],N[Dir1],MaxC): if ope3=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: if nops({degree(ope1,N[Dir1]),degree(ope2,N[Dir1]),degree(ope3,N[Dir1])})<>1 then RETURN(FAIL): fi: if nops({degree(numer(normal(ope1)),n[Dir1]), degree(numer(normal(ope2)),n[Dir1]), degree(numer(normal(ope3)),n[Dir1])})<>1 then RETURN(FAIL): fi: Deg:=degree(numer(normal(ope1)),n[Dir1]): Ord:=degree(ope1,N[Dir1]): gu:={[Pt1,ope1],[Pt2,ope2],[Pt3,ope3]}: Resh:=[seq(n[i1],i1=1..Dir1-1),seq(n[i1],i1=Dir1+1..k)]: for i from 4 to nops(lu) do Pt:=lu[i]: ku:=[seq(LaW(Steps,[op(1..Dir1-1,Pt),i1,op(Dir1..k-1,Pt)]), i1=1..MaxC^2+10)]: ope:=findrec(ku,Deg,Ord,n[Dir1],N[Dir1]): if ope<>FAIL then gu:=gu union {[Pt,ope]}: fi: od: GuessRatP(gu,Resh,MaxD,N[Dir1],Ord): end: #RecLaW2(Steps,Dir1,MaxC,MaxD,n,N): A general recurrence in direction #Dir1 for LaW(Steps,a), of max degree (in the non-active variables) MaxD #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #RecLaW2({[0,1],[1,0]},2,12,n,N); RecLaW2:=proc(Steps,Dir1,MaxC,MaxD,n,N) local L1,L2,neq,Resh,x,i1,k,a,b,i: k:=nops(Steps[1]); L1:=15: Resh:=[seq(x[i1],i1=1..k-1)]: neq:=max(seq(nops(GenRat(Resh,a,b,i,MaxD-i)[2]),i=0..MaxD)): for L2 from L1 while nops(Cube2(k-1,L1,L2))<=neq+10 do od: RecLaW1(Steps,Dir1,MaxC,MaxD,L1,L2,n,N): end: RecLaW:=proc(Steps,Dir1,MaxC,n,N) local MaxD1,MaxC1,gu,L: for L from 2 to MaxC do for MaxC1 from 2 to L do MaxD1:=L-MaxC1: gu:=RecLaW2(Steps,Dir1,MaxC1,MaxD1,n,N): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: #End Unrestriced walks #######End Parallel to the axes ################Diagonals ###Unrestriced Walks #RecLaW1D(Steps,MaxC,MaxD,L1,L2,n,N): A general recurrence in #diagonal direction #for LaW(Steps,a), using as data-set a cube [10,10+L]^(k-1) #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #RecLaW1D({[0,1],[1,0]},12,10,10,20,n,N); RecLaW1D:=proc(Steps,MaxC,MaxD,L1,L2,n,N) local gu,ope,lu,k,Pt1,Pt2,Pt3,ope1,ope2,ope3,i,gu1,gu2,gu3, Deg,Ord,Pt,i1,ku,Resh,j1,x: k:=nops(Steps[1]): lu:=Cube2(k-1,L1,L2): gu:={}: Pt1:=lu[1]: gu1:=[seq(LaW(Steps,[i,seq(Pt1[i1]+i,i1=1..nops(Pt1))]), i=1..MaxC^2+10)]: ope1:=Findrec(gu1,n,N,MaxC): if ope1=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt2:=lu[2]: gu2:=[seq(LaW(Steps,[i,seq(Pt2[i1]+i,i1=1..nops(Pt2))]), i=1..MaxC^2+10)]: ope2:=Findrec(gu2,n,N,MaxC): if ope2=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt3:=lu[3]: gu3:=[seq(LaW(Steps,[i,seq(Pt3[i1]+i,i1=1..nops(Pt3))]), i=1..MaxC^2+10)]: ope3:=Findrec(gu3,n,N,MaxC): if ope3=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: if nops({degree(ope1,N),degree(ope2,N),degree(ope3,N)})<>1 then RETURN(FAIL): fi: if nops({degree(numer(normal(ope1)),n), degree(numer(normal(ope2)),n), degree(numer(normal(ope3)),n)})<>1 then RETURN(FAIL): fi: Deg:=degree(numer(normal(ope1)),n): Ord:=degree(ope1,N): gu:={[Pt1,ope1],[Pt2,ope2],[Pt3,ope3]}: Resh:=[seq(x[i1],i1=2..k)]: for i from 4 to nops(lu) do Pt:=lu[i]: ku:=[seq(LaW(Steps,[i1,seq(i1+Pt[j1],j1=1..nops(Pt))]), i1=1..MaxC^2+10)]: ope:=findrec(ku,Deg,Ord,n,N): if ope<>FAIL then gu:=gu union {[Pt,ope]}: fi: od: ope:=GuessRatP(gu,Resh,MaxD,N,Ord): ope:=subs({n=n[1],seq(x[i1]=n[i1]-n[1],i1=2..k)},ope): end: #RecLaW2D(Steps,MaxC,MaxD,n,N): A general recurrence in direction #of the diag. # for LaW(Steps,a), of max degree (in the non-active variables) MaxD #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #RecLaW2D({[0,1],[1,0]},2,12,n,N); RecLaW2D:=proc(Steps,MaxC,MaxD,n,N) local L1,L2,neq,Resh,x,i1,k,a,b,i: k:=nops(Steps[1]); L1:=20: Resh:=[seq(x[i1],i1=1..k-1)]: neq:=max(seq(nops(GenRat(Resh,a,b,i,MaxD-i)[2]),i=0..MaxD)): for L2 from L1 while nops(Cube2(k-1,L1,L2))<=neq+15 do od: RecLaW1D(Steps,MaxC,MaxD,L1,L2,n,N): end: RecLaWD:=proc(Steps,MaxC,n,N) local MaxD1,MaxC1,gu,L: for L from 2 to MaxC do for MaxC1 from 2 to L do MaxD1:=L-MaxC1: gu:=RecLaW2D(Steps,MaxC1,MaxD1,n,N): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: #End Unrestriced walks #######End Diagonal ################Parallel to the axes general DataSet #FindGRL1L2(DataSet1,L1,L2,Dir1,MaxC,MaxD,n,N): A general recurrence in direction #Dir1 for a discrete function given by a cubical-data set, with domain #[0,L2]^k #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #FindGRL1L2(LaWCube({[0,1],[1,0]},40),10,40,2,12,10,n,N); FindGRL1L2:=proc(DataSet,L1,L2,Dir1,MaxC,MaxD,n,N) local gu,ope,lu,k,Pt1,Pt2,Pt3,ope1,ope2,ope3,i,gu1,gu2,gu3, Deg,Ord,Pt,i1,ku,Resh,T: if MaxC^2+10>=L2 then ERROR(`L2 must be bigger than MaxC^2+10, i.e.`, MaxC^2+10): fi: k:=nops(DataSet[1][1]): for i from 1 to nops(DataSet) do T[DataSet[i][1]]:=DataSet[i][2]: od: lu:={seq(DataSet[i][1],i=1..nops(DataSet))}: if lu<>Cube1(k,L2) then ERROR(`Domain of Data set does not match`,[0,L2]^k): fi: lu:=Cube2(k-1,L1,L2): gu:={}: Pt1:=lu[1]: gu1:=[seq(T[[op(1..Dir1-1,Pt1),i,op(Dir1..k-1,Pt1)]], i=1..MaxC^2+10)]: ope1:=Findrec(gu1,n[Dir1],N[Dir1],MaxC): if ope1=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt2:=lu[2]: gu2:=[seq(T[[op(1..Dir1-1,Pt2),i,op(Dir1..k-1,Pt2)]], i=1..MaxC^2+10)]: ope2:=Findrec(gu2,n[Dir1],N[Dir1],MaxC): if ope2=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt3:=lu[3]: gu3:=[seq(T[[op(1..Dir1-1,Pt3),i,op(Dir1..k-1,Pt3)]], i=1..MaxC^2+10)]: ope3:=Findrec(gu3,n[Dir1],N[Dir1],MaxC): if ope3=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: if nops({degree(ope1,N[Dir1]),degree(ope2,N[Dir1]),degree(ope3,N[Dir1])})<>1 then RETURN(FAIL): fi: if nops({degree(numer(normal(ope1)),n[Dir1]), degree(numer(normal(ope2)),n[Dir1]), degree(numer(normal(ope3)),n[Dir1])})<>1 then RETURN(FAIL): fi: Deg:=degree(numer(normal(ope1)),n[Dir1]): Ord:=degree(ope1,N[Dir1]): gu:={[Pt1,ope1],[Pt2,ope2],[Pt3,ope3]}: Resh:=[seq(n[i1],i1=1..Dir1-1),seq(n[i1],i1=Dir1+1..k)]: for i from 4 to nops(lu) do Pt:=lu[i]: ku:=[seq(T[[op(1..Dir1-1,Pt),i1,op(Dir1..k-1,Pt)]], i1=1..MaxC^2+10)]: ope:=findrec(ku,Deg,Ord,n[Dir1],N[Dir1]): if ope<>FAIL then gu:=gu union {[Pt,ope]}: fi: od: GuessRatP(gu,Resh,MaxD,N[Dir1],Ord): end: #FindGR(DataSet1,L2,Dir1,MaxC,MaxD,n,N): A general recurrence in direction #Dir1 for a discrete function given by a cubical-data set, with domain #[0,L2]^k #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #FindGR(LaWCube({[0,1],[1,0]},20),20,2,12,10,n,N); FindGR:=proc(DataSet,L2,Dir1,MaxC,MaxD,n,N) local L1,gu: if MaxC^2+10>=L2 then ERROR(`L2 must be bigger than MaxC^2+10, i.e.`, MaxC^2+10): fi: for L1 from 5 to L2 do gu:=FindGRL1L2(DataSet,L1,L2,Dir1,MaxC,MaxD,n,N): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: ##########end Parallel to the axes general DataSet ################Diagonal general DataSet #FindGRdL1L2(DataSet1,L1,L2,MaxC,MaxD,n,N): A general recurrence in direction #Dir1 for a discrete function given by a cubical-data set, with domain #[0,L2]^k #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #FindGRdL1L2(LaWCube({[0,1],[1,0]},40),10,40,2,12,10,n,N); FindGRdL1L2:=proc(DataSet,L1,L2,MaxC,MaxD,n,N) local gu,ope,lu,k,Pt1,Pt2,Pt3,ope1,ope2,ope3,i,gu1,gu2,gu3, Deg,Ord,Pt,i1,ku,Resh,T,x,i2: if MaxC^2+10>=L2 then ERROR(`L2 must be bigger than MaxC^2+10, i.e.`, MaxC^2+10): fi: k:=nops(DataSet[1][1]): for i from 1 to nops(DataSet) do T[DataSet[i][1]]:=DataSet[i][2]: od: lu:={seq(DataSet[i][1],i=1..nops(DataSet))}: if lu<>Cube1(k,L2) then ERROR(`Domain of Data set does not match`,[0,L2]^k): fi: lu:=Cube2(k-1,L1,L2): gu:={}: Pt1:=lu[1]: gu1:=[seq(T[[i, seq(i+Pt1[i1],i1=1..nops(Pt1))]], i=1..MaxC^2+10)]: ope1:=Findrec(gu1,n,N,MaxC): if ope1=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt2:=lu[2]: gu2:=[seq(T[[i, seq(i+Pt2[i1],i1=1..nops(Pt2))]], i=1..MaxC^2+10)]: ope2:=Findrec(gu2,n,N,MaxC): if ope2=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: Pt3:=lu[3]: gu3:=[seq(T[[i, seq(i+Pt3[i1],i1=1..nops(Pt3))]], i=1..MaxC^2+10)]: ope3:=Findrec(gu3,n,N,MaxC): if ope3=FAIL then #print(`Nothing with complexity <=`,MaxC): RETURN(FAIL): fi: if nops({degree(ope1,N),degree(ope2,N),degree(ope3,N)})<>1 then RETURN(FAIL): fi: if nops({degree(numer(normal(ope1)),n), degree(numer(normal(ope2)),n), degree(numer(normal(ope3)),n)})<>1 then RETURN(FAIL): fi: Deg:=degree(numer(normal(ope1)),n): Ord:=degree(ope1,N): gu:={[Pt1,ope1],[Pt2,ope2],[Pt3,ope3]}: Resh:=[seq(x[i1],i1=2..k)]: for i from 4 to nops(lu) do Pt:=lu[i]: ku:=[seq(T[[i1, seq(i1+Pt[i2],i2=1..nops(Pt))]], i1=1..MaxC^2+10)]: ope:=findrec(ku,Deg,Ord,n,N): if ope<>FAIL then gu:=gu union {[Pt,ope]}: fi: od: ope:=GuessRatP(gu,Resh,MaxD,N,Ord): ope:=subs({n=n[1],seq(x[i1]=n[i1]-n[1],i1=2..k)},ope): end: #FindGRd(DataSet1,L2,MaxC,MaxD,n,N): A general recurrence in DIAGONAL direction #for a discrete function given by a cubical-data set, with domain #[0,L2]^k #and maximum complexity MaxC. #n is the (indexed symbol for the variable) #For example, try: #FindGRd(LaWCube({[0,1],[1,0]},20),20,12,10,n,N); FindGRd:=proc(DataSet,L2,MaxC,MaxD,n,N) local L1,gu: if MaxC^2+10>=L2 then ERROR(`L2 must be bigger than MaxC^2+10, i.e.`, MaxC^2+10): fi: for L1 from 5 to L2 do gu:=FindGRdL1L2(DataSet,L1,L2,MaxC,MaxD,n,N): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: ##########end Diagonal general DataSet