######################################################################
##DRDS.txt: Save this file as  DRDS.txt                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read DRDS.txt<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by George Spahn and                                        #
#Doron Zeilberger, Rutgers University ,                              #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Spring 2023
 
print(`Created:  Spring  2023`):
print(` This is  DRDS.txt `):
print(`It is one of the packages that accompany the article `):
print(`Experimenting with Discrete Rational Dynamical Systems`):
print(`by George Spahn and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):


print(`---------------------------------------`):
 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the difference equations procedures done directly, type ezraD();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

Digits:=50:
with(combinat):
with(linalg):
with(Optimization):

with(numtheory):

ezraS:=proc()

if args=NULL then
 print(` The story procedures are: FindInvStory, GLaBv, GLynBsV, GSFPv, GSFPvSR,  PaperC, PaperD2, PaperD2sr, PaperD3, PaperD3sr, PaperD4sr, ProveFPdV  `):

 print(``):

else
ezra(args):
fi:

end:


 


ezraD1:=proc()
 
if args=NULL then
 print(` The supporting procedures for difference equations are:  NiceRec, NormOrbD, PerRecs, PerRecsS,`):

 print(``):

else
ezra(args):
fi:

end:

ezraD:=proc()
 
if args=NULL then
 print(` The  procedures dealing directly with difference equations, not via Targem, are:  FindInv, FindPerRecs, FindPerRecsA, FPd, GLaB, GLynB, GLynBs, InvLa, InvLy, PerRecsSS,PerSolsD, OrbD, ProveFPd, RRDE`):

 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()
 
if args=NULL then
 print(` The supporting procedures are: CheckInv, Coeffs,  DisOrb,  IsPer1, Iter, Jac,  LadasDB,  MyExactMax1, MyExactMax2`):
  print(` NumMax, NumMax1, Norm2, OurMax, PerOrb,PerOrb1, PerOrbs,  PerOrbsSp, RandPol,RanPt,  RandRF,	Recon, Targem, TryOut, TryOut1 `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: CrPol, CoGSFP, FP, FPsimp, FPp,  GSFP, LSFP, NormOrb, Orb, RRT `):
 print(` `):



elif nops([args])=1 and op(1,[args])=CheckInv then
print(`CheckInv(F,x,n,Inv): Checks that the invariant Inv for the recurrence given by F. Try:`):
print(`CheckInv((p+x[n])/x[n-1],x,InvLy(p,x,n));`):

elif nops([args])=1 and op(1,[args])=CrPol then
print(`CrPol(T,x,a,K): The crucial polynomial in a[1],.., a[nops(T)] whose non-negativity would imply rigorously that`):
print(`the unique locally stable equilibrium point in the positive orthant is indeed a global attractor. Try:`):
print(`CrPol([x[2],(1+x[2])/(1+x[1])],x,a,4);`):

elif nops([args])=1 and op(1,[args])=Coeffs then
print(`Coeffs(P,X): Given a polynomial in the list of variables X, outputs the set of coefficients. Try:`):
print(`Coeffs((x+y)^3,[x,y]);`):

elif nops([args])=1 and op(1,[args])=CoGSFP then
print(`CoGSFP(F,x,K1,K2): inputs a transformation F, in dimension k, say, (where k=nops(F)) phrsed in terms of x[1],x[2],..,x[k], and positive integer parameters K1,K2,A, outputs a conjectured globally stable fixed point of F`):
print(`in floating-point, but taking K1 random points in [0,1000]^k looking at the orbit to K1, and seeing if it seems to go to the same fixed points all the time. Try:`):
print(`CoGSFP([(x[1]+1)/(x[1]+6)],x,100,10);`):

elif nops([args])=1 and op(1,[args])=DisOrb then
print(`DisOrb(F,x,x0,N,pt): inputs  a mapping F (a list of length, k, say) descibing a mapping of k-dimensional space`):
print(`phrased in terms of the variables x[1], ..., x[k] (x is a symbol), and an initial point X0 (of length k) of mumbers`):
print(`a positive integer N, and a fixed point pt, `):
print(`outputs the sequence of consecutive ratios of the `):
print(` distances to the fixed point pt for the  N members of the orbit.`):
print(`It also returns the last index where the ratio is 1. Try: `):
print(`DisOrb([5/2*x[1]*(1-x[1])],x,[1/2],10,[3/5]);`):

elif nops([args])=1 and op(1,[args])=FindInv then
print(`FindInv(F,x,n,k,d): Given a recurrence F phrased in terms of x[n], or order k, tries to find a polynomial P(x[n],x[n-1],..,x[n-k+1]) such that`):
print(`P(x[n],x[n-1],..,x[n-k+1])/(x[n]*x[n-1]*...*x[n-k+1]) is an invariant. It outputs the invariant. Try:`):
print(`FindInv((p+x[n])/x[n-1],x,n,2,3);`):
print(`FindInv((p+x[n]+x[n-1])/x[n-2],x,n,3,4);`):

elif nops([args])=1 and op(1,[args])=FindInvStory then
print(`FindInvStory(K): The invariants of the rational recurrences`):
print(`x[n+1]=(p+x[n]+x[n-1]+...+x[n-k+2])/x[n-k+1]. Try:`):

print(`FindInvStory(4);`):

elif nops([args])=1 and op(1,[args])=FindPerRecs then
print(`FindPerRecs(K,P,x,n): A list of length K whose ith entry is the list of pairs [p,F] where p is a period and F is the corresponding set of recurrence of period p. Try:`):
print(`FindPerRecs(4,20,x,n);`):


elif nops([args])=1 and op(1,[args])=FindPerRecsA then
print(`FindPerRecsA(K,x,n): A list of length K whose ith entry is the list of pairs [p,F] where p is a period and F is the corresponding set of recurrence of period p, p<=3*k. Try:`):
print(`FindPerRecsA(4,x,n);`):

elif nops([args])=1 and op(1,[args])=FP then
print(`FP(F,x): The fixed points of the transformation F, asuming they are all numeric (i.e. it is the generic case where there are finitely many solutions`). 
print(`It uses  Groebner bases. Try:`):
print(`FP([5/2*x[1]*(1-x[1])],x);`):
print(`FP([x[2],(x[1]+x[2])/2],x);`):



elif nops([args])=1 and op(1,[args])=FPd then
print(`FPd(F,x,n,k): inputs a recurrence F or order k, expressesd as x[n+1]=F where F is a rational function of x[n],..., x[n-k+1], outputs the`):
print(`stable fixed point (if it exists) that is positive, let's call it pt, followed by the recurrence obeyed by x[n]-pt, that hopefully`):
print(`goes to zero. Try:`):
print(`FPd((1+x[n])/(1+x[n-1]),x,n,2);`):
print(`F:=RRDE(x,n,3,1,10): FPd(F,x,n,3);`):

elif nops([args])=1 and op(1,[args])=FPp then
print(`FPp(F,x): The fixed points of the transformation F with all postive coordinates. Try:`):
print(`FPp([5/2*x[1]*(1-x[1])],x);`):
print(`FPp([x[2],(x[1]+x[2])/2],x);`):

elif nops([args])=1 and op(1,[args])=FPsimp then
print(`FPsimp(F,x): The fixed points of the transformation F just using solve. Try:`):
print(`FPsimp([x[2],(x[1]+x[2])/2],x);`):





elif nops([args])=1 and op(1,[args])=GLaB then
print(`GLaB(p,a,b,c,X,Y,Z): The polynomial in X,Y,Z, such that three consecutive terms of the orbit  generalized Ladas recurrence x[n+1]=(p+x[n]+x[n-1])/x[n-2]. with initial conditions`):
print(`x[1]=a, x[2]=b,x[3]=c must satisfy.`):
print(`It also returns the smallest and largest members of the first 20000 terms`):
print(`Try:`):
print(`GLaB(2,1,1,1,X,Y,Z);`):

elif nops([args])=1 and op(1,[args])=GLaBv then
print(`GLaBv(): A paper about the boundedness of the Generalized Ladas equation x[n+1]=(p+x[n]+x[n-1])/x[n-2]. Try:`):
print(`GLaBv();`):


elif nops([args])=1 and op(1,[args])=GLynB then
print(`GLynB(p,a,b,X,Y): The theortical  minimum and maximum of the orbit of the generalized Lynnes recurrence x[n+1]=(p+x[n])/x[n-1]. with initial conditions`):
print(`x[1]=a, x[2]=b. It returns a list of length 4 consisiting of`):
print(`(i) The expression in X and Y such that if you replace X and Y by x[n-1] and x[n] respectively, in the orbit started at [a,b] you get 0`):
print(`(ii) the quartic polynomial in X whose middle roots are the theoretical min and max of any orbit stating with [a,b]`):
print(`(iii) The floating-point values the these middle roots`):
print(`(iv) The actual min and max among the first 20000 terms `):
print(`Try:`):
print(`GLynB(2,1,1,X,Y);`):

elif nops([args])=1 and op(1,[args])=GLynBs then
print(`GLynBs(p,a,b,X): inputs SYMBOLIC (or numeric) p,a,b, outputs The quartic polynomial in Z whose middle roots are the lower and upper bounds of the members of the orbit of`):
print(`the generalized Lynnes equation`):
print(`x[n+1]=(p+x[n])/x[n-1] and initial conditions x[1]=a, x[2]=b. `):
print(`GLynBs(p,a,b,X);`):

elif nops([args])=1 and op(1,[args])=GLynBsV then
print(`GLynBsV(): A verbose form of GLynBs(p,a,b,X) (q.v.), Try:`):
print(`GLynBsV();`):


elif nops([args])=1 and op(1,[args])=GSFP then
print(`GSFP(F,x,K): inputs a transformation F in k dimensions, say, given as a list of length k with an expressions in x[1], ..., x[k]`):
print(`and the symbol x, and a positive integer K, outputs the  unique globally stable fixed point, and an integer k <=K`):
print(`such that if L is the distance-squared of the orbit from the fixed points, then provably`):
print(`L[i]-L[i+k]>=0`):
print(`if there is no such k<=K it returns [pt,FAIL]`):
print(`Try:`):
print(`GSFP([x[2],(1+x[2])/(1+x[1])],x,6);`):
print(`GSFP([x[2],(4+x[2])/(1+x[1])],x,6);`):
print(`F:=Targem(NiceRec(x,n,[1],[2],20),x,n,2); GSFP(F,x,4);`):



elif nops([args])=1 and op(1,[args])=GSFPv then
print(`GSFPv(F,x,K): Verbose version of GSFP(F,x,K) (q.v.)`):
print(`Try:`):
print(`GSFPv([x[2],(1+x[2])/(1+x[1])],x,4);`):
print(`GSFPv([x[2],(4+x[2])/(1+x[1])],x,4);`):
print(`T:=Targem(NiceRec(x,n,[1],[2],12),x,n,2); GSFPv(T,x,4);`):


elif nops([args])=1 and op(1,[args])=GSFPvSR then
print(`GSFPvSR(F,x,K): Semi-rigorous version of  GSFPv(F,x,K) (q.v.)`):
print(`Try:`):
print(`GSFPvSR([x[2],(1+x[2])/(1+x[1])],x,5);`):
print(`GSFPvSR([x[2],(4+x[2])/(1+x[1])],x,5);`):
print(`T:=Targem(NiceRec(x,n,[1],[2],12),x,n,2); GSFPvSR(T,x,4);`):
print(`F:=RRT(x,3,1,10): GSFPvSR(F,x,6);`):


elif nops([args])=1 and op(1,[args])=InvLa then
print(`InvLa(p,x,n): The invariant (eq. (5.12) in Kulenovic/Ladas) for the generalized Ladas eq. x[n+1]=(p+x[n]+x[n-1])/x[n-2]. Try:`):
print(`InvLa(p,x,n);`):

elif nops([args])=1 and op(1,[args])=InvLy then
print(`InvLy(p,x,n): The invariant (eq. (5.12) in Kulenovic/Ladas) for the generalized Lynnes eq. x[n+1]=(p+x[n])/x[n-1]. It also verifies that it is indeed correct. Try:`):
print(`InvLy(p,x,n);`):

elif nops([args])=1 and op(1,[args])=IsPer1 then
print(`IsPer1(F,x,n,k,p): Is the recurrence F of order k, in terms of x[n],...x[n-k+1] periodic of period p?. Try:`):
print(`IsPer1(1/x[n],x,n,1,2);`):

elif nops([args])=1 and op(1,[args])=Iter then
print(`Iter(T,x,m): inputs a transformation T in terms of x[1], ..., x[k] where k is nops(T), and  a positive integer m, outputs`):
print(`its T^m. Try:`):
print(`Iter([31/10*x[1]*(1-x[1])],x,2);`):

elif nops([args])=1 and op(1,[args])=Jac then
print(`Jac(F,x): The Jacobian matrix of the transformation F in x[1], ..., x[k], where k=nops(F). Try:`):
print(`Jac([(1+x[1]+x[2])/(5+2*x[1]+x[2]),(7+3*x[1]+x[2])/(5+12*x[1]+x[2])],x);`):

elif nops([args])=1 and op(1,[args])=LadasDB then
print(`LadasDB(x,n): The  data base of transformations gotten by recurrences`):
print(`Expressing x[n+1] in terms of x[n], x[n-1]. Try:`):
print(`LadasDB(x,n);`):


elif nops([args])=1 and op(1,[args])=LSFP then
print(`LSFP(F,x): The Locally Stable fixed points of the transformation F with all postive coordinates. `):
print(`It also returns the set of eigenvalues for each point. Try:`):
print(`LSFP([5/2*x[1]*(1-x[1])],x);`):

elif nops([args])=1 and op(1,[args])=MyExactMax1 then
print(`MyExactMax1(f,x,a,b): Home-made maximum.  Inputs an expression of f,  in the variable x, and numbers  a and b`):
print(`outputs the absolute maximum in the interval a<=x<=b. Try:`):
print(`MyExactMax1(x*(1-x),x,0,1);`):

elif nops([args])=1 and op(1,[args])=MyExactMax2 then
print(`MyExactMax2(f,x,y,a): inputs a  polynomial f in x and y finds the location and value of the maximum of f`):
print(`in [0,a]^2. Try: `):
print(`MyExactMax2(6-(x-1)^2-(y-5)^2,x,y,10); `):



elif nops([args])=1 and op(1,[args])=NiceRec then
print(`NiceRec(x,n,L1,L2,N): a nice non-linear recurrence with equ. point 1, of order nops(L1)+1 (where nops(L1)=nops(L2)) `):
print(`and L1 and L2 are increasing. Try:`):
print(`NiceRec(x,n,[3],[7],20);`):

elif nops([args])=1 and op(1,[args])=NormOrb then
print(`NormOrb(F,x,x0,pt,N): Given a transformation F phrased in terms of x[1],...,x[nops(F)], an initial condition x0 (the same length of F) a point pt (of the same length of F and x0), and a positive`):
print(`integer N, outputs the first N entires of the Euclidean norm of the orbit of F minut pt. Try:`):
print(`NormOrb([x[2],(1+x[2])/(1+x[1])],x,[2,3],[1,1],10);`):

elif nops([args])=1 and op(1,[args])=NormOrbD then
print(`NormOrbD(F,x,n,INI,K): Given a recurrence of order k=nops(INI), x[n+1]=F(x[n],...,x[n-k+1]), initial conditions INI, such that 0 is a fixed point of F`):
print(`finds the first K terms of the quantities L[k*i+1]^2+...+L[k*i+k]^2. Try:`):
print(`NormOrbD((x[n]-x[n-1])/(2+x[n-1]),x,n,[3.,5.],10);`):

elif nops([args])=1 and op(1,[args])=NumMax then
print(`NumMax(f,x,k,b,h): Inputs an expression f in x[1], ..., x[k] (x is a symbol and k is a positive integer) and a small pos. number h,`):
print(`the resoultion, outputs the maximum in [0,b]^k with resolution h. Try:`):
print(`NumMax(x[1]+x[2]+x[1]*x[2],x,2,1,0.1);`):

elif nops([args])=1 and op(1,[args])=NumMax1 then
print(`NumMax1(f,x,a,b,h): the  max of the function f of x from 0 to b with resolution h. Try:`):
print(`NumMax1(x*(1-x),x,1,0.1);`):


elif nops([args])=1 and op(1,[args])=Norm2 then
print(`Norm2(x): The Euclidean norm of x. Try:`):
print(`Norm2([3,4]);`):


elif nops([args])=1 and op(1,[args])=Orb then
print(`Orb(F,x,x0,N): inputs  a mapping F (a list of length, k, say) descibing a mapping of k-dimensional space`):
print(`phrased in terms of the variables x[1], ..., x[k] (x is a symbol), an initial point X0 (of length k) of numbers`):
print(`and a positive integer N, outputs the orbit of length N. Try:`):
print(`Orb([5/2*x[1]*(1-x[1])],x,[1/2],10);`):
print(`Orb([x[2],(x[1]+x[2])/2],x,[1/2,1/2],10);`):



elif nops([args])=1 and op(1,[args])=PaperC then
print(`PaperC(x,n,d,A,K): An article with K Conjectured about Global Asymptotic Stability for n-dimensional rational transformation of [0,infinity]^n into itself`):
print(`by picking random raional transformations with positive coeficients with numerator and denominators of degree d and coefficients from 1 to A`):
print(`tolerating up to k iterations.`):
print(`Try:`):
print(`PaperC(x,2,1,20,10);`):


elif nops([args])=1 and op(1,[args])=PaperD then
print(`DEPERTATED DO NOT USE `):
print(`PaperD(A,k,K): An article with K theorems about the limits of randomly chosen rational recurrences of order k, `):
print(`with posivie coefficients from 1 to A (both numerator and denominator). Finding the`):
print(`numbers that they converge to, regardless of initial conditions and proving`):
print(`rigorously (and sometimes semi-rigorously) that they indeed  converge to the same number. Try:`):
print(`PaperD(20,2,10);`):


elif nops([args])=1 and op(1,[args])=PaperD2 then
print(`PaperD2(A,K1,K2): A paper with (at most K2) rigorously proved theorem about the limits of second-order rational recurrences. K1 is the complexity parameter. Try:`):
print(`PaperD2(30,4,10);`):




elif nops([args])=1 and op(1,[args])=PaperD2sr then
print(`PaperD2sr(A,K1,K2): A paper with (at most K2) Semi-Rigorously proved theorem about the limits of second-order rational recurrences. K1 is the complexity parameter. Try:`):
print(`PaperD2sr(30,10,4);`):

elif nops([args])=1 and op(1,[args])=PaperD3 then
print(`PaperD3(A,K1,K2): A paper with (at most K2) rigorously proved theorem about the limits of third-order rational recurrences. K1 is the complexity parameter. Try:`):
print(`WARNING: TAKES FOR EVER ON MY COMPUTER, DO NOT USE!, TRY PaperD3sr(A,K1,K2) instread `):
print(`Try the semi-rigorous version, PaperD3(A,K1,K2)`):
print(`PaperD3(30,10,2);`):


elif nops([args])=1 and op(1,[args])=PaperD3sr then
print(`PaperD3sr(A,K1,K2): A paper with (at most K2) semi-rigorously proved theorem about the limits of third-order rational recurrences. K1 is the complexity parameter. Try:`):
print(`PaperD3sr(30,5,2);`):


elif nops([args])=1 and op(1,[args])=PaperD4sr then
print(`PaperD4sr(A,K1,K2): A paper with (at most K2) semi-rigorously proved theorem about the limits of fourth-order rational recurrences. K1 is the complexity parameter. Try:`):
print(`PaperD4sr(30,5,2);`):

elif nops([args])=1 and op(1,[args])=PerOrb then
print(`PerOrb(L,K,eps): Inputs an orbit L, a positive integer k, and a small positive number eps, finds`):
print(`a period of order <=K, if it exists, up to eps closeness, or else returns FAIL. Try:`):
print(`L:=Orb([3.1*x[1]*(1-x[1])],x,[0.6],200): PerOrb(L,5,10^(-10));`):

elif nops([args])=1 and op(1,[args])=PerOrb1 then
print(`PerOrb1(L,k,eps): Inputs an orbit L, a positive integer k, and a small positive number eps, finds`):
print(`a periodic of order k, if it exists, up to eps closeness, or else returns FAIL. Try:`):
print(`L:=Orb([3.1*x[1]*(1-x[1])],x,[0.6],200): PerOrb1(L,2,10^(-10));`):

elif nops([args])=1 and op(1,[args])=PerOrbs then
print(`PerOrbs(F,x,A,N,eps,K,R): inputs a mapping F from k-dim space to itself (where k=nops(F)) phrased in terms of`):
print(`x[1], ..., x[k], a pos. number A (for picking random points), and a large positive integer N, outputs all the limiting periodic`):
print(`orbits of period<=K using orbits of length N, for R random points. If any of them`):
print(`returns FAIL, then it returns FAIL. Try:`):
print(`PerOrbs([5/2*x[1]*(1-x[1])],x,1,400,10^(-10),4,20);`):
print(`PerOrbs([x[2],(x[1]+x[2])/2],x,1,400,10^(-10),4,10);`):
print(`PerOrbs([(11*x[1]+x[2])/(5+x[1]+8*x[2]),(x[1]+x[2])/(2+x[1]+x[2])],x,1,400,10^(-10),4,30);`):
print(` PerOrbs(Targem(LadasDB(x,n)[1],x,n,3),x,10,400,10^(-8),10,30);`):



elif nops([args])=1 and op(1,[args])=PerRecs then
print(`PerRecs(k,p,x,n,a,b): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k)+a[0])/add(b[i]*x[n+1-i],i=1..k)+b[0]).`):
print(`It returns a list of length 2, the first consisting of the non-linear recurrences, the second of the linear ones.Try: `):
print(`PerRecs(1,2,x,n,a,b); PerRecs(2,4,x,n,a,b);PerRecs(2,5,x,n,a,b);PerRecs(2,5,x,n,a,b); PerRecs(3,6,x,n,a,b); PerRecs(2,6,x,n,a,b);PerRecs(3,8,x,n,a,b);`):

elif nops([args])=1 and op(1,[args])=PerRecsS then
print(`PerRecsS(k,p,x,n,a): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k-1)+a[0])/x[n-k+1]). Try:`):
print(`PerRecsS(2,5,x,n,a); `):


elif nops([args])=1 and op(1,[args])=PerRecsSS then
print(`PerRecsSS(k,p,x,n): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k-1)+a[0])/x[n-k+1]). `):
print(`where a[0],a[1],...,a[k-1] are in {0,1}.`):
print(`Try:`):
print(`PerRecsSS(2,5,x,n); `):

elif nops([args])=1 and op(1,[args])=PerSolsD then
print(`PerSolsD(F,x,n,k,p): Inputs a k-th order difference equation in the form x[n+1]=F(x[n],...,x[n-k+1]), for some expression F, and`):
print(`finds all solutions of period p. When p=1, it is the fixed point. Try:`):
print(`F:=RRDE(x,n,3,1,10);PerSolsD(F,x,n,3,1);`):

elif nops([args])=1 and op(1,[args])=PerOrbsSp then
print(`PerOrbsSp(F,x):  PerOrbsSp(F,x,A,N,eps,K,R): PerOrbs(F,x,20,400,10^(-10),10,100). Try: `):
print(`PerOrbsSp([(11*x[1]+1)/(7*x[1]+4)],x);  `):

elif nops([args])=1 and op(1,[args])=OrbD then
print(`OrbD(F,x,n,INI,N): Given a difference equation F or order k phrased as x[n+1]=F(x[n],...,x[n-k+1]), where R is an expression, initial condtions INI (a list of length k), and a positive integer`):
print(`N, finds the list of the first N+k members. Try:`):
print(`OrbD(x[n]+x[n-1],x,n,[1,1],10);`):

elif nops([args])=1 and op(1,[args])=ProveFPd then
print(`ProveFPd(F,x,n,k): Finds the stable fixed point of the recurrence of order k given by F, and proves it rigorously. Try: `):
print(`ProveFPd((1+x[n])/(1+x[n-1]),x,n,2);`):


elif nops([args])=1 and op(1,[args])=ProveFPdV then
print(`ProveFPdV(F,x,n,k): verbose form of ProvePDd(F,x,n,k) (q.v.). Try: `):
print(`ProveFPdV((1+x[n])/(1+x[n-1]),x,n,2);`):

elif nops([args])=1 and op(1,[args])=OurMax then
print(`OurMaxize(f,x): Given a polynomial f in the list of variables x, finds its maximum over the positive orthant. Try`):
print(`OurMax(3-x^2-y^2,[x,y]);`):

elif nops([args])=1 and op(1,[args])=RandPol then
print(`RandPol(x,n,d,A): A random polynomial in x[1], ..., x[n] of total degree <=d with coefficients from 1 to A. Try:`):
print(`RandPol(x,2,1,10);`):

elif nops([args])=1 and op(1,[args])=RanPt then
print(`RanPt(k,A): A random point in [0,A]^k: Try:`):
print(`RanPt(3,100):`):


elif nops([args])=1 and op(1,[args])=RandRF then
print(`RandRF(x,n,d,A): A random rational function in x[1], ..., x[n] whose top and bottom are of total degree <=d with coefficients from 1 to A. Try:`):
print(`RandRF(x,2,1,10);`):

elif nops([args])=1 and op(1,[args])=Recon then
print(`Recon(F,x,A,N): Given a transformation from [0,infinity]^k (where k=nops(F)) phrased in terms of x[1], ..., x[k], and a letter`):
print(`x, finds out whether there is ONE locally stable fixed points, and by exploring`):
print(`N random points in [0,A]^k finds the the largest iterate before becoming a shrinking orbit. `):
print(`It returns the unique locally stable fixed point (if it exists) followed by that iterate`):
print(`Try:`):
print(`Recon([5/2*x[1]*(1-x[1])],x,10,100);`):

elif nops([args])=1 and op(1,[args])=RRDE then
print(`RRDE(x,n,k,d,A): a random rational difference equation of order k, and degree d, with positive coefficients <=A. Try:`):
print(`RRDE(x,n,3,1,10);`):

elif nops([args])=1 and op(1,[args])=RRT then
print(`RRT(x,n,d,A): A random rational  transformation from R^n to R^n, phrased in terms of x[1], ..., x[n]`):
print(`where each component has top and bottom of total degree <=d, with coefficients from 1 to A. Try:`):
print(`RRT(x,2,1,10);`):

elif nops([args])=1 and op(1,[args])=Targem then
print(`Targem(F,x,n,k): Given a k-th order recurrene in the format x[n+1]=F(x[n], ..., x[n-k+1])`):
print(`finds the corresponding transformation from R^k to R^k such that [x[n-k+1], ..., x[n]] goes to`):
print(`[x[n-k+2],..., x[n+1]]=[x[n-k+2],..., x[n+1]]. Try:`):
print(`Targem(x[n-1]/(x[n-1]+x[n-2]),x,n,3);`):

elif nops([args])=1 and op(1,[args])=TryOut then
print(`TryOut(T,x,K): Given a transformation T phrased in terms of x[1],x[2],..,x[nops(T)], and aa positive integer K, conjectures the smallest integer k<=K`):
print(`such that`):
print(`the NormOrb[i]-NormOrb[i+k] is  positive. Or RETURNS FAIL. Try:`):
print(`TryOut([x[2],(1+x[2])/(1+x[1])],x,10);`):
print(`T:=RRT(x,2,1,20); TryOut(T,x,10);`):
print(`T:=RRT(x,3,1,20); TryOut(T,x,10);`):
print(`T:=Targem(NiceRec(x,n,[1],[2],30),x,n,2): TryOut(T,x,5);`):
print(`T:=Targem(NiceRec(x,n,[1,2],[2,3],40),x,n,3): TryOut(T,x,6);`):

elif nops([args])=1 and op(1,[args])=TryOut1 then
print(`TryOut1(T,x,k): Given a transformation T phrased in terms of x[1],x[2],..,x[nops(T)], and aa positive integer k, investigates whether`):
print(`there is a unique locally stable fixed point, and then whether in the NormOrb for 100 randoml chosen initial conditions`):
print(`the NormOrb[i]-NormOrb[i+k] is  positive. It also returns the point. Try:`):
print(`TryOut1([x[2],(1+x[2])/(1+x[1])],x,2);`):
print(`TryOut1([x[2],(1+x[2])/(1+x[1])],x,3);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


#Orb(F,x,x0,N): inputs  a mapping F (a list of length, k, say) descibing a mapping of k-dimensional space
#phrased in terms of the variables x[1], ..., x[k] (x is a symbol), an initial point X0 (of length k) of mumbers
#and a positive integer N, outputs the orbit of length N. Try:
#Orb([5/2*x[1]*(1-x[1])],x,[1/2],10);
#Orb([x[2],(x[1]+x[2])/2],x,[1/2,1/2],10);
Orb:=proc(F,x,x0,N) local i,gu,k,i1,lu:

if not (type(F,list) and type(x,symbol) and type(x0,list) and nops(F)=nops(x0) and type(N,integer) and N>0) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

k:=nops(F):

gu:=[x0]:

for i from 2 to N do
 lu:=gu[i-1]:
 lu:=expand(subs({seq(x[i1]=lu[i1],i1=1..k)},F)):
 gu:=[op(gu),lu]:
od:

gu:

end:


#Norm2(x): The Euclidean distance between x and y. Try:
#Norm2(x);
Norm2:=proc(x) local i:
sqrt(add(x[i]^2,i=1..nops(x))):
end:

#PerOrb1(L,k,eps): Inputs an orbit L, a positive integer k, and a small positive number eps, finds
#the orbit is periodic with period k up to eps closeness, if it exists, or returns FAIL. It also returns the orbit. Try:
#L:=Orb([3.1*x[1]*(1-x[1])],x,[0.6],100); PerOrb1(L,2,10^(-10));
PerOrb1:=proc(L,k,eps) local N,L1,M,katan,i1,vu:
N:=nops(L):
L1:=[op(trunc(N/2)..N,L)]:
vu:=max(seq(evalf(Norm2(L1[i1]-L1[i1+k])),i1=1..nops(L1)-k)):


if  vu<eps then
   M:=evalf([op(nops(L1)-k+1..nops(L1),L1)]):
   katan:=min(seq(evalf(Norm2(M[i1])),i1=1..nops(M))):
   for i1 from 1 while Norm2(M[i1])>katan do od:
   M:=[op(i1..nops(M),M),op(1..i1-1,M)]:
   RETURN(M):
fi:
FAIL:

end:

#PerOrb(L,K,eps): Inputs an orbit L, a positive integer K, and a small positive number eps, finds a period of length<=K,
#up to eps closeness. or returns FAIL. Try:
#L:=Orb([3.1*x[1]*(1-x[1])],x,[0.6],100); PerOrb(L,8,10^(-10));

PerOrb:=proc(L,K,eps) local k,gu:

for k from 1 to K do
 gu:=PerOrb1(L,k,eps):

  if gu<>FAIL then
   RETURN(gu):
  fi:
od:

FAIL:
end:

#RanPt(k,A): A random point in [0,A]^k: Try:
#RanPt(3,100):
RanPt:=proc(k,A) local i,ra:
ra:=rand(10*A..1000*A):
evalf([seq(ra()/1000,i=1..k)]):
end:



#PerOrbs(F,x,A,N,eps,K,R): inputs a mapping F of R^k to R^k (where k=nops(F)) phrased in terms of
#x[1], ..., x[k], a pos. number A, and a large positive integer N, outputs all the limiting periodic
#orbits of period<=K using orbits of length N, for R random points. If any of them
#returns FAIL, then it returns FAIL. Try:
#PerOrbs([5/2*x[1]*(1-x[1])],x,1,400,10^(-10),4,20);
#PerOrbs([x[2],(x[1]+x[2])/2],x,1,400,10^(-10),4,10);
PerOrbs:=proc(F,x,A,N,eps,K,R) local d,gu,x0,mu,r:

d:=nops(F):

gu:={}:

for r from 1 to R do
 x0:=RanPt(d,A):

  mu:=Orb(F,x,x0,N):
  mu:=PerOrb(mu,K,eps):
  if mu=FAIL then
   RETURN(FAIL):
  else
   gu:=gu union {evalf(mu,10)}:
  fi:
od:

gu:

end:

#PerOrbsS(F,x): same as PerOrbs(F,x,1,400,10^(-10),10,100):
PerOrbsSp:=proc(F,x):
PerOrbs(F,x,20,400,10^(-10),10,100):
end:
  
#LadasDB(x,n): The  data base of transformations gotten by recurrences
#Expressing x[n+1] in terms of x[n], x[n-1]. Try:
#LadasDB(x,n);
LadasDB:=proc(x,n):
[
x[n-1]/(x[n-1]+x[n-2]),
(x[n]+x[n-2])/x[n-1],
(1+x[n-2])/x[n],
(1+x[n])/(x[n-1]+x[n-2]),
(1/p^2+x[n])/(p*x[n-1]+x[n-2]),
(1+x[n-1])/(1+x[n-2]),
(x[n]+x[n-1])/(x[n]+x[n-2]),
(x[n-1]+x[n-2])/(x[n]+x[n-2]),
(x[n-1]+x[n-2])/x[n-2],
(1+x[n-1])/(1+x[n]+x[n-2]),
(1+x[n-1]+x[n-2])/x[n-2],
(1+x[n])/x[n-1],
(p+x[n])/x[n-1],
(1+x[n]+x[n-1])/x[n-2],
(p+x[n]+x[n-1])/x[n-2]
]:
end:



#Targem(F,x,n,k): Given a k-th order recurrene in the format x[n+1]=F(x[n], ..., x[n-k+1])
#finds the corresponding transformation from R^k to R^k such that [x[n-k+1], ..., x[n]] goes to
#[x[n-k+2],..., x[n+1]]=[x[n-k+2],..., x[n+1]]. Try:
#Targem(x[n-1]/(x[n-1]+x[n-2]),x,n,3);
Targem:=proc(F,x,n,k) local i:
[seq(x[i],i=2..k),
subs({seq(x[n-i]=x[k-i],i=0..k-1)},F)]:
end:




#FP(F,x): The fixed points of the transformation F. Try:
#FP([5/2*x[1]*(1-x[1])],x);
#FP([x[2],(x[1]+x[2])/2],x);
FP:=proc(F,x) local lu,eq,gu,n,i,i1,X,j,ku,mu:

n:=nops(F):

if n=1 then
 lu:=[solve({numer(x[1]-F[1])},{x[1]})]:
 RETURN({seq([subs(lu[i],x[1])],i=1..nops(lu))}):
fi:

X:=[seq(x[i1],i1=1..n)]:

eq:={seq(numer(x[i]-F[i]),i=1..n)}:


gu:=Groebner[Basis](eq,plex(seq(x[n+1-i],i=1..n))):




#if nops(gu)=1 then
# RETURN([[seq(x[i],i=1..n-1),solve(gu[1],x[n])]]):
#fi:

if n=2 then
if not degree(gu[2],x[2])=1  then
  RETURN(FAIL):
fi:


else
for i from 2 to n do
 if not (degree(gu[i],x[i])=1 and {seq(degree(gu[i],x[j]),j=2..i-1), seq(degree(gu[i],x[j]),j=i+1..n)}={0}) then
  RETURN(FAIL):
fi:
od:

fi:

mu:=[x[1],seq(solve(gu[i],x[i]),i=2..n)]:

ku:=[solve(gu[1],x[1])]:


[seq(subs(x[1]=ku[i1],mu),i1=1..nops(ku))]:

end:


#Jac(F,x): The Jacobian matrix of the transformation F in x[1], ..., x[k], where k=nops(F). Try:
#Jac([(1+x[1]+x[2])/(5+2*x[1]+x[2]),(7+3*[1]+x[2])/(5+12*x[1]+x[2])],x);
Jac:=proc(F,x) local k,i,j:
k:=nops(F):

normal([seq([seq(diff(F[i],x[j]),j=1..k)],i=1..k)]):

end:




#FPp(F,x): The fixed points of the transformation F with all postive coordinates. Try:
#FPp([5/2*x[1]*(1-x[1])],x);
#FPp([x[2],(x[1]+x[2])/2],x);
FPp:=proc(F,x) local lu,i,gu,muam,muam1,j:

lu:=FP(F,x):

gu:={}:

for i from 1 to nops(lu) do

 muam:=lu[i]:
  muam1:=evalf(muam):


 if {seq(type(muam1[j],numeric),j=1..nops(muam1))}={true} then

  if max(seq(abs(coeff(muam1[j],I,1)),j=1..nops(muam1)))<1/10^10 and   min(seq(coeff(muam1[j],I,0),j=1..nops(muam1)))>0  then
   gu:=gu union {muam}: 
  fi:

fi:

od:
gu:

end:





#LSFP(F,x): The Locally Stable fixed points of the transformation F with all postive coordinates. 
#It also returns the set of eigenvalues for each point. Try:
#LSFP([5/2*x[1]*(1-x[1])],x);
LSFP:=proc(F,x) local n,lu,gu,M,vu,M1,i1,mu,i,z,f:
option remember:
n:=nops(F):

lu:=FPp(F,x) :

M:=Jac(F,x):

gu:={}:

for i from 1 to nops(lu) do
mu:=lu[i]:
M1:=subs({seq(x[i1]=mu[i1],i1=1..n)},M):

f:=charpoly(M1,z):

vu:=[solve(f,z)]:

if max(seq(abs(evalf(vu[i1])),i1=1..nops(vu)))<1 then
 gu:=gu union {[mu,evalf(vu)]}:
fi:
od:
gu:

end:

#Iter(T,x,m): inputs a transformation T in terms of x[1], ..., x[k] where k is nops(T), and  a positive integer m, outputs
#its T^m. Try:
#Iter([31/10*x[1]*(1-x[1]),x,2);

Iter:=proc(T,x,m) local gu,i,k:
option remember:
k:=nops(T):
if m=0 then
 RETURN([seq(x[i],i=1..k)]):
elif m=1 then
 RETURN(T):
else
 gu:=Iter(T,x,m-1):
 RETURN(normal(subs({seq(x[i]=T[i],i=1..k)},gu))):
fi:
end:
 





#DisOrb(F,x,x0,N,pt): inputs  a mapping F (a list of length, k, say) descibing a mapping of k-dimensional space
#phrased in terms of the variables x[1], ..., x[k] (x is a symbol), and an initial point X0 (of length k) of mumbers
#a positive integer N, and a fixed point pt, 
#outputs the sequence of the ratios of the consecutive distances to the fixed point pt for the  N members of the orbit. Try:
#and a positive integer N, outputs the orbit of length N. 
#It also returns the last index where the ratio is 1
#Try:
#DisOrb([5/2*x[1]*(1-x[1])],x,[1/2],10,[3,5]);
DisOrb:=proc(F,x,x0,N,pt) local k,gu,i,mu:


k:=nops(F):

if nops(pt)<>k then
print(F, `and `, pt, `should be of the same length `):
 RETURN(FAIL):
fi:

if simplify(subs({seq(x[i]=pt[i],i=1..k)},F)-pt)<>[0$k] then
 RETURN(FAIL):
fi:


gu:=Orb(F,x,x0,N,pt):

gu:=evalf([seq(Norm2(gu[i]-pt),i=1..nops(gu))]):

if member(0,gu) then
 RETURN(gu):
fi:

mu:=[]:
for i from 1 to nops(gu)-1 do
 if abs(gu[i])>10^(-10) then
  mu:=[op(mu),gu[i+1]/gu[i]]:
 fi:
od:

k:=0:

for i from 1 to nops(mu) do
 if mu[i]>1 then
   k:=i:
 fi:
od:

[mu,k]:
end:



#NumMax1(f,x,b,h): the  max of the function f of x from 0 to b with resolution h. Try:
#NumMax1(x*(1-x),x,1,0.1);
NumMax1:=proc(f,x,b,h) local i,N:
N:=trunc(b/h):
max(seq(subs(x=i*h,f),i=0..N)):

end:


#NumMax(f,x,k,b,h): Inputs an expression f in x[1], ..., x[k] (x is a symbol and k is a positive integer) and a small pos. number h,
#the resoultion, outputs the maximum in [0,b]^k with resolution h. Try:
#NumMax(x[1]+x[2]+x[1]*x[2],x,2,1,0.1);
NumMax:=proc(f,x,k,b,h) local i,N:

if k=1 then
 RETURN(NumMax1(f,x[1],b,h)):
else
N:=trunc(b/h):
max(seq(NumMax(subs(x[k]=i*h,f),x,k-1,b,h),i=0..N)):
fi:

end:

#MyExactMax1Float(f,x,a,b): Home-made maximum.  Inputs an expression of f,  in the variable x, and numbers  a and b
#outputs the absolute maximum in the interval a<=x<=b. Try:
#MyExactMax1Float(x*(1-x),x,0,1);
MyExactMax1Float:=proc(f,x,a,b) local f1,gu,aluf,si,i,mu:

f1:=diff(f,x):

gu:=[solve(f1,x)]:

mu:=[]:

for i from 1 to nops(gu) do
 if evalf(gu[i])>=a and evalf(gu[i])<=b then
   mu:=[op(mu),gu[i]]:
 fi:
od:

mu:=[a,op(mu),b]:

aluf:=mu[1]:

si:=evalf(subs(x=mu[1],f)):

for i from 2 to nops(mu) do
 
  if evalf(subs(x=mu[i],f))>si then
   aluf:=mu[i]:
   si:=evalf(subs(x=mu[i],f)):
  fi:
od:

[aluf,si]:

end:


#MyExactMax1(f,x,a,b): Home-made maximum.  Inputs an expression of f,  in the variable x, and numbers  a and b
#outputs the absolute maximum in the interval a<=x<=b. Try:
#MyExactMax1(x*(1-x),x,0,1);
MyExactMax1:=proc(f,x,a,b) local f1,gu,aluf,si,i,mu:

f1:=diff(f,x):

gu:=[solve(f1,x)]:

mu:=[]:

for i from 1 to nops(gu) do
 if evalf(gu[i])>=a and evalf(gu[i])<=b then
   mu:=[op(mu),gu[i]]:
 fi:
od:

mu:=[a,op(mu),b]:

aluf:=mu[1]:

si:=subs(x=mu[1],f):

for i from 2 to nops(mu) do
 
  if evalf(subs(x=mu[i],f))>evalf(si) then
   aluf:=mu[i]:
   si:=subs(x=mu[i],f):
  fi:
od:

[aluf,si]:

end:

#MyExactMax2Float(f,x,y): inputs a  polynomial f in x and y finds the location and value of the maximum of f
#in [0,a]^2. Try:
#MyExactMax2Float(6-(x-1)^2-(y-5)^2,x,y,10);
MyExactMax2Float:=proc(f,x,y,a) local mu,eq,gu,aluf,i,si,ku:
mu:=[[MyExactMax1Float(subs(y=0,f),x,0,a)[1],0],[0,MyExactMax1Float(subs(x=0,f),y,0,a)[1]]]:

eq:={diff(f,x),diff(f,y)}:

gu:={solve(eq,{x,y})}:

gu:={seq(subs(gu[i],[x,y]),i=1..nops(gu))}:


ku:={}:

for i from 1 to nops(gu) do
 if evalf(gu[i][1])>=0 and evalf(gu[i][2])>=0 then
  ku:=ku union {gu[i]}:
 fi:
od:

mu:=[op(mu),op(ku)]:


aluf:=mu[1]:

si:=evalf(subs({x=aluf[1],y=aluf[2]},f)):

for i from 2 to nops(mu) do
 if evalf(subs({x=mu[i][1],y=mu[i][2]},f))>si then
   si:=evalf(subs({x=mu[i][1],y=mu[i][2]},f)):
   aluf:=mu[i]:
 fi:
od:

aluf,si:
   
end:




#MyExactMax2(f,x,y): inputs a  polynomial f in x and y finds the location and value of the maximum of f
#in [0,a]^2. Try:
#MyExactMax2(6-(x-1)^2-(y-5)^2,x,y,10);
MyExactMax2:=proc(f,x,y,a) local mu,eq,gu,aluf,i,si,ku:
mu:=[[MyExactMax1(subs(y=0,f),x,0,a)[1],0],[0,MyExactMax1(subs(x=0,f),y,0,a)[1]]]:

eq:={diff(f,x),diff(f,y)}:

gu:={solve(eq,{x,y})}:

gu:={seq(subs(gu[i],[x,y]),i=1..nops(gu))}:


ku:={}:

for i from 1 to nops(gu) do
 if evalf(gu[i][1])>=0 and evalf(gu[i][2])>=0 then
  ku:=ku union {gu[i]}:
 fi:
od:

mu:=[op(mu),op(ku)]:


aluf:=mu[1]:

si:=simplify(subs({x=aluf[1],y=aluf[2]},f)):

for i from 2 to nops(mu) do
 if evalf(subs({x=mu[i][1],y=mu[i][2]},f))>evalf(si) then
   si:=simplify(subs({x=mu[i][1],y=mu[i][2]},f)):
   aluf:=mu[i]:
 fi:
od:

aluf,si:
   
end:


#RandPol(x,n,d,A): A random polynomial in x[1], ..., x[n] of total degree <=d with coefficients from 1 to A. Try:
#RandPol(x,2,1,10);
RandPol:=proc(x,n,d,A) local ra,i,j:
ra:=rand(1..A):

if n=1 then
 RETURN(add(ra()*x[1]^i,i=0..d)):
fi:

expand(add(x[n]^j*RandPol(x,n-1,d-j,A),j=0..d)):

end:


#RandRF(x,n,d,A): A random rational function in x[1], ..., x[n] of total degree <=d with coefficients from 1 to A. Try:
#RandRF(x,2,1,10);
RandRF:=proc(x,n,d,A):
RandPol(x,n,d,A)/RandPol(x,n,d,A):
end:


#RRT(x,n,d,A): A random rational transformation from R^n to R^n, in terms of x[1] , ..., x[n] with integer coefficients from 1 to A.  
#with total degree 1 in tops and bottoms. Try:
#RRT(2,x,1,10);
RRT:=proc(x,n,d,A) local i:
[seq(RandRF(x,n,d,A),i=1..n)]:
end:


#Recon(F,x,A,N): Given a transformation from [0,infinity]^k (where k=nops(F)) phrased in terms of x[1], ..., x[k], and a letter
#x, finds out whether there is ONE locally stable fixed points, and by exploring
#N random points in [0,A]^k finds the the largest iterate before becoming a shrinking orbit. 
#It returns the unique locally stable fixed point (if it exists) followed by the value of the iterate
#Try:
#Recon([5/2*x[1]*(1-x[1])],x,10,100);
Recon:=proc(F,x,A,N) local gu,pt,x0,i,hal,gadol,k,j:
gu:=LSFP(F,x):

if gu={} then
 print(`There are no locally stable fixed points`):
 RETURN(FAIL):
fi:

if nops(gu)>1 then
 print(`There are more than one locally stable fixed points. Here there are`):
 print({seq(gu[i][1],i=1..nops(gu))}):
 RETURN(FAIL):
fi:

pt:=gu[1][1]:

k:=nops(F):

gadol:=0:

for j from 1 to k do

for i from 1 to N do

x0:=evalf(RanPt(k-1,A)):

x0:=[op(1..j-1,x0),0.02,op(j..k-1,x0)]:


hal:=DisOrb(F,x,x0,30,pt)[2]:
 if hal>gadol then
   gadol:=hal:
 fi:

od:

od:


for i from 1 to N do

x0:=evalf(RanPt(nops(F),A)):


hal:=DisOrb(F,x,x0,30,pt)[2]:
 if hal>gadol then
   gadol:=hal:
 fi:

od:

[pt,gadol]:
end:





	




#FPsimp(F,x): The fixed points of the transformation F just using solve. Try:
#FPsimp([x[2],(x[1]+x[2])/2],x);
FPsimp:=proc(F,x) local eq,n,i1,X,var,var1:

n:=nops(F):

X:=[seq(x[i1],i1=1..n)]:

eq:={seq(numer(x[i1]-F[i1]),i1=1..n)}:

var:=solve(eq,{op(X)}):

{seq(normal(subs(var1,X)),var1 in var)}:
end:




#CoGSFP(F,x,K1,K2): inputs a transformation F, in dimension k, say, (where k=nops(F)) phrsed in terms of x[1],x[2],..,x[k], and positive integer parameters K1,K2,A, outputs a conjectured globally stable fixed point of F
#in floating-point, but taking K1 random points in [0,1000]^k looking at the orbit to K1, and seeing if it seems to go to the same fixed points all the time. Try:
#CoGSFP([(x[1]+1)/(x[1]+6)],x,100,10);
CoGSFP:=proc(F,x,K1,K2) local k,ra,x0,L,di,pt,i,i1:
k:=nops(F):
ra:=rand(0..10^5):
x0:=evalf([seq(ra(),i=1..k)]):
L:=Orb(F,x,x0,K1):
di:=L[-1]-L[-2]:
if max(seq(abs(di[i1]),i1=1..nops(di)))>1/10^(Digits-5) then
 RETURN(FAIL):

else
pt:=L[-1]:
fi:

for i from 1 to K2 do

x0:=evalf([seq(ra(),i1=1..k)]):
L:=Orb(F,x,x0,K1):
di:=L[-1]-L[-2]:
if max(seq(abs(di[i1]),i1=1..nops(di)))>1/10^(Digits-5) then
 RETURN(FAIL):
fi:

di:=L[-1]-pt:

if max(seq(abs(di[i1]),i1=1..nops(di)))>1/10^(Digits-5) then
 RETURN(FAIL):
fi:
od:
pt:
end:






#RRDE(x,n,k,d,A): a random rational difference equation of order k, and degree d, with positive coefficients <=A. Try:
#RRDE(x,n,3,1,10);
RRDE:=proc(x,n,k,d,A) local gu,i:
gu:=RandRF(x,k,d,A):
gu:=subs({seq(x[i]=x[n+1-i],i=1..k)},gu):
gu:
end:

#OrbD(F,x,n,INI,N): Given a difference equation F or order k phrased as x[n+1]=F(x[n],...,x[n-k+1]), where R is an expression, initial condtions INI (a list of length k), and a positive integer
#N, finds the list of the first N+k members. Try:
#OrbD(x[n]+x[n-1],x,n,[1,1],10);
OrbD:=proc(F,x,n,INI,N) local gu,k,kha,i,i1:
k:=nops(INI):
gu:=INI:

for i from 1 to N do
 kha:=normal(subs({seq(x[n-i1]=gu[nops(gu)-i1],i1=0..k-1)},F)):
 gu:=[op(gu),kha]:
od:

gu:

end:


#OrbDslow(F,x,n,INI,N): Given a difference equation F or order k phrased as x[n+1]=F(x[n],...,x[n-k+1]), where R is an expression, initial condtions INI (a list of length k), and a positive integer
#N, finds the list of the first N+k members. Try:
#OrbDslow(x[n]+x[n-1],x,n,[1,1],10);
OrbDslow:=proc(F,x,n,INI,N) local gu,k,kha,i,i1:
k:=nops(INI):
gu:=INI:

for i from 1 to N do
if
 subs({seq(x[n-i1]=gu[nops(gu)-i1],i1=0..k-1)},denom(normal(F)))=0 then
 RETURN(FAIL):
else
 kha:=normal(subs({seq(x[n-i1]=gu[nops(gu)-i1],i1=0..k-1)},F)):
 gu:=[op(gu),kha]:
fi:
od:

gu:

end:

#PerSolsD(F,x,n,k,p): Inputs a k-th order difference equation in the form x[n+1]=F(x[n],...,x[n-k+1]), for some expression F, and
#finds all solutions of period p. When p=1, it is the fixed point. Try:
#F:=RRDE(x,n,3,1,10);PerSolsD(F,x,n,3,1);
PerSolsD:=proc(F,x,n,k,p) local y,INI,i,eq,var,var1,gu,lu,lu1,mu,mu1:
INI:=[seq(y[i],i=1..k)]:
gu:=OrbD(F,x,n,INI,k+p):

eq:=numer(normal({seq(gu[i]-gu[i+p],i=1..k)})):


var:={op(INI)}:


var1:=[solve(eq,var)]:

lu:={}:

for i from 1 to nops(var1) do
 lu:=lu union {subs(var1[i],INI)}:
od:


mu:={}:
for lu1 in lu do
mu1:=OrbDslow(F,x,n,lu1,p-k):
if mu1<>FAIL then
 mu1:=[op(1..p,mu1)]:
 mu:=mu union {mu1}:
fi:
od:
end:


#Coeffs(P,X): Given a polynomial in the list of variables X, outputs the set of coefficients. Try:
#Coeffs((x+y)^3,[x,y]);
Coeffs:=proc(P,X) local P1,X1,x,lu,lu1,i:
P1:=expand(P):
x:=X[nops(X)]:
if nops(X)=1 then
RETURN({seq(coeff(P1,x,i),i=0..degree(P1,x))}):
fi:
X1:=[op(1..nops(X)-1,X)]:
lu:={seq(coeff(P1,x,i),i=0..degree(P1,x))}:
{seq(op(Coeffs(lu1,X1)),lu1 in lu)}:
end:


#IsPer1(F,x,n,k,p): Is the recurrence F of order k, in terms of x[n],...x[n-k+1] periodic of period p?. Try:
#IsPer1(1/x[n],x,n,1,2);
IsPer1:=proc(F,x,n,k,p) local INI,i,X,L,i1:
INI:=[seq(X[i],i=1..k)]:
L:=normal(OrbD(F,x,n,INI,p)):
evalb({seq(numer(normal(L[i1+p]-L[i1])),i1=1..k)}={0} ):
end:



#IsLin(F,x,n,k): Is the recurrence F of order k, in terms of x[n],..,x[n-k+1] linear?
IsLin:=proc(F,x,n,k) local B,i:
B:=denom(normal(F)):
if degree(B,{seq(x[n+1-i],i=1..k)})=0 then
 true:
else
 false:
fi:
end:

#PerRecs(k,p,x,n,a,b): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k)+a[0])/add(b[i]*x[n+1-i],i=1..k)+b[0]). Try:
#PerRecs(1,2,x,n,a,b); PerRecs(2,4,x,n,a,b);PerRecs(2,5,x,n,a,b);PerRecs(2,5,x,n,a,b); PerRecs(3,6,x,n,a,b); PerRecs(2,6,a,b);PerRecs(3,8,x,n,a,b);
PerRecs:=proc(k,p,x,n,a,b) local i,F,INI,L,vars,Eqs,eqs,vars1,vars11,hal,ka,ka1,GU1,GU2:
INI:=[seq(x[i],i=1..k)]:
F:=(add(a[i]*x[n+1-i],i=1..k)+a[0])/(add(b[i]*x[n+1-i],i=1..k)+b[0]):
L:=normal(OrbD(F,x,n,INI,p)):
vars:={seq(a[i],i=0..k),seq(b[i],i=0..k)}:

Eqs:={seq(numer(normal(L[i+p]-L[i])),i=1..k)}:

eqs:={seq(op(Coeffs(Eqs[i],INI)),i=1..k)}:

vars1:=[solve(eqs,vars)]:


GU1:={}:
GU2:={}:

for i from 1 to nops(vars1) do
vars11:=vars1[i]:

if expand(subs(vars11,numer(F)))<>0 and expand(subs(vars11,denom(F)))<>0 then
hal:=normal(subs(vars11,F)):
if IsPer1(hal,x,n,k,p) then

ka:=divisors(p) minus {p}:


if not member(true,{seq(IsPer1(hal,x,n,k,ka1), ka1 in ka)}) then
if IsLin(hal,x,n,k) then
 GU2:=GU2 union {hal}:
else
 GU1:=GU1 union {hal}:
fi:
fi:
 fi:
fi:
od:
[GU1,GU2]:



end:

#PerRecsS(k,p,x,n,a): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k-1)+a[0])/x[n-k+1]). Try:
#PerRecsS(2,5,x,n,a); 
PerRecsS:=proc(k,p,x,n,a) local i,F,INI,L,vars,Eqs,eqs,vars1,vars11,hal,ka,ka1,GU1,GU2:
INI:=[seq(x[i],i=1..k)]:
F:=(add(a[i]*x[n+1-i],i=1..k-1)+a[0])/x[n-k+1]:
L:=normal(OrbD(F,x,n,INI,p)):
vars:={seq(a[i],i=0..k-1)}:

Eqs:={seq(numer(normal(L[i+p]-L[i])),i=1..k)}:

eqs:={seq(op(Coeffs(Eqs[i],INI)),i=1..k)}:

vars1:=[solve(eqs,vars)]:


GU1:={}:
GU2:={}:

for i from 1 to nops(vars1) do
vars11:=vars1[i]:

if expand(subs(vars11,numer(F)))<>0 and expand(subs(vars11,denom(F)))<>0 then
hal:=normal(subs(vars11,F)):
if IsPer1(hal,x,n,k,p) then

ka:=divisors(p) minus {p}:


if not member(true,{seq(IsPer1(hal,x,n,k,ka1), ka1 in ka)}) then
if IsLin(hal,x,n,k) then
 GU2:=GU2 union {hal}:
else
 GU1:=GU1 union {hal}:
fi:
fi:
 fi:
fi:
od:
[GU1,GU2]:



end:







#PerRecsSS1(k,p,x,n): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k-1)+a[0])/x[n-k+1]). 
#where a[0],a[1],...,a[k-1] are in {0,1}.
#Try:
#PerRecsSS1(2,5,x,n,a); 
PerRecsSS1:=proc(k,p,x,n) local INI,S,gu,s,F,s1,L,L1,ra,ka,ka1,i:
option remember:
ra:=rand(1..20):
INI:=[seq(x[i],i=1..k)]:
S:=powerset({1,seq(x[n+1-i],i=1..k-1)}) minus {{}}:

gu:={}:

for s in S do
F:=add(s1, s1 in s)/x[n-k+1]:
L1:=OrbD(F,x,n,[seq(ra(),i=1..k)],p):

if {seq(normal(L1[i+p]-L1[i]),i=1..k)}={0} then

L:=normal(OrbD(F,x,n,[seq(x[i],i=1..k)],p)):

if {seq(normal(L[i+p]-L[i]),i=1..k)}={0} then
gu:=gu union {F}:
fi:
fi:
od:

gu:

end:


#PerRecsSS(k,p,x,n): all the rational difference equations of order k with period p, of the form x[n+1]=(add(a[i]*x[n+1-i],i=1..k-1)+a[0])/x[n-k+1]). 
#where a[0],a[1],...,a[k-1] are in {0,1}.
#Try:
#PerRecsSS(2,5,x,n,a); 
PerRecsSS:=proc(k,p,x,n) local ka,p1:
ka:=divisors(p) minus {p}:
PerRecsSS1(k,p,x,n) minus {seq(op(PerRecsSS1(k,p1,x,n)),p1 in ka)}:
end:

#End Difference Equations directly



#InvLy(p,x,n): The invariant (eq. (5.12) in Kulenovic/Ladas) for the generalized Lynnes eq. x[n+1]=(p+x[n])/x[n-1]. Try:
#InvLy(p,x,n);
InvLy:=proc(p,x,n) local gu:
gu:=(p+x[n-1]+x[n])*(1+x[n-1])*(1+x[n])/(x[n-1]*x[n]):
if normal(subs(x[n+1]=(p+x[n])/x[n-1],subs(n=n+1,gu)) -gu)<>0 then
 RETURN(FAIL):
fi:
gu:
end:



#CheckInv(F,x,n,Inv): Checks that the invariant Inv for the recurrence given by F. Try:
#CheckInv((p+x[n])/x[n-1],x,InvLy(p,x,n));
CheckInv:=proc(F,x,n,Inv):
evalb( normal(subs(x[n+1]=F,subs(n=n+1,Inv)) -Inv)=0):

end:


#Vecs(k,d): The set of all vectors [a1,...,ak] of non-neg. integers that add-up to <=d. Try:
#Vecs(3,3);
Vecs:=proc(k,d) local gu,d1,gu1,i,gu11:
option remember:
if k=1 then
 RETURN({seq([i],i=0..d)}):
fi:

gu:={}:

for i from 0 to d do
 gu1:=Vecs(k-1,d-i):
 gu:=gu union {seq([i,op(gu11)],gu11 in gu1)}:
od:
gu:
end:

#GenPol(x,n,k,d,c): A generic polynomial in x[n],...,x[n-k+1] of degree d, followed by its set of coefficients, c[i] indexed starting at c[k]. Try:
#GenPol(x,2,2,c,1);
GenPol:=proc(x,n,k,d,c) local gu,v,i:
gu:=Vecs(k,d):
add(c[op(v)]*mul(x[n+1-i]^v[i],i=1..k), v in gu),{seq(c[op(v)], v in gu)}:
end:

#FindInv(F,x,n,k,d): Given a recurrence F phrased in terms of x[n], or order k, tries to find a polynomial P(x[n],x[n-1],..,x[n-k+1]) such that
#P(x[n],x[n-1],..,x[n-k+1])/(x[n]*x[n-1]*...*x[n-k+1]) is an invariant. It outputs the invariant. Try
#FindInv((p+x[n])/x[n-1],x,n,2,3);
FindInv:=proc(F,x,n,k,d) local gu,c,P,var,lu,i,eq,var1,var11, yof:
gu:=GenPol(x,n,k,d,c):
P:=gu[1]/mul(x[n+1-i],i=1..k):
var:=gu[2]:
lu:=numer(normal(subs(x[n+1]=F,subs(n=n+1,P)-P))):
eq:=Coeffs(lu,[seq(x[n+1-i],i=1..k)]):

var1:=solve(eq,var): 

yof:={}:

for var11 in var1 do
if op(1,var11)=op(2,var11) then
 yof:=yof union {op(1,var11)}:
fi:
od:


P:=normal(subs(var1,P)):

P:={seq(factor(coeff(P,i,1)), i in yof)}:

P minus {1}:
end:



#FindInvG(F,x,n,k,d,g): Given a recurrence F phrased in terms of x[n], or order k, tries to find a polynomial P(x[n],x[n-1],..,x[n-k+1]) such that
#P(x[n],x[n-1],..,x[n-k+1])/(x[n]*x[n-1]*...*x[n-k+1])^g is an invariant. It outputs the invariant. Try
#FindInvG((p+x[n])/x[n-1],x,n,2,3,1);
FindInvG:=proc(F,x,n,k,d,g) local gu,c,P,var,lu,i,eq,var1,var11, yof:
gu:=GenPol(x,n,k,d,c):
P:=gu[1]/mul(x[n+1-i],i=1..k)^g:
var:=gu[2]:
lu:=numer(normal(subs(x[n+1]=F,subs(n=n+1,P)-P))):
eq:=Coeffs(lu,[seq(x[n+1-i],i=1..k)]):

var1:=solve(eq,var): 

yof:={}:

for var11 in var1 do
if op(1,var11)=op(2,var11) then
 yof:=yof union {op(1,var11)}:
fi:
od:


P:=normal(subs(var1,P)):

P:={seq(factor(coeff(P,i,1)), i in yof)}:

P:
end:

#OurMax(f,x): Given a polynomial f in the list of variables x, finds its maximum over the positive orthant. Try
#OurMax(3-x^2-y^2,[x,y]);
OurMax:=proc(f,x) local f1,i,eq,var,lu,IntM,M,x1,ka:
f1:=evalf(f):
if nops(x)=1 then
RETURN(maximize(f1,x[1]=0..infinity)):
fi:

eq:={seq(diff(f1,x[i]),i=1..nops(x))}:
var:={op(x)}:
lu:=fsolve(eq,var):

if not type(lu,set) then
 RETURN(FAIL):
fi:

if max(op(subs(lu,x)))>=0 then
IntM:=subs(lu,f1):
else
 IntM:=-infinity:
fi:

M:=IntM:

for i from 1 to nops(x) do
 x1:=[op(1..i-1,x),op(i+1..nops(x),x)]:
ka:=OurMax(subs(x[i]=0,f1),x1):
 if ka=FAIL then
 RETURN(FAIL):
else
 M:=max(M,ka):
fi:
od:
M:
end:


##start dealing directly with difference equations, not via the translation to tranformations in R^k done via Targem.




#PaperC(x,n,d,A,K): An article with K Conjectured about Global Asymptotic Stability for n-dimensional rational transformation of [0,infinity]^n into itself
#by picking random raional transformations with positive coeficients with numerator and denominators of degree d and coefficients from 1 to A
#tolerating up to k iterations.
#Try:
#PaperC(x,2,1,20,10);
PaperC:=proc(x,n,d,A,K) local co,T,S,lu,t0,i,R:
t0:=time():
co:=0:


S:={}:

print(K, ` Conjectured Gloabal Asymptotically Stable Fixed Points of Certain Transormations from the positive orthant of`, R^n, `to itself `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
while co<K do

T:=RRT(x,n,d,A):


  lu:=CoGSFP(T,x,1000,100):
 if lu<>FAIL then
 co:=co+1:
  S:=S union {T}:
print(``):
print(`----------------------------------`):
print(``):
print(`Conjecture Number`, co):
print(``):
print(`The transformation`, [seq(x[i],i=1..n)], `goes to `, T, `has the following global attractor in the positive orthant.`):
print(``):
lprint(lu):
print(``):
print(`----------------------------------`):
print(``):

 fi:


od:
print(``):
print(`-------------------------`):
print(``):
print(`This ends this papeer that took`, time()-t0, `seconds to produce. `):
print(``):

end:



#FindPerRecs(K,P,x,n): A list of length K whose ith entry is the list of pairs [p,F] where p is a period and F is the corresponding set of recurrence of period p. Try:
#FindPerRecs(4,20,x,n);
FindPerRecs:=proc(K,P,x,n) local gu,k,p,gu1:
gu:=[]:

for k from 1 to K do
gu1:=[]:

 for p from 1 to P do
 if PerRecsSS(k,p,x,n)<>{} then
  gu1:=[op(gu1),[p, PerRecsSS(k,p,x,n)]]:
 fi:
 od:
gu:=[op(gu),gu1]:
od:
gu:
end:


#FindPerRecsA(K,x,n): A list of length K whose ith entry is the list of pairs [p,F] where p is a period <=3k and F is the corresponding set of recurrence of period p. Try:
#FindPerRecs(4,20,x,n);
FindPerRecsA:=proc(K,x,n) local gu,k,p,gu1:
gu:=[]:

for k from 1 to K do
gu1:=[]:

 for p from 1 to 3*k do
 if PerRecsSS(k,p,x,n)<>{} then
  gu1:=[op(gu1),[p, PerRecsSS(k,p,x,n)]]:
 fi:
 od:
gu:=[op(gu),gu1]:
od:
gu:
end:



#FPd(F,x,n,k): inputs a recurrence F or order k, expressesd as x[n+1]=F where F is a rational function of x[n],..., x[n-k+1], outputs the
#stable fixed point (if it exists) that is positive, let's call it pt, followed by the recurrence obeyed by x[n]-pt, that hopefully
#goes to zero. Try:
#FPd((1+x[n])/(1+x[n-1]),x,n,2);
FPd:=proc(F,x,n,k) local y,i,lu,sol,sol1,pt,F1,z:
lu:=normal(y-subs({seq(x[n+1-i]=y,i=1..k)},F)):
sol:=[solve(lu,y)]:
sol:
sol1:=[]:

for i from 1 to nops(sol) do
 if coeff(sol[i],I,1)=0 and evalf(coeff(sol[i],I,0))>0 then
 sol1:=[op(sol1),sol[i]]:
 fi:
od:

if nops(sol1)<>1 then
 RETURN(FAIL):
fi:

pt:=sol1[1]:

F1:=subs(z=x,normal(subs({seq(x[n+1-i]=z[n+1-i]+pt,i=1..k)},F)-pt)):
[pt,F1]:
end:


#NormOrbD(F,x,n,INI,K): Given a recurrence of order k, x[n+1]=F(x[n],...,x[n-k+1]), initial conditions INI, such that 0 is a fixed point of F
#finds the first K terms of the quantities L[k*i+1]^2+...+L[k*i+k]^2. Try
#NormOrbD((x[n]-x[n-1])/(2+x[n-1]),x,n,k,[3.,5.],10);
NormOrbD:=proc(F,x,n,INI,K) local k,L,i,j:
k:=nops(INI):

if simplify(subs({seq(x[n+1-i]=0,i=1..k)},F))<>0 then
 RETURN(FAIL):
fi:

L:=OrbD(F,x,n,INI,3*K):

normal([seq(add(L[i*k+j]^2,j=1..k),i=0..K-1)]):
end:


#ProveSFPd(F,x,n,k): Finds the stable fixed point of the recurrence of order k given by F, and proves it rigorously. TryL
#ProveSFPd((1+x[n])/(1+x[n-1]),x,n,2)
ProveFPd:=proc(F,x,n,k) local pt,F1,a,i,L,lu:
option remember:
pt:=FPd(F,x,n,k):
if pt=FAIL then
 RETURN(FAIL):
fi:
F1:=pt[2]:
pt:=pt[1]:
L:=NormOrbD(F1,x,n,[seq(a[i],i=1..k)],2):
lu:=numer(normal(1-L[2]/L[1])):

 if min(op(Coeffs(evalf(lu),[seq(a[i],i=1..k)])))>=0 then
[pt,1]:
elif Optimization[Minimize](evalf(lu),assume=nonnegative)[1]>=0 then
 [pt,1/2]:
else
[pt,0]:
fi:


end:



#ProveSFPdV(F,x,n,k): Verbose form of ProveFPd(F,x,n,k) (q.v.). Try:
#ProveSFPdV((1+x[n])/(1+x[n-1]),x,n,2)
ProveFPdV:=proc(F,x,n,k) local pt,F1,a,i,L,lu,z,gu,L1:
option remember:
gu:=ProveFPd(F,x,n,k):
if gu=FAIL or gu[2]=0 then
 RETURN(FAIL):
fi:

pt:=gu[1]:

print(`Every solution of the recurrence`, x[n+1]=F, `with positive initial conditions converges to`, pt, `that in floating-point is`, evalf(pt,10)):

F1:=FPd(F,x,n,k)[2]:

print(`We have to prove that for arbitrary initial conditions, the sequence x[n] converges to`, pt):
print(``):
print(`Let's make the change of variable `):
print(``):
print(x[n]=pt+z[n]):
print(``):
F1:=subs(x=z,F1):
print(``):
print(`then z[n] satisfies the recurrence`):
print(``):
print(z[n+1]=F1):
print(``):
print(`and in decimals: `):
print(``):
print(z[n+1]=evalf(F1,10)):
print(``):

print(`We have to prove that z[n] always tends to 0`):
print(``):
print(` Let `, [seq(z[i],i=1..k)], `be `, k , `consecutive terms of the solution sequence `):
print(``):
L1:=normal(OrbD(F1,z,n,[seq(z[i],i=1..k)],k)):
print(``):
print(`then the next`, k, `terms let's call them `, [seq(z[i],i=k+1..2*k)] ):
print(``):
print([op(k+1..2*k,L1)]):
print(``):
print(`We now claim that for all `, seq(z[i],i=1..k)):
print(``):
print(1-add(z[i]^2,i=k+1..2*k)/add(z[i]^2,i=1..k)>0):
print(``):

L:=NormOrbD(F1,z,n,[seq(z[i],i=1..k)],2):
lu:=normal(1-L[2]/L[1]):
print(`plugging it in, the right side is`, lu):
print(``):
print(`We have to prove that it is always non-negative.`):
print(``):
print(`Since the denominator is obviously positive, we have to prove that the numerator`):
print(``):
lprint(numer(lu)):
print(``):
print(`is non-negative `):
print(``):
print(`in decimals it is`):
lprint(evalf(numer(lu))):

if gu[2]=1 then
print(`Since all the coefficients are positive (check!), this is obviously the case. QED.`):
else
print(`This is a routine multivariable calculus exercise, that Maple kindly does for us. Alas, it uses numerics, so we only have a semi-rigorous proof, that however`):
print(`can be made fully rigorous given more time.`):
print(``):
fi:


end:





#PaperDold(A,k,K): 
#An article with K theorems about the limits of randomly chosen rational recurrences of order k, 
#with posivie coefficients from 1 to A (both numerator and denominator). Finding the
#numbers that they converge to, regardless of initial conditions and proving
#rigorously (and sometimes semi-rigorously) that they indeed  converge to the same number. Try:
#PaperDold(20,2,10);
PaperDold:=proc(A,k,K) local co,F,t0,n,x,gu:

t0:=time():
co:=0:

print(K, `theorems regarding the Global limits of certain randomly-chosen rational recurrences of order`, k, `with positive coefficientgs. `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
while co<K do
F:=RRDE(x,n,k,1,A):
gu:=ProveFPd(F,x,n,k):
if gu[2]=1 or gu[2]=1/2 then
 co:=co+1:
print(``):
print(`------------------------------------------`):
print(``):
print(`Theorem number`, co):
print(``):
ProveFPdV(F,x,n,k):
print(``):
fi:
od:
print(`---------------------------------------------------------`):
print(``):
print(`This ends this paper that took`,  time()-t0, `seconds to make. `):
print(``):
end:



#PaperD2(A,K1,K2): A paper with (at most K2) rigorously proved theorem about the limits of second-order rational recurrences. K1 is the complexity parameter. Try:
#PaperD2(30,10,4);
PaperD2:=proc(A,K1,K2) local co,F,n,x,a,b,T,gu:
co:=0:

print(K2, `Rigorously Proved Theorems about the global limits of certain second-order rational recurreces with arbitrary posiitive initial conditions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for a from 2 to trunc(A/2)  do
 for b from 3 to trunc(A/2)-3  while co<K2 do
  F:=NiceRec(x,n,[a],[b],A):
  T:=Targem(F,x,n,2):
  gu:=TryOut(T,x,K1):
  if gu<>FAIL then
   co:=co+1:
 print(``):
print(`Theorem number`, co, `: For any positive initial conditions x[1], x[2], the recurrence `):
print(``):
print(x[n+1]=F ):
print(``):
print(`converges to`, gu[1][1]):
print(``):
print(`Proof: Introducing the corresponding transformation from [0,infinity]^2 into [0,infinity]^2`):
print(`[x[1],x[2]] goes to`, T):
print(``):
print(`We have to prove that its orbit always converges to`, gu[1], `no matter where you start. `):
print(``):
print(`So let's prove the following`):
print(``):
GSFPv(T,x,gu[2]):
print(``):
fi:
od:
od:



end:


#PaperD2sr(A,K1,K2): A paper with (at most K2) rigorously proved theorem about the limits of second-order rational recurrences. K1 is the complexity parameter. Try:
#PaperD2sr(30,10,4);
PaperD2sr:=proc(A,K1,K2) local co,F,n,x,a,b,T,gu:
co:=0:

print(`Semi-Rigorously Proved Theorems about the global limits of certain second-order rational recurreces with arbitrary posiitive initial conditions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for a from 1 to trunc(A/2)  do
 for b from 1 to trunc(A/2)-3  while co<K2 do
  F:=NiceRec(x,n,[a],[b],A):
  T:=Targem(F,x,n,2):
  gu:=TryOut(T,x,K1):
  if gu<>FAIL then
   co:=co+1:
 print(``):
print(`Theorem number`, co, `: For any positive initial conditions x[1], x[2], the recurrence `):
print(``):
print(x[n+1]=F ):
print(``):
print(`converges to`, gu[1][1]):
print(``):
print(`Proof: Introducing the corresponding transformation from [0,infinity]^2 into [0,infinity]^2`):
print(`[x[1],x[2]] goes to`, T):
print(``):
print(`We have to prove that its orbit always converges to`, gu[1], `no matter where you start. `):
print(``):
print(`So let's semi-rigorously prove the following`):
print(``):
print(`Lemma: `):
print(``):
GSFPvSR(T,x,gu[2]):
print(``):
fi:
od:
od:



end:






#PaperD3(A,K1,K2): A paper with (at most K2) rigorously proved theorem about the limits of  third-order rational recurrences. K1 is the complexity parameter. Try:
#PaperD3(30,10,4);
PaperD3:=proc(A,K1,K2) local co,F,n,x,a1,a2,b1,b2,T,gu:
co:=0:

print(`WARNING: TAKES FOR EVER ON MY COMPUTER, DO NOT USE!, TRY PaperD3sr(A,K1,K2) instread `):

print(`Rigorously Proved Theorem about the global limits of certain second-order rational recurreces with arbitrary posiitive initial conditions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for a1 from 1 to trunc(A/2)  do
 for a2 from 1 to trunc((A-2-a1)/2) do
 for b1 from 2 to trunc(A/2)  do
  for  b2 from 3 to trunc((A-2-b1)/2) while co<K2 do
  F:=NiceRec(x,n,[a1,a2],[b1,b2],A):

  T:=Targem(F,x,n,3):
  gu:=TryOut(T,x,K1):
  if gu<>FAIL then
   co:=co+1:
 print(``):
print(`Theorem number`, co, `: For any positive initial conditions x[1], x[2], x[3], the recurrence `):
print(``):
print(x[n+1]=F ):
print(``):
print(`converges to`, gu[1][1]):
print(``):
print(`Proof: Introducing the corresponding transformation from [0,infinity]^3 into [0,infinity]^3`):
print(`[x[1],x[2],x[3]] goes to`, T):
print(``):
print(`We have to prove that its orbit always converges to`, gu[1], `no matter where you start. `):
print(``):
print(`So let's prove the following`):
print(``):
GSFPv(T,x,gu[2]):

print(``):
fi:
od:
od:
od:
od:



end:





#GSFPvSR(F,x,K): A semi-rigorous version of GSFPv(F,x,K)
#and the symbol x, outputs the  unique globally stable fixed point, and an integer k <=K
#that tells you the number of iterations needed until it enters a shrinking regime. 
#if there is no such k<=K it returns [pt,FAIL]
#Try:
#GSFPvSR([x[2],(1+x[2])/(1+x[1])],x,10);


GSFPvSR:=proc(F,x,K) local gu,pt,lu,k,i,a,L,R,L1,i1:
gu:=TryOut(F,x,K):
if gu=FAIL then
 RETURN(FAIL):
fi:


pt:=gu[1]:


if not (type(pt[1],integer) or type(pt[1],fraction)) then
 pt:=evalf(pt):
fi:

k:=gu[2]:

L1:=normal(Orb(F,x,[seq(a[i],i=1..nops(F))],k+1)):

L:=normal(NormOrb(F,x,[seq(a[i],i=1..nops(F))],pt,k+1)):

lu:=numer(normal(L[1]-(101/100)*L[1+k])):


print(`Lemma`):
print(` The transformation`,  [seq(x[i],i=1..nops(F))] , `goes to `, F, ` has the Global Attractor`, pt, `for the positive orthant in`, R^nops(F)):
print(``):

print(`Semi-Rigorous Proof:`):
print(``):
print(`Starting at any point in the positive orthant`, [seq(a[i],i=1..nops(F))], `the first`, k+1, `terms of the orbit are `):
print(``):
lprint(L1):
print(``):
print(`Let L be the distance-squared of the above members of the orbit to the fixed point`, pt, `Then L is `):
print(``):
lprint(L):
print(``):
print(`We claim that the`, k+1, `-th entry  times (101/100) is always less than the first entry`):
print(``):
print(`The difference between the distance-square to the point of the originating point and (101/100) times that distance`, k+1, `-th point is `):
print(``):
lu:=normal(L[1]-(101/100)*L[1+k]):
print(``):

if (type(pt[1],integer) or type(pt[1],fraction)) then
 lprint(lu):
else
 lprint(evalf(lu,10)):
fi:

print(``):

print(`The denominator is obviously positive`):
print(``):
print(`We have to prove that the numerator is non-negative in the positive orthant`):
print(``):
print(`The numerator is the polynomial`):
print(``):

lu:=numer(lu):



if (type(pt[1],integer) or type(pt[1],fraction)) then
 lprint(lu):
else
 lprint(evalf(lu,10)):
fi:


print(``):
print(`We claim that this polynomial is always non-negative in the positive orthant of`, R^nops(F)):
print(``):
print(`This is a routine exercise in multi-variable calculus, that we prefer not to do rigorously,since it would take too long, and checked it numerically for many random values`):
print(`and since we know that there exists a proof, we did not bother. Semi-Rigorous QED.`):
print(``):

end:

#PaperD3sr(A,K1,K2): A paper with (at most K2) semi-rigorously proved theorem about the limits of  third-order rational recurrences. K1 is the complexity parameter. Try:
#PaperD3sr(30,10,4);
PaperD3sr:=proc(A,K1,K2) local co,F,n,x,a1,a2,b1,b2,T,gu:
co:=0:

print(K2, `Semi-Rigorously Proved Theorems about the global limits of certain Third-order rational recurreces with arbitrary posiitive initial conditions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for a1 from 1 to trunc(A/2)  do
 for a2 from a1+1 to trunc((A-2-a1)/2) do
 for b1 from 2 to trunc(A/2)  do
  for  b2 from b1+1 to trunc((A-2-b1)/2) while co<K2 do
  F:=NiceRec(x,n,[a1,a2],[b1,b2],A):

  T:=Targem(F,x,n,3):
  gu:=TryOut(T,x,K1):
  if gu<>FAIL then
   co:=co+1:
 print(``):
print(`Theorem number`, co,`: `):
print(` For any positive initial conditions x[1], x[2], x[3], the recurrence `):
print(``):
print(x[n+1]=F ):
print(``):
print(`converges to`, gu[1][1]):
print(``):
print(`Proof: Introducing the corresponding transformation from [0,infinity]^3 into [0,infinity]^3`):
print(`[x[1],x[2],x[3]] goes to`, T):
print(``):
print(`We have to prove that its orbit always converges to`, gu[1], `no matter where you start `):
print(``):
print(`So let's (semi-rigorousl) prove the following`):
print(``):
GSFPvSR(T,x,gu[2]):

print(``):
fi:
od:
od:
od:
od:



end:


#NormOrb(F,x,x0,pt,N): Given a transformation F phrased in terms of x[1],...,x[nops(F)], an initial condition x0 (the same length of F) a point pt (of the same length of F and x0), and a positive
#integer N, outputs the first N entires of the Euclidean norm of the orbit of F minut pt. Try:
#NormOrb([x[2],(1+x[2])/(1+x[1])],x,[2,3],[1,1],10);
NormOrb:=proc(F,x,x0,pt,N) local gu,i:
gu:=Orb(F,x,x0,N):
[seq(Norm2(gu[i]-pt)^2,i=1..nops(gu))]:
end:




#GSFP(F,x,K): inputs a transformation F, given as a list of length k, say, with an expressions in x[1],...,x[k]
#and the symbol x, outputs the  unique globally stable fixed point, and an integer k <=K
#that tells you the number of iterations needed until it enters a shrinking regime. 
#if there is no such k<=K it returns [pt,FAIL]
#Try:
#GSFP([x[2],(1+x[2])/(1+x[1])],x,10);

GSFP:=proc(F,x,K) local gu,pt,lu,k,i,a,L:
option remember:
gu:=LSFP(F,x):


if gu={} then
# print(`There are no locally stable fixed points`):
 RETURN(FAIL):
fi:

if nops(gu)>1 then
 print(`There are more than one locally stable fixed points. Here there are`):
 print({seq(gu[i][1],i=1..nops(gu))}):
 RETURN(FAIL):
fi:

gu:=TryOut(F,x,K):

if gu=FAIL then
 RETURN(FAIL):
fi:


pt:=gu[1]:

k:=gu[2]:

L:=normal(NormOrb(F,x,[seq(a[i],i=1..nops(F))],pt,k+1)):
lu:=numer(normal(L[1]-(101/100)*L[1+k])):


if evalf(simplify(minimize(lu,seq(a[i]=0..infinity,i=1..nops(F)))))>=-10^(-10) then
 [pt,k]:
else
[pt,FAIL]:
fi:




end:




#CrPol(T,x,a,K): The crucial polynomial in a[1],.., a[nops(T)] whose non-negativity would imply rigorously that
#the unique locally stable equilibrium point in the positive orthant is indeed a global attractor. It also returns the point and the k<=K. Try:
#CrPol([x[2],(1+x[2])/(1+x[1])],x,a,4);
CrPol:=proc(T,x,a,K) local gu,k,L,i,pt:
gu:=TryOut(T,x,K):
if gu=FAIL then
 RETURN(FAIL):
fi:
k:=gu[2]:
pt:=gu[1]:

L:=NormOrb(T,x,[seq(a[i],i=1..nops(T))],pt,k+1):
[expand(numer(normal(L[1]-(101/100)*L[k+1]))),pt,k]:
end:

#GSFPv(F,x,K): A verbose form of GSFP(F,x,K)
#and the symbol x, outputs the  unique globally stable fixed point, and an integer k <=K
#that tells you the number of iterations needed until it enters a shrinking regime. 
#if there is no such k<=K it returns [pt,FAIL]
#Try:
#GSFPv([x[2],(1+x[2])/(1+x[1])],x,10);

GSFPv:=proc(F,x,K) local gu,pt,lu,k,i,a,L,R,L1,ku:
option remember:
gu:=GSFP(F,x,K):


if gu=FAIL then
 RETURN(FAIL):
fi:

if gu[2]=FAIL then
 RETURN(FAIL):
fi:


pt:=gu[1]:


k:=gu[2]:

L1:=normal(Orb(F,x,[seq(a[i],i=1..nops(F))],k+1)):

L:=normal(NormOrb(F,x,[seq(a[i],i=1..nops(F))],pt,k+1)):

lu:=numer(normal(L[1]-(101/100)*L[1+k])):


print(`Proof that the transformation`, [seq(x[i],i=1..nops(F))] , `goes to `, F, ` has the Global Attractor`, pt, `for the positive orthant in`, R^nops(F)):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Starting at any point in the positive orthant`, [seq(a[i],i=1..nops(F))], `the first`, k+1, `terms of the orbit are `):
print(``):
lprint(L1):
print(``):
print(`Let L be the distance-squared of the above members of the orbit to the fixed point`, pt, `Then L is `):
print(``):
lprint(L):
print(``):
print(`We claim that the first  entry  minus  101/100 times the `, k+1, `th entry is always non-negative `):
print(``):
print(`The difference between the distance-square to the point of the originating point and the distance to`, k+1, `-th point is, times 101/100 is`):
print(``):
lu:=normal(L[1]-(101/100)*L[1+k]):
print(``):
lprint(lu):
print(``):

print(`The denominator is obviously positive`):
print(``):
print(`We have to prove that the numerator is non-negative in the positive orthant`):
print(``):
print(`The numerator is the polynomial`):
print(``):
lu:=numer(lu):
if (type(pt[1],fraction) or type(pt[1],integer)) then
lprint(lu):
else
lprint(evalf(lu,10)):
fi:

print(``):
ku:=minimize(lu,a[1]=0..infinity,a[2]=0..infinity,location); 
print(``):
print(`The minimum is indeed`, ku[1], `and this takes place at`, subs(ku[2][1][1],[seq(a[i],i=1..nops(F))])  ):
print(``):

end:



#TryOut1(T,x,k): Given a transformation T phrased in terms of x[1],x[2],..,x[nops(T)], and aa positive integer k, investigates whether
#there is a unique locally stable fixed point, and then whether in the NormOrb for 100 randoml chosen initial conditions
#the NormOrb[i]-NormOrb[i+k] is  positive. Try:
#TryOut1([x[2],(1+x[2])/(1+x[1])],x,2);
TryOut1:=proc(T,x,k) local i,L,ra,gu,x0,pt,i1,ku,ku1,M,lu,yo:


gu:=LSFP(T,x):


if gu={} then
 print(`There are no locally stable fixed points`):
 RETURN(FAIL):
fi:

if nops(gu)>1 then
 print(`There are more than one locally stable fixed points. Here there are`):
 print({seq(gu[i][1],i=1..nops(gu))}):
 RETURN(FAIL):
fi:

pt:=gu[1][1]:




if  min(op(evalf(pt)))<0 then
 RETURN(FAIL):
fi:


if nops(T)=1 then
L:=evalf(NormOrb(T,x,evalf(pt)+[0.1],evalf(pt),100)):
else
L:=evalf(NormOrb(T,x,evalf(pt)+[0$(nops(T)-1),0.1],evalf(pt),100)):
fi:




if min(seq(L[i1]-(101/100)*L[i1+k],i1=1..nops(L)-k))<0 then
  RETURN([pt,false]):
fi:

ra:=rand(1..10000):

for i from 1 to 300 do
 x0:=evalf([seq(ra()/1000,i1=1..nops(T))]):

L:=evalf(NormOrb(T,x,x0,evalf(pt),100)):


if min(seq(L[i1]-(101/100)*L[i1+k],i1=1..nops(L)-k))<0 then
  RETURN([pt,false]):
fi:

od:



M:=normal(NormOrb(T,x,[seq(a[i],i=1..nops(T))],pt,k+1)):

lu:=numer(normal(M[1]-(101/100)*M[1+k])):



ku:=Kufsa([10$nops(T)]):

 yo:=min(seq(subs({seq(a[i1]=ku1[i1],i1=1..nops(T))},lu),ku1 in ku)):



if evalf(yo)<0 then
  RETURN([pt,false]):
fi:


[pt,true]:
end:


#TryOut(T,x,K): Given a transformation T phrased in terms of x[1],x[2],..,x[nops(T)], and aa positive integer K, conjectures the smallest integer k<=K
#such that
#the NormOrb[i]-NormOrb[i+k] is  positive. Or RETURNS FAIL. Try:
#TryOut([x[2],(1+x[2])/(1+x[1])],x,10);
TryOut:=proc(T,x,K) local k,gu:

for k from 1 to K do

gu:=TryOut1(T,x,k):

if gu=FAIL then
  RETURN(FAIL):
fi:
if gu[2] then
 RETURN([gu[1],k]):
fi:
od:

FAIL:
end:



#Kufsa(L): The discrete box with dimensions L, Try: Kufsa([2,2,2]);
Kufsa:=proc(L) local i,L1,x,gu,gu1:
if L=[] then
 RETURN({[]}):
fi:
L1:=[op(1..nops(L)-1,L)]:
gu:=Kufsa(L1):
{seq(seq([op(gu1),x],gu1 in gu),x=0..L[nops(L)])}:

end:




#InvLa(p,x,n): The invariant (eq. (5.12) in Kulenovic/Ladas) for the generalized Ladas eq. x[n+1]=(p+x[n]+x[n-1])/x[n-2]. Try:
#InvLa(p,x,n);
InvLa:=proc(p,x,n) local gu:
gu:=(x[n]^2*x[n-2]^2+p*x[n-1]^2+x[n]^2*x[n-2]+x[n]*x[n-2]^2+x[n]*x[n-1]^2+x[n-2]*x[n-1]^2+x[n-1]^3+p*x[n-1]+x[n]*x[n-2]+x[n]*x[n-1]+x[n-2]*x[n-1]+x[n-1]^2)/x[n]/x[n-1]/x[n-2]:
if normal(subs(x[n+1]=(p+x[n]+x[n-1])/x[n-2],subs(n=n+1,gu)) -gu)<>0 then
 RETURN(FAIL):
fi:
gu:
end:





#GLynB(p,a,b,X,Y): The theortical  minimum and maximum of the orbit of the generalized Lynnes recurrence x[n+1]=(p+x[n])/x[n-1]. with initial conditions
#x[1]=a, x[2]=b. It returns the polynimial in X whose middle roots are these, followed by their decimal values
#It also returns the smallest and largest members of the first 20000 terms
#Try:
#GLynB(2,1,1,X,Y);
GLynB:=proc(p,a,b,X,Y) local x,n,gu,lu,mu,L:

gu:=InvLy(p,x,n):

gu:=numer(gu-subs({x[n-1]=a,x[n]=b},gu)):

lu:=subs(x[n]=X,discrim(gu,x[n-1])):

mu:=[fsolve(lu,X)]:
L:=OrbD((p+x[n])/x[n-1],x,n,evalf([a,b]),20000):

gu:=subs({x[n-1]=X,x[n]=Y},gu):

[gu, lu, evalf([mu[2],mu[3]],10), evalf([min(L),max(L)],10)]:


end:



#GLynBs(p,a,b,X): inputs SYMBOLIC (or numeric) p,a,b, outputs The quartic polynomial in Z whose middle roots are the lower and upper bounds of the members of the orbit of
#the generalized Lynnes equation
#x[n+1]=(p+x[n])/x[n-1] and initial conditions x[1]=a, x[2]=b. 
#GLynBs(p,a,b,X);
GLynBs:=proc(p,a,b,X) local x,n,gu:

gu:=InvLy(p,x,n):

gu:=numer(gu-subs({x[n-1]=a,x[n]=b},gu)):

subs(x[n]=X,discrim(gu,x[n-1])):

end:




#GLaB(p,a,b,c,X,Y,Z): The polynomial in X,Y,Z, such that three consecutive terms of the orbit  generalized Ladas recurrence x[n+1]=(p+x[n]+x[n-1])/x[n-2]. with initial conditions
#x[1]=a, x[2]=b,x[3]=c must satisfy.
#It also returns the smallest and largest members of the first 20000 terms
#Try:
#GLaB(2,1,1,X,Y,Z);
GLaB:=proc(p,a,b,c,X,Y,Z) local x,n,gu,L:

gu:=InvLa(p,x,n):

gu:=numer(gu-subs({x[n-2]=a,x[n-2]=b,x[n]=c},gu)):

gu:=subs({x[n-2]=X,x[n-1]=Y,x[n]=Z},gu):
L:=OrbD((p+x[n]+x[n-1])/x[n-2],x,n,evalf([a,b,c]),20000):

[gu,  evalf([min(L),max(L)],10)]:

end:



#GLynBsV(): A verbose form of GLynBs(p,a,b,X) (q.v.), Try:
#GLynBsV();
GLynBsV:=proc() local p,a,b,n,x,X,Y,gu,lu:
gu:=collect(GLynBs(p,a,b,X),X):

print(`Upper Bounds and Lower Bounds on the members of the orbit of the generalized Lynnes Equation x[n+1]=(p+x[n])/x[n-1]`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Let a,b,p be a positive real numbers, and consider the terms of the non-linear recurrence`):
print(``):
print(x[n+1]=(p+x[n])/x[n-1]):
print(``):
print(`and in Maple notation`):
print(``):
print(x[n+1]=(p+x[n])/x[n-1]):
print(``):
print(`Subject to the  initial conditions`):
print(``):
lprint(x[1]=a, x[2]=b):
print(``):
print(`Then x[n] is always between m and M where m and M are the middle two roots of quartic equation in X`):
print(``):
print(gu=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu=0):
print(``):
print(`Proof: `):
print(``):
lu:=InvLy(p,x,n):
print(``):
print(`It is readily seen, by pluging in the recurrence that for any two consecutive terms of the orbit,x[n-1],x[n]`):
print(``):
print(`we have the following invariant `):
print(``):
print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):
print(`Plugging in we n=2, we have that x[n-1]=X, and x[n]=Y must obey`):
print(``):
lu:=subs({x[n-1]=X,x[n]=Y},lu)-subs({x[n-1]=a,x[n]=b},lu):
print(``):
print(lu=0):
print(``):
print(`simplifying, we see that `):
print(``):
print(numer(lu)=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(numer(lu)=0):
print(``):

print(`Both lower bounds and upper bounds must be roots of the discriminant with respect to Y of the above quadratic in Y, that is the above-mentioned quartic equation in X, namely`):
print(``):
print(gu=0):
print(``):
print(`Let's illustrate the theorem with p=5, a=3,b=7`):
print(``):
lu:=GLynB(5,3,7,X,Y):
print(``):
print(`The above lower and upper bounds are`):
print(``):
lprint(lu[3]):
print(``):
print(`and the minimum and maximum in the first 20000 terms are`):
print(``):
lprint(lu[4]):
print(``):
print(`The difference is`):
print(``):
lprint(lu[3]-lu[4]):
print(``):
print(`Not Bad!`):
print(``):
print(`-------------------------------------`):

end:




#GLaBv(): A paper about the boundedness of the Generalized Ladas equation x[n+1]=(p+x[n]+x[n-1])/x[n-2]. Try:
#GLaBv();
GLaBv:=proc() local p,a,b,c,n,x,X,Y,Z,lu,lu1:

lu:=InvLa(p,x,n):


print(`Upper Bounds and Lower Bounds on the members of the orbit of the generalized Ladas Equation x[n+1]=(p+x[n]+x[n-1]))/x[n-2]`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Let a,b,c,p be a positive real numbers, and consider the terms of the non-linear recurrence`):
print(``):
print(x[n+1]=(p+x[n]+x[n-1])/x[n-2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(x[n+1]=(p+x[n]+x[n-1])/x[n-2]):
print(``):
print(`Subject to the  initial conditions`):
print(``):
lprint(x[1]=a, x[2]=b,x[3]=c):
print(``):
print(`Then x[n] is always between m and M are the maximum that X can take in the bounded surface`):
print(``):
lu1:=numer(lu-subs({x[n-2]=a,x[n-1]=b,x[n]=c},lu)):
print(``):
lu1:=subs({x[n-2]=X,x[n-1]=Y,x[n]=Z},lu1):
print(``):
print(lu1=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu1=0):
print(``):

print(``):
print(``):
print(`Proof: `):
print(``):
lu:=InvLa(p,x,n):
print(``):
print(`It is readily seen, by pluging in the recurrence that for any three consecutive terms of the orbit,x[n-2]=X,x[n-1]=Y,x[n]=Z`):
print(``):
print(`we have the following invariant `):
print(``):
print(lu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(lu):
print(``):
print(`Plugging in we n=3, we have that x[n-2]=X, x[n-1]=Y, and x[n]=Z must obey`):
print(``):
lu:=subs({x[n-2]=X,x[n-1]=Y,x[n]=Z},lu)-subs({x[n-2]=a,x[n-1]=b,x[n]=c},lu):
print(``):
print(lu=0):
print(``):
print(`simplifying, we see that `):
print(``):
print(factor(numer(lu))=0):
print(``):
print(`Both lower bounds and upper bounds must be the smallest and largest values that X can take on the above bounded surface.`):
print(``):
print(`Let's illustrate the theorem with p=5, a=3,b=7,c=11`):
print(``):
lu:=GLaB(5,3,7,11,X,Y,Z):
print(``):
print(`The smallest and largest member among the first 20000 members are`):
print(``):
lprint(lu[2]):
print(``):
print(`-------------------------------------`):

end:


#NiceRec(x,n,L1,L2,N): a nice non-linear recurrence with equ. point 1, of order nops(L1)+1 (where nops(L1)=nops(L2)) 
#and L1 and L2 are increasing. Try:
#NiceRec(x,n,[3],[7],20);
NiceRec:=proc(x,n,L1,L2,N) local i1,k:

if not (type(L1,list) and type(L2,list) and nops(L1)=nops(L2) and min(seq(L1[i1+1]-L1[i1],i1=1..nops(L1)-1))>=0 and min(seq(L2[i1+1]-L2[i1],i1=1..nops(L2)-1))>=0 and
L1[nops(L1)]<=N-convert(L1,`+`)-1 and L2[nops(L2)]<=N-convert(L2,`+`)-1 ) then
 RETURN(FAIL)
fi:

k:=nops(L1)+1:

(1+add(L1[i1]*x[n-k+i1],i1=1..k-1)+(N-convert(L1,`+`))*x[n])/(1+add(L2[i1]*x[n-k+i1],i1=1..k-1)+(N-convert(L2,`+`))*x[n]):

end:



#PaperD4sr(A,K1,K2): A paper with (at most K2) semi-rigorously proved theorem about the limits of  fourth-order rational recurrences. K1 is the complexity parameter. Try:
#PaperD4sr(30,5,2);
PaperD4sr:=proc(A,K1,K2) local co,F,n,x,a1,a2,a3,b1,b2,b3,T,gu:
co:=0:

#KAKI
print(K2, `Semi-Rigorously Proved Theorems about the global limits of certain Fourth-order rational recurreces with arbitrary posiitive initial conditions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for a1 from 1 to trunc(A/2)  do
for a2 from a1+1 to trunc(A/2)  do
 for a3 from a2+1 to trunc((A-2-a1-a2)/2) do
 for b1 from 2 to trunc(A/2)  do
  for  b2 from b1+1 to trunc((A-2-b1)/2)  do
  for  b3 from b2+1 to trunc((A-2-b1-b2)/2) while co<K2 do
  F:=NiceRec(x,n,[a1,a2,a3],[b1,b2,b3],A):

  T:=Targem(F,x,n,4):
  gu:=TryOut(T,x,K1):
  if gu<>FAIL then
   co:=co+1:
 print(``):
print(`Theorem number`, co,`: `):
print(` For any positive initial conditions x[1], x[2], x[3], x[4], the recurrence `):
print(``):
print(x[n+1]=F ):
print(``):
print(`converges to`, gu[1][1]):
print(``):
print(`Proof: Introducing the corresponding transformation from [0,infinity]^4 into [0,infinity]^4`):
print(`[x[1],x[2],x[3],x[4]] goes s to`, T):
print(``):
print(`We have to prove that its orbit always converges to`, gu[1], `no matter where you start `):
print(``):
print(`So let's (semi-rigorousl) prove the following`):
print(``):
GSFPvSR(T,x,gu[2]):

print(``):
fi:
od:
od:
od:
od:
od:
od:



end:




#FindInvStory(K): The invariants of the rational recurrences
#x[n+1]=(p+x[n]+x[n-1]+...+x[n-k+1])/x[n-k+1]. Try:
#FindInvStory(4);
FindInvStory:=proc(K) local k,F,gu,x,n,i1,p,L:

print(``):
print(`Invariants for recurrences of the form x[n+1]=(p+x[n]+...+x[n-k+1])/x[n-k] for k from 1 to`, K):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for k from 1 to K do
print(``):
F:=(p+add(x[n-i1],i1=0..k-2))/x[n-k+1]:
print(``):
gu:=FindInv(F,x,n,k,k+1):
if gu={} then
 RETURN(FAIL):
fi:
 gu:=gu[1]:

print(``):
print(`-------------------------------`):
print(``):

print(`Theorem Number:`, k, `For p positive, the recurrence `):
print(``):
lprint(x[n+1]=F):
print(``):
print(`has the following invariant`):
print(``):
print(gu):
print(``):
print(`and in Maple format: `):
print(``):
lprint(gu):
print(``):
print(`Proof: Routine`):
print(``):
print(`Corollary: For any positive initial conditions the orbit is bounded`):
print(``):
print(`Let's illustate it with the case p=2 and initial conditions`, seq(x[i1]=i1+2,i1=1..k)):
print(``):
L:=OrbD(subs(p=2,F),x,n,  evalf([seq(i1+2,i1=1..k)]),1000):

print(``):
print(`The smallest entry among the first 1000 terms is`, min(L)):
print(``):
print(`while the largest is`, max(L)):
print(``):
od:

print(``):
print(`--------------------------`):
print(``):
end: