###################################################################### ##LADAS: Save this file as LADAS # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read LADAS: # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Dec. 2008-March 2009 with(linalg): print(`Created: Dec. 23, 2008-March 2009`): print(`Version of March 10, 2009`): print(` This is LADAS `): print(`It is a Maple package that accompanies the article `): print(`Teaching the Computer how to Discover(!) and then Prove(!!)`): print(`(all by Itself(!!!)) Analogs of Collatz's Notorious 3x+1 Conjecture`): print(`by Doron Zeilberger.`): print(`This article is 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/tokhniot/LADAS .`): print(`For a list of the main procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(`For a list of the supporting procedures type ezra1();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(`For a list of the supporting procedures dealing with`): print(`regular expressions, type ezraRE();,`): print(`but there is no specific help yet for these procedures`): print(``): with(combinat): ezraCheck:=proc(): if args=NULL then print(` The checking procedures are: CheckHafelSre, Concretize1 `): else ezra(args): fi: end: ezra1:=proc(): if args=NULL then print(` The supporting procedures are: AllRules, BatuakhRa, BruteSolve `): print(` Canon, CheckRule, CondPair, `): print(` DiscoverAtoms, DiscoverAtomsf, DressUp `): print(`EfoBa, FindMates, FindPreScheme1, FindPreSchemeG,`): print(` FindPreScheme1G, gDiscoverAtomsf,`): print(` GuessST, HafelS, HafelSw , HafelSwn, `): print(` HopeForOrbits, IsCont, IsFL, IsPrim, KavuaScheme `): print(`MakeTrans,MakeTransR, MakomRishon, MaxTraj1, `): print(` MitatSdom1,MitatSdom2, Novea1, Novea1V`): print(`Novea2, Novea2V, ParseWord,ParseWord1, Per1, `): print(` PtorAB, score1, Sipur1, Sipur1V, Squeeze1, `): print(` sSymbTrajs, Starters, StateScheme, sYeladim `): print(`SymbTraj1, TestABij, `): print(` Tsamek, VecToRule, VerifySchemeE, WeedOut, Yeladim, YeladimAB`): print(`YeladimABf `): else ezra(args): fi: end: ezra:=proc(): if args=NULL then print(` The main procedures are: BatuakhRaSymb `): print(` Children, DAP, DAPf,DAPfA, DAPsf `): print(`DAPsurv, DAPsurvf, DAPsurvY, DAPsurvYf, `): print(` DG1, FindPreScheme, GRE, GREplus, `): print(` gDAP, gDAPsurvf, gGRE, gGREplus,`): print(`gSymbTrajsABf, gYeladimABf, Hafel, HafelSre, `): print(` ImproveScheme, MakeWebBook, MakeScheme, MakeSchemeVerbose `): print(` MaxTrajs, MaxTrajsa0,Novea, Orbits,`): print(` Per, PerMod, Pers, PersMod, Persa0, randPers, Sipur, SymbTrajs, `): print(` SymbTrajsAB,SymbTrajsABf,SymbTrajsC. Tovim `): print(` Traj, VerifyScheme, VerifySchemeV, YeladimRE `): elif nops([args])=1 and op(1,[args])=AllRules then print(`AllRules(a,b): all (generalized) Amleh-Grove-Kent-Ladas-style rules, phrased in terms of symbols a and b`): print(`Try: AllRules(a,b); `): elif nops([args])=1 and op(1,[args])=BatuakhRa then print(`BatuakhRa(gadol,b1,c1,a,b,A,B): Given three`): print(`linear expression, gadol, b1,c1, phrased`): print(`phrased in terms of symbols a and b and`): print(`pos. integers A and B, finds out whether`): print(`it disobeys the necessary conditions`): print(`abs(b1)<=abs(gadol) , abs(c1)<=abs(gadol) `): print(`abs(gadol)<=A*abs(b1)+B*abs(c1)`): print(`it requires that the coeff. of a and b gadol`): print(`have the same sign.`): print(`For example, try:`): print(`BatuakhRa(a+b,a/4,b/4,a,b,1,2) ;`): elif nops([args])=1 and op(1,[args])=BatuakhRaSymb then print(`BatuakhRaSymb(gadol,b1,c1,a,b,A,B,m,r,K1,K2): Given `): print(`a linear expression, gadol, in terms of symbols a and b`): print(` and linear expressions b1,c1, phrased in terms of a and b`): print(`and symbolic integers m[1],m[2], ..., m[r] checks whether`): print(`for the specializations m[1], ... m[r] between K1 and K2`): print(`pos. integers A and B, finds out whether`): print(`it disobeys the necessary conditions`): print(`abs(b1)<=abs(gadol) , abs(c1)<=abs(gadol) `): print(`abs(gadol)<=A*abs(b1)+B*abs(c1)`): print(`it requires that the coeff. of a and b of gadol`): print(`have the same sign.`): print(`For example, try:`): print(`BatuakhRaSymb(a+b,a/2^m[1],b/2^m[1],a,b,1,2,m,1,2,3) ;`): elif nops([args])=1 and op(1,[args])=BatuakhRaOld then print(`BatuakhRaOld(L,a,b,A,B): Given a trajectory L`): print(`phrased in terms of symbols a and b and`): print(`pos. integers A and B, finds out whether`): print(`it disobeys the necessary conditions`): print(`abs(L[i])<=abs(L[3]) and `): print(`abs(L[3])<=A*abs(L[nops(L)-1])+B*abs(L[nops(L)])`): print(`For example, try:`): print(`BatuakhRaOld([a,b,a+b,a/4,b/4],a,b,1,2) ;`): elif nops([args])=1 and op(1,[args])=BruteSolve then print(`BruteSolve(L,L1,a,b,m,K): by brute force finds a`): print(`realizations of the symbolic list L phrased in terms`): print(`of a and b such that plugging-in specific integers`): print(`a0 b0 for and b will give a sequence of integers`): print(`corresponding such that L==L1 mod m, where a and`): print(`b range from -K to K. For example, try: `): print(`BruteSolve([a,b,a/2,a/2-b/4],[2,2,1,1],a,b,2,10);`): elif nops([args])=1 and op(1,[args])=Canon then print(`Canon(L): the canonical form of a period L`): print(`For example, try: Canon([3,6,3,6]);`): elif nops([args])=1 and op(1,[args])=CheckHafelSre then print(`CheckHafelSre(F,a,b,R,k,r): `): print(`Checks HafelSre(F,a,b,R); with r random values`): print(`with multipliticies between 1 and k For example, try:`): print(`CheckHafelSre(E2,a,b,[[1,2,1],[[1,2],n]],10,4);`): elif nops([args])=1 and op(1,[args])=CheckRule then print(`CheckRule(F,a,b): checks whether the rule F is a good one`): print(`For example, try: CheckRule([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b);`): elif nops([args])=1 and op(1,[args])=Children then print(`Children(F,a,b,t): Given a rule F, phrased in terms of a and`): print(`b, and a trajectory t, finds all its legal extensions`): print(`by one step. For example, try:`): print(`Children(F3,a,b,[[a,b,a-b],[1,2,1]]);`): elif nops([args])=1 and op(1,[args])=Concretize1 then print(`Concretize1(R,p,N): given a regular expression R`): print(`in terms of a list of parameters p and a list of`): print(`numeric pos. integers N of the same size,`): print(`realizes the regular expression R`): print(`For example, try:`): print(`Concretize1([[1,2],[[2,1,2],a]],[a],[2]);`): elif nops([args])=1 and op(1,[args])=CondPair then print(`CondPair(Z1,Z2,a,b): given pairs Z1, Z2 of linear combinations`): print(`of a and b finds the condition that there is hope for`): print(`there to be a non-zero solution to Z1=Z2`): print(`For example, try:`): print(`CondPair([c1*a+c2*b,c3*a+c4*b],[a,b],a,b);`): elif nops([args])=1 and op(1,[args])=DAP then print(`DAP(F,a,b,k,M): Given a rule F, phrased in terms`): print(`of a and b defining a recurrence`): print(`x[n]=F(x[n-1],x[n-2]),`): print(`finds a boundedness scheme (if it exists)`): print(`of depth <=k, followed by all the periodic orbits`): print(`of 3<=length <=M`): print(`followed by a proof that these are all of them.`): print(`For example, try:`): print(`DAP(E4,a,b,10,6);`): elif nops([args])=1 and op(1,[args])=DAPf then print(`DAPf(F,a,b,k,M): fast version of DAP(F,a,b,k,M); `): print(`Given a rule F, phrased in terms`): print(`of a and b defining a recurrence`): print(`x[n]=F(x[n-1],x[n-2]),`): print(`finds a boundedness scheme (if it exists)`): print(`of depth <=k, followed by all the periodic orbits`): print(`of 3<=length <=M`): print(`followed by a proof that these are all of them.`): print(`For example, try:`): print(`DAPf(E4,a,b,10,6);`): elif nops([args])=1 and op(1,[args])=DAPfA then print(`DAPfA(F,a,b,k,M): appreviate form of DAPf(F,a,b,k,M); `): print(`For example, try:`): print(`DAPfA(E2,a,b,10,20);`): elif nops([args])=1 and op(1,[args])=DAPsf then print(`DAPsf(F,a,b,k,M1,M2): A combination of DAP and DAPf`): print(`does like DAP (i.e. slowly, but stricter) until M1`): print(`and between M1+1 and M2, does like DAPf`): print(`DAPsf(E4,a,b,10,14,30);`): elif nops([args])=1 and op(1,[args])=DAPsurv then print(`DAPsurv(F,a,b,k,M,nosaf): All the trajectories' profiles`): print(`with the DAP (see entry) convention that survive nosaf `): print(`number of steps. For example, try:`): print(`DAPsurv(E4,a,b,10,20,2);`): elif nops([args])=1 and op(1,[args])=DAPsurvf then print(`DAPsurvf(F,a,b,k,M,nosaf): (fast version)`): print(`All the trajectories' profiles`): print(`with the DAP (see entry) convention that survive nosaf `): print(`number of steps. For example, try:`): print(`DAPsurvf(E4,a,b,10,20,2);`): elif nops([args])=1 and op(1,[args])=DAPsurvY then print(`DAPsurvY(F,a,b,k,M,nosaf): like DAPsurv, but nice`): print(`For example, try:`): print(`DAPsurvY(E6,a,b,10,20,2);`): elif nops([args])=1 and op(1,[args])=DAPsurvYf then print(`DAPsurvYf(F,a,b,k,M,nosaf): (fast version) like DAPsurvf, but nice`): print(`For example, try:`): print(`DAPsurvYf(E6,a,b,10,20,2);`): elif nops([args])=1 and op(1,[args])=DiscoverAtoms then print(`DiscoverAtoms(F,a,b,k,M,nosaf): All the `): print(`possible atoms in DAPsurv(F,a,b,k,M,nosaf) (see entry).`): print(`For example, try:`): print(`DiscoverAtoms(E4,a,b,10,20,2);`): elif nops([args])=1 and op(1,[args])=DiscoverAtomsf then print(`DiscoverAtomsf(F,a,b,k,M,nosaf): (fast version) All the `): print(`possible atoms in DAPsurvf(F,a,b,k,M,nosaf) (see entry).`): print(`For example, try:`): print(`DiscoverAtomsf(E4,a,b,10,20,2);`): elif nops([args])=1 and op(1,[args])=DG1 then print(`DG1(F,a,b): The directed graph corresponding to the second-order`): print(`recurrence given by rule F. For example, try:`): print(`DG1([[(-a+b)/2,-a-b],[-a-b,(-a+b)/2]],a,b);`): elif nops([args])=1 and op(1,[args])=DressUp then print(`DressUp(SymT,a,b,m,K): Given a symbolic trajectory`): print(`SymT in terms of a and b followed by the mod m description`): print(`find all numerical manifestations that follow the template`): print(`with the last entries less than K.`): print(`For example, try:`): print(`DressUp([[a,b,a-b],[1,2,1]],a,b,2,4);`): elif nops([args])=1 and op(1,[args])=EfoBa then print(`EfoBa(F,a,b,i,j): if you apply rule F `): print(`(phrased in terms of variables a and b) to inputs`): print(`such that a=i mod m and b=j mod m, what can`): print(`the mod of the output be?. For example, try:`): print(`EfoBa(F2,a,b,1,1);`): elif nops([args])=1 and op(1,[args])=FindMates then print(`FindMates(i,j,m,c): Given an integer m, and integers i,j`): print(`between 1 and m inclusive, finds all the d's`): print(`such that i*c+j*d==0 (mod m)`): print(`For example, try;`): print(`FindMates(1,1,2,1);`): elif nops([args])=1 and op(1,[args])=GRE then print(`GRE(F,a,b,k,M,nosaf,m): like DAPsurvf, but guesses`): print(`a set of regular expressions covering the language`): print(`of trajectories given by rule F `): print(`For example, try:`): print(`GRE(E2,a,b,10,20,7,m);`) elif nops([args])=1 and op(1,[args])=GREplus then print(`GREplus(F,a,b,k,M,nosaf,m): like GRE, but also`): print(`with suffixes. In other words, it guesses`): print(`a set of regular expressions covering the language`): print(`of trajectories given by rule F `): print(`For example, try:`): print(`GREplus(E2,a,b,10,30,3,m);`) elif nops([args])=1 and op(1,[args])=GuessST then print(`GuessST(F,a,b,k,M,nosaf): guesses symbolic trajectories`): print(`for rule F. For example,try:`): print(`GuessST(E3,a,b,10,40,4);`): elif nops([args])=1 and op(1,[args])=gDAPf then print(`gDAPf(F,a,b,S,M): Given a rule F, phrased in terms`): print(`of a and b defining a recurrence`): print(`x[n]=F(x[n-1],x[n-2]),`): print(`and a prescheme S`): print(`finds all legal trajectories`): print(`followed by all the periodic orbits`): print(`of 3<=length <=M`): print(`followed by a proof that these are all of them.`): print(`For example, try:`): print(`gDAPf(E2,a,b,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]} ,10);`): elif nops([args])=1 and op(1,[args])=gDAPsurvf then print(`gDAPsurvf(F,a,b,S,M,nosaf): (fast version)`): print(`All the trajectories' profiles`): print(`with the gDAP (see entry) convention that survive nosaf `): print(`number of steps. For example, try:`): print(`gDAPsurvf(E2,a,b,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]}, 20,2);`): elif nops([args])=1 and op(1,[args])=gGREplus then print(`gGREplus(F,a,b,S,M,nosaf,m): like gGRE, but also`): print(`with suffixes. In other words, it guesses`): print(`a set of regular expressions covering the language`): print(`of trajectories given by rule F `): print(`For example, try:`): print(`gGREplus(E2,a,b,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]},30,3,m);`) elif nops([args])=1 and op(1,[args])=gDiscoverAtomsf then print(`gDiscoverAtomsf(F,a,b,S,M,nosaf): (fast version) All the `): print(`possible atoms in gDAPsurvf(F,a,b,S,M,nosaf) (see entry).`): print(`For example, try:`): print(`gDiscoverAtomsf(E4,a,b,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]},20,2);`): elif nops([args])=1 and op(1,[args])=gGRE then print(`gGRE(F,a,b,k,M,nosaf,m): like gDAPsurvf, but guesses`): print(`a set of regular expressions covering the language`): print(`of trajectories given by rule F `): print(`For example, try:`): print(`gGRE(E2,a,b, {[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]} ,20,7,m);`) elif nops([args])=1 and op(1,[args])=gMakeScheme then print(`gMakeScheme(F,a,b,K,orekh,Ma): tries to find a`): print(`scheme for the rule F`): print(`For example, try:`): print(`gMakeScheme(E4,a,b,30,30,4);`): elif nops([args])=1 and op(1,[args])=gYeladimABf then print(`gYeladimABf(F,a,b,t,PS): Given a rule F, phrased in terms of a and`): print(`b, and a trajectory t, and a prescheme PS`): print(`finds all its legal extensions`): print(`by one step. For example, try:`): print(`gYeladimABf(E3,a,b,[[a,-b,a+b],[1,2,1]], { [[1,1],[1,2]], [[1,2],[1,2]],[[2,1],[1,2]] } );`): elif nops([args])=1 and op(1,[args])=Hafel then print(`Hafel(F,a,b,a0,b0): Given a rule F (of type m=nops(F))`): print(`where F[i][j] is the affine linear expression`): print(`in a and b that applies if a=i(mod m) and b=j(mod(m))`): print(`applied to a0,b0. For example try,`): print(`Hafel([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3);`): elif nops([args])=1 and op(1,[args])=HafelS then print(`HafelS(F,a,b,Zug,i,j): Given a rule F (of type m=nops(F))`): print(`where F[i][j] is the linear expression`): print(`in a and b that applies if a=i(mod m) and b=j(mod(m))`): print(`applied to symbolic Zug, assuming that`): print(` Zug[1]=i(mod m) and Zug[2]=j(mod(m))`): print(` For example try,`): print(`HafelS([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,[a,b],2,1);`): elif nops([args])=1 and op(1,[args])=HafelSw then print(`HafelSw(F,a,b,Zug,w): Given a rule F (of type m=nops(F))`): print(`where F[i][j] is the linear expression`): print(`in a and b that applies if a=i(mod m) and b=j(mod(m))`): print(`applied to the current window`): print(`[a,b], assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m))`): print(`keeps doing it with the parities given by the word in`): print(`{1,2, ..., m}*, w,`): print(` For example try,`): print(`HafelSw(E2,a,b,[a,b],[1,1,2]);`): elif nops([args])=1 and op(1,[args])=HafelSwn then print(`HafelSwn(F,a,b,Zug,w,n): Given a rule F (of type m=nops(F))`): print(`where F[i][j] is the linear expression`): print(`in a and b that applies if a=i(mod m) and b=j(mod(m))`): print(`applied to the current window`): print(`[a,b], assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m))`): print(`keeps doing it with the parities given by the word in`): print(`{1,2, ..., m}*, w, iterated n times, for symbolic n`): print(` For example try,`): print(`HafelSwn(E2,a,b,[a,b],[1,1,2],n);`): elif nops([args])=1 and op(1,[args])=HafelSre then print(`HafelSre(F,a,b,R): Given a rule F (of type m=nops(F))`): print(`where F[i][j] is the linear expression`): print(`in a and b that applies if a=i(mod m) and b=j(mod(m))`): print(`applied to the current window`): print(`[a,b], assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m))`): print(`keeps doing it with the parities given by the RE R`): print(`For example, try:`): print(`HafelSre(E2,a,b,[[2,1,1],[[1],n]]);`): elif nops([args])=1 and op(1,[args])=HopeForOrbits then print(`HopeForOrbits(F,a,b,k,M,nosaf,m): looks at the`): print(`condition, phrased in terms of m[1],m[2]'s etc.`): print(`that there is a non-trivial orbit. For example try:`): print(`HopeForOrbits(E2,a,b,10,20,7,m);`) : elif nops([args])=1 and op(1,[args])=ImproveScheme then print(`ImproveScheme(Sch,F,a,b): Given a scheme for a `): print(`second-order recurrence x[n]=F(x[n-1],x[n-2])`): print(`where F is phrased in terms of a and b,`): print(`improves it by adding one more condition`): print(`For example, try`): print(`ImproveScheme(SchF1V,F1,a,b);`): elif nops([args])=1 and op(1,[args])=IsCont then print(`IsCont(w,gw,m,t): is the concrete word w a continuation`): print(`of a specialization of the generalized word m`): print(`phrased in terms of the index variable m`): print(`output is the smallest possible suffix of length <=t or FAIL`): print(`For example try: IsCont([1,2,1,2,1,2,1,2,2,2],[[[1,2],m[1]]],m,4);`): elif nops([args])=1 and op(1,[args])=IsFL then print(`IsFL(G,v): Given a directed graph G, and a vertex,`): print(`decides whether it is friend-less (except itself)`): print(`For example, try:`): print(`IsFL({[1,{1,2}],[2,{2}]},1);`): elif nops([args])=1 and op(1,[args])=IsPrim then print(`IsPrim(P): Is the period P primitive?`): print(`For example try: `): print(`IsPrim([3,-1,2,3,-1,2]);`): elif nops([args])=1 and op(1,[args])=KavuaScheme then print(`KavuaScheme(Sch,a,b): An upper bound for a scheme Sch`): print(`phrased in terms of a and b. For example, try:`): print(`KavuaScheme(MakeScheme(E2,a,b,10)[1],a,b);`): elif nops([args])=1 and op(1,[args])=MakeWebBook then print(`MakeWebBook(kvu,a,b, K, M, K1,orekh,Ma,nosaf,m,x,n,shem1,shem2):`): print(`makes a web-book from the set of previously computed`): print(`good books. Each rule has two outputs, cat(shem1,i) for`): print(`the scheme and cat(shem2, i) for the rest (Sipur1V)`): print(`For example, try:`): print(`MakeWebBook({E2,E3},a,b,100, 20, 15, 15, 4, 4,m,x,n,Sch,Sip);`): elif nops([args])=1 and op(1,[args])=MakeScheme then print(`MakeScheme(F,a,b,K): Given a second-order recurrence`): print(`x[n]=F(x[n-1],x[n-2]) given in terms of rule F`): print(`phrased in terms of a and b, tries to`): print(`make a scheme with up-to K improvements.`): print(`If it succeeds, it returns true, followed by`); print(`the successful scheme, if it fails, it returns`): print(`false followed by the unsuccessful last attempt`): print(`(with K additions)`): print(`For example, try:`): print(`MakeScheme([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,10);`): elif nops([args])=1 and op(1,[args])=MakeSchemeVerbose then print(`MakeSchemeVerbose(F,a,b,K): Verbose`): print(`version of MakeScheme(F,a,b,K) (see entry)`): print(`For example, try:`): print(`MakeSchemeVerbose(E2,a,b,10):`): elif nops([args])=1 and op(1,[args])=MakeTrans then print(`MakeTrans(C1,C2,a,b): given two lists of lists of size m, say,`): print(`makes the transformation with these coefficients`): print(`For example, try: MakeTrans([[1,2],[2,1]],[[1,-2],[-2,1]],a,b);`): elif nops([args])=1 and op(1,[args])=MakeTransR then print(`MakeTransR(C1,a,b): given a list of lists of size m, say,`): print(`makes a random transformation with C1 describing the`): print(`coefficients of x[n-1].`): print(`For example, try: MakeTransR([[1,2],[2,1]],a,b);`): elif nops([args])=1 and op(1,[args])=MakomRishon then print(`MakomRishon(L): the largest location in a circular list L.`): print(`For example, try:`): print(`MakomRishon([1,3,-1,2]);`): elif nops([args])=1 and op(1,[args])=MaxTrajs then print(`MaxTrajs(F,a,b,M,K): all the triples`): print(` [abs(a0),abs(b0),MaxTraj1(F,a,b,a0,b0,K)]`): print(`where the last component is bigger than abs(a0)+abs(b0)`): print(`for a0 and b0 relatively prime between -M and M.`): print(`For example, try:`): print(`MaxTrajs([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,5,100);`); elif nops([args])=1 and op(1,[args])=MaxTrajsa0 then print(`MaxTrajsa0(F,a,b,a0,M,K): all the triples`): print(` [abs(a0),abs(b0),MaxTraj1(F,a,b,a0,b0,K)]`): print(`where the last component is bigger than abs(a0)+abs(b0)`): print(`for a0 given as (fourth) input`): print(`and b0 relatively prime to a0, between -M and M.`): print(`For example, try:`): print(`MaxTrajsa0([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,1,5,100);`); elif nops([args])=1 and op(1,[args])=MaxTraj1 then print(`MaxTraj1(F,a,b,a0,b0,K): the max. magnitude of elements`): print(`of the trajectory of length K starting`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`MaxTraj1([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100);`): elif nops([args])=1 and op(1,[args])=MitatSdom then print(`MitatSdom(w,S): Given a word w and a set of templates`): print(`and double templates`): print(`expresses it in terms of [v1,[T1,m1],[T2,m2],v2] etc. where`): print(`v1 is a prefix and v2 is a suffix. For example, try`): print(`MitatSdom([1,2,1,2,1,2,1,1,1,1,1,1],{[[2,1],[1]]});`): elif nops([args])=1 and op(1,[args])=MitatSdom1 then print(`MitatSdom1(w,T): Given a word w and a single template T`): print(`expresses it in terms of [v1,[T,m],v2] where`): print(`v1 is a prefix and v2 is a suffix. For example, try`): print(`MitatSdom1([1,2,1,2,1,2,1,2],[2,1]);`): elif nops([args])=1 and op(1,[args])=MitatSdom2 then print(`MitatSdom2(w,T): Given a word w and a double template T=[T1,T2]`): print(`expresses it in terms of [v1,[T1,m1],v2,[T2,m2],v3] where`): print(`v1 is a prefix , v2 is inside, and v3 is a suffix. For example, try`): print(`MitatSdom2([1,2,1,2,1,2,1,2,1,1,1,1,1],[[2,1],[1]]);`): elif nops([args])=1 and op(1,[args])=Novea then print(`Novea(S,L,a,b): Given a set of linear forms S in (a,b)`): print(`and a single linear form L in (a,b), decides whether`): print(`|L|<=A (for some A) is implied by the inequalities`): print(`|S1|<=A, |S2|<=A, for members of S, taken one-at-a-time`): print(`or two-at-a-time. For example,try:`): print(`Novea({a,b,a-b},b-(9/10)*a,a,b);`): elif nops([args])=1 and op(1,[args])=Novea1 then print(`Novea1(L1,L2): Given Linear forms L1, L2`): print(`decides whether abs(L1/L2)=1. For example, try: `): print(`Novea1((a+b)/2,(-a-b)/2);`): elif nops([args])=1 and op(1,[args])=Novea2 then print(`Novea2(L1,L2,L3,a,b): Given Linear forms L1, L2,L3 in`): print(`(a,b) expresses L3 as a linear combination of L1 and`): print(`L2 and decides whether abs(L1)<=A, abs(L2)<=A implies`): print(`abs(L3)<=A. For example, try`): print(`Novea2(a,b,(a+b)/2,a,b);`): elif nops([args])=1 and op(1,[args])=Novea2V then print(`Verbose form of Novea2(L1,L2,L3,a,b,A)`): print(`Novea2V(L1,L2,L3,a,b,A): Given Linear forms L1, L2,L3 in`): print(`(a,b) expresses L3 as a linear combination of L1 and`): print(`L2 and decides whether abs(L1)<=A, abs(L2)<=A implies`): print(`abs(L3)<=A. For example, try`): print(`Novea2V(a,b,(a+b)/2,a,b,A);`): elif nops([args])=1 and op(1,[args])=Orbits then print(`Orbits(F,a,b,k): all the orbits of period k of the recurrence`): print(`x[n]=F(x[n-1],x[n-2]), where F is phrased in terms of a and b.`): print(`For example, try: `): print(`Orbits(F3,a,b,3);`): elif nops([args])=1 and op(1,[args])=ParseWord then print(`ParseWord(w,t1,t2): Given a word w, represents it as`): print(`v1 w1^a*w2^b....v2, with nops(v1)<=t1`): print(`and nops(v2)<=t2. For example, try`): print(`ParseWord([1,2,1,1,2,1,1,2,1,1,1,1,1,1],0,0);`): elif nops([args])=1 and op(1,[args])=ParseWord1 then print(`ParseWord1(w): Given a word w, represents it as`): print(`w1^a*w2^b.... For example, try`): print(`ParseWord1([1,2,1,1,2,1,1,2,1,1,1,1,1,1]);`): elif nops([args])=1 and op(1,[args])=Per then print(`Per(F,a,b,a0,b0,K): the ultimate period`): print(`if it is found, or FAIL, trying the trajectory`): print(`up to length K starting`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`Per([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100);`): elif nops([args])=1 and op(1,[args])=PerMod then print(`PerMod(F,a,b,a0,b0,K): the ultimate period`): print(`if it is found of the mod m of trajectories`): print(`, or FAIL, trying the trajectory`): print(`up to length K starting`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`PerMod([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100);`): elif nops([args])=1 and op(1,[args])=Per1 then print(`Per1(L): Given a list L,decides whether it ends`): print(`with a period (repeat), and returns it.`): print(`If it fails,it returns FAIL.`): print(`For example, try: `): print(`Per1([3,4,1,3,2,1,3,2]);`): elif nops([args])=1 and op(1,[args])=Pers then print(`Pers(F,a,b,K,M): the ultimate periods`): print(`for trajectories`): print(`up to length K starting for all initial conditions`): print(`(a0,b0) between -M and M`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`Pers([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,100,5);`): elif nops([args])=1 and op(1,[args])=PersMod then print(`PersMod(F,a,b,K,M): the ultimate periods`): print(`for trajectories mod m`): print(`up to length K starting for all initial conditions`): print(`(a0,b0) between -M and M`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`PersMod([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,100,5);`): elif nops([args])=1 and op(1,[args])=Persa0 then print(`Persa0(F,a,b,a0,K,M): the ultimate period`): print(`if it is found, or FAIL, trying the trajectory`): print(`up to length K starting for all initial conditions`): print(`(a0(fixed),b0) and b0 between -M and M`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`Persa0([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,1,100,5);`): elif nops([args])=1 and op(1,[args])=FindPreScheme then print(`FindPreScheme(F,a,b,K,orekh,Ma): tries to find a`): print(`prescheme for the rule F`): print(`For example, try:`): print(`FindPreScheme(E2,a,b,30,30,4);`): elif nops([args])=1 and op(1,[args])=FindPreSchemeG then print(`FindPreSchemeG(F,a,b,K,orekh,Ma,PreLim): tries to find a`): print(`prescheme for the rule F, after letting it run for PreLim steps`): print(`For example, try:`): print(`FindPreSchemeG(E2,a,b,30,30,4,10);`): elif nops([args])=1 and op(1,[args])=FindPreScheme1 then print(`FindPreScheme1(F,a,b,i,j,K,orekh,Ma): tries to find a`): print(`prescheme for the rule F for x[n-1] = i mod m, x[n]=j mod m`): print(`For example, try:`): print(`FindPreScheme1(E2,a,b,1,1,30,30,2);`): elif nops([args])=1 and op(1,[args])=FindPreScheme1G then print(`FindPreScheme1G(F,a,b,i,j,K,orekh,Ma,PreLim): tries to find a`): print(`prescheme for the rule F for x[n-1] = i mod m, x[n]=j mod m`): print(`after letting it run for PreLim steps.`): print(`For example, try:`): print(`FindPreScheme1G(E2,a,b,1,1,30,30,2,10);`): elif nops([args])=1 and op(1,[args])=PtorAB then print(`PtorAB(L,a,b,A,B): Given a list of linear expressions in`): print(`a and b, finds the solution of `): print(`abs(L[1])=3,`): print(`finds all good trajectories of length k, `): print(`(where the third element is the largest in absolute value)`): print(`starting with a and b. `): print(`For example, try:`): print(`SymbTrajs(F1,a,b,4);`): elif nops([args])=1 and op(1,[args])=gSymbTrajsABf then print(`gSymbTrajsABf(F,a,b,k,PS) `): print(`Given a rule F phrased in terms of a and b`): print(`and an integer k>=3,and a prescheme PS`): print(`finds all good trajectories of length k, `): print(`(where the third element is the largest in absolute value)`): print(`and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!)`): print(`starting with a and b. `): print(`For example, try:`): print(`gSymbTrajsABf(E2,a,b,4,{[[1,1],[1,1]], [[1,2],[1,1]], [[2,1],[1,1]]});`): elif nops([args])=1 and op(1,[args])=SymbTrajsAB then print(`SymbTrajsAB(F,a,b,k,A,B): Given a rule F phrased in terms of a and b`): print(`and an integer k>=3,`): print(`finds all good trajectories of length k, `): print(`(where the third element is the largest in absolute value)`): print(`and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!)`): print(`starting with a and b. `): print(`For example, try:`): print(`SymbTrajsAB(E2,a,b,4,1,1);`): elif nops([args])=1 and op(1,[args])=SymbTrajsABf then print(`SymbTrajsABf(F,a,b,k,A,B): fast version of SymbTrajsAf(F,a,b,k,A,B)`): print(`in other words `): print(`Given a rule F phrased in terms of a and b`): print(`and an integer k>=3,`): print(`finds all good trajectories of length k, `): print(`(where the third element is the largest in absolute value)`): print(`and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!)`): print(`starting with a and b. `): print(`For example, try:`): print(`SymbTrajsABf(E2,a,b,4,1,1);`): elif nops([args])=1 and op(1,[args])=SymbTrajsC then print(`SymbTrajsC(F,a,b,k): Given a rule F phrased in terms of a and b`): print(`and an integer k>=2,`): print(`finds all good trajectories of length k, `): print(`starting with a and b. `): print(`For example, try:`): print(`SymbTrajsC(F1,a,b,4);`): elif nops([args])=1 and op(1,[args])=TestABij then print(`TestABij(F,a,b,A,B,i,j,K,orekh): tests whether it is true for`): print(`starting values in [-K,K] that abs(x[n])<=A*abs(x[-1])+B*abs(x[-2])`): print(`whenever x[n-1] mod m=i and x[n] mod m=j`): print(`up to trajectories of length orekh. For example, try:`): print(`TestABij(E2,a,b,1,1,1,1,10,30);`): elif nops([args])=1 and op(1,[args])=Tovim then print(`Tovim(a,b, K, M, K1,orekh,Ma,nosaf,m): `): print(`trying out all possible 256 rules, and finds the good ones`): print(`For example, try`): print(`Tovim(a,b,200, 20, 15, 15, 4, 7 ,m);`): elif nops([args])=1 and op(1,[args])=Traj then print(`Traj(F,a,b,a0,b0,K): the trajectory of length K starting`): print(`at x[1]=a0, x[2]=b0, of the recurrence`): print(`x[n]=F(x[n-1],x[n]) given in terms of the rule F`): print(`phrased in terms of the variables a and b.`): print(`For example, try:`): print(`Traj([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100);`): elif nops([args])=1 and op(1,[args])=Tsamek then print(`Shrinks the lists by dividing them by their igcd`): elif nops([args])=1 and op(1,[args])=VecToRule then print(`VecToRule(v,a,b): inputs a -1,1 vector of length 8 and outputs the appropriate rule of type mod 2`): print(`For example, try:`): print(`VecToRule([1,-1,1,-1,1,-1,1,-1],a,b):`): elif nops([args])=1 and op(1,[args])=VerifyScheme then print(`VerifyScheme(Sch,F,a,b): Given a scheme for a `): print(`second-order recurrence x[n]=F(x[n-1],x[n-2])`): print(`where F is phrased in terms of a and b,`): print(`checks whether it is a valid scheme. For example, try`): print(`VerifyScheme(SchF1,F1,a,b);`): elif nops([args])=1 and op(1,[args])=VerifySchemeV then print(`VerifySchemeV(Sch,F,a,b,NAME): a verbose form of `): print(`VerifyScheme(Sch,F,a,b): in other words, given a scheme for a `): print(`second-order recurrence x[n]=F(x[n-1],x[n-2])`): print(`where F is phrased in terms of a and b,`): print(`checks whether it is a valid scheme. NAME is the`): print(` name of the statement. For example, try`): print(`VerifySchemeV(MakeScheme(E2,a,b,10),E2,a,b,2);`): elif nops([args])=1 and op(1,[args])=VerifySchemeE then print(`VerifySchemeE(Sch,F,a,b,M,K): Empirically verifies the`): print(`inequalities implied by the scheme for initial values`): print(`between -M and M looking at trajectories of length K.`): print(`For example, try:`): print(`VerifySchemeE(MakeScheme(E2,a,b,10),E2,a,b,10,50);`): elif nops([args])=1 and op(1,[args])=WeedOut then print(`WeedOut(G): Given a directed graph, kicks out`): print(`all friend-less vertices (and of course their outgoing edges)`): print(`For example, try:`): print(`WeedOut(DG1(E1,a,b));`): elif nops([args])=1 and op(1,[args])=Yeladim then print(`Yeladim(F,a,b,t): Given a rule F, phrased in terms of a and`): print(`b, and a trajectory t, finds all its legal extensions`): print(`in a period where the maximal element is the third`): print(`by one step. For example, try:`): print(`Yeladim(F3,a,b,[[a,b,a-b],[1,2,1]]);`): elif nops([args])=1 and op(1,[args])=YeladimAB then print(`YeladimAB(F,a,b,t,A,B): Given a rule F, phrased in terms of a and`): print(`b, and a trajectory t, finds all its legal extensions`): print(`in a period where the maximal element is the third`): print(`and |x[3]|<=A|x[n-1]|+B|x[n]|`): print(`by one step. For example, try:`): print(`YeladimAB(F3,a,b,[[a,b,a-b],[1,2,1]],1,1);`): elif nops([args])=1 and op(1,[args])=YeladimABf then print(`YeladimABf(F,a,b,t,A,B): fast version of YeladimAB(F,a,b,t,A,B)`): print(`Given a rule F, phrased in terms of a and`): print(`b, and a trajectory t, finds all its legal extensions`): print(`in a period where the maximal element is the third`): print(`and |x[3]|<=A|x[n-1]|+B|x[n]|`): print(`by one step. For example, try:`): print(`YeladimABf(F3,a,b,[[a,b,a-b],[1,2,1]],1,1);`): elif nops([args])=1 and op(1,[args])=YeladimRE then print(`YeladimRE(F,a,b,Ini,A,B,R,m,r,w,K1,K2): Given `): print(`a rule F phrased in terms of a and b and`): print(`a regular expression R phrased in terms of m`): print(`and a word w finds all the words of length one more`): print(`that continues [op(R),op(w)] and do `): print(`not violate the necessary conditions`): print(`as checked by BatuakhRaSymb (with parameters r, K1, K2).`): print(`For example, try:`): print(`YeladimRE(E7,a,b,[a,-b,a+b],1,2,[[[2,1,1,1],m[1]]],m,1,[],1,10);`): else print(`There is no ezra for`,args): fi: end: ####precomputed good rules Tov1:= [[[-1/2*a-1/2*b, -a-b], [-a-b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, -a-b], [-a-b, -1/ 2*a-1/2*b]], [[1/2*a+1/2*b, -a-b], [-a-b, -1/2*a-1/2*b]], [[-1/2*a+1/2*b, a-b], [-a-b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, a-b], [-a-b, -1/2*a-1/2*b]], [[-1/2*a+1/2 *b, -a+b], [-a-b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, -a+b], [-a-b, -1/2*a-1/2*b]], [[-1/2*a-1/2*b, a+b], [-a-b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, a+b], [-a-b, -1/2*a -1/2*b]], [[-1/2*a-1/2*b, -a-b], [a-b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, -a-b], [a -b, -1/2*a-1/2*b]], [[-1/2*a+1/2*b, a-b], [a-b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, a-b], [a-b, -1/2*a-1/2*b]], [[-1/2*a+1/2*b, -a+b], [a-b, -1/2*a-1/2*b]], [[1/2* a+1/2*b, -a+b], [a-b, -1/2*a-1/2*b]], [[-1/2*a-1/2*b, a+b], [a-b, -1/2*a-1/2*b] ], [[1/2*a-1/2*b, a+b], [a-b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, a+b], [a-b, -1/2*a -1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a+b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, -a-b], [ -a+b, -1/2*a-1/2*b]], [[-1/2*a+1/2*b, a-b], [-a+b, -1/2*a-1/2*b]], [[1/2*a+1/2* b, a-b], [-a+b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, -a+b], [-a+b, -1/2*a-1/2*b]], [[ -1/2*a+1/2*b, -a+b], [-a+b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, -a+b], [-a+b, -1/2*a -1/2*b]], [[-1/2*a-1/2*b, a+b], [-a+b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, a+b], [-a +b, -1/2*a-1/2*b]], [[-1/2*a-1/2*b, -a-b], [a+b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, -a-b], [a+b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, a-b], [a+b, -1/2*a-1/2*b]], [[-1/2* a+1/2*b, a-b], [a+b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, a-b], [a+b, -1/2*a-1/2*b]], [[-1/2*a+1/2*b, -a+b], [a+b, -1/2*a-1/2*b]], [[1/2*a+1/2*b, -a+b], [a+b, -1/2*a -1/2*b]], [[-1/2*a-1/2*b, a+b], [a+b, -1/2*a-1/2*b]], [[1/2*a-1/2*b, a+b], [a+b , -1/2*a-1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a-b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, - a-b], [-a-b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, -a-b], [-a-b, 1/2*a-1/2*b]], [[-1/2* a+1/2*b, a-b], [-a-b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, a-b], [-a-b, 1/2*a-1/2*b]], [[-1/2*a+1/2*b, -a+b], [-a-b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, -a+b], [-a-b, 1/2*a -1/2*b]], [[-1/2*a-1/2*b, a+b], [-a-b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, a+b], [-a- b, 1/2*a-1/2*b]], [[-1/2*a-1/2*b, -a-b], [a-b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, -a -b], [a-b, 1/2*a-1/2*b]], [[-1/2*a+1/2*b, a-b], [a-b, 1/2*a-1/2*b]], [[1/2*a+1/ 2*b, a-b], [a-b, 1/2*a-1/2*b]], [[-1/2*a+1/2*b, -a+b], [a-b, 1/2*a-1/2*b]], [[1 /2*a+1/2*b, -a+b], [a-b, 1/2*a-1/2*b]], [[-1/2*a-1/2*b, a+b], [a-b, 1/2*a-1/2*b ]], [[1/2*a-1/2*b, a+b], [a-b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, a+b], [a-b, 1/2*a-\ 1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a+b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, -a-b], [-a +b, 1/2*a-1/2*b]], [[-1/2*a+1/2*b, a-b], [-a+b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, a -b], [-a+b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, -a+b], [-a+b, 1/2*a-1/2*b]], [[-1/2*a +1/2*b, -a+b], [-a+b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, -a+b], [-a+b, 1/2*a-1/2*b]] , [[-1/2*a-1/2*b, a+b], [-a+b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, a+b], [-a+b, 1/2*a -1/2*b]], [[-1/2*a-1/2*b, -a-b], [a+b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, -a-b], [a+ b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, a-b], [a+b, 1/2*a-1/2*b]], [[-1/2*a+1/2*b, a-b ], [a+b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, a-b], [a+b, 1/2*a-1/2*b]], [[-1/2*a+1/2* b, -a+b], [a+b, 1/2*a-1/2*b]], [[1/2*a+1/2*b, -a+b], [a+b, 1/2*a-1/2*b]], [[-1/ 2*a-1/2*b, a+b], [a+b, 1/2*a-1/2*b]], [[1/2*a-1/2*b, a+b], [a+b, 1/2*a-1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a-b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, -a-b], [-a-b, -1/2 *a+1/2*b]], [[1/2*a+1/2*b, -a-b], [-a-b, -1/2*a+1/2*b]], [[-1/2*a+1/2*b, a-b], [-a-b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, a-b], [-a-b, -1/2*a+1/2*b]], [[-1/2*a+1/2 *b, -a+b], [-a-b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, -a+b], [-a-b, -1/2*a+1/2*b]], [[-1/2*a-1/2*b, a+b], [-a-b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, a+b], [-a-b, -1/2*a +1/2*b]], [[-1/2*a-1/2*b, -a-b], [a-b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, -a-b], [a -b, -1/2*a+1/2*b]], [[-1/2*a+1/2*b, a-b], [a-b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, a-b], [a-b, -1/2*a+1/2*b]], [[-1/2*a+1/2*b, -a+b], [a-b, -1/2*a+1/2*b]], [[1/2* a+1/2*b, -a+b], [a-b, -1/2*a+1/2*b]], [[-1/2*a-1/2*b, a+b], [a-b, -1/2*a+1/2*b] ], [[1/2*a-1/2*b, a+b], [a-b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, a+b], [a-b, -1/2*a +1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a+b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, -a-b], [ -a+b, -1/2*a+1/2*b]], [[-1/2*a+1/2*b, a-b], [-a+b, -1/2*a+1/2*b]], [[1/2*a+1/2* b, a-b], [-a+b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, -a+b], [-a+b, -1/2*a+1/2*b]], [[ -1/2*a+1/2*b, -a+b], [-a+b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, -a+b], [-a+b, -1/2*a +1/2*b]], [[-1/2*a-1/2*b, a+b], [-a+b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, a+b], [-a +b, -1/2*a+1/2*b]], [[-1/2*a-1/2*b, -a-b], [a+b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, -a-b], [a+b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, a-b], [a+b, -1/2*a+1/2*b]], [[-1/2* a+1/2*b, a-b], [a+b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, a-b], [a+b, -1/2*a+1/2*b]], [[-1/2*a+1/2*b, -a+b], [a+b, -1/2*a+1/2*b]], [[1/2*a+1/2*b, -a+b], [a+b, -1/2*a +1/2*b]], [[-1/2*a-1/2*b, a+b], [a+b, -1/2*a+1/2*b]], [[1/2*a-1/2*b, a+b], [a+b , -1/2*a+1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a-b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, - a-b], [-a-b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, -a-b], [-a-b, 1/2*a+1/2*b]], [[-1/2* a+1/2*b, a-b], [-a-b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, a-b], [-a-b, 1/2*a+1/2*b]], [[-1/2*a+1/2*b, -a+b], [-a-b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, -a+b], [-a-b, 1/2*a +1/2*b]], [[-1/2*a-1/2*b, a+b], [-a-b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, a+b], [-a- b, 1/2*a+1/2*b]], [[-1/2*a-1/2*b, -a-b], [a-b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, -a -b], [a-b, 1/2*a+1/2*b]], [[-1/2*a+1/2*b, a-b], [a-b, 1/2*a+1/2*b]], [[1/2*a+1/ 2*b, a-b], [a-b, 1/2*a+1/2*b]], [[-1/2*a+1/2*b, -a+b], [a-b, 1/2*a+1/2*b]], [[1 /2*a+1/2*b, -a+b], [a-b, 1/2*a+1/2*b]], [[-1/2*a-1/2*b, a+b], [a-b, 1/2*a+1/2*b ]], [[1/2*a-1/2*b, a+b], [a-b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, a+b], [a-b, 1/2*a+ 1/2*b]], [[-1/2*a-1/2*b, -a-b], [-a+b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, -a-b], [-a +b, 1/2*a+1/2*b]], [[-1/2*a+1/2*b, a-b], [-a+b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, a -b], [-a+b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, -a+b], [-a+b, 1/2*a+1/2*b]], [[-1/2*a +1/2*b, -a+b], [-a+b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, -a+b], [-a+b, 1/2*a+1/2*b]] , [[-1/2*a-1/2*b, a+b], [-a+b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, a+b], [-a+b, 1/2*a +1/2*b]], [[-1/2*a-1/2*b, -a-b], [a+b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, -a-b], [a+ b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, a-b], [a+b, 1/2*a+1/2*b]], [[-1/2*a+1/2*b, a-b ], [a+b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, a-b], [a+b, 1/2*a+1/2*b]], [[-1/2*a+1/2* b, -a+b], [a+b, 1/2*a+1/2*b]], [[1/2*a+1/2*b, -a+b], [a+b, 1/2*a+1/2*b]], [[-1/ 2*a-1/2*b, a+b], [a+b, 1/2*a+1/2*b]], [[1/2*a-1/2*b, a+b], [a+b, 1/2*a+1/2*b]]] : ####end precomputed good rules ####sample functions of the F kind F1:=[[(a+b)/2,a-b],[a-b,(a+b)/2]]: SchF1:= {[[1, 1], {[1, 1], [1, 2]}, {a,b}], [[1, 2], {[2, 1]},{a,b,a-b,a-2*b}], [[2, 1], {[1, 1]},{a,b,a-b}], [[2, 2], {[2, 1], [2, 2]},{a,b} ]}: SchF1V:= {[[1, 1], {[1, 1], [1, 2]}, {a,b}], [[1, 2], {[2, 1]},{a,b}], [[2, 1], {[1, 1]},{a,b}], [[2, 2], {[2, 1], [2, 2]},{a,b} ]}: F1a:=[[(-a+b)/2,a-b],[a-b,(a+b)/2]]: SchF1a:= {[[1, 1], {[1, 1], [1, 2]}, {a,b}], [[1, 2], {[2, 1]},{a,b,a-b,a-2*b}], [[2, 1], {[1, 1]},{a,b,a-b}], [[2, 2], {[2, 1], [2, 2]},{a,b} ]}: F2:=[[(a+b)/2,-a+b],[-a+b,(a+b)/2]]: F3:=[[(-a+b)/2,-a-b],[-a-b,(-a+b)/2]]: F4:=[[2/3*a+1/3*b, -5/3*a-2/3*b, a], [a, 2/3*a-2/3*b, a+2/3*b], [2/3*a+b, 4/3*a-b, 1/3*b]]: F5:=[[2/3*a+1/3*b, 1/3*a+1/3*b, a], [a, -1/3*a+1/3*b, 2/3*b], [-1/3*a+b, 1/3*a-b, -1/3*b]]: ####end sample functions of the F kind ####Sample functions of the E kind (following Amleh et. al.) E1:=[[(a+b)/2,a+b],[a+b,(a+b)/2]]: E1d:=[[(a-b)/2,a-b],[a-b,(a-b)/2]]: E2:=[[(a+b)/2,-a+b],[-a+b,(a+b)/2]]: E2d:=[[(a-b)/2,-a-b],[-a-b,(a-b)/2]]: E3:=[[(a+b)/2,a-b],[a-b,(a+b)/2]]: E3d:=[[(a-b)/2,a+b],[a+b,(a-b)/2]]: E4:=[[(a+b)/2,-a-b],[-a-b,(a+b)/2]]: E4d:=[[(a-b)/2,-a+b],[-a+b,(a-b)/2]]: E5:=[[(-a+b)/2,a+b],[a+b,(-a+b)/2]]: E5d:=[[(-a-b)/2,a-b],[a-b,(-a-b)/2]]: E6:=[[(-a+b)/2,-a+b],[-a+b,(-a+b)/2]]: E6d:=[[(-a-b)/2,-a-b],[-a-b,(-a-b)/2]]: E7:=[[(-a+b)/2,a-b],[a-b,(-a+b)/2]]: E7d:=[[(-a-b)/2,a+b],[a+b,(-a-b)/2]]: E8:=[[(-a+b)/2,-a-b],[-a-b,(-a+b)/2]]: E8d:=[[(-a-b)/2,-a+b],[-a+b,(-a-b)/2]]: #Hafel(F,a,b,a0,b0): Given a rule F (of type m=nops(F)) #where F[i][j] is the affine linear expression #in a and b that applies if a=i(mod m) and b=j(mod(m)) #applied to a0,b0. For example try, #Hafel([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3); Hafel:=proc(F,a,b,a0,b0) local m,i,j: m:=nops(F): if not type(F,list) or {seq(nops(F[i]),i=1..m)}<>{m} then print(`bad input`): RETURN(FAIL): fi: i:=a0 mod m: j:=b0 mod m: if i=0 then i:=m: fi: if j=0 then j:=m: fi: subs({a=a0,b=b0},F[i][j]): end: #CheckRule(F,a,b): checks whether the rule F is a good one #For example, try: CheckRule([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b); CheckRule:=proc(F,a,b) local m,i,j: option remember: m:=nops(F): if not type(F,list) or {seq(nops(F[i]),i=1..m)}<>{m} then RETURN(false): fi: for i from 1 to m do for j from 1 to m do if not type(subs({a=i,b=j},F[i][j]),integer) then RETURN(false): fi: od: od: true: end: #Traj(F,a,b,a0,b0,K): the trajectory of length K starting #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #Traj([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100); Traj:=proc(F,a,b,a0,b0,K) local gu,a1,b1,b2,i: if not CheckRule(F,a,b) then RETURN(FAIL): fi: gu:=[a0,b0]: a1:=a0: b1:=b0: for i from 3 to K do b2:=Hafel(F,a,b,a1,b1): a1:=b1: b1:=b2: gu:=[op(gu),b1]: od: gu: end: #Per1(L): Given a list L,decides whether it ends #with a period (repeat), and returns it. #If it fails,it returns FAIL. #For example, try: #Per1([3,4,1,3,2,1,3,2]); Per1:=proc(L) local i,m: m:=nops(L): for i from 2 to trunc(m/2) do if [op(m-2*i+1..m-i,L)]=[op(m-i+1..m,L)] then RETURN([op(m-i+1..m,L)]): fi: od: FAIL: end: #Per(F,a,b,a0,b0,K): the ultimate period #if it is found, or FAIL, trying the trajectory #up to length K starting #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #Per([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100); Per:=proc(F,a,b,a0,b0,K) local gu,a1,b1,b2,i,lu: if not CheckRule(F,a,b) then RETURN(FAIL): fi: gu:=[a0,b0]: a1:=a0: b1:=b0: for i from 3 to K do b2:=Hafel(F,a,b,a1,b1): a1:=b1: b1:=b2: gu:=[op(gu),b1]: lu:=Per1(gu): if lu<>FAIL then RETURN(lu): fi: od: FAIL: end: #Pers(F,a,b,K,M): the ultimate period #if it is found, or FAIL, trying the trajectory #up to length K starting for all initial conditions #(a0,b0) between -M and M #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #Pers([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,100,5); Pers:=proc(F,a,b,K,M) local gu,a0,b0: if not CheckRule(F,a,b) then RETURN(FAIL): fi: gu:={}: for a0 from -M to M do for b0 from -M to M do gu:=gu union {Canon(Per(F,a,b,a0,b0,K))}: od: od: gu: end: #MakomRishon(L): the largest location in a circular list L #For example, try: #MakomRishon([1,3,-1,2]); MakomRishon:=proc(L) local aluf,si1,si2,i: if nops(L)=1 then RETURN(1): fi: aluf:=1: si1:=abs(L[1]): si2:=abs(L[2]): for i from 2 to nops(L)-1 do if abs(L[i])>si1 then aluf:=i: si1:=abs(L[i]): si2:=abs(L[i+1]): elif abs(L[i])=si1 and abs(L[i+1])>si2 then aluf:=i: si2:=abs(L[i+1]): fi: od: for i from nops(L) to nops(L) do if abs(L[i])>si1 then aluf:=i: si1:=abs(L[i]): elif abs(L[i])=si1 and abs(L[1])>si2 then aluf:=i: si2:=abs(L[1]): fi: od: aluf: end: #Canon(L): the canonical form of a period L #For example, try: Canon([3,6,3,6]); Canon:=proc(L) local i,L1,L2,m,j: if L=FAIL then RETURN(FAIL): fi: if nops(L)=1 then RETURN([1]): fi: L1:=L: m:=gcd(L1[1],L1[2]): if m=0 then RETURN(L1): fi: L2:=[seq(L1[i]/m,i=1..nops(L1))]: if {seq(type(L2[i],integer),i=1..nops(L2))}={true} then L1:=L2: fi: i:=MakomRishon(L1): if L1[i]<0 then L1:=[seq(-L1[j],j=1..nops(L1))]: fi: [op(i..nops(L1),L1),op(1..i-1,L1)]: end: #MaxTraj1(F,a,b,a0,b0,K): the max. magnitude of elements #of the trajectory of length K starting #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #MaxTraj1([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100); MaxTraj1:=proc(F,a,b,a0,b0,K) local gu,i: gu:=Traj(F,a,b,a0,b0,K): max(seq(abs(gu[i]),i=3..nops(gu))): end: #MaxTrajs(F,a,b,M,K): all the triples [abs(a0),abs(b0),MaxTraj1(F,a,b,a0,b0,K)] #where the last component is bigger than abs(a0)+abs(b0) #for a0 and b0 relatively prime between -M and M #For example, try: #MaxTrajs([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,5,100); MaxTrajs:=proc(F,a,b,M,K) local gu,a0,b0,c0,ka: gu:={}: for a0 from -M to M do for b0 from -M to M do if gcd(a0,b0)=1 then ka:=abs(a0)+abs(b0): c0:=MaxTraj1(F,a,b,a0,b0,K): if c0>ka then gu:=gu union {[abs(a0),abs(b0),c0]}: fi: fi: od: od: gu: end: #Persa0(F,a,b,a0,K,M): the ultimate period #if it is found, or FAIL, trying the trajectory #up to length K starting with a0 and initial conditions #x[2]=b0, between -M and M #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #Persa0([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,1,5); Persa0:=proc(F,a,b,a0,K,M) local gu,b0: if not CheckRule(F,a,b) then RETURN(FAIL): fi: gu:={}: for b0 from -M to M do if gcd(a0,b0)=1 then gu:=gu union {Canon(Per(F,a,b,a0,b0,K))}: fi: od: gu: end: #MaxTrajsa0(F,a,b,a0,M,K): #all the triples [abs(a0),abs(b0),MaxTraj1(F,a,b,a0,b0,K)] #where the last component is bigger than abs(a0)+abs(b0) #for a0 given as (fourth) input #and b0 relatively prime to a0 between -M and M #For example, try: #MaxTrajsa0([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,1,5,100); MaxTrajsa0:=proc(F,a,b,a0,M,K) local gu,b0,c0,ka,champ,shi: gu:={}: shi:=1: champ:=-M: for b0 from -M to M do if gcd(a0,b0)=1 then ka:=abs(a0)+abs(b0): c0:=MaxTraj1(F,a,b,a0,b0,K): if c0>ka then gu:=gu union {[abs(a0),abs(b0),c0]}: if c0/abs(b0)>shi then shi:=c0/abs(b0): champ:=b0: fi: fi: fi: od: gu,evalf(shi),champ: end: #MakeTrans(C1,C2,a,b): given two lists of lists of size m, say, #makes the transformation with these coefficients #For example, try: MakeTrans([[1,2],[2,1]],[[1,-2],[-2,1]],a,b); MakeTrans:=proc(C1,C2,a,b) local gu,gu1,m,i,j: m:=nops(C1): if nops(C2)<>m then RETURN(FAIL): fi: if {seq(nops(C1[i]),i=1..nops(C1))}<>{m} or {seq(nops(C2[i]),i=1..nops(C2))}<>{m} then RETURN(FAIL): fi: gu:=[]: for i from 1 to m do gu1:=[]: for j from 1 to m do gu1:=[op(gu1),normal((C1[i][j]*a+C2[i][j]*b)/m)]: od: gu:=[op(gu),gu1]: od: gu: end: #FindMates(i,j,m,c): Given an integer m, and integers i,j #between 1 and m inclusive, finds all the d's #such that i*c+j*d==0 (mod m) #For example, try; #FindMates(1,1,2,1); FindMates:=proc(i,j,m,c) local d,gu: gu:={}: for d from -m to m do if i*c+j*d mod m=0 then gu:=gu union {d}: fi: od: gu: end: #MakeTransR(C1,a,b): given a list of lists of size m, say, #makes a random transformation with C1 describing the #coefficients of x[n-1]. #For example, try: MakeTransR([[1,2],[2,1]],a,b); MakeTransR:=proc(C1,a,b) local gu,gu1,m,i,j,lu,d: m:=nops(C1): if {seq(nops(C1[i]),i=1..nops(C1))}<>{m} then RETURN(FAIL): fi: gu:=[]: for i from 1 to m do gu1:=[]: for j from 1 to m do lu:=FindMates(i,j,m,C1[i][j]): if lu={} then print([i,j,m,C1[i][j]]): RETURN(FAIL): fi: d:=lu[rand(1..nops(lu))()]; gu1:=[op(gu1),d]: od: gu:=[op(gu),gu1]: od: MakeTrans(C1,gu,a,b): end: Tsamek1:=proc(L1) local m,i: if L1=FAIL then RETURN(FAIL): fi: m:=igcd(op(L1)): if m<>0 then [seq(L1[i]/m,i=1..nops(L1))]: else L1: fi: end: Tsamek:=proc(S) local s: {seq(Tsamek1(s), s in S)}: end: #EfoBa(F,a,b,i,j): if you apply rule F #(phrased in terms of variables a and b) to inputs #such that a=i mod m and b=j mod m, what can #the mod of the output be?. For example, try: #EfoBa(F2,a,b,1,1); EfoBa:=proc(F,a,b,i,j) local m,A,B,i1,gu: m:=nops(F): gu:=expand(subs({a=m*A+i,b=m*B+j},F[i][j])) mod m: if type(gu,integer) then if gu=0 then gu:=m: fi: RETURN({gu}): else RETURN({seq(i1,i1=1..m)}): fi: end: #DG1(F,a,b): The directed graph corresponding to the second-order #recurrence given by rule F. For example, try: #DG1(F2,a,b) DG1:=proc(F,a,b) local m,gu,i,j,s,lu: option remember: m:=nops(F): gu:={}: for i from 1 to m do for j from 1 to m do lu:=EfoBa(F,a,b,i,j): gu:=gu union {[[i,j],{seq([j,s],s in lu)} ]}: od: od: gu: end: #Novea2(L1,L2,L3,a,b): Given Linear forms L1, L2,L3 in #(a,b) expresses L3 as a linear combination of L1 and #L2 and decides whether abs(L1)<=A, abs(L2)<=A implies #abs(L3)<=A. For example, try #Novea2(a,b,(a+b)/2,a,b); Novea2:=proc(L1,L2,L3,a,b) local c1,c2,gu,var,eq: var:={c1,c2}: gu:=expand(L3-c1*L1-c2*L2): eq:={coeff(gu,a,1),coeff(gu,b,1)}: var:=solve(eq,var): if var=NULL then RETURN(false): fi: c1:=subs(var,c1): c2:=subs(var,c2): if not (type(c1,numeric) and type(c2,numeric)) then RETURN(false): fi: if abs(c1)+abs(c2)<=1 then RETURN(true): else RETURN(false): fi: end: #Novea1(L1,L2): Given Linear forms L1,L2 in #is L2=+-L1. For example, try: #Novea1((a+b)/2,-(a+b)/2); Novea1:=proc(L1,L2) local gu: gu:=normal(L1/L2): if abs(gu)=1 then true: else false: fi: end: #Novea(S,L,a,b): Given a set of linear forms S in (a,b) #and a single linear form L in (a,b), decides whether #|L|<=A (for some A) is implied by the inequalities #|S1|<=A, |S2|<=A, for members of S, taken one-at-a-time #or two-at-a-time. For example,try: #Novea({a,b,a-b},b-(9/10)*a,a,b); Novea:=proc(S,L,a,b) local i,j: for i from 1 to nops(S) do if Novea1(S[i],L) then RETURN(true): fi: od: for i from 1 to nops(S) do for j from i+1 to nops(S) do if Novea2(S[i],S[j],L,a,b) then RETURN(true): fi: od: od: false: end: #SchF1:= #{[[1, 1], {[1, 1], [1, 2]}, {a,b}], #[[1, 2], {[2, 1]},{a,b,a-b,a-2*b}], #[[2, 1], {[1, 1]},{a,b,a-b}], #[[2, 2], {[2, 1], [2, 2]},{a,b} ]}: #VerifyScheme1(Sch,F,a,b): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #checks whether it is a valid scheme. For example, try #VerifyScheme1(SchF1,F1,a,b); VerifyScheme1:=proc(Sch,F,a,b) local i,s,v,N,IE,n,Ta,IE1,lu: for i from 1 to nops(Sch) do Ta[Sch[i][1]]:=Sch[i][3]: od: for i from 1 to nops(Sch) do lu:=Sch[i]: v:=lu[1]: N:=lu[2]: IE:=lu[3]: for n in N do IE1:=Ta[n]: IE1:=subs({a=b,b=F[v[1]][v[2]]},IE1): if not {seq(Novea(IE,s,a,b), s in IE1)}={true} then RETURN(false): fi: od: od: true: end: #VerifyScheme(Sch,F,a,b): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #checks whether it is a valid scheme. For example, try #VerifyScheme(SchF1,F1,a,b); VerifyScheme:=proc(Sch,F,a,b): VerifyScheme1(Sch[1],F,a,b): end: #Sch:= #{[[1, 1], {[1, 1], [1, 2]}, {a,b}], #[[1, 2], {[2, 1]},{a,b}], #[[2, 1], {[1, 1]},{a,b}], #[[2, 2], {[2, 1], [2, 2]},{a,b} ]}: #ImproveScheme(Sch,F,a,b): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #improves it by adding one more condition #For example, try #ImproveScheme(SchF1V,F1,a,b); ImproveScheme:=proc(Sch,F,a,b) local i,s,v,N,IE,n,Ta,IE1,lu,Tn,i1: for i from 1 to nops(Sch) do Ta[Sch[i][1]]:=Sch[i][3]: Tn[Sch[i][1]]:=Sch[i][2]: od: for i from 1 to nops(Sch) do lu:=Sch[i]: v:=lu[1]: N:=lu[2]: IE:=lu[3]: for n in N do IE1:=Ta[n]: IE1:=subs({a=b,b=F[v[1]][v[2]]},IE1): for s in IE1 do if not Novea(IE,s,a,b) then IE:=IE union {s}: Ta[lu[1]]:= Ta[lu[1]] union {s}: RETURN( { seq([Sch[i1][1],Tn[Sch[i1][1]],Ta[Sch[i1][1]] ],i1=1..nops(Sch)) } ): fi: od: od: od: Sch: end: #MakeScheme(F,a,b,K): Given a second-order recurrence #x[n]=F(x[n-1],x[n-2]) given in terms of rule F #phrased in terms of a and b, tries to #make a scheme with up-to K improvements. #If it succeeds, it returns true, followed by #the successful scheme, if it fails, it returns #false followed by the unsuccessful last attempt #(with K additions) #For example, try: #MakeScheme(F1,a,b,10): MakeSchemeOld:=proc(F,a,b,K) local gu,Sch,Sch1,kha,i: gu:=DG1(F,a,b): Sch:={seq([op(gu[i]),{a,b}],i=1..nops(gu))}: gu:=VerifyScheme1(Sch,F,a,b): if gu then RETURN(Sch): fi: Sch1:=ImproveScheme(Sch,F,a,b): for i from 1 to K while Sch<>Sch1 do kha:=ImproveScheme(Sch1,F,a,b): Sch:=Sch1: Sch1:=kha: od: if VerifyScheme1(Sch1,F,a,b) then RETURN([Sch1,KavuaScheme(Sch1,a,b)]): else RETURN(FAIL): fi: end: #SymbTraj1(F,a,b,v): Given a rule F phrased in terms of a and b #and a word v in the alphabet {1,...,m} (where m:=nops(F)) #finds the trajectory of length nops(v)+1, starting with #a and b and assuming that the i-th term is v[i] mod m. #For example, try: #SymbTraj1(F1,a,b,[1,2,1,2]); SymbTraj1:=proc(F,a,b,v) local i,gu: if nops(v)<2 then RETURN(FAIL): fi: gu:=[a,b]: for i from 1 to nops(v)-1 do gu:=[op(gu),subs({a=gu[nops(gu)-1],b=gu[nops(gu)]},F[v[i]][v[i+1]])]: od: gu: end: #Starters(F,a,b): Given a rule F, phrased in terms of a and b #outputs all the possible symbolic starters, where the third #element is larger than the first two #(lists of length 3). #For example, try: #Starters(F1,a,b); Starters:=proc(F,a,b) local G,i,gu,V,T,lu,v,ku,k1: G:=DG1(F,a,b): V:={seq(G[i][1],i=1..nops(G))}: for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: for v in V do lu:=[[a,b,F[v[1]][v[2]]],v]: if solve({abs(lu[1][3])>abs(lu[1][2]),abs(lu[1][3])>abs(lu[1][1])},{a,b}) <>NULL then ku:=T[v]: gu:=gu union {seq([lu[1],[op(lu[2]),k1[2]]],k1 in ku)}: fi: od: gu: end: #Yeladim(F,a,b,t): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #Yeladim(F1,a,b,[[a,b,a-b],[1,2,1]]); Yeladim:=proc(F,a,b,t) local gu,st,G,a1,b1,c1,ku,k1,T,i: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): if solve({abs(c1)NULL then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #SymbTrajs(F,a,b,k): Given a rule F phrased in terms of a and b #and an integer k>=3, #finds all good trajectories of length k, #(where the third element is the largest in absolute value) #starting with a and b #For example, try: #SymbTrajs(F1,a,b,4); SymbTrajs:=proc(F,a,b,k) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Starters(F,a,b)): fi: mu:=SymbTrajs(F,a,b,k-1): gu:={}: for m in mu do gu:=gu union Yeladim(F,a,b,m): od: gu: end: #sYeladim(F,a,b,t): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal strong extensions #by one step. For example, try: #sYeladim(F1,a,b,[[a,b,a-b],[1,2,1]]); sYeladim:=proc(F,a,b,t) local gu,st,G,a1,b1,c1,ku,k1,T,i,i1: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): if solve({abs(c1)NULL then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #sSymbTrajs(F,a,b,k): Given a rule F phrased in terms of a and b #and an integer k>=3, #finds all good trajectories (in the strong sense) of length k, #(where the third element is the largest in absolute value) #starting with a and b #For example, try: #sSymbTrajs(F1,a,b,4); sSymbTrajs:=proc(F,a,b,k) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Starters(F,a,b)): fi: mu:=sSymbTrajs(F,a,b,k-1): gu:={}: for m in mu do gu:=gu union sYeladim(F,a,b,m): od: gu: end: #IsPrim(P): Is the period P primitive? #For example try: #IsPrim([3,-1,2,3,-1,2]); IsPrim:=proc(P) local gu,k,g,i1: k:=nops(P): gu:=numtheory[divisors](k) minus {k}: for g in gu do if {seq(evalb(P[i1]=P[i1+g]),i1=1..k-g)}={true} then RETURN(false): fi: od: true: end: #Orbits(F,a,b,k): all the orbits of period k of the recurrence #x[n]=F(x[n-1],x[n-2]), where F is phrased in terms of a and b #For example, try: #Orbits(F3,a,b,3); Orbits:=proc(F,a,b,k) local gu,eq,mu,mu1,mu2,var,i,m,i1: m:=nops(F): mu:=SymbTrajs(F,a,b,k+2): gu:={}: for i from 1 to nops(mu) do mu1:=mu[i][1]: mu2:=mu[i][2]: eq:={mu1[k+1]-a,mu1[k+2]-b}: var:=solve(eq,{a,b}): if var<>NULL then mu1:=subs(var,mu1): if mu1<>[0$(k+2)] then if mu2[1]=mu2[k+1] and mu2[2]=mu2[k+2] then mu1:=[op(1..k,mu1)]: mu2:=[op(1..k,mu2)]: if not ({seq(type(mu1[i1],numeric),i1=1..nops(mu1))}={true} and {seq(type(mu1[i1],integer),i1=1..nops(mu1))}<>{true}) and IsPrim(mu1) then if not {seq(type(mu1[i1],integer),i1=1..nops(mu1))}={true} or {seq(mu1[i1]-mu2[i1] mod m,i1=1..nops(mu1))}<>{0} then gu:=gu union {[mu1,mu2]}: fi: fi: fi: fi: fi: od: gu: end: #Children(F,a,b,t): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #Children(F1,a,b,[[a,b,a-b],[1,2,1]]); Children:=proc(F,a,b,t) local st,G,a1,b1,c1,ku,k1,T,i: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : end: #SymbTrajsC(F,a,b,k): Given a rule F phrased in terms of a and b #and an integer k>=2, #finds all good trajectories of length k, #starting with a and b #For example, try: #SymbTrajsC(F1,a,b,4); SymbTrajsC:=proc(F,a,b,k) local gu,mu,m: option remember: if k<2 then RETURN(FAIL): fi: if k=2 then RETURN({[[a,b],[1,1]],[[a,b],[1,2]],[[a,b],[2,1]],[[a,b],[2,2]]}): fi: mu:=SymbTrajsC(F,a,b,k-1): gu:={}: for m in mu do gu:=gu union Children(F,a,b,m): od: gu: end: #PersMod(F,a,b,a0,b0,K,M): the ultimate period #if it is found, or FAIL, trying the mod m trajectory #up to length K starting for all initial conditions #(a0,b0) between -M and M #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #PersMod([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,100,5); PersMod:=proc(F,a,b,K,M) local gu,a0,b0: if not CheckRule(F,a,b) then RETURN(FAIL): fi: gu:={}: for a0 from -M to M do for b0 from -M to M do gu:=gu union {Canon(PerMod(F,a,b,a0,b0,K))}: od: od: gu: end: #PerMod(F,a,b,a0,b0,K): the ultimate period #if it is found for the mod m, or FAIL, trying the trajectory #up to length K starting #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #PerMod([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,2,3,100); PerMod:=proc(F,a,b,a0,b0,K) local gu,a1,b1,b2,i,lu,m,i1: if not CheckRule(F,a,b) then RETURN(FAIL): fi: m:=nops(F): gu:=[a0,b0]: a1:=a0: b1:=b0: for i from 3 to K do b2:=Hafel(F,a,b,a1,b1): a1:=b1: b1:=b2: gu:=[op(gu),b1]: lu:=Per1Mod([seq(gu[i1] mod m,i1=1..nops(gu))]): if lu<>FAIL then RETURN(lu): fi: od: FAIL: end: #Per1Mod(L): Given a list L,decides whether it ends #with a period (repeat), and returns it. #If it fails,it returns FAIL. #For example, try: #Per1([3,4,1,3,2,1,3,2]); Per1Mod:=proc(L) local i,m,j1: m:=nops(L): for i from 2 to trunc(m/6) do if nops({seq([op(m-j1*i+1..m-(j1-1)*i,L)],j1=1..5)})=1 then RETURN([op(m-i+1..m,L)]): fi: od: FAIL: end: #KavuaScheme(Sch,a,b): An upper bound for a scheme Sch #phrased in terms of a and b. For example, try: #KavuaScheme(MakeScheme(E2,a,b,10)[1],a,b); KavuaScheme:=proc(Sch,a,b) local i,j1: [ max(seq(seq(abs(coeff(Sch[i][3][j1],a,1)), j1=1..nops(Sch[i][3])),i=1..nops(Sch))), max(seq(seq(abs(coeff(Sch[i][3][j1],b,1)), j1=1..nops(Sch[i][3])),i=1..nops(Sch)))]: end: #VerifySchemeE(Sch,F,a,b,M,K): Empirically verifies the #inequalities implied by the scheme for initial values #between -M and M looking at trajectories of length K #For example, try: #VerifySchemeE(MakeScheme(E2,a,b,10),E2,a,b,10,50); VerifySchemeE:=proc(Sch1,F,a,b,M,K) local mu,A,B,a0,b0,t,i1: if Sch1=FAIL then RETURN(FAIL): fi: mu:=Sch1[2]: A:=mu[1]: B:=mu[2]: for a0 from -M to M do for b0 from -M to M do t:=Traj(F,a,b,a0,b0,K): if max(seq(abs(t[i1]),i1=1..nops(t)))>A*abs(a0)+B*abs(b0) then RETURN(false, [a0,b0]): fi: od: od: true: end: #YeladimAB(F,a,b,t,A,B): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #YeladimAB(F1,a,b,[[a,b,a-b],[1,2,1]],2,1); YeladimAB:=proc(F,a,b,t,A,B) local gu,st,G,a1,b1,c1,ku,k1,T,i,t1: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): t1:=[op(t[1]),c1]: if PtorABf(t1,a,b,A,B)<>[] then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #SymbTrajsAB(F,a,b,k,A,B): Given a rule F phrased in terms of a and b #and an integer k>=3, #finds all good trajectories of length k, #(where the third element is the largest in absolute value) #and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!) #starting with a and b #For example, try: #SymbTrajsAB(E2,a,b,4,1,1); SymbTrajsAB:=proc(F,a,b,k,A,B) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Starters(F,a,b)): fi: mu:=SymbTrajsAB(F,a,b,k-1,A,B): gu:={}: for m in mu do gu:=gu union YeladimAB(F,a,b,m,A,B): od: gu: end: #DAP(F,a,b,k,M): Given a rule F, phrased in terms #of a and b defining a recurrence #x[n]=F(x[n-1],x[n-2]), #finds a boundedness scheme (if it exists) #of depth <=k, followed by all the periodic orbits #of 3<=length <=M #followed by a proof that these are all of them. #For example, try: #DAP(E4,a,b,10,6); DAP:=proc(F,a,b,k,M) local gu,mu,S,i,Su,g1,g1a,A,B,gu1,var: option remeber: S:=MakeScheme(F,a,b,k): if S=FAIL then print(`No boundedness scheme of depth <= `, k, ` found. `) : RETURN(FAIL): fi: mu:=S[2]: S:=S[1]: A:=mu[1]: B:=mu[2]: gu:={}: for i from 3 to M do gu1:=SymbTrajsAB(F,a,b,i,A,B): Su[i]:=gu1: for g1 in gu1 do var:=solve({g1[1][1]=g1[1][i-1],g1[1][2]=g1[1][i]},{a,b}): g1a:=subs(var,g1[1]): g1a:=[op(1..i-2,g1a)]: if g1a<>[0$i] and g1[2][1]=g1[2][i-1] and g1[2][2]=g1[2][i] and IsPrim(g1a) then gu:=gu union {[ g1a,[op(1..i-2,g1[2])] ]}: fi: od: od: gu,[seq(Squeeze1(Su[i],nops(F)),i=3..M)]: end: #DAPsurv(F,a,b,k,M,nosaf): All the trajectories' profiles #with the DAP (see entry) convention that survive nosaf number of steps. #For example, try: #DAPsurv(E4,a,b,10,6,2); DAPsurv:=proc(F,a,b,k,M,nosaf) local gu, mu,lu,lu1,i: option remember: mu:=DAP(F,a,b,k,M+nosaf)[2]: gu:=[]: for i from nosaf+3 to nosaf+M do lu:=mu[i-2]: lu:={seq(lu1[2],lu1 in lu)}: lu:={seq([op(1..nops(lu1)-nosaf,lu1)],lu1 in lu)}: gu:=[op(gu),lu]: od: gu: end: #PtorAB(L,a,b,A,B): Given a list of linear expressions in #a and b, finds the solution of #abs(L[1])4 and L[nops(L)-1]=L[3] then kha:=kha union {abs(L[nops(L)])FAIL then RETURN([op(1..lu[2]-1,gu)]): fi: od: FAIL: end: #DressUp(SymT,a,b,m,K): Given a symbolic trajectory #SymT in terms of a and b followed by the mod m description #find all numerical manifestations that follow the templace #with the last entries less than K. #For example, try: #DressUp([[a,b,a-b],[1,2,1]],a,b,2,4); DressUp:=proc(SymT,a,b,m,K) local mu1,mu2,i,j,mu10,i1,gu: mu1:=SymT[1]: mu2:=SymT[2]: gu:={}: for i from -K to K do for j from -K to K do mu10:=subs({a=i,b=j},mu1): if {seq(type(mu10[i1],integer),i1=1..nops(mu10))}={true} then if {seq(evalb((mu10[i1]-mu2[i1]) mod m =0 ),i1=1..nops(mu2))}={true} and max(seq(abs(mu10[i1]),i1=1..nops(mu10)))=abs(mu10[3]) then gu:=gu union {mu10}: fi: fi: od: od: gu: end: #MakeSchemeVerbose(F,a,b,K): Verbose #version of MakeScheme(F,a,b,K) (see entry) #For example, try: #MakeSchemeVerbose(E2,a,b,10): MakeSchemeVerbose:=proc(F,a,b,K) local T,mu,m,x,n,i,j,S,s,A,B,lu,lu1,S1,s1: m:=nops(F): S:=MakeScheme(F,a,b,K): S1:={seq(S[1][i][1],i=1..nops(S[1]))}: if S=FAIL then RETURN(FAIL): fi: print(`Consider the following recurrence epressing x[n+1]` ): print(`in terms of x[n-1] and x[n] according to their mod `, m): print(`behavior, as follows`): for s1 in S1 do i:=s1[1]: j:=s1[2]: print(`if x[n-1], x[n] are, mod`, m, [i mod m, j mod m], ` respectively `): print(`then `, x[n+1]=subs({a=x[n-1],b=x[n]},F[i][j]) ): print(): od: mu:=S[2]: A:=mu[1]: B:=mu[2]: S:=S[1]: for s in S do T[s[1]]:=s[3]: od: print(`The following inequalities hold (according to the cases)`): for s1 in S1 do i:=s1[1]: j:=s1[2]: lu:=T[[i,j]]: print(`if x[n-1], x[n] are, mod`, m, [i mod m, j mod m], ` respectively `): print(` then `): for lu1 in lu do print(abs(subs({a=x[n-1],b=x[n]},lu1)) <= A*abs(x[-1])+B*abs(x[0]) ): od: od: end: #ParseWord1(w): Given a word w, represents it as #w1^a*w2^b.... For example, try #ParseWord1([1,2,1,1,2,1,1,2,1,1,1,1,1,1]); ParseWord1:=proc(w) local j,v,w1,lu,i1,j1: v:=Per1(w): if w=[] then RETURN([]): fi: if v=FAIL then RETURN(FAIL): fi: if nops(v)=2 and v[1]=v[2] then v:=[v[1]]: fi: for j from 2 to trunc(nops(w)/nops(v)) while [op(nops(w)-j*nops(v)+1..nops(w),w)]= [seq(seq(v[i1],i1=1..nops(v)),j1=1..j)] do od: j:=j-1: w1:=[op(1..nops(w)-j*nops(v),w)]: lu:=ParseWord1(w1): if lu=FAIL then RETURN(FAIL): fi: [op(lu),[v,j]]: end: #ParseWord(w,t1,t2): Given a word w, represents it as #v1w1^a*w2^b....v2 For example, try #with nops(v1)<=t1 and nops(v2)<=t2. For example try: #ParseWord([1,2,1,1,2,1,1,2,1,1,1,1,1,1],0,0); ParseWord:=proc(w,t1,t2) local s1,s2,w1,ku,v1,v2: for s1 from 0 to t1 do for s2 from 0 to t2 do w1:=[op(s1+1..nops(w)-s2,w)]: v1:=[op(1..s1,w)]: v2:=[op(nops(w)-s2+1..nops(w),w)]: ku:=ParseWord1(w1): if ku<>FAIL then RETURN([v1,ku,v2]): fi: od: od: FAIL: end: #Canon1(w): the canonical form of the word w Canon1:=proc(w) local i: i:=MakomRishon(w): [op(i..nops(w),w),op(1..i-1,w)]: end: #DiscoverAtoms(F,a,b,k,M,nosaf): All the #possible atoms in DAPsurv(F,a,b,k,M,nosaf) (see entry). #For example, try: #DiscoverAtoms(E4,a,b,10,6,2); DiscoverAtoms:=proc(F,a,b,k,M,nosaf) local gu,mu,m,gadol,lu,i,g: option remember: gu:=DAP(F,a,b,k,M)[1]: gadol:=max(seq(nops(g), g in gu)): if gadol>M then print(`Make M much bigger`): RETURN(FAIL): fi: mu:=DAPsurv(F,a,b,k,M,nosaf): mu:=mu[nops(mu)]: mu:={seq(ParseWord(m,4,4),m in mu)} minus {FAIL}: mu:={seq(m[2],m in mu)}: gadol:=max(seq(nops(mu[i]),i=1..nops(mu))): if gadol>2 then print(`Not yet implemented`): RETURN(FAIL): fi: lu:={}: for i from 1 to nops(mu) do if nops(mu[i])=1 then lu:=lu union {[ Canon1(mu[i][1][1]) ]}: elif nops(mu[i])=2 then lu:=lu union {[ Canon1(mu[i][1][1]),Canon1(mu[i][2][1]) ]}: fi: od: lu: lu: end: #MitatSdom1(w,T): Given a word w and a single template T #expresses it in terms of [v1,[T,m],v2] where #v1 is a prefix and v2 is a suffix. For example, try #MitatSdom1([1,2,1,2,1,2,1,2],[2,1]); MitatSdom1:=proc(w,T) local i,v1,w1,t,j: t:=nops(T): for i from 1 to nops(w)-t while [op(i..i+t-1,w)]<>T do od: if i>nops(w)-2*t then RETURN(FAIL): fi: v1:=[op(1..i-1,w)]: w1:=[op(i..nops(w),w)]: if [op(1..t,w1)]<>T then RETURN(FAIL): fi: j:=1: w1:=[op(t+1..nops(w1),w1)]: while nops(w1)>=t and [op(1..t,w1)]=T do j:=j+1: w1:=[op(t+1..nops(w1),w1)]: od: [v1,[T,j],w1]: end: #MitatSdom2(w,T): Given a word w and a double template T=[T1,T2] #expresses it in terms of [v1,[T1,m1],[T2,m2],v2] where #v1 is a prefix and v2 is a suffix. For example, try #MitatSdom2([1,2,1,2,1,2,1,1,1,1,1,1],[[2,1],[1]]); MitatSdom2:=proc(w,T) local i,v1,w1,T1,T2,t1,j,lu: if nops(T)<>2 then RETURN(FAIL): fi: T1:=T[1]: T2:=T[2]: t1:=nops(T1): for i from 1 to nops(w)-t1 while [op(i..i+t1-1,w)]<>T1 do od: if i>nops(w)-2*t1 then RETURN(FAIL): fi: v1:=[op(1..i-1,w)]: w1:=[op(i..nops(w),w)]: if [op(1..t1,w1)]<>T1 then RETURN(FAIL): fi: j:=1: w1:=[op(t1+1..nops(w1),w1)]: while nops(w1)>=t1 and [op(1..t1,w1)]=T1 do j:=j+1: w1:=[op(t1+1..nops(w1),w1)]: od: lu:=MitatSdom1(w1,T2): if lu=FAIL then RETURN(FAIL): fi: [v1,[T1,j],op(lu)]: end: #score1(w): Given a word as a regular expression finds #its score. For example, try score1([[1,2],[[1,2],5],[[2,3,1],2],[3,1,2]]); score1:=proc(w) local lu,i: lu:=0: for i from 1 to nops(w) do if nops(w[i])=2 and type(w[i][1],list) then lu:=lu+w[i][2]*nops(w[i][1]): fi: od: lu: end: #MitatSdom(w,S): Given a word w and a set of templates #and double templates #expresses it in terms of [v1,[T1,m1],[T2,m2],v2] etc. where #v1 is a prefix and v2 is a suffix. For example, try #MitatSdom([1,2,1,2,1,2,1,1,1,1,1,1],{[[2,1],[1]]}); MitatSdom:=proc(w,S) local gu,T,gu1,aluf,si,i,aluf1: gu:={}: for T in S do if nops(T)=1 then gu1:=MitatSdom1(w,T[1]): if gu1<>FAIL then gu:=gu union {gu1}: fi: elif nops(T)=2 then gu1:=MitatSdom2(w,T): if gu1<>FAIL then gu:=gu union {gu1}: fi: fi: od: if gu={} then RETURN(FAIL): fi: aluf:=gu[1]: si:=score1(gu[1]): for i from 2 to nops(gu) do if score1(gu[i])>si then aluf:=gu[i]: si:=score1(aluf): fi: od: aluf: aluf1:=[]: for i from 1 to nops(aluf) do if aluf[i]<>[] then aluf1:=[op(aluf1),aluf[i]]: fi: od: aluf1: end: #DAPsurvY(F,a,b,k,M,nosaf): #Like DAPsurv(F,a,b,k,M,nosaf) (see entry) but the output is nice #For example, try: #DAPsurvY(E4,a,b,10,20,2); DAPsurvY:=proc(F,a,b,k,M,nosaf) local gu,lu,i,w: gu:=DAPsurv(F,a,b,k,M,nosaf): if gu=FAIL then RETURN(FAIL): fi: lu:=DiscoverAtoms(F,a,b,k,M,nosaf): if lu=FAIL or lu={} then RETURN(FAIL): fi: [seq({seq(MitatSdom(w,lu), w in gu[i])},i=1..nops(gu))]: end: #DAPsurvYf(F,a,b,k,M,nosaf): #Like DAPsurvf(F,a,b,k,M,nosaf) (see entry) but the output is nice #For example, try: #DAPsurvYf(E4,a,b,10,20,2); DAPsurvYf:=proc(F,a,b,k,M,nosaf) local gu,lu,i,w: gu:=DAPsurvf(F,a,b,k,M,nosaf): if gu=FAIL then RETURN(FAIL): fi: lu:=DiscoverAtomsf(F,a,b,k,M,nosaf): if lu=FAIL or lu={} then RETURN(FAIL): fi: [seq({seq(MitatSdom(w,lu), w in gu[i])},i=1..nops(gu))]: end: #(a and b are generic positive integers) #YeladimABp(F,a,b,t,A,B): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #YeladimABp(E3,a,b,[[a,-b,a+b],[1,2,1]],1,2); YeladimABp:=proc(F,a,b,t,A,B) local gu,st,G,a1,b1,c1,ku,k1,T,i: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): if solve ({abs(c1)<=abs(t[1][3]),abs(t[1][3])<=A*abs(b1)+B*abs(c1), a>0, b>0},{a,b}) <>NULL then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #Startersp(F,a,b): Given a rule F, phrased in terms of a and b #outputs all the possible symbolic starters, where the third #element is larger than the first two #(lists of length 3). #For example, try: #Startersp(F1,a,b); Startersp:=proc(F,a,b) local G,i,gu,V,T,lu,v,ku,k1: G:=DG1(F,a,b): V:={seq(G[i][1],i=1..nops(G))}: for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: for v in V do lu:=[[a,b,F[v[1]][v[2]]],v]: if solve({abs(lu[1][3])>abs(lu[1][2]),abs(lu[1][3])>abs(lu[1][1]), a>0,b>0},{a,b}) <>NULL then ku:=T[v]: gu:=gu union {seq([lu[1],[op(lu[2]),k1[2]]],k1 in ku)}: fi: if solve({abs(lu[1][3])>abs(lu[1][2]),abs(lu[1][3])>abs(lu[1][1]), a>0,b<0},{a,b}) <>NULL then ku:=T[v]: gu:=gu union {seq([subs(b=-b,lu[1]),[op(lu[2]),k1[2]]],k1 in ku)}: fi: od: gu: end: #SymbTrajsABp(F,a,b,k,A,B): Given a rule F phrased in terms of a and b #and an integer k>=3, #finds all good trajectories of length k, #(where the third element is the largest in absolute value) #and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!) #starting with a and b #For example, try: #SymbTrajsABp(E3,a,b,4,1,2); SymbTrajsABp:=proc(F,a,b,k,A,B) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Startersp(F,a,b)): fi: mu:=SymbTrajsABp(F,a,b,k-1,A,B): gu:={}: for m in mu do gu:=gu union YeladimABp(F,a,b,m,A,B): od: gu: end: #DAPp(F,a,b,k,M): Given a rule F, phrased in terms #of a and b defining a recurrence #x[n]=F(x[n-1],x[n-2]), #finds a boundedness scheme (if it exists) #of depth <=k, followed by all the periodic orbits #of 3<=length <=M #followed by a proof that these are all of them. #For example, try: #DAPp(E4,a,b,10,6); DAPp:=proc(F,a,b,k,M) local gu,mu,S,i,Su,g1,g1a,A,B,gu1,var: option remeber: S:=MakeScheme(F,a,b,k): if S=FAIL then print(`No boundedness scheme of depth <= `, k, ` found. `) : RETURN(FAIL): fi: mu:=S[2]: S:=S[1]: A:=mu[1]: B:=mu[2]: gu:={}: for i from 3 to M do gu1:=SymbTrajsABp(F,a,b,i,A,B): Su[i]:=gu1: for g1 in gu1 do var:=solve({g1[1][1]=g1[1][i-1],g1[1][2]=g1[1][i]},{a,b}): g1a:=subs(var,g1[1]): g1a:=[op(1..i-2,g1a)]: if g1a<>[0$i] and g1[2][1]=g1[2][i-1] and g1[2][2]=g1[2][i] and IsPrim(g1a) then gu:=gu union {[ g1a,[op(1..i-2,g1[2])] ]}: fi: od: od: gu,[seq(Squeeze1(Su[i],nops(F)),i=3..M)]: end: #DAPsurvp(F,a,b,k,M,nosaf): All the trajectories' profiles #with the DAP (see entry) convention that survive nosaf number of steps. #For example, try: #DAPsurvp(E4,a,b,10,6,2); DAPsurvp:=proc(F,a,b,k,M,nosaf) local gu, mu,lu,lu1,i: option remember: mu:=DAPp(F,a,b,k,M+nosaf)[2]: gu:=[]: for i from nosaf+3 to nosaf+M do lu:=mu[i-2]: lu:={seq(lu1[2],lu1 in lu)}: lu:={seq([op(1..nops(lu1)-nosaf,lu1)],lu1 in lu)}: gu:=[op(gu),lu]: od: gu: end: #DiscoverAtomsp(F,a,b,k,M,nosaf): All the #possible atoms in DAPsurv(F,a,b,k,M,nosaf) (see entry). #For example, try: #DiscoverAtomsp(E4,a,b,10,6,2); DiscoverAtomsp:=proc(F,a,b,k,M,nosaf) local gu,mu,m,gadol,lu,i,g: option remember: gu:=DAPp(F,a,b,k,M)[1]: gadol:=max(seq(nops(g), g in gu)): if gadol>M then print(`Make M much bigger`): RETURN(FAIL): fi: mu:=DAPsurvp(F,a,b,k,M,nosaf): mu:=mu[nops(mu)]: mu:={seq(ParseWord(m,4,4),m in mu)} minus {FAIL}: mu:={seq(m[2],m in mu)}: gadol:=max(seq(nops(mu[i]),i=1..nops(mu))): if gadol>2 then print(`Not yet implemented`): RETURN(FAIL): fi: lu:={}: for i from 1 to nops(mu) do if nops(mu[i])=1 then lu:=lu union {[ Canon1(mu[i][1][1]) ]}: elif nops(mu[i])=2 then lu:=lu union {[ Canon1(mu[i][1][1]),Canon1(mu[i][2][1]) ]}: fi: od: lu: lu: end: #DAPsurvYp(F,a,b,k,M,nosaf): #Like DAPsurv(F,a,b,k,M,nosaf) (see entry) but the output is nice #For example, try: #DAPsurvYp(E4,a,b,10,20,2); DAPsurvYp:=proc(F,a,b,k,M,nosaf) local gu,lu,i,w: gu:=DAPsurvp(F,a,b,k,M,nosaf): if gu=FAIL then RETURN(FAIL): fi: lu:=DiscoverAtomsp(F,a,b,k,M,nosaf): if lu=FAIL or lu={} then RETURN(FAIL): fi: [seq({seq(MitatSdom(w,lu), w in gu[i])},i=1..nops(gu))]: end: #BatuakhRaOld(L,a,b,A,B): Given a trajectory L #phrased in terms of symbols a and b and #pos. integers A and B, finds out whether #it disobeys the necessary conditions #abs(L[i])<=abs(L[3]) and #abs(L[3])<=A*abs(L[nops(L)-1])+B*abs(L[nops(L)]) #For example, try: #BatuakhRaOld([a,b,a+b,a/4,b/4],a,b,1,2) ; BatuakhRaOld:=proc(L,a,b,A,B) local gadol,b1,c1,c11,lu: gadol:=L[3]: b1:=L[nops(L)-1]: c1:=L[nops(L)]: if coeff(gadol,a,1)*coeff(gadol,b,1)<0 then RETURN(false): fi: if coeff(gadol,a,1)<0 or coeff(gadol,b,1)<0 then gadol:=expand(-gadol): fi: if not coeff(c1,a,1)*coeff(c1,b,1)<0 then c11:=abs(coeff(c1,a,1))*a+abs(coeff(c1,b,1))*b: if coeff(c11,a,1)>=coeff(gadol,a,1) and coeff(c11,b,1)>=coeff(gadol,b,1) and c11<>gadol then RETURN(true): fi: fi: lu:=expand(A*(abs(coeff(b1,a,1))*a+abs(coeff(b1,b,1))*b )+ B*(abs(coeff(c1,a,1))*a+abs(coeff(c1,b,1))*b )): if coeff(lu,a,1)<=coeff(gadol,a,1) and coeff(lu,b,1)<=coeff(gadol,b,1) and lu<>gadol then RETURN(true): fi: false: end: #(a and b are generic positive integers) #YeladimABfOld(F,a,b,t,A,B): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #YeladimABfOld(E3,a,b,[[a,-b,a+b],[1,2,1]],1,2); YeladimABfOld:=proc(F,a,b,t,A,B) local gu,st,G,a1,b1,c1,ku,k1,T,i: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): if not BatuakhRa([op(t[1]),c1],a,b,A,B) then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #(a and b are generic positive integers) #YeladimABfOld(F,a,b,t,A,B): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #YeladimABfOld(E3,a,b,[[a,-b,a+b],[1,2,1]],1,2); YeladimABfOld:=proc(F,a,b,t,A,B) local gu,st,G,a1,b1,c1,ku,k1,T,i: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): if not BatuakhRaOld([op(t[1]),c1],a,b,A,B) then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #SymbTrajsABf(F,a,b,k,A,B): Given a rule F phrased in terms of a and b #and an integer k>=3, #finds all good trajectories of length k, #(where the third element is the largest in absolute value) #and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!) #starting with a and b #For example, try: #SymbTrajsABf(E3,a,b,4,1,2); SymbTrajsABf:=proc(F,a,b,k,A,B) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Startersp(F,a,b)): fi: mu:=SymbTrajsABf(F,a,b,k-1,A,B): gu:={}: for m in mu do gu:=gu union YeladimABf(F,a,b,m,A,B): od: gu: end: #DAPf(F,a,b,k,M): Given a rule F, phrased in terms #of a and b defining a recurrence #x[n]=F(x[n-1],x[n-2]), #finds a boundedness scheme (if it exists) #of depth <=k, followed by all the periodic orbits #of 3<=length <=M #followed by a proof that these are all of them. #For example, try: #DAPf(E4,a,b,10,6); DAPf:=proc(F,a,b,k,M) local gu,mu,S,i,Su,g1,g1a,A,B,gu1,var: option remeber: S:=MakeScheme(F,a,b,k): if S=FAIL then print(`No boundedness scheme of depth <= `, k, ` found. `) : RETURN(FAIL): fi: mu:=S[2]: S:=S[1]: A:=mu[1]: B:=mu[2]: gu:={}: for i from 3 to M do gu1:=SymbTrajsABf(F,a,b,i,A,B): Su[i]:=gu1: for g1 in gu1 do var:=solve({g1[1][1]=g1[1][i-1],g1[1][2]=g1[1][i]},{a,b}): g1a:=subs(var,g1[1]): g1a:=[op(1..i-2,g1a)]: if g1a<>[0$i] and g1[2][1]=g1[2][i-1] and g1[2][2]=g1[2][i] and IsPrim(g1a) then gu:=gu union {[ g1a,[op(1..i-2,g1[2])] ]}: fi: od: od: gu,[seq(Squeeze1(Su[i],nops(F)),i=3..M)]: end: #DAPsurvf(F,a,b,k,M,nosaf): All the trajectories' profiles #with the DAP (see entry) convention that survive nosaf number of steps. #For example, try: #DAPsurvf(E4,a,b,10,6,2); DAPsurvf:=proc(F,a,b,k,M,nosaf) local gu, mu,lu,lu1,i: option remember: mu:=DAPf(F,a,b,k,M+nosaf): if mu=FAIL then RETURN(FAIL): fi: mu:=mu[2]: gu:=[]: for i from nosaf+3 to nosaf+M do lu:=mu[i-2]: lu:={seq(lu1[2],lu1 in lu)}: lu:={seq([op(1..nops(lu1)-nosaf,lu1)],lu1 in lu)}: gu:=[op(gu),lu]: od: gu: end: #DiscoverAtomsf(F,a,b,k,M,nosaf): All the #possible atoms in DAPsurvf(F,a,b,k,M,nosaf) (see entry). #For example, try: #DiscoverAtomsf(E4,a,b,10,6,2); DiscoverAtomsf:=proc(F,a,b,k,M,nosaf) local gu,mu,m,gadol,lu,i,g: option remember: gu:=DAPf(F,a,b,k,M)[1]: gadol:=max(seq(nops(g), g in gu)): if gadol>M then print(`Make M much bigger`): RETURN(FAIL): fi: mu:=DAPsurvf(F,a,b,k,M,nosaf): mu:=mu[nops(mu)]: mu:={seq(ParseWord(m,4,4),m in mu)} minus {FAIL}: mu:={seq(m[2],m in mu)}: gadol:=max(seq(nops(mu[i]),i=1..nops(mu))): if gadol>2 then print(`Not yet implemented`): RETURN(FAIL): fi: lu:={}: for i from 1 to nops(mu) do if nops(mu[i])=1 then lu:=lu union {[ Canon1(mu[i][1][1]) ]}: elif nops(mu[i])=2 then lu:=lu union {[ Canon1(mu[i][1][1]),Canon1(mu[i][2][1]) ]}: fi: od: lu: lu: end: #SymbTrajsABsf(F,a,b,k,A,B,M1): Given a rule F phrased in terms of a and b #and an integer k>=3, #finds all good trajectories of length k, if k<=M1 in the #slow (strict) sense, if k>M1 in the weak (fast) sense #(where the third element is the largest in absolute value) #and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!) #starting with a and b #For example, try: #SymbTrajsABsf(E3,a,b,4,1,2,10); SymbTrajsABsf:=proc(F,a,b,k,A,B,M1) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Startersp(F,a,b)): fi: if k<=M1 then RETURN(SymbTrajsAB(F,a,b,k,A,B)): fi: mu:=SymbTrajsABsf(F,a,b,k-1,A,B,M1): gu:={}: for m in mu do gu:=gu union YeladimABf(F,a,b,m,A,B): od: gu: end: #DAPsf(F,a,b,k,M1,M2): A combination of DAP and DAPf #does like DAP (i.e. slowly, but stricter) until M1 #and between M1+1 and M2, does like DAPf #DAPsf(E4,a,b,10,14,30); DAPsf:=proc(F,a,b,k,M1,M2) local gu,mu,S,i,Su,g1,g1a,A,B,gu1,var: option remeber: S:=MakeScheme(F,a,b,k): if S=FAIL then print(`No boundedness scheme of depth <= `, k, ` found. `) : RETURN(FAIL): fi: mu:=S[2]: S:=S[1]: A:=mu[1]: B:=mu[2]: gu:={}: for i from 3 to M1 do gu1:=SymbTrajsAB(F,a,b,i,A,B): Su[i]:=gu1: for g1 in gu1 do var:=solve({g1[1][1]=g1[1][i-1],g1[1][2]=g1[1][i]},{a,b}): g1a:=subs(var,g1[1]): g1a:=[op(1..i-2,g1a)]: if g1a<>[0$i] and g1[2][1]=g1[2][i-1] and g1[2][2]=g1[2][i] and IsPrim(g1a) then gu:=gu union {[ g1a,[op(1..i-2,g1[2])] ]}: fi: od: od: for i from M1+1 to M2 do gu1:=SymbTrajsABsf(F,a,b,i,A,B,M1): Su[i]:=gu1: for g1 in gu1 do var:=solve({g1[1][1]=g1[1][i-1],g1[1][2]=g1[1][i]},{a,b}): g1a:=subs(var,g1[1]): g1a:=[op(1..i-2,g1a)]: if g1a<>[0$i] and g1[2][1]=g1[2][i-1] and g1[2][2]=g1[2][i] and IsPrim(g1a) then gu:=gu union {[ g1a,[op(1..i-2,g1[2])] ]}: fi: od: od: gu,[seq(Squeeze1(Su[i],nops(F)),i=3..M2)]: end: #BlowUp1(w,m): Given a word w finds w repeated as many times #as possible, of length m BlowUp1:=proc(w,m) local k,v,i, a: k:=nops(w): v:=[]: a:=trunc(m/k): for i from 1 to a do v:=[op(v),op(w)]: od: {[op(v),op(1..m-a*k,w)]}: end: #BlowUp2(w1,w2,m): Given two words w1 and w2 finds #the set of all words of the form w1^a*w2^b #as possible up to length m BlowUp2:=proc(w1,w2,m) local k,i, a,gu: k:=nops(w1): gu:={}: a:=trunc(m/k): for i from 0 to a do gu:=gu union {[op(BlowUp1(w1,k*i)[1]),op(BlowUp1(w2,m-k*i)[1])]}: od: gu: end: #BlowUp(Lw,m): Combines BlowUp1 and BlowUp2 BlowUp:=proc(Lw,m) : if nops(Lw)<1 or nops(Lw)>2 then print(`Not yet implemented`): RETURN(FAIL): fi: if nops(Lw)=1 then RETURN(BlowUp1(Lw[1],m)): else RETURN(BlowUp2(op(Lw),m) union BlowUp1(Lw[1],m)): fi: end: #GuessST(F,a,b,k,M,nosaf): guesses symbolic trajectories #for rule F. For example,try: #GuessST(E3,a,b,10,40,4); GuessST:=proc(F,a,b,k,M,nosaf) local mu,gadol,i,m,T,t,lu,gu1,gu,akha,w, akha1,ot1: mu:=DiscoverAtomsf(F,a,b,k,M,nosaf): gadol:=max(seq(nops(mu[i]),i=1..nops(mu))): if gadol>2 then print(`Not yet implemented`): RETURN(FAIL): fi: for i from 1 to 2 do T[i]:={}: od: for m in mu do T[nops(m)]:=T[nops(m)] union {m}: od: if T[2]={} then lu:=T[1]: else T[1]:={seq(t[1],t in T[1])} minus {seq(op(t),t in T[2])}: if T[1]<>{} then RETURN(FAIL): fi: lu:=T[2]: fi: if nops(lu)>1 then print(lu): print(`Not yet implemented`): RETURN(FAIL): fi: lu:=lu[1]: gu:=DAPsurvf(F,a,b,k,M,nosaf): gu1:=[seq(gu[i] minus BlowUp(lu,i+2),i=1..nops(gu))]: lu,gu1: if gu1[nops(gu1)]={} then RETURN([lu,gu1]): fi: i:=nops(gu1): akha:=gu1[i]: ot1:={seq(w[1],w in akha)}: if nops(ot1)>1 then RETURN(FAIL): fi: ot1:=ot1[1]: akha1:={seq([op(2..nops(w),w)], w in akha)}: if akha1 minus BlowUp(lu,i+1)={} then RETURN([lu,ot1,gu1]): fi: FAIL: end: #IsFL(G,v): Given a directed graph G, and a vertex, #decides whether it is friend-less (except itself) #For example, try: #IsFL({[1,{1,2}],[2,{2}]},1); IsFL:=proc(G,v) local V,T,i,V1,v1: V:={seq(G[i][1],i=1..nops(G))}: for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: V1:=V minus {v}: for v1 in V1 do if member(v,T[v1]) then RETURN(false): fi: od: true: end: #WeedOut(G): Given a directed graph, kicks out #all friend-less vertices (and of course their outgoing edges) #For example, try: #WeedOut(DG1(E1,a,b)); WeedOut:=proc(G) local g,G1: G1:=G: for g in G do if IsFL(G,g[1]) then G1:=G1 minus {g}: fi: od: G1: end: #MakeScheme(F,a,b,K): Given a second-order recurrence #x[n]=F(x[n-1],x[n-2]) given in terms of rule F #phrased in terms of a and b, tries to #make a scheme with up-to K improvements. #If it succeeds, it returns true, followed by #the successful scheme, if it fails, it returns #false followed by the unsuccessful last attempt #(with K additions) #For example, try: #MakeScheme(F1,a,b,10): MakeScheme:=proc(F,a,b,K) local gu,Sch,Sch1,kha,i: option remember: gu:=WeedOut(DG1(F,a,b)): Sch:={seq([op(gu[i]),{a,b}],i=1..nops(gu))}: gu:=VerifyScheme1(Sch,F,a,b): if gu then RETURN(Sch): fi: Sch1:=ImproveScheme(Sch,F,a,b): for i from 1 to K while Sch<>Sch1 do kha:=ImproveScheme(Sch1,F,a,b): Sch:=Sch1: Sch1:=kha: od: if VerifyScheme1(Sch1,F,a,b) then RETURN([Sch1,KavuaScheme(Sch1,a,b)]): else RETURN(FAIL): fi: end: ###procedures on sets of words detecting regular expressions##### ezraRE:=proc(): print(`the available procedures are`): print(`WtF(w), GatherAtoms(S,k), Khaverim(w):`): print(`Grab1(w,A), Parse11(w,A), Parse1(w,A), Parse(w,A)`): print(`Hakhlil(w,A,m), Comps1(L,n), SpellOut(gw,m,n), Pars(w,A) `): print(`Kase1(S,A,m), Kase(S,A,m), Nake(gw)`): print(`FRE(LS,m) , Meyutar(gw1,S,m,K)`): print(`OneStepClean(S,m,K)`): end: #FtWold(L): Given a word in frequency notation #[[u1,a1],[u2,a2], ..., [un,an]] spells it out. #For example try: #FtWold([[[2,3],2],[[4,5,6],3]]); FtWold:=proc(L) local i,j,w: w:=[]: for i from 1 to nops(L) do for j from 1 to L[i][2] do w:=[op(w),op(L[i][1])]: od: od: w: end: #IsPP(w,k): Is the word w a perfect k-th power #For example, try: IsPP([1,2,1,2,1,2],3); IsPP:=proc(w,k) local n,a,i,j: n:=nops(w): a:=n/k: if not type(a,integer) then RETURN(false): fi: if not nops({seq([seq(w[i+a*j],i=1..a)],j=0..k-1)})=1 then RETURN(false) fi: true: end: #LargestPower(w): Given a word w, finds the largest integer #that it is a perfect power with(numtheory): LargestPower:=proc(w) local n,gu,i,ha: n:=nops(w): gu:=sort(convert(divisors(n),list)): ha:=1: for i from 1 to nops(gu) do if IsPP(w,gu[i]) then ha:=gu[i]: fi: od: ha: end: #LPsw(w): given a word w finds the subword that is the #largest power, and returns the stuff before it, #itself expressed as a power and the stuff after it. #For example try: LPsw([1,1,2,1,2,1,2,1,1,1]): LPsw:=proc(w) local aluf,si,i,j,w1,si1,aluf1: si:=1: aluf:={}: for i from 1 to nops(w) do for j from i+1 to nops(w) do w1:=[op(i..j,w)]: if LargestPower(w1)>si then si:=LargestPower(w1): aluf:={[i,j]}: elif LargestPower(w1)=si then aluf:=aluf union {[i,j]}: fi: od: od: aluf1:=aluf[1]: si1:=aluf1[2]-aluf1[1]: if nops(aluf)>1 then for i from 2 to nops(aluf) do if aluf[i][2]-aluf[i][1]>si1 then si1:=aluf[i][2]-aluf[i][1]: aluf1:=aluf[i]: fi: od: fi: aluf1,si: end: #WtF(w): converts w to frequency notation #For example, try: WtF([1,2,1,2,1,2]); WtF:=proc(w) local w1,w2,w3,aluf,si,a,i,j: if w=[] then RETURN([]): fi: if nops(w)=1 then RETURN([[w,1]]): fi: aluf:=LPsw(w): si:=aluf[2]: i:=aluf[1][1]: j:=aluf[1][2]: w1:=[op(1..i-1,w)]: w2:=[op(i..j,w)]: w3:=[op(j+1..nops(w),w)]: a:=nops(w2)/si: w2:=[op(1..a,w2)]: [op(WtF(w1)),[w2,si],op(WtF(w3))]: end: #Khaverim(w): all the cyclic mates of the word w Khaverim:=proc(w) local i: {seq([op(i..nops(w),w),op(1..i-1,w)],i=2..nops(w))}: end: #GatherAtoms(S,k): given a set of words S #gathers those that repeat at least k times. #For example, try: #GatherAtoms({[1,1,1,1,1]},3); GatherAtoms:=proc(S,k) local gu,s,lu,i,g,gu1,T: gu:={}: for s in S do lu:=WtF(s): for i from 1 to nops(lu) do if lu[i][2]>=k then gu:=gu union {lu[i][1]}: fi: od: od: gu1:={}: while gu<>{} do g:=gu[1]: gu1:=gu1 union {g}: gu:=(gu minus Khaverim(g)) minus {g} : od: for i from 1 to max(seq(nops(g), g in gu1)) do T[i]:={}: od: lu:={}: for g in gu1 do T[nops(g)]:=T[nops(g)] union {g}: lu:=lu union {nops(g)}: od: lu:=convert(lu,list): lu:=sort(lu): lu:=[seq(lu[nops(lu)-i+1],i=1..nops(lu))]: [seq(op(T[lu[i]]),i=1..nops(lu))]: end: #Grab1(w,A): Given a word w, and a set of atoms #expresses w as a regular expression in the A's #For example try: Grab1([1,2,1,2,1,2],{[1,2]}); Grab1:=proc(w,A) local i,a,k,w1: for i from 1 to nops(A) do a:=A[i]: if nops(w)>=nops(a) and [op(1..nops(a),w)]=a then k:=1: w1:=[op(nops(a)+1..nops(w),w)]: while nops(w1)>=nops(a) and [op(1..nops(a),w1)]=a do k:=k+1: w1:=[op(nops(a)+1..nops(w1),w1)]: od: RETURN([a,k],w1): fi: od: FAIL: end: #For example try: Parse1([1,2,1,2,1,2,1,1,1,1,1],[[1,2],[1]]); Parse11:=proc(w,A) local w1,gu,lu,mu: gu:=Grab1(w,A): if gu=FAIL then RETURN(FAIL): fi: lu:=gu[1]: w1:=gu[2]: if w1=[] then RETURN([lu]): else mu:=Parse11(w1,A): if mu=FAIL then RETURN(FAIL): else RETURN([lu,op(mu)]): fi: fi: end: Parse1:=proc(w,A) local j,w1,lu: lu:=Parse11(w,A): if lu<>FAIL then RETURN(lu): fi: for j from 1 to nops(w)/2 do w1:=[op(1..nops(w)-j,w)]: lu:=Parse11(w1,A): if lu<>FAIL then RETURN([op(lu),[op(nops(w)-j+1..nops(w),w)]]): fi: od: FAIL: end: Parse:=proc(w,A) local i,w1,lu: lu:=Parse1(w,A): if lu<>FAIL then RETURN(lu): fi: for i from 1 to nops(w)/2 do w1:=[op(i+1..nops(w),w)]: lu:=Parse1(w1,A): if lu<>FAIL then RETURN([[op(1..i,w)],op(lu)]): fi: od: FAIL: end: #Hakhlil(w,A,m): given a word, and a list of atoms #finds out whether it is can be generalized #and if it is returns a regular expression #in terms of m, (followed by the list of #cardinalities of the featured atoms followed #by the list of cardinalities of the intermediate things #otherwise it returns FAIL. #For example try: #Hakhlil([1,2,1,2,1,2,1,2],[[1,2]],m); Hakhlil:=proc(w,A,m) local lu,i,c,gu,ku1,ku2: lu:=Pars(w,A): if lu=FAIL then RETURN(FAIL): fi: gu:=[]: for i from 1 to nops(lu) do if nops(lu[i])=2 and type(lu[i][1],list) and not member(lu[i][1],A) then RETURN(FAIL): fi: if nops(lu[i])=2 and member(lu[i][1],A) and type(lu[i][1],list) and lu[i][2]<2 then RETURN(FAIL): fi: od: c:=0: gu:=[]: ku1:=[]: ku2:=[]: for i from 1 to nops(lu) do if type(lu[i],list) and nops(lu[i])=2 and type(lu[i][1],list) then ku1:=[op(ku1),nops(lu[i][1])]: c:=c+1: gu:=[op(gu),[lu[i][1],m[c]] ]: else gu:=[op(gu),lu[i] ]: ku2:=[op(ku2),nops(lu[i])]: fi: od: [gu,ku1,ku2]: end: #Comps1(L,n): Given a list of pos. integers L and #an integer n finds all vectors v of size #nops(L) such that L[1]*v[1]+ ...+L[nops(L)]*v[nops(L)]=n Comps1:=proc(L,n) local c,gu,gu1,g1,i: option remember: c:=nops(L): if n=0 then RETURN({[0$c]}): fi: if nops(L)=1 then if type(n/L[1],integer) then RETURN({[n/L[1]]}): else RETURN({}): fi: fi: gu:={}: for i from 0 to trunc(n/L[c]) do gu1:=Comps1([op(1..c-1,L)],n-i*L[c]): gu:=gu union {seq([op(g1),i],g1 in gu1)}: od: gu: end: #SpellOut(gw,m,n): given a generalized word [GW,ku1,ku2] #phrased in terms of the indexed variable m, and #a pos. integer n, finds all the words of length n that #fit that template SpellOut:=proc(gw,m,n) local w,ku1,n1,mu,m1,c,i1,gu: w:=gw[1]: ku1:=gw[2]: c:=nops(ku1): n1:=n-convert(gw[3],`+`): mu:=Comps1(ku1,n1): gu:={}: for m1 in mu do gu:=gu union { FtW(subs({seq(m[i1]=m1[i1],i1=1..c)},w)) }: od: gu: end: #FtW(L): Given a word in frequency notation #[[u1,a1],[u2,a2], ..., [un,an]] interspersed with #plain words, spells it out. #For example try: #FtW([ [1,3,4], [[2,3],2],[[4,5,6],3], [2,4,3]]); FtW:=proc(L) local i,j,w,lu1: w:=[]: for i from 1 to nops(L) do lu1:=L[i]: if not (nops(lu1)=2 and type(lu1[1],list) ) then w:=[op(w),op(lu1)]: else for j from 1 to L[i][2] do w:=[op(w),op(L[i][1])]: od: fi: od: w: end: #SpellOut(gw,m,n): given a generalized word [GW,ku1,ku2] #phrased in terms of the indexed variable m, and #a pos. integer n, finds all the words of length n that #fit that template SpellOut:=proc(gw,m,n) local w,ku1,n1,mu,m1,c,i1,gu: w:=gw[1]: ku1:=gw[2]: c:=nops(ku1): n1:=n-convert(gw[3],`+`): mu:=Comps1(ku1,n1): gu:={}: for m1 in mu do gu:=gu union { FtW(subs({seq(m[i1]=m1[i1],i1=1..c)},w)) }: od: gu: end: Pars:=proc(w,A) local w1,v,i,j: w1:=Parse(w,A): for i from 1 to nops(w1) while (nops(w1[i])=2 and type(w1[i][1],list) and w1[i][2]=1) do od: i:=i-1: v:=[seq(op(w1[j][1]),j=1..i)]: if v<>[] then [v,op(i+1..nops(w1),w1)]: else w1: fi: end: #Kase1(S,A,m): given a set of words S and a list of atoms A #tries to find one generalized word # covering of S by a set of regular expressions #in the atoms of A phrased in terms of the #indexed variable m. It also outputs the left-over #For example, try: #Kase1({[1,2,1,2,1,2]},[[1,2]],m); Kase1:=proc(S,A,m) local w,gw,gu: for w in S do gw:=Hakhlil(w,A,m): if gw<>FAIL then gu:=S minus SpellOut(gw,m,nops(w)): fi: RETURN(gw,gu): od: FAIL: end: #Kase(S,A,m): given a set of words S and a list of atoms A #tries to find one generalized word # covering of S by a set of regular expressions #in the atoms of A phrased in terms of the #indexed variable m. It also outputs the left-over #For example, try: #Kase({[1,2,1,2,1,2]},[[1,2]],m); Kase:=proc(S,A,m) local gu,S1,lu: gu:={}: S1:=S: lu:=Kase1(S1,A,m): while lu<>FAIL and lu[1]<>FAIL do gu:=gu union {lu[1]}: S1:=lu[2]: lu:=Kase1(S1,A,m): if S1={} then RETURN(gu,{}): fi: od: gu,S1: end: #Nake(gw): cleans up the generalized word gw Nake:=proc(gw) local a,b,b1: a:=gw[nops(gw)]: if nops(a)=2 and type(a[1],list) then RETURN(gw): fi: b:=gw[nops(gw)-1]: b1:=b[1]: if nops(a)>nops(b1) then RETURN(FAIL): fi: if [op(1..nops(a),b1)]=a then RETURN([op(1..nops(gw)-1,gw)]): else RETURN(FAIL): fi: end: NAKE:=proc(S) local s: {seq(Nake(s),s in S)}: end: #Kas(LS,A,m): given a list of sets of words S and a list of atoms A #tries to find one generalized word # covering of S by a set of regular expressions #in the atoms of A phrased in terms of the #indexed variable m. It also outputs the left-over #For example, try: #Kas([{[1,2,1,2,1,2]},{[1,2,1,2,1,2,1]}],[[1,2]],m); Kas:=proc(LS,A,m) local gu1,gu2,gu3,gu4,lu1,lu2,lu3,lu4,lu11,gu: if nops(LS)<=4 then print(`Make LS at least of length 4`): RETURN(FAIL): fi: gu1:=LS[nops(LS)]: gu2:=LS[nops(LS)-1]: gu3:=LS[nops(LS)-2]: gu4:=LS[nops(LS)-3]: lu1:=Kase(gu1,A,m): lu2:=Kase(gu2,A,m): lu3:=Kase(gu3,A,m): lu4:=Kase(gu4,A,m): if {lu1[2],lu2[2],lu3[2],lu4[2]}<>{{}} then RETURN(FAIL): fi: lu1:=lu1[1]:lu2:=lu2[1]:lu3:=lu3[1]:lu4:=lu4[1]: lu1:={seq(Nake(lu11[1]),lu11 in lu1)}: lu2:={seq(Nake(lu11[1]),lu11 in lu2)}: lu3:={seq(Nake(lu11[1]),lu11 in lu3)}: lu4:={seq(Nake(lu11[1]),lu11 in lu4)}: gu:=lu1 union lu2 union lu3 union lu4: if member(FAIL,gu) then RETURN(FAIL): else RETURN(gu): fi: end: #Halbesh(gw): dresses up the generalized word gw Halbesh:=proc(gw) local ku1,ku2,i: ku1:=[]: if not (nops(gw[1])=2 and type(gw[1][1],list)) then ku2:=[nops(gw[1])]: else ku2:=[]: fi: for i from 1 to nops(gw) do if (nops(gw[i])=2 and type(gw[i][1],list)) then ku1:=[op(ku1),nops(gw[i][1])]: fi: od: [gw,ku1,ku2]: end: #IsImplied1(gw1,gw2,m,K): is the generalized word gw1 implied #by the generalized word gw2? using words of up to length K? #For example, try: #IsImplied1([[1,2,1],[[1,2],m[1]] ],[[[1,2,1],m[1]],[[1,2],m[2]]] ,m,30); IsImplied1:=proc(gw1,gw2,m,K) local i: evalb( {seq( SpellOut(Halbesh(gw1),m,i) minus SpellOut(Halbesh(gw2),m,i), i=20..K)}={{}}): end: #Meyutar(gw1,S2,m,K): is the generalized word gw1 implied #by the rest of the generalized word gw2? using words of up to length K? #For example, try: #Meyutar([[1,2,1],[[1,2],m[1]] ],{[[[1,2,1],m[1]],[[1,2],m[2]]]} ,m,30); Meyutar:=proc(gw1,S,m,K) local S1,gw2: S1:=S minus {gw1}: for gw2 in S1 do if IsImplied1(gw1,gw2,m,K) then RETURN(true): fi: od: false: end: OneStepClean:=proc(S,m,K) local s: for s in S do if Meyutar(s,S,m,K) then RETURN(S minus {s}): fi: od: S: end: #Finds regular expressions that cover the list of sets LS #For example, try #FRE1(DAPsurvf(E2,a,b,10,32,4),m); FRE1:=proc(LS,m) local A,gu,gu1,gu2: A:=GatherAtoms(LS[nops(LS)],3): if A=[] then RETURN(FAIL): fi: gu:=Kas(LS,A,m) : if gu=FAIL then RETURN(FAIL): fi: gu1:=OneStepClean(gu,m,30): while gu1<>gu do gu2:=OneStepClean(gu1,m,30): gu:=gu1: gu1:=gu2: od: gu: end: AllPref1:=proc(w) local i: {seq([op(1..i,w)],i=0..nops(w))}:end: AllPref:=proc(S) local s: {seq(op(AllPref1(s)), s in S)}:end: #Finds regular expressions that cover the list of sets LS #For example, try #FRE(DAPsurvf(E2,a,b,10,32,4),m); FRE:=proc(LS,m) local gu,L,bita,L1,gw,i: gu:=FRE1(LS,m): if gu=FAIL then RETURN(gu): fi: bita:=max(seq(seq(nops(gw[i][1]),i=1..nops(gw)),gw in gu))+1: if bita>nops(LS) then RETURN(FAIL): fi: L:= { seq(seq(op(SpellOut(Halbesh(gw),m,i)), gw in gu),i=1..nops(LS))}: L:=AllPref(L): L1:={seq(op(LS[i]),i=1..nops(LS)-bita)}: if L1 minus L<>{} then RETURN(FAIL): else RETURN(gu): fi: end: ###End procedures on sets of words detecting regular expressions##### #GRE(F,a,b,k,M,nosaf,m): like DAPsurvf, but guesses #a set of regular expressions covering the language #of trajectories given by rule F #For example, try: #GRE(E2,a,b,10,20,7,m);`): GRE:=proc(F,a,b,k,M,nosaf,m) local gu: gu:=DAPsurvf(F,a,b,k,M,nosaf): if gu=FAIL then RETURN(FAIL): fi: FRE(gu,m): end: #Begin applying words #HafelS(F,a,b,Zug,i,j): Given a rule F (of type m=nops(F)) #where F[i][j] is the affine linear expression #in a and b that applies if a=i(mod m) and b=j(mod(m)) #applied to symbolic Zug, assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m)) # For example try, #HafelS([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,[a,b],2,1); HafelS:=proc(F,a,b,Zug,i,j) local m,i1: m:=nops(F): if not type(F,list) or {seq(nops(F[i1]),i1=1..m)}<>{m} then print(`bad input`): RETURN(FAIL): fi: [Zug[2],subs({a=Zug[1],b=Zug[2]},F[i][j])]: end: #HafelSw(F,a,b,Zug,w): Given a rule F (of type m=nops(F)) #where F[i][j] is the affine linear expression #in a and b that applies if a=i(mod m) and b=j(mod(m)) #applied to the current window #[a,b], assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m)) #keeps doing it with the parities given by the word in #{1,2, ..., m}* w. # For example try, #HafelSw([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,[a,b],[2,1,2]); HafelSw:=proc(F,a,b,Zug,w) local gu ,i: gu:=Zug: for i from 1 to nops(w)-1 do gu:=HafelS(F,a,b,gu,w[i],w[i+1]): od: gu: end: #HafelSwn1(F,a,b,Zug,w,n): Given a rule F (of type m=nops(F)) #where F[i][j] is the linear expression #in a and b that applies if a=i(mod m) and b=j(mod(m)) #applied to the current window #[a,b], assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m)) #keeps doing it with the parities given by the word in #{1,2, ..., m}*, w, applied n times and w[1], for symbolic n # For example try, #HafelSwn1([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,[a,b],[2,1,2],n); HafelSwn1:=proc(F,a,b,Zug,w,n) local gu ,A,w1,C1,C2,mua,mub,lu,mu: w1:=[op(w),w[1]]: gu:=HafelSw(F,a,b,Zug,w1): A:=[[coeff(gu[1],a,1),coeff(gu[1],b,1)], [coeff(gu[2],a,1),coeff(gu[2],b,1)]]: lu:={eigenvalues(A)}: if nops(lu)=1 then mu:=(C1+C2*n)*lu[1]^n: else mu:=C1*lu[1]^n+C2*lu[2]^n: fi: mua:=subs(solve({subs(n=0,mu)=Zug[1],subs(n=1,mu)=gu[1]},{C1,C2}),mu): mub:=subs(solve({subs(n=0,mu)=Zug[2],subs(n=1,mu)=gu[2]},{C1,C2}),mu): [mua,mub]: end: #HafelSwn(F,a,b,Zug,w,n): Given a rule F (of type m=nops(F)) #where F[i][j] is the linear expression #in a and b that applies if a=i(mod m) and b=j(mod(m)) #applied to the current window #[a,b], assuming that Zug[1]=i(mod m) and Zug[2]=j(mod(m)) #keeps doing it with the parities given by the word in #{1,2, ..., m}*, w, applied n times # For example try, #HafelSwn([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,[a,b],[2,1,2],n); HafelSwn:=proc(F,a,b,Zug,w,n) local lu: lu:=HafelSwn1(F,a,b,Zug,w,n): lu:=subs(n=n-1,lu): expand(HafelSw(F,a,b,lu,w)): end: #HafelSreOld(F,a,b,R): The regular expression applied to [a,b] #using the rule F. For example try #HafelSreOld(E2,a,b,[[1,2,1],[[1,2],n]]); HafelSreOld:=proc(F,a,b,R) local gu,w1: if nops(R)=1 and not (nops(R[1])=2 and type(R[1][1],list)) then RETURN(HafelSw(F,a,b,R[1])): fi: if not (nops(R[1])=2 and type(R[1][1],list)) then w1:=[op(R[1]),R[2][1][1]]: gu:=HafelSw(F,a,b,[a,b],w1): gu:=HafelSwn(F,a,b,gu,R[2][1],R[2][2]): fi: gu: end: #HafelSre1NoGood(F,a,b,R): The regular expression applied to [a,b] #using the rule F. For example try #HafelSre1NoGood(E2,a,b,[[1,2,1],[[1,2],n]]); HafelSre1NoGood:=proc(F,a,b,R) local gu,w1,n,C1,C2,C3,C4,A, lu, eq1,var1, eq2,var2,w,mu1,mu2: if nops(R)=1 and not (nops(R[1])=2 and type(R[1][1],list)) then RETURN(HafelSw(F,a,b,[a,b],R[1])): fi: if nops(R)=1 and nops(R[1])=2 and type(R[1][1],list) then RETURN(HafelSwn(F,a,b,[a,b],op(R[1]))): fi: if nops(R)=2 and not (nops(R[1])=2 and type(R[1][1],list)) then w:=R[2][1]: n:=R[2][2]: w1:=[op(w),w[1]]: gu:=HafelSw(F,a,b,[a,b],w1): A:=[[coeff(gu[1],a,1),coeff(gu[1],b,1)], [coeff(gu[2],a,1),coeff(gu[2],b,1)]]: lu:={eigenvalues(A)}: if nops(lu)=1 then mu1:=(C1+C2*n)*lu[1]^n: mu2:=(C3+C4*n)*lu[1]^n: else mu1:=C1*lu[1]^n+C2*lu[2]^n: mu2:=C3*lu[1]^n+C4*lu[2]^n: fi: eq1:={subs(n=1,mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w)])[1], subs(n=2,mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w)])[1]}: var1:={C1,C2}: mu1:=simplify(subs(solve(eq1,var1),mu1)): eq2:={subs(n=1,mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w)])[2], subs(n=2,mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w)])[2]}: var2:={C3,C4}: mu2:=simplify(subs(solve(eq2,var2),mu2)): RETURN([mu1,mu2]): fi: end: #HafelSreOld(F,a,b,R): The regular expression applied to [a,b] #using the rule F. For example try #HafelSre(E2,a,b,[[1,2,1],[[1,2],n]]); HafelSreOld:=proc(F,a,b,R) local gu,w1,n,C1,C2,C3,C4,A, lu, eq1,var1, eq2,var2,mu1,mu2,w,zug, C11, C12, C21, C22, D11, D12, D21, D22, Ax, gux, lux, m, x, x1: if nops(R)=1 and not (nops(R[1])=2 and type(R[1][1],list)) then RETURN(HafelSw(F,a,b,[a,b],R[1])): fi: if nops(R)=1 and nops(R[1])=2 and type(R[1][1],list) then RETURN(HafelSwn(F,a,b,[a,b],op(R[1]))): fi: if nops(R)=2 and not (nops(R[1])=2 and type(R[1][1],list)) then w:=R[2][1]: n:=R[2][2]: w1:=[op(w),w[1]]: gu:=HafelSw(F,a,b,[a,b],w1): A:=[[coeff(gu[1],a,1),coeff(gu[1],b,1)], [coeff(gu[2],a,1),coeff(gu[2],b,1)]]: lu:={eigenvalues(A)}: if nops(lu)=1 then mu1:=(C1+C2*n)*lu[1]^n: mu2:=(C3+C4*n)*lu[1]^n: else mu1:=C1*lu[1]^n+C2*lu[2]^n: mu2:=C3*lu[1]^n+C4*lu[2]^n: fi: eq1:={subs(n=1,mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w)])[1], subs(n=2,mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w)])[1]}: var1:={C1,C2}: mu1:=simplify(subs(solve(eq1,var1),mu1)): eq2:={subs(n=1,mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w)])[2], subs(n=2,mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w)])[2]}: var2:={C3,C4}: mu2:=simplify(subs(solve(eq2,var2),mu2)): RETURN([mu1,mu2]): fi: if nops(R)=2 and nops(R[1])=2 and type(R[1][1],list) and nops(R[2])=2 and type(R[2][1],list) then RETURN(HafelSre(F,a,b,[[],op(R)])): fi: if nops(R)=3 and not (nops(R[1])=2 and type(R[1][1],list)) and (nops(R[2])=2 and type(R[2][1],list)) and not (nops(R[1])=2 and type(R[1][1],list)) then zug:=HafelSre(F,a,b,[op(1..2,R)]): RETURN(HafelSw(F,a,b,zug,[R[2][1][nops(R[2][1])], op(R[3])])): fi: if nops(R)=3 and not (nops(R[1])=2 and type(R[1][1],list)) then w:=R[2][1]: n:=R[2][2]: x:=R[3][1]: m:=R[3][2]: w1:=[op(w),w[1]]: x1:=[op(x),x[1]]: gu:=HafelSw(F,a,b,[a,b],w1): gux:=HafelSw(F,a,b,[a,b],x1): A:=[[coeff(gu[1],a,1),coeff(gu[1],b,1)], [coeff(gu[2],a,1),coeff(gu[2],b,1)]]: Ax:=[[coeff(gux[1],a,1),coeff(gux[1],b,1)], [coeff(gux[2],a,1),coeff(gux[2],b,1)]]: lu:={eigenvalues(A)}: lux:={eigenvalues(Ax)}: if nops(lu)=2 and nops(lux)=2 then mu1:=C11*lu[1]^n*lux[1]^m+C12*lu[1]^n*lux[2]^m+ C21*lu[2]^n*lux[1]^m+C22*lu[2]^n*lux[2]^m: mu2:=D11*lu[1]^n*lux[1]^m+D12*lu[1]^n*lux[2]^m+ D21*lu[2]^n*lux[1]^m+D22*lu[2]^n*lux[2]^m: elif nops(lu)=1 and nops(lux)=2 then mu1:=C11*lu[1]^n*lux[1]^m+C12*lu[1]^n*lux[2]^m+ C21*lu[1]^n*n*lux[1]^m+C22*lu[1]^n*n*lux[2]^m: mu2:=D11*lu[1]^n*lux[1]^m+D12*lu[1]^n*lux[2]^m+ D21*lu[1]^n*n*lux[1]^m+D22*lu[1]^n*n*lux[2]^m: elif nops(lu)=2 and nops(lux)=1 then mu1:=C11*lu[1]^n*lux[1]^m+C12*lu[1]^n*lux[1]^m*m+ C21*lu[2]^n*lux[1]^m+C22*lu[2]^n*lux[1]^m*m: mu2:=D11*lu[1]^n*lux[1]^m+D12*lu[1]^n*lux[1]^m*m+ D21*lu[2]^n*lux[1]^m+D22*lu[2]^n*lux[1]^m*m: elif nops(lu)=1 and nops(lux)=1 then mu1:=C11*lu[1]^n*lux[1]^m+C12*lu[1]^n*lux[1]^m*m+ C21*lu[1]^n*ln*lux[1]^m+C22*lu[1]^n*lux[2]^m*n*m: mu2:=D11*lu[1]^n*lux[1]^m+D12*lu[1]^n*lux[1]^m*m+ D21*lu[1]^n*n*lux[1]^m+D22*lu[1]^n*lux[1]^m*n*m: else RETURN(FAIL): fi: eq1:={ subs({n=1,m=1},mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(x)])[1], subs({n=1,m=2},mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(x),op(x)])[1], subs({n=2,m=1},mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w),op(x)])[1], subs({n=2,m=2},mu1)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w),op(x),op(x)]) [1]}: var1:={C11,C12,C21,C22}: mu1:=simplify(subs(solve(eq1,var1),mu1)): eq1:={ subs({n=1,m=1},mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(x)])[2], subs({n=1,m=2},mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(x),op(x)])[2], subs({n=2,m=1},mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w),op(x)])[2], subs({n=2,m=2},mu2)=HafelSw(F,a,b,[a,b],[op(R[1]),op(w),op(w),op(x),op(x)]) [2]}: var1:={D11,D12,D21,D22}: mu2:=simplify(subs(solve(eq1,var1),mu2)): RETURN([mu1,mu2]): fi: if nops(R)=3 and nops(R[1])=2 and type(R[1][1],list) and nops(R[2])=2 and type(R[2][1],list) and not (nops(R[3])=2 and type(R[3][1],list)) then zug:=HafelSre(F,a,b,[op(1..2,R)]): RETURN(HafelSw(F,a,b,zug,[R[2][1][nops(R[2][1])], op(R[3])])): fi: if nops(R)=2 and nops(R[1])=2 and type(R[1][1],list) and not (nops(R[2])=2 and type(R[2][1],list)) then zug:=HafelSre(F,a,b,[R[1]]): RETURN(HafelSw(F,a,b,zug,[R[1][1][nops(R[1][1])], op(R[2])])): fi: print(`Not yet implemented`): RETURN(FAIL): end: #BatuakhRa(gadol,b1,c1,a,b,A,B): Given three #linear expression, gadol, b1,c1, phrased #phrased in terms of symbols a and b and #pos. integers A and B, finds out whether #it disobeys the necessary conditions #abs(b1)<=abs(gadol) , abs(c1)<=abs(gadol) #abs(gadol)<=A*abs(b1)+B*abs(c1) #it requires that the coeff. of a and b of gadol #have the same sign. #For example, try: #BatuakhRa(a+b,a/4,b/4,a,b,1,2) ; BatuakhRa:=proc(Gadol,b1,c1,a,b,A,B) local c11,lu,gadol: gadol:=Gadol: if coeff(gadol,a,1)*coeff(gadol,b,1)<0 then RETURN(false): fi: if coeff(gadol,a,1)<0 or coeff(gadol,b,1)<0 then gadol:=expand(-gadol): fi: if not coeff(c1,a,1)*coeff(c1,b,1)<0 then c11:=abs(coeff(c1,a,1))*a+abs(coeff(c1,b,1))*b: if coeff(c11,a,1)>=coeff(gadol,a,1) and coeff(c11,b,1)>=coeff(gadol,b,1) and c11<>gadol then RETURN(true): fi: fi: lu:=expand(A*(abs(coeff(b1,a,1))*a+abs(coeff(b1,b,1))*b )+ B*(abs(coeff(c1,a,1))*a+abs(coeff(c1,b,1))*b )): if coeff(lu,a,1)<=coeff(gadol,a,1) and coeff(lu,b,1)<=coeff(gadol,b,1) and lu<>gadol then RETURN(true): fi: false: end: #(a and b are generic positive integers) #YeladimABf(F,a,b,t,A,B): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #YeladimABfOld(E3,a,b,[[a,-b,a+b],[1,2,1]],1,2); YeladimABf:=proc(F,a,b,t,A,B) local gu,st,G,a1,b1,c1,ku,k1,T,i: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): if not BatuakhRa(t[1][3],b1,c1,a,b,A,B) then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #BatuakhRaSymb(gadol,b1,c1,a,b,A,B,m,r,K1,K2): Given #a linear expression, gadol, in terms of symbols a and b # and linear expressions b1,c1, phrased in terms of a and b #and symbolic integers m[1],m[2], ..., m[r] checks whether #for the specializations m[1], ... m[r] between K1 and K2 #pos. integers A and B, finds out whether #it disobeys the necessary conditions #abs(b1)<=abs(gadol) , abs(c1)<=abs(gadol) #abs(gadol)<=A*abs(b1)+B*abs(c1) #it requires that the coeff. of a and b of gadol #have the same sign. #For example, try: #BatuakhRaSymb(a+b,a/2^m[1],b/2^m[1],a,b,1,2,m,1,2,3) ; BatuakhRaSymb:=proc(Gadol,b1,c1,a,b,A,B,m,r,K1,K2) local i,j,b1a,c1a: if r>2 then print(`Not yet implemented`): fi: if r=1 then for i from K1 to K2 do b1a:=subs({m[1]=i},b1): c1a:=subs({m[1]=i},c1): if not BatuakhRa(Gadol,b1a,c1a,a,b,A,B) then RETURN(false): fi: od: RETURN(true): fi: if r=2 then for i from K1 to K2 do for j from K1 to K2 do b1a:=subs({m[1]=i,m[2]=j},b1): c1a:=subs({m[1]=i,m[2]=j},c1): if not BatuakhRa(Gadol,b1a,c1a,a,b,A,B) then RETURN(false): fi: od: od: RETURN(true): fi: end: #YeladimRE(F,a,b,Ini,A,B,R,m,r,w,K1,K2): Given #a rule F phrased in terms of a and b and #a regular expression R phrased in terms of m #and a word w finds all the words of length one more #that continues [op(R),op(w)] and do not violate the necessary conditions #as checked by BatuakhRaSymb (with parameters r, K1, K2). #For example, try: #YeladimRE(E7,a,b,[a,-b,a+b],1,2,[[[2,1,1,1],m[1]]],m,1,[],1,10); YeladimRE:=proc(F,a,b,Ini,A,B,R,m,r,w,K1,K2) local k,gu,lu,w1,akha,i, lu1,Gadol: if nops(Ini)<>3 then print(`Fourth argument should be a list of length 3`): RETURN(FAIL): fi: Gadol:=Ini[3]: akha:=R[nops(R)][1]: akha:=akha[nops(akha)]: k:=nops(F): lu:=HafelSre(F,a,b,R): lu:=subs({a=Ini[1],b=Ini[2]},lu): gu:={}: for i from 1 to k do w1:=[akha,op(w),i]: lu1:=HafelSw(F,a,b,lu,w1): if not BatuakhRaSymb(Gadol,op(lu1),a,b,A,B,m,r,K1,K2) then gu:=gu union {[op(w),i]}: fi: od: gu: end: #SpellOutE(gw,m,n1): spells out the generalized word gw #express version. For example, try: #SpellOutE([ [[2,1],m[1]], [[1,2,2],m[2]] ],m,10); SpellOutE:=proc(gw,m,n1) local ku,i,w1,gw1,lu,lu1: if nops(gw[1])=2 and type(gw[1][1],list) then ku:=[seq(nops(gw[i][1]),i=1..nops(gw))]: RETURN(SpellOut([gw,ku,[]],m,n1)): else w1:=gw[1]: gw1:=[op(2..nops(gw),gw)]: lu:=SpellOutE(gw1,m,n1-nops(w1)): RETURN({seq([op(w1),op(lu1)], lu1 in lu)}): fi: end: #DAPfA(F,a,b,k,M): a list of all the set of legal words for rule #F of length up to M, starting at length 3 DAPfA:=proc(F,a,b,k,M) local gu,i,j,w: gu:=DAPf(F,a,b,k,M+1): gu:=[seq({seq(gu[2][i][j][2],j=1..nops(gu[2][i]))},i=1..nops(gu[2]))]: [seq({seq([op(1..nops(w)-1,w)],w in gu[i+1])},i=1..nops(gu)-1)]: end: #gDAPfA(F,a,b,S,M): a list of all the set of legal words for rule #F of length up to M, starting at length 3 gDAPfA:=proc(F,a,b,S,M) local gu,i,j,w: gu:=gDAPf(F,a,b,S,M+1): gu:=[seq({seq(gu[2][i][j][2],j=1..nops(gu[2][i]))},i=1..nops(gu[2]))]: [seq({seq([op(1..nops(w)-1,w)],w in gu[i+1])},i=1..nops(gu)-1)]: end: #IsCont(w,gw,m,t): is the concrete word w a continuation #of a specialization of the generalized word m #phrased in terms of the index variable m #output is the smallest possible suffix of length <=t #if it fails it returms FAIL #For example try: IsCont([1,2,1,2,1,2,1,2,2,2],[[[1,2],m[1]]],m,4); IsCont:=proc(w,gw,m,t) local w1,j: for j from 0 to min(t,nops(w)-1) do w1:=[op(1..nops(w)-j,w)]: if member(w1,SpellOutE(gw,m,nops(w1))) then RETURN([op(nops(w)-j+1..nops(w),w)]): fi: od: FAIL: end: #GREplus(F,a,b,k,M,nosaf,m): like GRE, but #also with suffixes #GREplus(E2,a,b,10,20,7,m);`): GREplus:=proc(F,a,b,k,M,nosaf,m) local mu,lu,gu,s,S,lu1,ku,i,mu1: lu:=GRE(F,a,b,k,M,nosaf,m): if lu=FAIL then RETURN(FAIL): fi: mu:=DAPfA(F,a,b,k,M): mu1:=[op(1..2,DAPf(F,a,b,k,3)[2][1][1][1])]: gu:={}: for i from max(M-3*nosaf,10) to nops(mu) do S:=mu[i]: for s in S do for lu1 in lu do ku:=IsCont(s,lu1,m,nosaf+3): if ku<>FAIL then if ku<>[] then gu:=gu union {[op(lu1),ku]}: else gu:=gu union {lu1}: fi: fi: od: od: od: gu,mu1: end: #HopeForOrbits(F,a,b,k,M,nosaf,m): looks at the #condition, phrased in terms of m[1],m[2]'s etc. #that there is a non-trivial orbit. For example try: #HopeForOrbits(E2,a,b,10,20,7,m); HopeForOrbits:=proc(F,a,b,k,M,nosaf,m) local gu,mu1,g,eq,ku,S,lu: gu:=GREplus(F,a,b,k,M,nosaf,m): mu1:=gu[2]: gu:=gu[1]: S:={}: for g in gu do lu:=HafelSre(F,a,b,g): eq:=expand([lu[1]-mu1[1],lu[2]-mu1[2]]): ku:=[[coeff(eq[1],a,1),coeff(eq[1],b,1)],[coeff(eq[2],a,1),coeff(eq[2],b,1)]]: ku:=expand(ku[1][1]*ku[2][2]-ku[1][2]*ku[2][1]): S:=S union {[g,ku]}: od: S: end: #####Empirically finding a scheme #TestABij(F,a,b,A,B,i,j,K,orekh): tests whether it is true for #starting values in [-K,K] that abs(x[n])<=A*abs(x[-1])+B*abs(x[-2]) #whenever x[n-1] mod m=i and x[n-1] mod m=j #up to trajectories of length orekh. For example, try: #TestABij(E2,a,b,1,1,1,1,10,30); TestABij:=proc(F,a,b,A,B,i,j,K,orekh) local a0,b0,t,i1,m: m:=nops(F): for a0 from -K to K do for b0 from -K to K do if gcd(a0,b0) mod m<>0 then t:=Traj(F,a,b,a0,b0,orekh): for i1 from 3 to nops(t)-1 do if t[i1-1]-i mod m=0 and t[i1]-j mod m=0 and not abs(t[i1])<=A*abs(t[1])+B*abs(t[2]) then RETURN(false): fi: od: fi: od: od: true: end: #FindPreScheme1(F,a,b,i,j,K,orekh,Ma): tries to find a #prescheme for the rule F #For example, try: #FindPreScheme1(E2,a,b,1,1,30,30,2); FindPreScheme1:=proc(F,a,b,i,j,K,orekh,Ma) local A,B,M: for M from 1 to Ma do for A from 1 to M-1 do B:=M-A: if TestABij(F,a,b,A,B,i,j,K,orekh) then RETURN([A,B]): fi: od: od: FAIL: end: #FindPreScheme(F,a,b,K,orekh,Ma): tries to find a #prescheme for the rule F #For example, try: #FindPreScheme(E2,a,b,30,30,2); FindPreScheme:=proc(F,a,b,K,orekh,Ma) local gu,i,j,m,gu1: m:=nops(F): gu:={}: for i from 1 to m do for j from 1 to m do if not (i=m and j=m) then gu1:=FindPreScheme1(F,a,b,i,j,K,orekh,Ma): if gu1=FAIL then RETURN(FAIL): else gu:=gu union {[[i,j],gu1]}: fi: fi: od: od: gu: end: #TestABijG(F,a,b,A,B,i,j,K,orekh,PreLim): tests whether it is true for #starting values in [-K,K] that abs(x[n])<=A*abs(x[-1])+B*abs(x[-2]) #whenever x[n-1] mod m=i and x[n-1] mod m=j #up to trajectories of length orekh, starting at PreLim. For example, try: #TestABijG(E2,a,b,1,1,1,1,100,30,10); TestABijG:=proc(F,a,b,A,B,i,j,K,orekh,PreLim) local a0,b0,t,i1,m: m:=nops(F): for a0 from -K to K do for b0 from -K to K do if gcd(a0,b0) mod m<>0 then t:=Traj(F,a,b,a0,b0,orekh): for i1 from PreLim to nops(t)-1 do if t[i1-1]-i mod m=0 and t[i1]-j mod m=0 and not abs(t[i1])<=A*abs(t[1])+B*abs(t[2]) then RETURN(false): fi: od: fi: od: od: true: end: #FindPreScheme1G(F,a,b,i,j,K,orekh,Ma,PreLim): tries to find a #prescheme for the rule F letting it run first PreLim moves #For example, try: #FindPreScheme1G(E2,a,b,1,1,30,30,2); FindPreScheme1G:=proc(F,a,b,i,j,K,orekh,Ma,PreLim) local A,B,M: for M from 1 to Ma do for A from 1 to M-1 do B:=M-A: if TestABijG(F,a,b,A,B,i,j,K,orekh,PreLim) then RETURN([A,B]): fi: od: od: FAIL: end: #FindPreSchemeG(F,a,b,K,orekh,Ma,PreLim): tries to find a #prescheme for the rule F (starting at PreLim) #For example, try: #FindPreSchemeG(E2,a,b,30,30,2,10); FindPreSchemeG:=proc(F,a,b,K,orekh,Ma,PreLim) local gu,i,j,m,gu1: m:=nops(F): gu:={}: for i from 1 to m do for j from 1 to m do if not (i=m and j=m) then gu1:=FindPreScheme1G(F,a,b,i,j,K,orekh,Ma,PreLim): if gu1=FAIL then RETURN(FAIL): else gu:=gu union {[[i,j],gu1]}: fi: fi: od: od: gu: end: ####End Empirically finding a scheme ######generalized DAPf #gSymbTrajsABf(F,a,b,k,PS): Given a rule F phrased in terms of a and b #and an integer k>=3, as well as a pre-scheme #finds all good trajectories of length k, #(where the third element is the largest in absolute value) #and |x[3]|<=A|x[n-1]|+B|x[n]| (cyclically!) #if x[n-1] mod m=i, x[n] mod m=j when [[i,j],[A,B]] is a member of PS #starting with a and b #For example, try: #gSymbTrajsABf(E2,a,b,4,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]} ); gSymbTrajsABf:=proc(F,a,b,k,PS) local gu,mu,m: option remember: if k<3 then RETURN(FAIL): fi: if k=3 then RETURN(Startersp(F,a,b)): fi: mu:=gSymbTrajsABf(F,a,b,k-1,PS): gu:={}: for m in mu do gu:=gu union gYeladimABf(F,a,b,m,PS): od: gu: end: #gYeladimABf(F,a,b,t,PS): Given a rule F, phrased in terms of a and #b, and a trajectory t, finds all its legal extensions #by one step. For example, try: #gYeladimABf(E3,a,b,[[a,-b,a+b],[1,2,1]],{[[1,1],[1,1]], [[1,2],[1,1]],[[2,1],[1,1]]} ); gYeladimABf:=proc(F,a,b,t,PS) local gu,st,G,a1,b1,c1,ku,k1,T,i,A,B,TA,ps,AB: G:=DG1(F,a,b): for i from 1 to nops(G) do T[G[i][1]]:=G[i][2]: od: for ps in PS do TA[ps[1]]:=ps[2]: od: gu:={}: st:=[t[2][nops(t[2])-1],t[2][nops(t[2])]]: ku:=T[st]: a1:=t[1][nops(t[1])-1]: b1:=t[1][nops(t[1])]: c1:=subs({a=a1,b=b1},F[st[1]][st[2]]): AB:=TA[st]: A:=AB[1]: B:=AB[2]: if not BatuakhRa(t[1][3],b1,c1,a,b,A,B) then gu:=gu union {seq([[op(t[1]),c1],[op(t[2]),k1[2]]],k1 in ku)} : fi: gu: end: #gDAPf(F,a,b,S,M): Given a rule F, phrased in terms #of a and b defining a recurrence #x[n]=F(x[n-1],x[n-2]), #and a prescheme S #of depth <=k, followed by all the periodic orbits #of 3<=length <=M #followed by a proof that these are all of them. #For example, try: #gDAPf(E4,a,b,10,6); gDAPf:=proc(F,a,b,S,M) local gu,i,Su,g1,g1a, gu1,var: option remeber: gu:={}: for i from 3 to M do gu1:=gSymbTrajsABf(F,a,b,i,S): Su[i]:=gu1: for g1 in gu1 do var:=solve({g1[1][1]=g1[1][i-1],g1[1][2]=g1[1][i]},{a,b}): g1a:=subs(var,g1[1]): g1a:=[op(1..i-2,g1a)]: if g1a<>[0$i] and g1[2][1]=g1[2][i-1] and g1[2][2]=g1[2][i] and IsPrim(g1a) then gu:=gu union {[ g1a,[op(1..i-2,g1[2])] ]}: fi: od: od: gu,[seq(Squeeze1(Su[i],nops(F)),i=3..M)]: end: #gDAPsurvf(F,a,b,S,M,nosaf): All the trajectories' profiles #with the gDAPf (see entry) convention that survive nosaf number of steps. #For example, try: #gDAPsurvf(E2,a,b,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]},6,2); gDAPsurvf:=proc(F,a,b,S,M,nosaf) local gu, mu,lu,lu1,i: option remember: mu:=gDAPf(F,a,b,S,M+nosaf): if mu=FAIL then RETURN(FAIL): fi: mu:=mu[2]: gu:=[]: for i from nosaf+3 to nosaf+M do lu:=mu[i-2]: lu:={seq(lu1[2],lu1 in lu)}: lu:={seq([op(1..nops(lu1)-nosaf,lu1)],lu1 in lu)}: gu:=[op(gu),lu]: od: gu: end: #gDiscoverAtomsf(F,a,b,S,M,nosaf): All the #possible atoms in gDAPsurvf(F,a,b,S,M,nosaf) (see entry). #For example, try: #gDiscoverAtomsf(E4,a,b, {[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]},20,2); gDiscoverAtomsf:=proc(F,a,b,S,M,nosaf) local gu,mu,m,gadol,lu,i,g: option remember: gu:=gDAPf(F,a,b,S,M)[1]: gadol:=max(seq(nops(g), g in gu)): if gadol>M then print(`Make M much bigger`): RETURN(FAIL): fi: mu:=gDAPsurvf(F,a,b,S,M,nosaf): mu:=mu[nops(mu)]: mu:={seq(ParseWord(m,4,4),m in mu)} minus {FAIL}: mu:={seq(m[2],m in mu)}: gadol:=max(seq(nops(mu[i]),i=1..nops(mu))): if gadol>2 then print(`Not yet implemented`): RETURN(FAIL): fi: lu:={}: for i from 1 to nops(mu) do if nops(mu[i])=1 then lu:=lu union {[ Canon1(mu[i][1][1]) ]}: elif nops(mu[i])=2 then lu:=lu union {[ Canon1(mu[i][1][1]),Canon1(mu[i][2][1]) ]}: fi: od: lu: lu: end: #gGRE(F,a,b,S,M,nosaf,m): like gDAPsurvf, but guesses #a set of regular expressions covering the language #of trajectories given by rule F #For example, try: #gGRE(E2,a,b, {[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]},20,7,m); gGRE:=proc(F,a,b,S,M,nosaf,m) local gu: gu:=gDAPsurvf(F,a,b,S,M,nosaf): if gu=FAIL then RETURN(FAIL): fi: FRE(gu,m): end: #gGREplus(F,a,b,S,M,nosaf,m): like GRE, but #also with suffixes #gGREplus(E2,a,b,{[[1,1],[1,1]],[[1,2],[1,1]],[[2,1],[1,1]]},20,7,m);`): gGREplus:=proc(F,a,b,P,M,nosaf,m) local mu,lu,gu,s,S,lu1,ku,i,mu1: lu:=gGRE(F,a,b,P,M,nosaf,m): if lu=FAIL then RETURN(FAIL): fi: mu:=gDAPfA(F,a,b,P,M): mu1:=[op(1..2,gDAPf(F,a,b,P,3)[2][1][1][1])]: gu:={}: for i from max(M-3*nosaf,10) to nops(mu) do S:=mu[i]: for s in S do for lu1 in lu do ku:=IsCont(s,lu1,m,nosaf+3): if ku<>FAIL then if ku<>[] then gu:=gu union {[op(lu1),ku]}: else gu:=gu union {lu1}: fi: fi: od: od: od: gu,mu1: end: ####start gMakeScheme #gMakeScheme(F,a,b,K): Given a second-order recurrence #x[n]=F(x[n-1],x[n-2]) given in terms of rule F #phrased in terms of a and b, tries to #make a scheme with up-to K improvements. #If it succeeds, it returns true, followed by #the successful scheme, if it fails, it returns #false followed by the unsuccessful last attempt #(with K additions) #For example, try: #gMakeScheme(E4,a,b,30,30,4): gMakeScheme:=proc(F,a,b,K,orekh,Ma) local gu,Sch,Sch1,kha,i: option remember: Sch:=FindPreScheme(F,a,b,K,orekh,Ma): RETURN(Sch): gu:=gVerifyScheme1(Sch,F,a,b): if gu then RETURN(Sch): fi: Sch1:=ImproveScheme(Sch,F,a,b): for i from 1 to K while Sch<>Sch1 do kha:=ImproveScheme(Sch1,F,a,b): Sch:=Sch1: Sch1:=kha: od: if VerifyScheme1(Sch1,F,a,b) then RETURN([Sch1,KavuaScheme(Sch1,a,b)]): else RETURN(FAIL): fi: end: #gVerifyScheme1(Sch,F,a,b): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #checks whether it is a valid scheme. For example, try #gVerifyScheme1(SchF1,F1,a,b); gVerifyScheme1:=proc(Sch,F,a,b) local i,s,v,N,IE,n,Ta,IE1,lu: for i from 1 to nops(Sch) do Ta[Sch[i][1]]:=Sch[i][3]: od: RETURN(op(Ta)): for i from 1 to nops(Sch) do lu:=Sch[i]: v:=lu[1]: N:=lu[2]: IE:=lu[3]: for n in N do IE1:=Ta[n]: IE1:=subs({a=b,b=F[v[1]][v[2]]},IE1): if not {seq(Novea(IE,s,a,b), s in IE1)}={true} then RETURN(false): fi: od: od: true: end: ####VerifySchemeVerbose #NoveaVerbose(S,L,a,b,A): Given a set of linear forms S in (a,b) #and a single linear form L in (a,b), decides whether #|L|<=A (for some A) is implied by the inequalities #|S1|<=A, |S2|<=A, for members of S, taken one-at-a-time #or two-at-a-time. For example,try: #NoveaVerbose({a,b,a-b},b-(9/10)*a,a,b,A); NoveaVerbose:=proc(S,L,a,b,A) local i,j: for i from 1 to nops(S) do if Novea1(S[i],L) then Novea1V(S[i],L,A): fi: od: for i from 1 to nops(S) do for j from i+1 to nops(S) do if Novea2(S[i],S[j],L,a,b) then Novea2V(S[i],S[j],L,a,b,A): fi: od: od: end: #Novea1V(L1,L2,A): Given Linear forms L1,L2 in #is L2=+-L1. For example, try: #Novea1V((a+b)/2,-(a+b)/2,A); Novea1V:=proc(L1,L2,A) local gu,M1,M2: gu:=normal(L1/L2): if abs(gu)=1 then print(` let M1 be`, L1, `and M2 be `, L2): print(`The absolute values of`, M1, M2, `are the same`): print(`so obviously, since |M1|<=A, it follows that |M2|<=A`): else false: fi: end: #Novea2V(L1,L2,L3,a,b): Given Linear forms L1, L2,L3 in #(a,b) expresses L3 as a linear combination of L1 and #L2 and decides whether abs(L1)<=A, abs(L2)<=A implies #abs(L3)<=A. For example, try #Novea2V(a,b,(a+b)/2,a,b); Novea2V:=proc(L1,L2,L3,a,b) local c1,c2,gu,var,eq,M1,M2,M3: var:={c1,c2}: gu:=expand(L3-c1*L1-c2*L2): eq:={coeff(gu,a,1),coeff(gu,b,1)}: var:=solve(eq,var): if var=NULL then RETURN(false): fi: c1:=subs(var,c1): c2:=subs(var,c2): if not (type(c1,numeric) and type(c2,numeric)) then RETURN(false): fi: if abs(c1)+abs(c2)<=1 then print(`Let M1 be`, L1, `Let M2 be `, L2, `and let M3 be `, L3): print(` Since `): print(M3=c1*M1+c2*M2): print(` we have by the triangle inequality`): print(`since the sum of `, abs(c1), `and `, abs(c2), `is <= 1 `): print(`that if |M1| and |M2| are both <=A`): print(`it follows that`): print(`|M3| <= A`): else RETURN(false): fi: end: #VerifySchemeV(Sch,F,a,b,A,B,Shem): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #checks whether it is a valid scheme. Shem is the name of #the lemma. #For example, try #VerifySchemeV(SchF1,F1,a,b,6); VerifySchemeV:=proc(Sch,F,a,b,Shem) local A,B: A:=Sch[2][1]: B:=Sch[2][2]: StateScheme(F,Sch,a,b,x,n,Shem): VerifyScheme1V(Sch[1],F,a,b,A,B,Shem): end: #VerifyScheme1Vold(Sch,F,a,b,A,B): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #checks whether it is a valid scheme. For example, try #VerifyScheme1Vold(SchF1,F1,a,b,A,B); VerifyScheme1Vold:=proc(Sch,F,a,b,A,B) local i,s,v,N,IE,n,Ta,IE1,lu,x: for i from 1 to nops(Sch) do Ta[Sch[i][1]]:=Sch[i][3]: od: for i from 1 to nops(Sch) do lu:=Sch[i]: v:=lu[1]: N:=lu[2]: IE:=lu[3]: print(`If right now we are in state`, v): print(`we have to prove that the linear expressions`, IE ): print(`are in absolute value <=`): print(A*abs(x[-1])+B*abs(x[0])): for n in N do print(`and we move to state`, n): IE1:=Ta[n]: print(`The linear expressions in a=x[n-1], b=x[n] that we know`): print(`by induction to be, in absolute value, <=`): print(A*abs(x[-1])+B*abs(x[0])): print(` are `, IE1): print(`applying the rule, in terms of the new a=x[n] and b=x[n+1]`): IE1:=subs({a=b,b=F[v[1]][v[2]]},IE1): print(IE1): for s in IE1 do NoveaVerbose(IE,s,a,b): od: od: od: end: #The whole story for rule F, phrased in terms of a and b #with parameters K and M (for Pers), orekh. and Ma for #FindPreScheme, nosaf and m for gGRE and gGREplus. #Shem is the name #For example, try: # Sipur1V(E2,a,b,100, 20, 15, 15, 4, 4,m,x,n,2); Sipur1V:=proc(F,a,b, K, M, K1,orekh,Ma,nosaf,m,x,n,Shem) local S,a0,m1,i,j,mu, s,gu,S2,R,ka,Z1: m1:=2: print(cat(`Theorem `, Shem,` (by Shalosh B. Ekhad)`) ): print(`Consider the following recurrence`): for i from 1 to 2 do for j from 1 to 2 do print(`if x[n-1], x[n] mod 2 equal, resp. `, i mod 2, j mod 2, `then`, x[n+1]=subs({a=x[n-1],b=x[n]},F[i][j])): od: od: for a0 from 1 to 3 do mu:=Persa0(F,a,b,a0,K,M): if member(FAIL,mu) then RETURN(FAIL): fi: od: mu:=Pers(F,a,b,K,M): if member(FAIL,mu) then print(`It seems to have trajectories that go to infinity`): RETURN(FAIL): else print(`based on empirical observation, we are lead to the following`): print(`Conjecture: Every trajectory of the rule F`): print(`with gcod(x[-1],x[0])=1 eventually`): print(`ends in one of the following orbits`): print(mu): fi: S:=FindPreScheme(F,a,b,K1,orekh,Ma): if S=FAIL then print(`Couldn't find a prescheme with the given paramaters.`): print(`Perhaps none exists`): RETURN(FAIL): fi: print(`We have the following bounds for the general term`, x[n]): if nops({seq(s[2], s in S)})=1 then print(abs(x[n])<=S[1][1][1]*abs(x[-1])+S[1][1][2]*abs(x[0])): S2:= MakeScheme(F,a,b,10): if S2<>FAIL then print(`This is proved (rigorously!) by the existence of the following scheme`): print( S2): else print(`I couldn't confirm it rigorously`): fi: else for s in S do print(`If x[n-1] mod`, m1, `equals `, s[1][1], `and `): print(`x[n] mod`, m1, ` equals`, s[1][2] , `then `): print(abs(x[n])<=s[2][1]*abs(x[-1])+s[2][2]*abs(x[0])): od: print(`This can presumbly be proved rigorously by an appropriate scheme`): print(`but it is not implemented yet`): fi: gu:=gGREplus(F,a,b,S,M,nosaf,m): ka:=gu[2]: gu:=gu[1]: print(`Without loss of generality we can have the first two elements as`,ka): print(`Because of the above inequalities, it is possible to prove`): print(`(rigorously!) that the only words in the alphabet, {0,1}`): print(`that a trajectory of the rule modulo 2, can have are`): print(gu mod 2): #print(`It is easy to find the symbolic trajectories and`): #print(` equating the first two elements with the last two elements`): #print(`yields a=0, b=0 at each case, except when you get the above orbits`): for i from 1 to nops(gu) do R:=gu[i]: print(`For the regular expression`, R): print(`applying the rule, we see, by induction`): print(`that the last two elements are`): Z1:=subs({a=ka[1],b=ka[2]}, HafelSre(F,a,b,R)): print(Z1 ): print(`By equating the first two elements with the last two elements`): print(`we see that a and b must be 0 (or one of the above-found orbits)`): print(`unless the following condition holds`): print(numer(factor(CondPair(Z1,ka,a,b)))=0): print(`and this can never happen, unless it (possibly) coincides with the`): print(`the above-found orbits (you prove it!)`): od: end: #The whole story for rule F, phrased in terms of a and b #with parameters K and M (for Pers), orekh. and Ma for #FindPreScheme, nosaf and m for gGRE and gGREplus. #For example, try: # Sipur1(E2,a,b,100, 20, 15, 15, 4, 7 ,m); Sipur1:=proc(F,a,b, K, M, K1,orekh,Ma,nosaf,m) local S,a0,mu,gu,lu: lu:=[F]: for a0 from 1 to 3 do mu:=Persa0(F,a,b,a0,K,M): if member(FAIL,mu) then RETURN(lu): fi: od: mu:=Pers(F,a,b,K,M): if member(FAIL,mu) then RETURN(lu): else lu:=[op(lu),mu]: fi: S:=FindPreScheme(F,a,b,K1,orekh,Ma): if S=FAIL then RETURN(lu): fi: lu:=[op(lu),S]: gu:=gGREplus(F,a,b,S,M,nosaf,m): if gu=FAIL then RETURN(lu): else lu:=[op(lu),[gu]]: fi: S:=MakeScheme(F,a,b,10): if S=FAIL then RETURN(lu): else lu:=[op(lu),S]: fi: lu: end: #VecToRule(v,a,b): inputs a -1,1 vector of length 8 and outputs the appropriate rule of type mod 2 #For example, try: #VecToRule([1,-1,1,-1,1,-1,1,-1],a,b): VecToRule:=proc(v,a,b): if nops(v)<>8 or ({op(v)} minus {-1,1})<>{} then RETURN(FAIL): fi: [[(v[1]*a+v[2]*b)/2,v[3]*a+v[4]*b],[v[5]*a+v[6]*b,(v[7]*a+v[8]*b)/2]]: end: #Vecs(L): all (-1,1) vectors of length L Vecs:=proc(L) local gu,g: option remember: if L=0 then RETURN([[]]): fi: gu:=Vecs(L-1): [seq([op(g),-1],g in gu), seq([op(g),1],g in gu)]: end: #AllRules(a,b): all (generalized) Amleh-Grove-Kent-Ladas-style rules AllRules:=proc(a,b) local gu,v: gu:=Vecs(8): [seq(VecToRule(v,a,b), v in gu)]: end: #Sipur(a,b, K, M, K1,orekh,Ma,nosaf,m): trying out all possible 256 rules, and carries the #analysis as far as possible, using Sipur1. For example, try: #Sipur(a,b,200, 20, 15, 15, 4, 7 ,m); Sipur:=proc(a,b, K, M, K1,orekh,Ma,nosaf,m) local R,gu,F,lu,P,i: R:=AllRules(a,b): gu:=[]: for F in R do gc(): lu:=Sipur1(F,a,b, K, M, K1,orekh,Ma,nosaf,m): gu:=[op(gu),lu]: od: for i from 1 to 5 do P[i]:=[]: od: for i from 1 to nops(gu) do P[nops(gu[i])]:=[op(P[nops(gu[i])]),gu[i]]: od: print(`Out of the 256 rules, tried with the parameters`): print([K, M, K1,orekh,Ma,nosaf]): print(`We have the following break-up`): print(`The following `, nops(P[1])): print(` rules do not seem to always stay bounded`): print(P[1]): print(`The following `, nops(P[2])): print(` seem to always stay bounded, and tend to the indicated`): print(`periodic ultimate periods (or their multiples), but`): print(`no prescheme, let alone scheme was found for them`): print(P[2]): print(`The following `, nops(P[3])): print(` seem to always stay bounded, and tend to the indicated`): print(`periodic ultimate periods (or their multiples), and`): print(`to have a prescheme, as indicated, but couldn't find`): print(`regular expressions obeying the prescheme`): print(P[3]): print(`The following `, nops(P[4])): print(` seem to always stay bounded, and tend to the indicated`): print(`periodic ultimate periods (or their multiples), and`): print(`to have a prescheme, as indicated, and possess`): print(`regular expressions obeying the prescheme`): print(`but I couldn't find an (ordinary) scheme to prove the`): print(` boundedness`): print(P[4]): print(`The following `, nops(P[5])): print(` provably always stay bounded, and tend to the indicated`): print(`periodic ultimate periods (or their multiples), and`): print(`to have a prescheme, as indicated, and possess`): print(`regular expressions obeying the prescheme`): print(`and I found an (ordinary) scheme to prove the`): print(`boundedness. In other words there are as good as can be`): print(P[5]): [seq(P[i],i=1..5)]: end: #CondPair(Z1,Z2,a,b): given pairs Z1, Z2 of linear combinations #of a and b finds the condition that there is hope for #there to be a non-zero solution to Z1=Z2 #For example, try: #CondPair([c1*a+c2*b,c3*a+c4*b],[a,b],a,b); CondPair:=proc(Z1,Z2,a,b) local A,B: A:=expand(Z1[1]-Z2[1]): B:=expand(Z1[2]-Z2[2]): expand(coeff(A,a,1)*coeff(B,b,1)-coeff(A,b,1)*coeff(B,a,1)): end: HaimMila:=proc(L): if nops(L)=2 and type(L[1],list) then RETURN(false): else RETURN(true): fi: end: #HafelSre(F,a,b,R): The regular expression applied to [a,b] #using the rule F. For example try #HafelSre(E2,a,b,[[1,2,1],[[1,2],n]]); HafelSre:=proc(F,a,b,R) local R1,kama,akha,Mik,Devek,Zug,lu,lifakha: if nops(R)=1 then R1:=R[1]: if HaimMila(R1) then RETURN(HafelSw(F,a,b,[a,b],R1)): else RETURN(HafelSwn(F,a,b,[a,b],op(R1))): fi: fi: kama:=nops(R): R1:=[op(1..nops(R)-1,R)]: Mik:=HafelSre(F,a,b,R1): akha:=R[kama]: lifakha:=R[kama-1]: if not HaimMila(lifakha) and HaimMila(akha) then Devek:=[lifakha[1][nops(lifakha[1])],akha[1]]: elif HaimMila(lifakha) and not HaimMila(akha) then Devek:=[lifakha[nops(lifakha)],akha[1][1] ]: elif not HaimMila(lifakha) and not HaimMila(akha) then Devek:=[lifakha[1][nops(lifakha[1])],akha[1][1]]: elif HaimMila(lifakha) and HaimMila(akha) then Devek:=[lifakha[nops(lifakha)],akha[1]]: fi: Zug:=HafelSw(F,a,b,Mik,Devek): if HaimMila(akha) then RETURN(HafelSw(F,a,b,Zug,akha)): else lu:=HafelSwn(F,a,b,[a,b],op(akha)): RETURN(subs({a=Zug[1],b=Zug[2]}, lu ) ): fi: end: #Concretize1(R,p,N): given a regular expression R #in terms of a list of parameters p and a list of #numeric pos. integers N of the same size, #realizes the regular expression R #For example, try: #Concretize1([[1,2],[[2,1,2],a]],[a],[2]); Concretize1:=proc(R,p,N) local w,i,R1,j,j1: if nops(p)<>nops(N) then RETURN(FAIL): fi: R1:=subs({seq(p[i]=N[i],i=1..nops(p))},R): w:=[]: for i from 1 to nops(R1) do if not HaimMila(R1[i]) then j:=R1[i][2]: for j1 from 1 to j do w:=[op(w),op(R1[i][1])]: od: else w:=[op(w),op(R1[i])]: fi: od: w: end: #CheckHafelSre(F,a,b,R,k,r): #Checks HafelSre(F,a,b,R); with r random values #with multipliticies between 1 and k For example, try: #CheckHafelSre(E2,a,b,[[1,2,1],[[1,2],n]],10,4); CheckHafelSre:=proc(F,a,b,R,k,r) local p,lu,i1,ra,N,i,w,ku: p:=[]: for i from 1 to nops(R) do if not HaimMila(R[i]) then p:=[op(p),R[i][2]]: fi: od: lu:=HafelSre(F,a,b,R): ra:=rand(1..k): ku:={}: for i from 1 to r do N:=[seq(ra(),i1=1..nops(p))]: w:=Concretize1(R,p,N): ku:= ku union { evalb( HafelSw(F,a,b,[a,b],w) =expand(subs({seq(p[i1]=N[i1],i1=1..nops(p) )},lu)))}: od: ku: end: #StateScheme(F,S,a,b,x,n,Shem): states the inequalities scheme #for the rule F, in human, verbose langhage, for example, try #StateScheme(E2,MakeScheme(E2,a,b,10),a,b,x,n,2) ; StateScheme:=proc(F,S,a,b,x,n,Shem) local i,j,S1,s1, s3,A,B: print(cat( `Lemma `, Shem, ` (by Shalosh B. Ekhad): `)): print(`Let F be the recurrence x[n+1]=F(x[n],x[n-1]) where`): for i from 1 to 2 do for j from 1 to 2 do print(`if x[n-1] and x[n] mod 2 are`, i mod 2, j mod 2, `respectively then`): print(x[n+1]=subs({a=x[n-1],b=x[n]},F[i][j])): od: od: A:=S[2][1]: B:=S[2][2]: print(`----------------------------------------------------`): print(`The following inequality holds: `, abs(x[n]) <= A*abs(x[-1])+ B*abs(x[0]) ): print(`---------------------------------------------------------`): print(`Proof (by S. B. Ekhad): To facilitate an inductive argument, we need`): print(`the following stronger statement.`): print(): print(`Stronger Lemma: `): S1:=S[1]: for s1 in S1 do print(` If x[n-1] and x[n] mod 2 are `, s1[1][1] mod 2, s1[1][2] mod 2 , `resp.`): print(`then the following inequalities are true`): for s3 in s1[3] do print(abs(subs({a=x[n-1], b=x[n]},s3) ) <=A*abs(x[-1])+B*abs(x[0]) ): od: print(`-------`): od: end: #VerifyScheme1V(Sch,F,a,b,A,B,Shem): Given a scheme for a #second-order recurrence x[n]=F(x[n-1],x[n-2]) #where F is phrased in terms of a and b, #checks whether it is a valid scheme. For example, try #VerifyScheme1V(SchF1,F1,a,b,A,B); VerifyScheme1V:=proc(Sch,F,a,b,A,B,Shem) local i,s,v,N,IE,n,Ta,ie1, lu,x, IE1,ie,nk,pu: for i from 1 to nops(Sch) do Ta[Sch[i][1]]:=Sch[i][3]: od: print(`..................`): print(`Proof of the Stronger Lemma:`): for i from 1 to nops(Sch) do lu:=Sch[i]: v:=lu[1]: N:=lu[2]: IE:=lu[3]: print(`If right now x[n-1] and x[n] mod 2 are, resp.`,v[1] mod 2, v[2] mod 2): print(`The induction hypothesis tells us that`): for ie in IE do print(abs(subs({a=x[n-1],b=x[n]},ie))<= A*abs(x[-1])+B*abs(x[0])): od: print(`--------------------------------------`): if nops(N)=1 then print(`After applying F, we arrive at the following state`): print(`x[n],x[n+1] mod 2 are, resp. `, N[1][1] mod 2, N[1][2] mod 2): else print(`After applying F, we arrive at one of the following states`): for nk in N do print(`x[n],x[n+1] mod 2 are, resp. `, nk[1] mod 2, nk[2] mod 2): od: fi: for nk in N do print(`If x[n] and x[n+1] mod 2 are, resp.`, nk[1] mod 2, nk[2] mod 2,`we have to prove `): IE1:=Ta[nk]: for ie1 in IE1 do print(abs(subs({a=x[n],b=x[n+1]},ie1))<= A*abs(x[-1])+B*abs(x[0])): od: IE1:=subs({a=b,b=F[v[1]][v[2]]},IE1): print(`in terms of x[n-1],x[n], we have to prove that`): for ie1 in IE1 do print(abs(subs({a=x[n-1],b=x[n]},ie1))<= A*abs(x[-1])+B*abs(x[0])): od: pu:=IE1: for ie1 in IE1 do if member(ie1,IE) or member(expand(-ie1),IE) then pu:=pu minus {ie1}: fi: od: if pu={} then print(`But these are all repeats of the already known inequalities`): break: fi: if pu<>IE1 then print(`removing those inequalities that we already know, this means that we have to prove`): for ie1 in pu do print(abs(subs({a=x[n-1],b=x[n]},ie1))<= A*abs(x[-1])+B*abs(x[0])): od: fi: IE1:=pu: print(`abbreviating `, a=x[n-1], b=x[n] ): for s in IE1 do NoveaVerbose(IE,s,a,b,A*abs(x[-1])+B*abs(x[0])): od: od: print(`This concludes the proof of the case where x[n-1] and x[n] mod 2 are,resp.`,v[1] mod 2, v[2] mod 2): print(`..........`): od: print(`Quod Erat Demonstrandum`): end: #Tovim(a,b, K, M, K1,orekh,Ma,nosaf,m): #trying out all possible 256 rules, and finds the good ones #For example, try #Tovim(a,b,200, 20, 15, 15, 4, 7 ,m); Tovim:=proc(a,b, K, M, K1,orekh,Ma,nosaf,m) local R,gu,F,lu: R:=AllRules(a,b): gu:=[]: for F in R do gc(): lu:=Sipur1(F,a,b, K, M, K1,orekh,Ma,nosaf,m): if nops(lu)=5 then gu:=[op(gu),F]: fi: od: gu: end: #MakeWebBook(kvu,a,b, K, M, K1,orekh,Ma,nosaf,m,x,n,shem1,shem2): #makes a web-book from the set of previously computed #good books. Each rule has two outputs, cat(shem1,i) for #the scheme and cat(shem2, i) for the rest (Sipur1V) #For example, try: #MakeWebBook({E2,E3},a,b,100, 20, 15, 15, 4, 4,m,x,n,Sch,Sip); MakeWebBook:=proc(kvu,a,b, K, M, K1,orekh,Ma,nosaf,m,x,n,shem1,shem2) local i,Sch,F: for i from 1 to nops(kvu) do gc(): F:=kvu[i]: writeto(cat(shem1,i)): Sch:=MakeScheme(F,a,b,10): VerifySchemeV(Sch,F,a,b,i): writeto(cat(shem2,i)): Sipur1V(F,a,b, K, M, K1,orekh,Ma,nosaf,m,x,n,i): od: print(`The whole thing took`, time(), `seconds of CPU time`): end: #randPers(F,a,b,a0,b0,K,M): the ultimate periods #if found, or FAIL, trying M trajectories with #random starting values between -K to K #at x[1]=a0, x[2]=b0, of the recurrence #x[n]=F(x[n-1],x[n]) given in terms of the rule F #phrased in terms of the variables a and b. #For example, try: #randPers([[(a+b)/2,a-b],[a-b,(a+b)/2]],a,b,100,5); randPers:=proc(F,a,b,K2,K1,M) local gu,a0,b0,ra,i,c0: if not CheckRule(F,a,b) then RETURN(FAIL): fi: gu:={}: ra:=rand(1..K2): for i from 1 to M do a0:=ra(): b0:=ra(): c0:=gcd(a0,b0): a0:=a0/c0: b0:=b0/c0: gu:=gu union {Canon(Per(F,a,b,a0,b0,K1))}: od: gu: end: