######################################################################
##RDS.txt: Save this file as  RDS.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read RDS<Enter>                                    #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Emilie Hogan (now Purvine)                              # 
#and Doron Zeilberger, Rutgers University ,                          #
######################################################################
 
#Created: May 21, 2010
#Version of May 2023 
print(`Created: May 21, 2010`):
print(` This is RDS.txt `):
print(`It is one of the packages that accompany the article `):
print(` A New Algorithm for Proving Global Asymptotic Stability of Rational Difference Equations`):
print(`by Emilie Hogan (now Purvine) and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

_EnvAllSolutions:=true:

#===========================================================================
ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: FindKline, IsKgood `):
 print(`  MaxCi, MaxF1, MaxF11,`):
 print(`MaxF1a, MaxFline, MaxFnaive, MaxYakhas,`):
 print(`  MinPol, MySolve2 `):

else
ezra(args):
fi:

end:

#===========================================================================
ezra:=proc()

if args=NULL then
 print(`The main procedures are: CGA, ConjK, ConjKnew `):
 print(` EfesT,   FindGoodks, FindGoodksR, Findk, FixFoint, Investigate `):
 print(`  Kline, Kpt, KptS, LSFP, LSFPs, MaxF, MinPol, ProveK, ProveKar`):
 print(` ProveKnr, `):
 print(`  Traj, Yakhas, YakhasS, Y1k, Y2k, Y3k `):



elif nops([args])=1 and op(1,[args])=CGA then
print(`CGA(R,x,y,K,L,N): Given a rational recurrence x[n+1]=R(x[n],x[n-1])`):
print(`in terms of a rational function R(x,y) conjectures`):
print(`a global attractor, if it exists.`):
print(`K and L are positive integers. We pick L randomly chosen`):
print(`initial points in [0,K]^2 and look at the trajectory N`):
print(`times and see if they all seem to converge to the global `):
print(`attractor .`):
print(`For example, try: `):
print(` CGA((24+3*y)/(1+x),x,y,10,20,300); `):


elif nops([args])=1 and op(1,[args])=ConjK then
print(`ConjK(R,x,y,Maxk): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]) finds a point that seems to be`):
print(`a global attractor and conjectures a k for it`):
print(`For example, try:`):
print(`ConjK((24+3*y)/(1+x),x,y,20);`):

elif nops([args])=1 and op(1,[args])=ConjKnew then
print(`ConjKnew(R,x,y,Maxk,L,T): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]) finds a point that seems to be`):
print(`a global attractor and conjectures a k for it`):
print(`1/10^L is the accuracy and T is the parameter for the number of`):
print(`initial points in the hill-climbing algorithm`):
print(`For example, try:`):
print(`ConjKnew((1+y)/(1+x),x,y,20,2,2);`):

elif nops([args])=1 and op(1,[args])=EfesT then
print(`EfesT(R,x,y): Given a function R(x,y) first finds out`):
print(`whether it is a locally-stable fixed point and then`):
print(`it returns the transformed transformation with 0 as fixed point`):
print(`For example, try:`):
print(`EfesT((24+3*y)/(1+x),x,y);`):


elif nops([args])=1 and op(1,[args])=FindGoodks then
print(`FindGoodks(R,x,y,Maxk,cu,KAMA): Given a rational function R(x,y)`):
print(`describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to`):
print(`itself, finds the fixed point, whether it is locally asympototically`):
print(`stable, and then finds all k<=Maxk that are shrinking exponent in the`):
print(`sense that every k iterations shrings the norm by <=cu.`):
print(`warning: this is only a guess! `):
print(`FindGoodks((1+y)/(1+x),x,y,10,9/10,1000);`):


elif nops([args])=1 and op(1,[args])=FindGoodksR then
print(`FindGoodksR(R,x,y,Maxk,cu,a,b,RangeA,RangeB,res) :`):
print(`Given a rational recurrence x[n]=R(x[n-1],x[n-2])`):
print(`phrased in terms of parameters a and b`):
print(`finds all good k's <=Maxk that are valid for all`):
print(`(a,b) in the rectangle in Parameter space `):
print(` Range A x RangeB with resolution res `):
print(` For example, try: `):
print(` FindGoodksR(Y2k(M+q,q,x,y),x,y,10,9/10,M,q,[1,3],[1,3],1) ; `):

elif nops([args])=1 and op(1,[args])=Findk then
print(`Findk(R,x,y,Maxk,cu): Given a rational function R(x,y)`):
print(`describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to`):
print(`itself, finds the fixed point, whether it is locally asympototically`):
print(`stable, and then finds the smallest k<=Maxk that is `):
print(` a shrinking exponent in the`):
print(`sense that every k iterations shrings the norm by <=cu.`):
print(`and proves it non-rigorously. `):
print(`if there is no such k<=Maxk it returns FAIL`);
print(`Findk((1+y)/(1+x),x,y,10,9/10);`):

elif nops([args])=1 and op(1,[args])=FindKold then
print(`FindKold(R,x,y,cu,Maxk,K,eps1,M): given a`):
print(` rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]), and a real number r1`):
print(`using K points between [-a,r1+a] and [r1+a,-a]`):
print(`retuns a table for the best k for each line`):
print(`x+y=r1 starting at r1=-2*a with intervals of eps1`):
print(`FindKold((24+3*y)/(1+x),x,y,.9,20,10,1,20);`):

elif nops([args])=1 and op(1,[args])=FindKline then
print(`FindkLine(R,x,y,r1,cu,Maxk,K): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]), and a real number r1`):
print(`and a pos. integer Maxk, finds a k <=Maxk`):
print(`that seems to be a good one for all x+y=r1`):
print(`using K points between [-a,r1+a] and [r1+a,-a]`):
print(`where a is the equi.pt. of R`):
print(`If it fails it returns FAIL`):
print(`Try:`):
print(`FindKline((24.+3*y)/(1.+x),x,y,3.,.9,10,100);`):

elif nops([args])=1 and op(1,[args])=FixPoint then
print(`FixPoint(R,x,y): Given a rational recurrence x[n+1]=R(x[n],x[n-1])`):
print(`in terms of a rational function R(x,y) finds its fixed points.`):
print(`For example, try:`):
print(`FixPoint((24+3*y)/(1+x),x,y);`):


elif nops([args])=1 and op(1,[args])=Investigate then
print(`Investigate(R,x,y,a,b,Ra,Rb,k): investigates`):
print(`whether a recurrence x[n]=R(x[n-1],x[n-2])`):
print(`given by a rational function R(x,y) and`):
print(`feauturing parameters a and b,`):
print(`is k-contracting (with contracting coeff.<cu) `):
print( `in Ra times Rb with`):
print(`the given fixed-point. For example, try:`):
print(`Investigate(Y1kN(a,b,x,y),x,y,a,b,[1,5],[1,5],2,0,9/10);`):

elif nops([args])=1 and op(1,[args])=IsKgood then
print(`IsKgood(R,x,y,r1,cu,k,K): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]), and a real number r1`):
print(`and a k, that seems to be a good one for all x+y=r1`):
print(`using K points between [-a,r1+a] and [r1+a,-a]`):
print(`decides whether this k seems to be a good one for the line x+y=r1`):
print(`Try:`):
print(`IsKgood((24+3*y)/(1+x),x,y,3,.9,2,10);`):

elif nops([args])=1 and op(1,[args])=Kline then
print(`Kline(R,x,y,r1,Maxk,cu,eps1): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]), and a line x+y=r1`):
print(`finds the smallest k <=Maxk such that iterating it k times`):
print(`followed by the ratio`):
print(`at pt shrinks the norm by at least cu. `):
print(`eps1 is the resolution. For example, try:`):
print(`Kline((24+3*y)/(1+x),x,y,3,20, .9,.1);`):

elif nops([args])=1 and op(1,[args])=Kpt then
print(`Kpt(R,x,y,pt,Maxk,cu): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]), and a point pt,`):
print(`finds the smallest k <= Maxk such that iterating it k times`):
print(`at pt shrinks the norm by at least cu`):
print(` followed by the ratio. For example, try:`):
print(`Kpt((24+3*y)/(1+x),x,y,[1,2],20,.9);`):

elif nops([args])=1 and op(1,[args])=KptS then
print(`KptS(R,x,y,pt,Maxk): given a rational transformation R(x,y)`):
print(`describing a second-order recurrence`):
print(`x[n]=R(x[n-1],x[n]), and a point pt,`):
print(`finds all the k from 1 to Maxk such that iterating it k times`):
print(`at pt shrinks the norm, , followed by the ratio. For example, try:`):
print(`KptS((24+3*y)/(1+x),x,y,[1,2],20);`):

elif nops([args])=1 and op(1,[args])=LSFP then
print(`LSFP(R,x,y): Given a rational recurrence x[n+1]=R(x[n],x[n-1])`):
print(`in terms of a rational function R(x,y) finds its `):
print(`locally stable fixed points`):
print(`For example, try:`):
print(`LSFP((24+3*y)/(1+x),x,y);`):


elif nops([args])=1 and op(1,[args])=LSFPs then
print(`LSFPs(R,x,y,para): Given a rational recurrence x[n+1]=R(x[n],x[n-1])`):
print(`in terms of a rational function R(x,y) finds its `):
print(`and a parameter set para, finds the fixed points`):
print(`and for each the conditions that (i) it is positive`):
print(`(ii) that it is locally stable.`):
print(`For example, try:`):
print(`LSFPs(Y2k(M,q,x,y),x,y,{M,q});`):



elif nops([args])=1 and op(1,[args])=MaxCi then
print(`MaxCi(F,x,y,a,R,K): The maximum location and value of the`);
print(`function F(x,y) on the part of the circle`):
print(`x^2+y^2=R^2, x,y>=-a with resolution K`):
print(`Try:`):
print(`MaxCi(x+y,x,y,1,2,10);`):


elif nops([args])=1 and op(1,[args])=MaxF then
print(`MaxF(F,x,y,a,b,L,T): The maximum location and value`):
print(`of F by T^2 hill-climbings at equally-spaced`):
print(`starting points up to resolution 1/10^L`):
print(`For example, try:`):
print(`MaxF(1/(1+x^2+y^2),x,y,1,2,5,2);`):

elif nops([args])=1 and op(1,[args])=MaxF1 then
print(`MaxF1(F,x,y,a,b,pt,L): The maximum location and value`):
print(`of F by hill-climbing starting at pt up to resolution 1/10^L .`):
print(`For example, try:`);
print(`MaxF1(1/(1+x^2+y^2),x,y,1,2,[3/2,3/2],10);`):

elif nops([args])=1 and op(1,[args])=MaxF1a then
print(`MaxF1a(F,x,y,a,eps,pt): The maximum location and value`):
print(`of F by hill-climbing starting at pt`):
print(`For example, try:`):
print(`MaxF1a(1/(1+x^2+y^2),x,y,1,2,.1,[1.5,1.5]);`):

elif nops([args])=1 and op(1,[args])=MaxF11 then
print(`MaxF11(F,x,y,a,b,eps,pt): The best neighbor`):
print(`at the point pt1 with resolution eps in [a,b]^2 .`):
print(`For example, try:`):
print(`MaxF11(1/(1+x^2+y^2),x,y,1,2,.1,[1.5,1.5]);`):

elif nops([args])=1 and op(1,[args])=MaxFline then
print(`MaxFline(F,x,y,R,a,eps): The maximum of the function F on the line`):
print(`x+y=R, for x>=-a, y>=-a with resolution epsilon.`):
print(`it reuturns the champion and the value`):
print(`For example, try:`):
print(`MaxFline(x+y,x,y,10,1,.1);`):


elif nops([args])=1 and op(1,[args])=MaxFnaive then
print(`MaxFnaive(F,x,y,a,M): The location and value of the `):
print(`maximum of the function F`);
print(` in the region -a<=x,y<=M take-away {[0,0]}`);
print(`For example, try:`);
print(`MaxFnaive(x*y,x,y,-10,10,1.);`):

elif nops([args])=1 and op(1,[args])=MaxYakhas then
print(`MaxYakhas(R,x,y,Rad1,k,K): The maximum of Yakhas-k`):
print(`on x^2+y^2=R^2, x,y>=-a with resolution K`):
print(`Try:`):
print(`MaxYakhas(Y2k(1,2,x,y),x,y,.4,2,10);`):

elif nops([args])=1 and op(1,[args])=MinPol then
print(`MinPol(P,c,d): Given a polynomial P in the variables c and d`):
print(`finds its minimum in the first quadrant c>=0,d>=0.`):
print(`For example, try:`):
print(` MinPol(c^2+d^2+1,c,d); `):

elif nops([args])=1 and op(1,[args])=MyMin2 then
print(`MyMin2(P,x,y): Given a polynomial P of x, y`):
print(`finds its minimum (in floating) in the positive quadrant`):
print(`For example, try:`):
print(`MyMin2(x^2+y^2+1,x,y);`):

elif nops([args])=1 and op(1,[args])=MySolve2 then
print(`MySolve2(P,Q,c,d): inputs two polynomials in (c,d)`):
print(`and outputs all the soltions in floats using`):
print(`the Sylvester matrix. For example, try:`):
print(`MySolve2(c^2-4*c*d+d^2,c^2+3*c*d+d^2,c,d);`):

elif nops([args])=1 and op(1,[args])=ProveK then
print(`ProveK(R,x,y,k,cu): Given a rational function R(x,y)`):
print(`describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to`):
print(`itself, finds the fixed point, and whether it is`):
print(` locally asympototically `):
print(`stable, and then proves that for the given k,R^k is `):
print(` a shrinking-norm map, with shrinking factor<cu `):
print(`For example, try`):
print(`ProveK((1+y)/(1+x),x,y,3,9/10);`):

elif nops([args])=1 and op(1,[args])=ProveKar then
print(`ProveKar(R,x,y,Maxk,cu): Given a rational function R(x,y)`):
print(`describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to`):
print(`itself, finds the fixed point, and whether it is`):
print(` locally asympototically `):
print(`stable, and then finds a k<=Maxk, such that R^k is `):
print(` a shrinking-norm map, with shrinking factor<cu `):
print(`Warning: this is almost rigorous.`):
print(`For example, try`):
print(`ProveKar((1+y)/(1+x),x,y,10,9/10);`):

elif nops([args])=1 and op(1,[args])=ProveKnr then
print(`ProveKnr(R,x,y,Maxk,cu,K,K1): Given a rational function R(x,y)`):
print(`describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to`):
print(`itself, finds the fixed point, and whether it is`):
print(` locally asympototically `):
print(`stable, and then finds a k<=Maxk, such that R^k is `):
print(` a shrinking-norm map, with shrinking factor<cu `):
print(`K and K1 are pos. integer, K larger than 10, K1 larger than 100`):
print(` (recommended). It looks at [0,K]^2 and tries K1 random values `):
print(`Warning: this is non-rigorous, but much faster`):
print(`For example, try`):
print(`ProveKnr((1+y)/(1+x),x,y,10,9/10,10,100);`):

elif nops([args])=1 and op(1,[args])=Traj then
print(`Traj(R,x,y,pt,N): the trajectory of the mapping (x,y)->(y,R(x,y))`):
print(`N times, starting with pt`):
print(`For example, try:`):
print(`Traj((24+3*y)/(1+x),x,y,[1,2],20);`):

elif nops([args])=1 and op(1,[args])=Yakhas then
print(`Yakhas(R,x,y,pt,a,k): the ratio that should be <1 `):
print(`w.r.t. to the trans. R`):
print(`and the fixed point a? For example, try:`):
print(` Yakhas((24+3*y)/(1+x),x,y,[1,2],6,5);`):

elif nops([args])=1 and op(1,[args])=YakhasS then
print(`YakhasS(R,x,y,pt,k): the ratio that should be <1 w.r.t. to the trans. R`):
print(`with symbolic parameters`):
print(`For example, try:`):
print(`YakhasS(Y2k(a,b,x,y),x,y,[c,d],5,(a+b-1)/2);`):

elif nops([args])=1 and op(1,[args])=Y1k then
print(`Y1k(p,q,x,y): (6.54) in K.-Ladas' book`):
print(` p and q are symbolic and p<q. For example, try: `):
print(`Try:Y1k(3,2,x,y);`):

elif nops([args])=1 and op(1,[args])=Y2k then
print(`Y2k(M,q,x,y): the Y2k transformation with`):
print(`parameters that guarantee a rational fixed-point`):
print(`M must be bigger than q `):
print(`Try:Y2k(5,3,x,y);`):

elif nops([args])=1 and op(1,[args])=Y3k then
print(`Y3k(c1,c2,d0,d1,p,x,y): the Y3k transformation with`):
print(`R(x,y)=(1+c1*x+c2*y)/(d0+d1*x+d2*y)`):
print(`with d2=(p^2-(c1+c2-d0)^2)/4-d1`):
print(`Try:Y3k(c1,c2,d0,d1,p,x,y);`):

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

#===========================================================================
#FixPoint(R,x,y): Given a rational recurrence x[n+1]=R(x[n],x[n-1])
#in terms of a rational function R(x,y) finds its fixed points
#For example, try:
#FixPoint((24+3*y)/(1+x),x,y);
FixPoint:=proc(R,x,y) local a:
a:=[solve(numer(x-subs({y=x},R)),{x})]:
[seq(subs(a[i],x),i=1..nops(a))]:
end:

#===========================================================================
#LSFP(R,x,y): Given a rational recurrence x[n+1]=R(x[n],x[n-1])
#in terms of a rational function R(x,y) finds its 
#locally stable fixed points
#For example, try:
#LSFP((24+3*y)/(1+x),x,y);
LSFP:=proc(R,x,y) local gu,p,q,g,lam,mu,lu,lu1:
option remember:
gu:=FixPoint(R,x,y):
mu:={}:

for g in gu do
 if evalf(g)>=0 then
  q:=subs({x=g,y=g},diff(R,x)):
  p:=subs({x=g,y=g},diff(R,y)):
  lu:=solve({lam^2-p*lam-q},{lam}):
 lu:={ seq(subs(lu1,lam),lu1 in lu)}:


   if {seq(evalb(evalf(abs(lu1))<1),lu1 in lu)}={true} then
    mu:=mu union {g}:
   fi:
 fi:
od:
mu:
   
end:

#===========================================================================
#Traj(R,x,y,pt,N): the trajectory of the mapping (x,y)->(y,R(x,y))
#N times, starting with pt
#For example, try:
#Traj((24+3*y)/(1+x),x,y,[1,2],20);

Traj:=proc(R,x,y,pt,N) local i,gu:
gu:=pt:
for i from 1 to N do
gu:=[op(gu),subs({x=gu[nops(gu)-1],y=gu[nops(gu)]},R)]:
od:
gu:

end:

#===========================================================================
#CGA(R,x,y,K,L,N): Given a rational recurrence x[n+1]=R(x[n],x[n-1])
#in terms of a rational function R(x,y) conjectures
#a global attractor, if it exists.
#K and L are positive integers. We pick L randomly chosen
#initial points in [0,K]^2 and look at the trajectory N
#times and see if they all seem to converge to the global
#attractor
#For example, try:
#CGA((24+3*y)/(1+x),x,y,10,20,100);
CGA:=proc(R,x,y,K,L,N) local gu,a,ra,i,traj1,pt,ku:

gu:=LSFP(R,x,y):

if nops(gu)=0 then
 print(`There are not even locally asympt. stable points`):
 RETURN(FAIL):
elif nops(gu)>1 then
 print(`There are more than one locally asympt. stable points`):
 RETURN(FAIL):
else
a:=evalf(gu[1]):
fi:

print(`a is`, a):

ra:=rand(0..K*1000):

for i from 1 to L do
 pt:=[ra()/1000., ra()/1000.]:
 traj1:=evalf(Traj(R,x,y,pt,N),10):
  #if abs(evalf(traj1[N])-evalf(a))>10^(-Digits+3) or abs(traj1[N-1]-a)>10^(-Digits+3) then
  ku:=abs(evalf(traj1[N])-evalf(a)):
 if ku>10^(-10)  then
  print(`ku is`, ku):
    RETURN(FAIL):
  fi:
od:

a:
end:

#===========================================================================
#Yakhas(R,x,y,pt,a,k): the ratio that should be <1 w.r.t. to the trans. R
#and the fixed point a?
#For example, try:
#Yakhas((24+3*y)/(1+x),x,y,[1,2],6,5);
Yakhas:=proc(R,x,y,pt,a,k) local t1,bu,i1:
option remember:
t1:=normal(Traj(R,x,y,pt,k)):
t1:=[seq(t1[i1]-a,i1=1..nops(t1))]:
bu:=(t1[k+1]^2+t1[k+2]^2)/(t1[1]^2+t1[2]^2):
normal(bu):

end:

#===========================================================================
#EfesT(R,x,y): Given a function R(x,y) first finds out
#whether it is a locally-stable fixed point and then
#it returns the transformed transformation with 0 as fixed point
#It also returns the fixed point
#For example, try:
#EfesT((24+3*y)/(1+x),x,y);
EfesT:=proc(R,x,y) local lu,a:
option remember:
lu:=LSFP(R,x,y):
if nops(lu)<>1 then
 RETURN(FAIL):
fi:
a:=lu[1]:
normal(subs({x=x+a,y=y+a},R)-a),a:
end:

#===========================================================================
#Y2kNew(M,q,x,y): the Y2k transformation with
#parameters that guarantee a rational fixed-point
#with 2*M>2*q+1
#Try:Y2k(5,3,x,y);
Y2kNew:=proc(M,q,x,y) local p,M1,q1:
M1:=2*M:
q1:=2*q+1:
p:=(M1-q1+1)*(M1+q1-1)/4:
(p+q1*y)/(1+x):
end:

#===========================================================================
#Y2k(M,q,x,y): the Y2k transformation with
#parameters that guarantee a rational fixed point
#Try:Y2k(5,3,x,y);
Y2k:=proc(M,q,x,y) local p:
p:=(M-q+1)*(M+q-1)/4:
(p+q*y)/(1+x):
end:

#===========================================================================
#Y1kOld(p,q,x,y): (6.54) in K./Ladas
#Try:Y1k(5,3,x,y);
Y1k:=proc(p,q,x,y):
(x+p*y)/(q+x):
end:

#===========================================================================
#Y1kN(p,q,x,y): (6.54) in K./Ladas
#Try:Y1kNew(5,3,x,y);
Y1kN:=proc(a,b,x,y) local p,q:
p:=a:
q:=a+b:
Y1k(p,q,x,y):
end:

#===========================================================================
#ConjK(R,x,y,MaxK): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]) finds a point that seems to be
#a global attractor and conjectures a k for it
#For example, try:
#ConjK((24+3*y)/(1+x),x,y,20);
ConjK:=proc(R,x,y,MaxK) local R1,gu,a,i,j,k,k1:
gu:=EfesT(R,x,y):

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

a:=gu[2]:
R1:=gu[1]:

for k from 1 to MaxK while 
max(seq(seq(Yakhas(R1,x,y,[-a+.01+i,-a+.01+j],0,k),i=1..5),j=1..5))>1 do od:

if k=MaxK+1 then
 RETURN(FAIL):
fi:

k1:=k:
gu:=
max(seq(seq(Yakhas(R1,x,y,[-a+.01+i,-a+.01+j],0,k1),i=1..30),j=1..30)):

if gu>1 then
for k from k1+1 to MaxK while 
max(seq(seq(Yakhas(R1,x,y,[-a+.01+i,-a+.01+j],0,k),i=1..30),j=1..30))>= 1 do od:
fi:

if k=MaxK+1 then
 RETURN(FAIL):
fi:

gu:=
max(seq(seq(Yakhas(R1,x,y,[-a+.01+i/10.,-a+.01+j/10.],0,k),i=1..40),j=1..40 )):

if gu>1 then
 RETUNR(FAIL):
fi:

k:

end:

#===========================================================================
#KptS(R,x,y,pt,Maxk): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]), and a point pt,
#finds all the k from 1 to Maxk such that iterating it k times
#followed by the ratio
#at pt shrinks the norm. For example, try:
#KptS((24+3*y)/(1+x),x,y,[1,2],20);
KptS:=proc(R,x,y,pt,MaxK) local R1,gu,k,mu:
gu:=EfesT(R,x,y):


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


R1:=gu[1]:

mu:={}:

for k from 1 to MaxK do

gu:=Yakhas(R1,x,y,pt,0,k):

if gu<1 then
 mu:=mu union {[k,gu]}:
fi:

od:

mu:

end:

#===========================================================================
#Kpt(R,x,y,pt,Maxk,cu): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]), and a point pt,
#finds the smallest k <=Maxk such that iterating it k times
#followed by the ratio
#at pt shrinks the norm by at least cu. For example, try:
#Kpt((24+3*y)/(1+x),x,y,[1,2],20);
Kpt:=proc(R,x,y,pt,MaxK,cu) local R1,gu,k,mu:

gu:=EfesT(R,x,y):


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


R1:=gu[1]:

mu:={}:

for k from 1 to MaxK do

gu:=Yakhas(R1,x,y,pt,0,k):

if gu<cu then
 RETURN([k,gu]):
fi:

od:

mu:

end:

#===========================================================================
#MaxCi(F,x,y,a,R,K): The maximum location and value of the
#function F(x,y) on the part of the circle
#x^2+y^2=R^2, x,y>=-a with resolution K
#Try:
#MaxCi(x+y,x,y,-1,2,10);
MaxCi:=proc(F,x,y,a,R,K) local theta,t1,t2,orekh,i:
if a/R<1 then
theta:=arcsin(evalf(a/R)):
t1:=-theta:
t2:=evalf(Pi/2+theta):
else
t1:=0: t2:=evalf(2*Pi):
fi:

orekh:=t2-t1:
max(seq(subs({x=R*cos(i/K*orekh+t1),y=R*sin(i/K*orekh+t1)},F),i=0..K)):
end:

#===========================================================================
#MaxYakhas(R,x,y,Rad1,K): The maximum of Yakhas
#on x^2+y^2=R^2, x,y>=-a with resolution K
#Try:
#MaxYakhas(Y2k(2,3,x,y),x,y,.4,2,10);
MaxYakhas:=proc(R,x,y,Rad1,k,K) local gu,a,R1,lu,lu1,pt,i:
gu:=EfesT(R,x,y):
a:=gu[2]:
R1:=gu[1]:

lu1:=
{seq([Rad1*cos(i*evalf(2*Pi)/K),Rad1*sin(i*evalf(2*Pi)/K)],
i=0..K-1)}:

lu:={}:
for pt in lu1 do
 if pt[1]>-a and  pt[2]>-a then
  lu:=lu union {pt}:
 fi:
od:

max(seq(Yakhas(R1,x,y,pt,0,k), pt in lu)):
end:

#===========================================================================
#MaxFline(F,x,y,R,a,eps): The maximum of the function F on the line
#x+y=R, for x>=-a, y>=-a with resolution epsilon
#For example, try:
#MaxFline(x+y,x,y,10,1,.1);
#it reuturns the champion and the value
MaxFline:=proc(F,x,y,R,a,eps) local K,i,si,champ,muam,halev:
K:=(R+2*a)/eps:
champ:=[-a,R+a]:
si:=subs({x=champ[1],y=champ[2]},F):

for i from 1 to K do
 muam:=[-a+i*eps,R+a-i*eps]:
 halev:=subs({x=muam[1],y=muam[2]},F):
  
if halev>si then
   si:=halev:
   champ:=muam:
fi:
od:
champ,si:

end:

#===========================================================================
#MaxFnaive(F,x,y,a,M): The location and value of the 
#maximum of the function F
# in the region -a<=x,y<=M take-away {[0,0]}
#For example, try:
#MaxFnaive(x*y,x,y,-10,10,1.);
MaxFnaive:=proc(F,x,y,a,M,eps) local i,j,K,champ,si,muam,halev:
K:=trunc((M+a)/eps):
champ:=[-a,-a]:
si:=subs({x=champ[1],y=champ[2]},F):

for i from 0 to K do
for j from 0 to K do
 muam:=[-a+i*eps,-a+j*eps]:
 halev:=subs({x=muam[1],y=muam[2]},F):
 if halev>si then
   champ:=muam:
   si:=halev:
 fi:
od:
od:
champ,si:
end:

#===========================================================================
#IsKgood(R,x,y,r1,cu,k,K): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]), and a real number r1
#and a k, that seems to be a good one for all x+y=r1
#using K points between [-a,r1+a] and [r1+a,-a]
#decides whether this k seems to be a good one for the line x+y=r1
#For example, try:
#IsKgood((24+3*y)/(1+x),x,y,3,.9,2,10);
IsKgood:=proc(R,x,y,r1,cu,k,K) 
local gu,a,eps,i1,R1:
gu:=EfesT(R,x,y):
if gu=FAIL then
  RETURN(FAIL):
fi:
a:=gu[2]:
R1:=gu[1]:

eps:=(r1+2*a)/K:

gu:=max(seq(Yakhas(R1,x,y,[r1+a-eps*i1,-a+eps*i1],0,k),i1=0..K)):

if gu<cu then
true:
else
false:
fi:


end:

#===========================================================================
#FindKline(R,x,y,r1,cu,Maxk,K): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]), and a real number r1
#using K points between [-a,r1+a] and [r1+a,-a]
#finds the smallest k seems to be a good one for the line x+y=r1
#For example, try:
#FindKline((24+3*y)/(1+x),x,y,3,.9,20,10);
FindKline:=proc(R,x,y,r1,cu,Maxk,K) 
local k:

for k from 1 to Maxk do
if IsKgood(R,x,y,r1,cu,k,K) then
 RETURN(k):
fi:
od:
FAIL:
end:

#===========================================================================
#FindKold(R,x,y,cu,Maxk,K,eps1,M): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]), and a real number r1
#using K points between [-a,r1+a] and [r1+a,-a]
#retuns a table for the best k for each line
#x+y=r1 starting at r1=-2*a with intervals of eps1
#FindKold((24+3*y)/(1+x),x,y,3,.9,20,10,1,20);
FindKold:=proc(R,x,y,cu,Maxk,K,eps1,M)  local r1,gu,i1,gu1,a,k,R1:

a:=EfesT(R,x,y)[2]:
R1:=EfesT(R,x,y)[1]:
gu:=[]:

for i1 from 0 to trunc((M+2*a)/eps1) do
  r1:=-2*a+eps1*i1:
 gu1:=FindKline(R,x,y,r1+0.01,cu,Maxk,K,eps1,M):
  if gu1=FAIL then
   RETURN(FAIL,gu,[R1,r1+0.01]):
  else
   gu:=[op(gu),[r1,gu1]]:
  fi:

od:

k:=max(seq(gu[i1][2],i1=1..nops(gu))):

if {seq(IsKgood(R,x,y,gu[i1][1]+0.01,cu,k,K),i1=1..nops(gu))}={true} then
  RETURN(k):
else
 RETURN(gu):
fi:
  

end:


###begin of stuff from MaxF###########

#===========================================================================
#MaxF11(F,x,y,a,b,eps,pt): The best neighbor
#at the point pt1 with resolution eps in [a,b]^2
#For example, try:
#MaxF11(1/(1+x^2+y^2),x,y,1,2,.1,[1.5,1.5]);
MaxF11:=proc(F,x,y,a,b,eps,pt) local Neis,Neis1, cha,rec,n,muam:

Neis:={[pt[1]+eps,pt[2]],[pt[1]-eps,pt[2]],
[pt[1],pt[2]+eps],[pt[1],pt[2]-eps]} minus {[0,0]}:

Neis1:={}:

for n in Neis do
 if n[1]>=a and n[1]<=b and n[2]>=a and n[2]<=b then
  Neis1:=Neis1 union {n}:
 fi:
od:

Neis:=Neis1:

cha:=pt:

rec:=subs({x=pt[1],y=pt[2]},F):

for n in Neis do
 muam:=subs({x=n[1],y=n[2]},F):
 if muam>rec then
   cha:=n:
   rec:=muam:
 fi:
od:

if rec<0 then
 ERROR(`Something is wrong due to floating-point, make DIGITS bigger`):
fi:
cha,rec:
end:

#===========================================================================
#MaxF1a(F,x,y,a,b,eps,pt): The maximum location and value
#of F by hill-climbing starting at pt
#For example, try:
#MaxF1a(1/(1+x^2+y^2),x,y,1,2,.1,[1.5,1.5]);
MaxF1a:=proc(F,x,y,a,b,eps,pt) local pt1,pt2,pt3:

pt1:=pt:

pt2:=MaxF11(F,x,y,a,b,eps,pt1)[1]:


while pt1<>pt2 do
 pt3:=MaxF11(F,x,y,a,b,eps,pt2)[1]:
pt1:=pt2:
pt2:=pt3:
od:

pt2,subs({x=pt2[1],y=pt2[2]},F):

end:

#===========================================================================
#MaxF1(F,x,y,a,b,pt,L): The maximum location and value
#of F by hill-climbing starting at pt up to resolution 1/10^L
#For example, try:
#MaxF1(1/(1+x^2+y^2),x,y,1,2,[3/2,3/2],10);
MaxF1:=proc(F,x,y,a,b,pt,L) local i,pt1:

pt1:=pt:
for i from 1 to L do
 pt1:=MaxF1a(F,x,y,a,b,1/10^i,pt1)[1]:
od:

pt1,subs({x=pt1[1],y=pt1[2]},F):

end:

#===========================================================================
#MaxF(F,x,y,a,b,L,T): The maximum location and value
#of F by T^2 hill-climbings at equally-spaced
#starting points up to resolution 1/10^L
#For example, try:
#MaxF(1/(1+x^2+y^2),x,y,1,2,5,2);
MaxF:=proc(F,x,y,a,b,L,T) local i,j,gu:

gu:=
{seq(seq([MaxF1(F,x,y,a,b,[a+(b-a)/T*i+1/11/T,a+(b-a)/T*j+1/11/T],L)],i=1..T),j=1..T)}:

max(seq(gu[i][2],i=1..nops(gu))):

end:

###end of stuff from MaxF###########

#===========================================================================
#ConjKnew(R,x,y,Maxk,L,T): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]) finds a point that seems to be
#a global attractor and conjectures a k for it
#1/10^L is the accuracy and T is the parameter for the number of
#initial points in the hill-climbing algorithm
#For example, try:
#ConjKnew((1+y)/(1+x),x,y,20,2,2);
ConjKnew:=proc(R,x,y,MaxK,L,T) local R1,a,gu,k,F,c,d,mu1,mu2,mu3
:
gu:=EfesT(R,x,y):
R1:=gu[1]:
a:=evalf(gu[2]):

for k from 1 to MaxK do

print(`k is `, k):
 F:=Yakhas(R1,x,y,[c,d],0,k):

 mu1:=MaxF(F,c,d,-a,-a+1.,2,2):



 if mu1>=1   then
    next:
 fi:

 mu2:=MaxF(F,c,d,-a,-a+1.,L,T):



 if mu2>=1   then
    next:
 fi:

 mu3:=MaxF(F,c,d,-a,-a+100.,L,T):



 if mu3<1   then
 RETURN(k,max(mu1,mu3,mu3)):
  fi:
od:
FAIL:
end:


#===========================================================================
#ConjKfast(R,x,y,Maxk): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]) finds a point that seems to be
#a global attractor and conjectures a k for it
#ConjKfast((1+y)/(1+x),x,y,20);
ConjKfast:=proc(R,x,y,MaxK) local a,gu,i,j:
:
gu:=EfesT(R,x,y):
a:=evalf(gu[2]):

max(
seq(seq(Kpt(R,x,y,[-a+0.01+i/10,-a+0.01+j/10],MaxK,1)[1],i=1..10),j=1..10),
seq(seq(Kpt(R,x,y,[-a+0.01+i/2,-a+0.01+j/2],MaxK,1)[1],i=2..10),j=2..10),
seq(seq(Kpt(R,x,y,[-a+0.01+i,-a+0.01+j],MaxK,1)[1],i=2..40),j=2..40)):


end:

#===========================================================================
#MinPol(P,c,d): Given a polynomial P in the variables c and d
#finds its minimum in the first quadrant c>=0,d>=0.
#For example, try:
#MinPol(c^2+d^2+1,c,d);
MinPol:=proc(P,c,d) local i, P1,ku,mu1,mu2,mu1a,mu2a,mu3:
P1:=subs(c=0,P):
mu1:=minimize(P1,d=0..infinity):
P1:=subs(d=0,P):
mu2:=minimize(P1,c=0..infinity):

P1:=subs(c=10000,P):
mu1a:=minimize(P1,d=0..infinity):
P1:=subs(d=10000,P):
mu2a:=minimize(P1,c=0..infinity):

mu3:=min(mu1,mu2,mu1a,mu2a):

ku:=[solve({diff(P,d),diff(P,c)}, {c,d})]:
for i from 1 to nops(ku) do
 if evalf(subs(ku[i],c))>=0 and evalf(subs(ku[i],d))>=0 then
  mu3:=min(mu3, evalf(subs(ku[i],P))):
 fi:
od:
mu3:


end:

#===========================================================================
#ProveKOld(R,x,y,Maxk,cu): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then finds a k<=Maxk that is shrinking exponent in the
#sense that every k iterations shrings the norm by <=cu.
#For example, try
#ProveKOld((1+y)/(1+x),x,y,10,9/10);
ProveKOld:=proc(R,x,y,Maxk,cu) local R1,a,c,d,Y,mu,k,i,j,lu:
R1:=EfesT(R,x,y);
if R1=FAIL then
 print(`Not locally asymp. stable`):
 RETURN(FAIL):
fi:
a:=R1[2]:


for k from 1 to Maxk do 
print(`trying k=`, k):
Y:=Yakhas(R,x,y,[c,d],a,k):
mu:=numer(cu-Y):

lu:=min(seq(seq(subs({c=0.11*i,d=0.11*j},mu),i=0..trunc(2*a/0.11)),
j=0..trunc(2*a/0.11))):

if lu<0 then
 next:
fi:

lu:=min(seq(seq(subs({c=1.11*i,d=1.11*j},mu),i=0..10),j=0..10)):

if lu<0 then
 next:
fi:

lu:=min(seq(seq(subs({c=11.1*i,d=11.1*j},mu),i=0..10),j=0..10)):

if lu<0 then
 next:
fi:

print(`There is hope for k=`, k, `trying to prove it rigorously`):
if evalf(MinPol(mu,c,d))>=0 then
print(`Yes!, k=`, k, `works `):
 RETURN(k):
fi:

od:

FAIL:

end:

#===========================================================================
#Kpt1(R,x,y,pt,Maxk,cu): given a rational transformation R(x,y)
#describing a second-order recurrence
#x[n]=R(x[n-1],x[n]), and a point pt,
#finds the smallest k <=Maxk such that iterating it k times
#followed by the ratio
#at pt shrinks the norm by at least cu. For example, try:
#Kpt1((24+3*y)/(1+x),x,y,[1,2],20);
Kpt1:=proc(R,x,y,pt,MaxK,cu) local R1,gu,k,mu,a:

a:=EfesT(R,x,y)[2]:


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


mu:={}:

for k from 1 to MaxK do

gu:=Yakhas(R,x,y,pt,a,k):

if gu<cu then
 RETURN([k,gu]):
fi:

od:

mu:

end:

#===========================================================================
#LSFPs(R,x,y,para): Given a rational recurrence x[n+1]=R(x[n],x[n-1])
#in terms of a rational function R(x,y) finds its 
#and a parameter set para, finds the fixed points
#and for each the conditions that (i) it is positive
#(ii) that it is locally stable
#For example, try:
#LSFPs(Y2k(M,q,x,y),x,y,{M,q});
LSFPs:=proc(R,x,y,para) local gu,p,q,g,lam,mu,lu,lu1:
option remember:
gu:=FixPoint(R,x,y):

mu:={}:

for g in gu do
  q:=subs({x=g,y=g},diff(R,x)):
  p:=subs({x=g,y=g},diff(R,y)):
  lu:=solve({lam^2-p*lam-q},{lam}):
 lu:={ seq(subs(lu1,lam),lu1 in lu)}:

    mu:=mu union {[g,lu]}:
od:
mu:
   
end:

#===========================================================================
#Y3k(c1,c2,d0,d1,p,x,y): the Y3k transformation with
#R(x,y)=(1+c1*x+c2*y)/(d0+d1*x+d2*y)
#with d2=(p^2-(c1+c2-d0)^2)/4-d1
#Try:Y3k(c1,c2,d0,d1,p,x,y);
Y3k:=proc(c1,c2,d0,d1,p,x,y) local d2:
d2:=(p^2-(c1+c2-d0)^2)/4-d1:
(1+c1*x+c2*y)/(d0+d1*x+d2*y):
end:

#===========================================================================
#YakhasS(R,x,y,pt,k,fp): the ratio that should be <1 w.r.t. to the trans. R
#with symbolic parameters and the indicated fixed-point
#For example, try:
#YakhasS(Y2k(a,b,x,y),x,y,[c,d],5,(a+b-1)/2);
YakhasS:=proc(R,x,y,pt,k,a) local t1,bu,i1:
option remember:

if normal(subs({x=a,y=a},R)-a)<>0 then
 print(a, ` is not a fixed point`):
 RETURN(FAIL):
fi:
t1:=normal(Traj(R,x,y,pt,k)):
t1:=[seq(t1[i1]-a,i1=1..nops(t1))]:
bu:=(t1[k+1]^2+t1[k+2]^2)/(t1[1]^2+t1[2]^2):
normal(bu):

end:

#===========================================================================
#ProveKnr(R,x,y,Maxk,cu,K,K1): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then finds a k<=Maxk that is shrinking exponent in the
#sense that every k iterations shrings the norm by <=cu.
#and proves it non-rigorously. K and K1 are positive integers
#should be at least 10, the larger the better
#(it looks at the region [0,K]^2 and picks K1 random points)
#It does it semi-rigorously!
#For example, try
#ProveKnr((1+y)/(1+x),x,y,10,9/10,10,100);
ProveKnr:=proc(R,x,y,Maxk,cu,K,K1) local R1,a,c,d,Y,mu,k,i,j,ka:
R1:=EfesT(R,x,y);
if R1=FAIL then
 print(`Not locally asymp. stable`):
 RETURN(FAIL):
fi:
a:=R1[2]:


for k from 1 to Maxk do 

print(`trying k=`, k):
Y:=Yakhas(R,x,y,[c,d],a,k):

ka:=max(seq(seq(subs({c=i+1/11,d=j+1/11},Y),i=0..5),j=0..5)):

if ka>=cu then
print(`The max ratio is>= `, evalf(ka)):
next:
fi:

ka:=max(seq(seq(subs({c=i+1/11,d=j+1/11},Y),i=0..10),j=0..10)):
if ka>=cu then
print(`The max ratio is>= `, evalf(ka)):
next:
else
print(`There is hope, the max. ratio is probably close to`,evalf(ka)):
fi:

mu:=numer(cu-Y):


if PPnr(mu,c,d,K,a,K1) then
 RETURN(k):
else
print(`k=`, k, `didn't work, we will keep going`):
 next:
fi:

od:

FAIL:

end:

#===========================================================================
#PPnr(f,x,y,K,a,K1): a non-rigorous proof that f(x,y) is positive
#in [0,infinity]^2 by using K as "far"
#a is the attractor such that f(a,a)=0
#K1 is the number of random points examined
PPnr:=proc(f,x,y,K,a,K1) local i,j,ra,pt,m:

if lcoeff(subs(x=0,f),y)<0 then
 RETURN(false);
fi:

if lcoeff(subs(y=0,f),x)<0 then
 RETURN(false);
fi:

if lcoeff(lcoeff(f,y),x)<0 then
 RETURN(false);
fi:


if lcoeff(lcoeff(f,x),y)<0 then
 RETURN(false);
fi:

if min(seq(subs({x=i,y=K},f),i=0..K))<K then
RETURN(false):
fi:

if min(seq(subs({y=i,x=K},f),i=0..K))<K then
RETURN(false):
fi:


m:=min(seq(seq(subs({x=i,y=j},f),i=0..K),j=0..K)):

if m<0 then
  RETURN(false):
fi:

if subs({x=a,y=a},diff(f,x))<>0 or subs({x=a,y=a},diff(f,y))<>0 then
RETURN(false):
fi:


ra:=rand(0..K*1000):


for i from 1 to K1 do
 pt:=[ra()/1000.,ra()/1000.]:
if subs({x=pt[1],y=pt[2]},f)<0 then
  RETURN(false):
fi:
od:
true:
end:

#===========================================================================
#ProveKnr(R,x,y,Maxk,cu,K,K1): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then finds a k<=Maxk that is shrinking exponent in the
#sense that every k iterations shrings the norm by <=cu.
#and proves it non-rigorously. K and K1 are positive integers
#should be at least 10, the larger the better
#(it looks at the region [0,K]^2 and picks K1 random points)
#It does it semi-rigorously!
#For example, try
#ProveKnr((1+y)/(1+x),x,y,10,9/10,10,100);
ProveKnrN:=proc(R,x,y,Maxk,cu,K,K1) local R1,a,c,d,Y,mu,k,i,j,lu:
R1:=EfesT(R,x,y);
if R1=FAIL then
 print(`Not locally asymp. stable`):
 RETURN(FAIL):
fi:
a:=R1[2]:


for k from 1 to Maxk do 

print(`trying k=`, k):
Y:=Yakhas(R,x,y,[c,d],a,k):

mu:=numer(cu-Y):

lu:=min(seq(seq(subs({c=0.11*i,d=0.11*j},mu),i=0..trunc(2*a/0.11)),
j=0..trunc(2*a/0.11))):

if lu<0 then
 next:
fi:

if PPnr(mu,c,d,K,a,K1) then
 RETURN(k):
else
print(`k=`, k, `didn't work, we will keep going`):
 next:
fi:

od:

FAIL:

end:

#===========================================================================
#Testk(R,x,y,k,cu,KAMA): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then tests whether k is OK, non-rigoroulsy
#Testk(1+y)/(1+x),x,y,10,9/10,KAMA);
Testk:=proc(R,x,y,k,cu,KAMA) local R1,a,a1,ma,i,j,ra,ma1,pt:
R1:=EfesT(R,x,y);
if R1=FAIL then
 print(`Not locally asymp. stable`):
 RETURN(FAIL):
fi:
a:=R1[2]:

if type(a,integer) then
 a1:=.5:
else
 a1:=0.:
fi:

ma:=max(seq(seq(Yakhas(R,x,y,evalf([a1+i,a1+j]),a,k),i=0..5),j=0..5)):

if ma>cu then
 RETURN(false):
fi:

 ma1:=max(
seq(seq(Yakhas(R,x,y,evalf([a+i/100,a+j/100]),a,k),i=1..10),j=1..10),
seq(seq(Yakhas(R,x,y,evalf([a-i/100,a+j/100]),a,k),i=1..10),j=1..10),
seq(seq(Yakhas(R,x,y,evalf([a+i/100,a-j/100]),a,k),i=1..10),j=1..10),
seq(seq(Yakhas(R,x,y,evalf([a-i/100,a-j/100]),a,k),i=1..10),j=1..10)
):

if ma1>cu then
RETURN(false):
fi:


ra:=rand(1..100000):
for i from 1 to KAMA do
 pt:=[ra()/1000.,ra()/1000.]:
 if Yakhas(R,x,y,pt,a,k)>cu then
    RETURN(false):
 fi:
od:

true:

end:

#===========================================================================
#FindGoodks(R,x,y,Maxk,cu,KAMA): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then finds all k<=Maxk that are shrinking exponent in the
#sense that every k iterations shrings the norm by <=cu.
#Warning: this is just a guess!
#and proves it non-rigorously. 
#FindGoodks((1+y)/(1+x),x,y,10,9/10,1000);
FindGoodks:=proc(R,x,y,Maxk,cu,KAMA) local k,gu:

gu:={}:
for k from 1 to Maxk do
 if Testk(R,x,y,k,cu,KAMA) then
   gu:=gu union {k}:
 fi:
od:

gu:

end:

#===========================================================================
#Findk(R,x,y,Maxk,cu): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then finds the smallest k<=Maxk that are shrinking exponent in the
#sense that every k iterations shrings the norm by <=cu.
#and proves it non-rigorously. 
#Findk((1+y)/(1+x),x,y,10,9/10);
Findk:=proc(R,x,y,Maxk,cu) local k,gu:

gu:={}:
for k from 1 to Maxk do
 if Testk(R,x,y,k,cu) then
 RETURN(k):
 fi:
od:

FAIL:

end:

#===========================================================================
#FindGoodksR(R,x,y,Maxk,cu,a,b,RangeA,RangeB,res) :
#Given a rational recurrence x[n]=R(x[n-1],x[n-2])
#phrased in terms of parameters a and b
#finds all good k's <=Maxk that are valid for all
#(a,b) in the rectangle in Parameter space
#Range A x RangeB with resolution res
#For example, try:
#FindGoodksR(Y2k(M+q,q,x,y),x,y,10,9/10,M,q,[1,3],[1,3],1) ;
FindGoodksR:=proc(R,x,y,Maxk,cu,a,b,RangeA,RangeB,res) 
local gu,a0,a1,b0,b1,A,B,i,j,gu1:
a0:=RangeA[1]: a1:=RangeA[2]:
b0:=RangeB[1]: b1:=RangeB[2]:

gu:=FindGoodks(subs({a=a0,b=b0},R),x,y,Maxk,cu);

if gu={} then
 RETURN({}):
fi:

for i from 0 to trunc((a1-a0)/res)+1 do
 for j from 0 to trunc((b1-b0)/res)+1 do
  A:=a0+i*res:  B:=b0+j*res:
 print(`[A,B}=`,[A,B]):
  gu1:=FindGoodks(subs({a=A,b=B},R),x,y,Maxk,cu);
  print(`gu1 is`, gu1):
  gu:=gu intersect gu1:
   print(`gu is`, gu):
  if gu={} then
   RETURN({}):
  fi:
 od:
od:

gu:
end:




##Start MyMin2
ez:=proc(): 
print(` FastDet1(A), rPxy(d,x,y,K1,K2), MySolve2(P,Q,x,y) `): 
print(` MyMin1(P,x), MyMin2(P,x,y) `):
end:

Digits:=100:
#===========================================================================
#rPxy(d,x,y,K1,K2)
rPxy:=proc(d,x,y,K1,K2) local i,j,ra:
ra:=rand(K1..K2):
add(add(ra()*x^i*y^j,i=0..d),j=0..d):
end:

#===========================================================================
FastDet1:= proc(A)
local var, deg, det, det1, n, i, j, p, pp, vals, a, b, len;
   n := {LinearAlgebra[Dimension](A)};
   if nops(n) <> 1 then return FAIL fi:
#"Matrix is not square." fi;
   n := op(1,n);
   var := indets(A);
   if nops(var) <> 1 then return FAIL: fi:
#"Matrix entries must be univariate polynomials." fi;
   var := op(1,var);
   deg := add(max(seq(degree(A[i,j]),i=1..n)),j=1..n);
   a := ceil(-deg/2); b := ceil(deg/2);
   det := 0; det1 := 1;
   p := 2^24; pp := 1;
   while det <> det1 do
     det1 := det;
     p := prevprime(p);
     vals := [seq(Det(subs(var=i,A)) mod p, i=a..b)];
     det := Interp([seq(i,i=a..b)], vals, var) mod p;
     det := [seq(coeff(det, var, i), i=0..deg)];
     if det1 <> 1 then
       len := max(nops(det), nops(det1));
       det := [op(det), seq(0, i=1..len-nops(det))];
       det1 := [op(det1), seq(0, i=1..len-nops(det1))];
       det := [seq(chrem([det[i],det1[i]], [p,pp]), i=1..len)];
     fi;
     pp := p*pp;
     det := mods(det, pp);
   od;
   return add(det[i+1]*var^i, i=0..deg);
end:

#===========================================================================
#MySolve2(P,Q,x,y): Given two polynomials P and Q in x and y
#finds all solutions of {P(x,y)=0,Q(x,y)=0}
#For example, try:
#MySolve2(x*y+1,2*x*y^2+3,x,y);
MySolve2:=proc(P,Q,x,y) local gux,P1,j,i,lu,mu,vu:
Digits:=100:
gux:=FastDet1(LinearAlgebra[SylvesterMatrix](P,Q,x)):

lu:=[solve(gux,y)]:
print(lu);

mu:={}:

for i from 1 to nops(lu) do
P1:=subs(y=lu[i],P):
vu:=[solve(P1,x)]:

 for j from 1 to nops(vu) do
  if abs(subs({y=evalf(lu[i]),x=evalf(vu[j])},Q))<10^(-10) then
   mu:=mu union {[vu[j],lu[i]]}:
  fi:
 od:
od:

mu:

end:

#===========================================================================
#MySolve2f(P,Q,x,y): Given two polynomials P and Q in x and y
#finds all solutions of {P(x,y)=0,Q(x,y)=0}
#For example, try:
#MySolve2f(x*y+1,2*x*y^2+3,x,y);
MySolve2f:=proc(P,Q,x,y) local gux,P1,j,i,lu,mu,vu:

gux:=FastDet1(LinearAlgebra[SylvesterMatrix](P,Q,x)):

lu:=[fsolve(gux,y)]:
#print(`FSOLVE(GUX,Y)`,lu);

mu:={}:

for i from 1 to nops(lu) do
if lu[i]>=0 then
P1:=subs(y=lu[i],P):

if P1<>0 then
vu:=[fsolve(P1,x)]:
#print(`FSOLVE(p1,X)`,[lu[i],vu]);

 for j from 1 to nops(vu) do
#  if abs(subs({y=lu[i],x=vu[j]},Q))<10^(-10) then
   if vu[j]>=0 then
   mu:=mu union {[vu[j],lu[i]]}:
   fi;
#  fi:
 od:
fi:

fi;
od:
#print(seq([mu[i],[subs({y=mu[i][2],x=mu[i][1]},P),subs({y=mu[i][2],x=mu[i][1]},Q)]],i=1..nops(mu)));
mu:

end:

#===========================================================================
#MyMin1(P,x): Given a polynomial P of x
#finds its minimum (in floating) in the positive real line
#For example, try:
#MyMin1(x^2+1,x,);
MyMin1:=proc(P,x) local P1,lu,i,mi,muam,cha:

if lcoeff(P,x)<0 then
 RETURN(-infinity):
fi:

P1:=diff(P,x):

lu:=[fsolve(P1,x)]:

cha:=0:
mi:=subs(x=0,P):

for i from 1 to nops(lu) do
 if lu[i]>0 then
 muam:=subs(x=lu[i],P):
 if muam<mi then
  cha:=lu[i]:
   mi:=muam:
 fi:
 fi:
od:
mi,cha:
end:

#===========================================================================
#MyMin2(P,x,y): Given a polynomial P of x, y
#finds its minimum (in floating) in the positive quadrant
#For example, try:
#MyMin2(x^2+y^2+1,x,y);
MyMin2:=proc(P,x,y) local lu,mi,muam,i,cha:
if lcoeff(lcoeff(P,x),y)<0 then
 RETURN(-infinity):
fi:
if lcoeff(lcoeff(P,y),x)<0 then
 RETURN(-infinity):
fi:

if lcoeff(subs(x=0,P),y)<0 then
 RETURN(-infinity):
fi:


if lcoeff(subs(y=0,P),x)<0 then
 RETURN(-infinity):
fi:

cha:=[0,0]:
mi:=subs({x=0,y=0},P):

muam:=MyMin1(subs(y=0,P),x):

if muam[1]<mi then
 mi:=muam[1]:
 cha:=[muam[2],0]:
fi:


muam:=MyMin1(subs(x=0,P),y):

if muam[1]<mi then
 mi:=muam[1]:
 cha:=[0,muam[2]]:
fi:


lu:=MySolve2f(diff(P,x),diff(P,y),x,y):


for i from 1 to nops(lu) do
 if lu[i][1]>=0 and lu[i][2]>=0 then
  muam:=subs({x=lu[i][1],y=lu[i][2]},P):
  if muam<mi then
   mi:=muam:
   cha:=lu[i]:
  fi:
 fi:
od:
mi,cha:


end:
##End MyMin2

#===========================================================================
#ProveK(R,x,y,k,cu): Given a rational function R(x,y)
#describing the mapping (x,y)->(y,R(x,y)) from [0,infinity]^2 to
#itself, finds the fixed point, whether it is locally asympototically
#stable, and then finds a k<=Maxk that is shrinking exponent in the
#sense that every k iterations shrings the norm by <=cu.
#For example, try
#ProveK((1+y)/(1+x),x,y,3,9/10);
ProveK:=proc(R,x,y,k,cu) local R1,a,c,d,Y,mu,kap:
option remember:

R1:=EfesT(R,x,y);

if R1=FAIL then
 print(`Not locally asymp. stable`):
 RETURN(false):
fi:

a:=R1[2]:

Y:=Yakhas(R,x,y,[c,d],a,k):

mu:=numer(cu-Y):

kap:=MyMin2(mu,c,d):

if kap=-infinity then
 RETURN(false):
fi:

if kap[1]>=0 then
	RETURN(true);
else
	RETURN(false);
fi;

#if kap[1]<10^(-20) and abs(kap[2][1]-a)<10^(-20) and abs(kap[2][2]-a)<10^(-20)
#   then
# RETURN(true):
#else
#  RETURN(false):
#fi:

end:

#===========================================================================
#Investigate(R,x,y,a,b,Ra,Rb,k): investigates
#whether a recurrence x[n]=R(x[n-1],x[n-2])
#given by a rational function R(x,y) and
#feauturing parameters a and b,
#is k-contracting in Ra times Rb with
#the given fixed-point. For example, try:
#Investigate(Y1kN(a,b,x,y),x,y,a,b,[1,5],[1,5],2,0,9/10);
Investigate:=proc(R,x,y,a,b,Ra,Rb,k,fp,cu) local i,j:
if normal(subs({x=fp,y=fp},R)-fp)<>0 then
 print(fp, `is not a fixed point `):
 RETURN(false):
fi:

if not ProveK(subs({a=Ra[1],b=Rb[1]},R), x,y,k,cu) then
 RETURN(false):
fi:

if not ProveK(subs({a=Ra[1],b=Rb[2]},R),x,y,k,cu) then
 RETURN(false):
fi:

if not ProveK(subs({a=Ra[2],b=Rb[1]},R),x,y,k,cu) then
 RETURN(false):
fi:

if not ProveK(subs({a=Ra[2],b=Rb[2]},R),x,y,k,cu) then
 RETURN(false):
fi:


for i from Ra[1] to Ra[2] do
for j from Rb[1] to Rb[2] do

if not ProveK(subs({a=i,b=j},R),x,y,k,cu) then
 RETURN(false):
fi:

od:
od:

true:

end: