###################################################################### ##COLLATZ: Save this file as COLLATZ To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read COLLATZ : # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: June 27, 2008 print(`Created: June 27, 2008`): print(` This is COLLATZ, a Maple package`): print(`to investigate (only empirically) (one dimensional) extensions of the 3x+1 conjecture`): print(`It is one two Maple packages that accompany 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(`and also available from his website`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`For a list of the procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): C1:=[(3*n+1)/2,n/2]: Z1:=[(n+2)*5/3,(n+1)/3,4*n/3]: Z2:=[(n+2)*4/3,(n+1)/3,4*n/3]: Z3:=[(n+2)*4/3,(n+1)/3,n/3]: with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are:`): print( ` Can1, ContOrbS, Hafel , HafelS, HafelSs, MaxMag, Merkhak,`): print(` MerkhakT, MostStuborn1, MostStuborn2 `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: `): print(` Akshanim, Orb, Orbs, OrbsS, PerOrbs, PerOrbsN, Survivors `): elif nops([args])=1 and op(1,[args])=Akshanim then print(`Akshanim(n0,L,n,K): the most stuborn starting values`): print(`for rule L up to n0 ending in a period, allowing`): print(`up to K iterations, followed by those that`): print(`require more than K iterations, followed by the record.`): print(`For example, try:`): print(`Akshanim(100,[(3*n+1)/2,n/2],n,2000) ;`): print(`For the Z rule, try:`): print(`Akshanim(100,[(n+2)*5/3,(n+1)/3,4*n/3],n,2000); `): elif nops([args])=1 and op(1,[args])=Can1 then print(`Can1(L): the canonical form of an orbit, with the smallest at the`): print(`beginning. For example, try: Can1([2,4,1,5]);`): elif nops([args])=1 and op(1,[args])=ContOrbS then print(`ContOrbS(A,L,n): given an orbit`): print(`(list of expressions of the`): print(`form an+b applies rule L to the last entry`): print(`and appends its iamge (if possible)`): print(`splits it into several`): print(`suborbits for which you can apply rule L`): print(`(phrased in terms of n)`): print(`For example, try:`): print(`ContOrbsS([n],[(3*n+1)/2,n/2],n); `): elif nops([args])=1 and op(1,[args])=Hafel then print(`Hafel(n0,L,n): given a rule L as a list of length m`): print(`of expressions in n such that if n mod m =i then`): print(`L[i] is applied, find its image. For example, for the`): print(`Collatz rule, applied to 10, try: `): print(`Hafel(10,[(3*n+1)/2,n/2],n); `): elif nops([args])=1 and op(1,[args])=HafelS then print(`HafelS(B,L,n): given a rule L as a list of length m`): print(`of expressions in n such that if n mod m =i then`): print(`L[i] is applied, find the image under that law when`): print(`it is applied on the affine-linear expression`): print(`of the form a1*n+a2, with a1 and a2 integers`): print(`and returns FAIL if it can't tell (i.e. not enough information)`): print(`For example, for the`): print(`Collatz rule, applied to 4*n+1 try: `): print(`HafelS(4*n+1,[(3*n+1)/2,n/2],n); `): print(`For the Z rule, try:`): print(`HafelS(3*n+1,[(n+2)*5/3,(n+1)/3,4*n/3],n); `): elif nops([args])=1 and op(1,[args])=HafelSs then print(`HafelSs(B,L,n): given an expression B of the`): print(`form an+b applies rule L if possible`): print(`and if not, splits the input into several`): print(`subclasses for which you can apply rule L`): print(`(phrased in terms of n)`): print(`For example, try:`): print(`HafelSs(n,[(3*n+1)/2,n/2],n);`): elif nops([args])=1 and op(1,[args])=MaxMag then print(`MaxMag(n0,L,n,K): the max. ratio of the largest element`): print(`in the orbit of n0 starting at L to n0, together with the champion`): print(`or FAIL`): print(`allowing up to K iterations`): print(`For example, try:`): print(`MaxMag(100,[(3*n+1)/2,n/2],n,1000);`): elif nops([args])=1 and op(1,[args])=Merkhak then print(`Merkhak(n0,L,n): the distance of integer n0 from 1`): print(`using Collatz-rule L phrased in terms of n. `): print(`Warning (it tacitly assumed that all orbits go to 1)`): print(`For example, try `): print(`Merkhak(17,[(3*n+1)/2,n/2],n);`): elif nops([args])=1 and op(1,[args])=MerkhakT then print(`MerkhakT(n0,L,n): the number of iterations`): print(`it takes to reduce n0`): print(`using Collatz-rule L phrased in terms of n. `): print(`Warning (it tacitly assumed that all orbits go to 1)`): print(`For example, try `): print(`MerkhakT(17,[(3*n+1)/2,n/2],n);`): elif nops([args])=1 and op(1,[args])=MostStuborn1 then print(`MostStuborn1(L,n,N0): Given a Collatz-style rule L,`): print(`phrased in terms of n, and a pos. integer N0`): print(`finds the integer between 1 and N0 that takes`): print(`longest to reach 1, followed by the record`): print(`(it tacitly assumes that all trajectories go to 1)`): print(`For example, try:`): print(`MostStuborn1([(3*n+1)/2,n/2],n,100);`): elif nops([args])=1 and op(1,[args])=MostStuborn2 then print(`MostStuborn2(L,n,N0): Given a Collatz-style rule L,`): print(`phrased in terms of n, and a pos. integer N0`): print(`finds the integer between 1 and N0 that takes`): print(`longest to get smaller, followed by the record`): print(`(it tacitly assumes that all trajectories go to 1)`): print(`For example, try:`): print(`MostStuborn1([(3*n+1)/2,n/2],n,100);`): elif nops([args])=1 and op(1,[args])=Orb then print(`Orb(n0,L,n,K): the orbit of rule L, phrased in terms of`): print(`n, starting at n0 of length up to K terms`): print(`the output is the beginning followed by a period, or`): print(`FAIL. For example, try:`): print(`Orb(10,[3*n+1,n/2],n,100); `): elif nops([args])=1 and op(1,[args])=Orbs then print(`Orbs(N0,L,n,K): Given a Collatz-like rule L phrased in terms of`): print(`n, and a pos. integer N0, and a large pos. integer K`): print(`finds the set of ultimate orbits obtained by running all initial`): print(`values up to N0 (using brute force), including FAIL if it`): print(`sems to run off to infinity. For example, try:`): print(`Orbs(1000,[3*n+1,n/2],n,100); `): print(`For the Z rule, try:`): print(`Orbs(1000,[(n+2)*5/3,(n+1)/3,4*n/3],n,100); `): elif nops([args])=1 and op(1,[args])=OrbsS then print(`OrbsS(L,n,k): Given a Collatz-type rule L, phrased in terms of`): print(`the variable n, lists all symbolic orbits of length k.`): print(`For example, for all the orbits of length 4 for the`): print(`Collatz rule, do:`): print(`OrbsS([(3*n+1)/2,n/2],n,4); `): print(`For the Z rule, try:`): print(`OrbsS([(n+2)*5/3,(n+1)/3,4*n/3],n,4); `): elif nops([args])=1 and op(1,[args])=PerOrbs then print(`PerOrbs(L,n,k): all the primitive`): print(` periodic orbits of the Collatz-type`): print(`rule L (phrased in terms of n) of length k. For example, try:`): print(`PerOrbs([(3*n+1)/2,n/2],n,4); `): print(`For the Z rule, try:`): print(`PerOrbs([(n+2)*5/3,(n+1)/3,4*n/3],n,4); `): elif nops([args])=1 and op(1,[args])=PerOrbsN then print(`PerOrbsN(L,n,k): all the primitive`): print(` periodic orbits of the Collatz-type`): print(`rule L (phrased in terms of n) of length k, including negatives`): print(`For example, try:`): print(`PerOrbsN([(3*n+1)/2,n/2],n,4); `): print(`For the Z rule, try:`): print(`PerOrbsN([(n+2)*5/3,(n+1)/3,4*n/3],n,4); `): elif nops([args])=1 and op(1,[args])=Survivors then print(`Survivors(L,n,k): Given a Collatz-type rule L, phrased in terms of`): print(`the variable n, lists all symbolic orbits of length k`): print(`that are still expanding after k generations`): print(`For example, for the`): print(`Collatz rule, do:`): print(`Survivors([3*n+1,n/2],n,4); `): print(`For the Z rule, try:`): print(`Survivors([(n+2)*5/3,(n+1)/3,4*n/3],n,4); `): else print(`There is no ezra for`,args): fi: end: #Hafel(n0,L,n): given a rule L as a list of length m #of expressions in n such that if n mod m =i then #L[i] is applied, find its image. For example, for the #Collatz rule, applied to 10, try: #Hafel(10,[3*n+1,n/2],n); Hafel:=proc(n0,L,n) local i,m: m:=nops(L): i:= n0 mod m: if i=0 then i:=m: fi: subs(n=n0,L[i]): end: #Orb(n0,L,n,K): the orbit of rule L, phrased in terms of #n, starting at n0 of length up to K terms #the output is the beginning followed by a period, or #FAIL. For example, try: #Orb(10,[3*n+1,n/2],n,100); Orb:=proc(n0,L,n,K) local i,j,cu,gu: cu:=n0: gu:=[n0]: for i from 1 to K do cu:=Hafel(cu,L,n): if member(cu,gu) then for j from nops(gu) by -1 to 1 while gu[j]<>cu do od: RETURN([ [op(1..j-1,gu)],[op(j..nops(gu),gu)] ]): else gu:=[op(gu),cu]: fi: od: [gu,FAIL]: end: #Can1(L): the canonical form of an orbit, with the smallest at the #beginning. For example, try: Can1([2,4,1,5]); Can1:=proc(L) local i,m: if L=FAIL then RETURN(L): fi: m:=min(op(L)): for i from 1 while L[i]<>m do od: [op(i..nops(L),L),op(1..i-1,L)]: end: #Orbs(N0,L,n,K): Given a Collatz-like rule L phrased in terms of #n, and a pos. integer N0, and a large pos. integer K #finds the set of ultimate orbits obtained by running all initial #values up to N0 (using brute force), including FAIL if it #sems to run off to infinity. For example, try: #Orbs(1000,[3*n+1,n/2],n,100); Orbs:=proc(N0,L,n,K) local n0: {seq(Can1(Orb(n0,L,n,K)[2]),n0=1..N0)}: end: #HafelS(B,L,n): given a rule L as a list of length m #of expressions in n such that if n mod m =i then #L[i] is applied, find the image under that law when #it is applied on the affine-linear expression #of the form a1*n+a2, with a1 and a2 integers #and returns FAIL if it can't tell (i.e. not enough information) #For example, for the #Collatz rule, applied to 4*n+1 try: #HafelS(4*n+1,[3*n+1,n/2],n); #For the Z rule, try: #HafelS(3*n+1,[(n+2)*5/3,(n+1)/3,4*n/3],n); HafelS:=proc(B,L,n) local i,m,a1,a2: m:=nops(L): a1:=coeff(B,n,1): a2:=coeff(B,n,0): if not type(a1/m,integer) then RETURN(FAIL): fi: i:= a2 mod m: if i=0 then i:=m: fi: subs(n=B,L[i]): end: #HafelSs(B,L,n): given an expression B of the #form an+b applies rule L if possible #and if not, splits the input into several #subclasses for which you can apply rule L #(phrased in terms of n) #For example, try: #HafelSs(n,[(3*n+1)/2,n/2],n); HafelSs:=proc(B,L,n) local i,m,a1,A,B1,mu: m:=nops(L): a1:=coeff(B,n,1): if a1 mod m=0 then RETURN({[B,HafelS(B,L,n)]}): fi: A:=m/gcd(a1,m): mu:=[]: for i from 1 to A-1 do B1:=subs(n=A*n+i,B): mu:=[op(mu),[B1,HafelS(B1,L,n)]]: od: B1:=subs(n=A*n,B): [op(mu),[B1,HafelS(B1,L,n)]]: end: #ContOrbS(A,L,n): given an orbit #(list of expressions of the #form an+b applies rule L to the last entry #and appends its iamge (if possible) #splits it into several #suborbits for which you can apply rule L #(phrased in terms of n) #For example, try: #ContOrbsS([n],[(3*n+1)/2,n/2],n); ContOrbsS:=proc(LiB,L,n) local i,m,a1,A,B1,mu,B: m:=nops(L): B:=LiB[nops(LiB)]: a1:=coeff(B,n,1): if a1 mod m=0 then RETURN({[op(LiB),HafelS(B,L,n)]}): fi: A:=m/gcd(a1,m): mu:={}: for i from 1 to A-1 do B1:=subs(n=A*n+i,B): mu:=mu union { [op(subs(n=A*n+i,LiB)),HafelS(B1,L,n)]}: od: B1:=subs(n=A*n,B): mu:=mu union { [op(subs(n=A*n,LiB)),HafelS(B1,L,n)]}: mu: end: #OrbsS(L,n,k): Given a Collatz-type rule L, phrased in terms of #the variable n, lists all symbolic orbits of length k #For example, for all the orbits of length 4 for the #Collatz rule, do: #OrbsS([3*n+1,n/2],n,4); #For the Z rule, try: #OrbsS([(n+2)*5/3,(n+1)/3,4*n/3],n,4); OrbsS:=proc(L,n,k) local gu,g: option remember: if k=1 then RETURN({[n]}): fi: gu:=OrbsS(L,n,k-1): {seq(op(ContOrbsS(g,L,n)), g in gu)}: end: with(numtheory): #IsPrim(L): is the periodic sequence L primitive? #For example, try: #IsPrim([4,1,3,4,1,3]); IsPrim:=proc(L) local m,gu,g,ka,j: m:=nops(L): gu:=divisors(m) minus {m}: for g in gu do ka:={seq([op(j*g+1..(j+1)*g,L)],j=0..m/g-1)}: if nops(ka)=1 then RETURN(false): fi: od: true: end: #PerOrbsOLD(L,n,k): all the periodic orbits of the Collatz-type #rule L (phrased in terms of n) of length k. For example, try: ##PerOrbsOLD([3*n+1,n/2],n,4); #For the Z rule, try: #PerOrbsOLD([(n+2)*5/3,(n+1)/3,4*n/3],n,4); PerOrbsOLD:=proc(L,n,k) local mu,gu,m,n0,m1,gu1,g1,g: mu:=OrbsS(L,n,k+1): gu:={}; for m in mu do if m[1]=m[k+1] then gu:=gu union {m}: fi: n0:=subs(solve({m[1]=m[k+1]},{n}),n): if type(n0, integer) and n0>=0 then m1:=subs(n=n0,m): if not m1=[0$nops(m1)] then gu:=gu union {[op(1..k,m1)]}: fi: fi: od: gu1:={seq(Can1(g), g in gu)}: gu:={}: for g1 in gu1 do if IsPrim(g1) then gu:=gu union {g1}: fi: od: gu: end: #Akshanim(n0,L,n,K): the most stuborn starting values #for rule L up to n0 ending in a period, allowing #up to K iterations, followed by those that #require more than K iterations, followed by the record. #For example, try: #Akshanim(100,[(3*n+1)/2,n/2],n,2000) ; Akshanim:=proc(n0,L,n,K) local gu1,gu2,i,lu,shi: gu1:={}: gu2:={}: shi:=0: for i from 1 to n0 do lu:=Orb(i,L,n,K) : if lu[2]=FAIL then gu2:=gu2 union {i}: else if nops(lu[1])>shi then gu1:={i}: shi:=nops(lu[1]): elif nops(lu[1])=shi then gu1:=gu1 union {i}: fi: fi: od: gu1,shi,gu2: end: #Survivors(L,n,k): Given a Collatz-type rule L, phrased in terms of #the variable n, lists all symbolic orbits of length k #that are still expanding after k generations #For example, for the #Collatz rule, do: #Survivors([3*n+1,n/2],n,4); #For the Z rule, try: #Survivors([(n+2)*5/3,(n+1)/3,4*n/3],n,4); Survivors:=proc(L,n,k) local gu,g,mu,g1,g11,a,b,A,B,ka,g11a,i: option remember: if k=1 then RETURN({[n]}): fi: gu:=Survivors(L,n,k-1): mu:={}: for g in gu do g1:=ContOrbsS(g,L,n): for g11 in g1 do a:= coeff(g11[1],n,1): if a<>0 or (a=0 and nops(g11)=nops(convert(g11,set))) then A:=coeff(g11[nops(g11)],n,1): b:= coeff(g11[1],n,0): B:=coeff(g11[nops(g11)],n,0): if A>=a and B>b then mu:=mu union {g11}: elif A>a and B>=b then mu:=mu union {g11}: elif Ab then ka:=trunc((B-b)/(a-A)): for i from 0 to ka do g11a:=subs(n=i, g11): if nops(g11a)=nops({op(g11a)}) then mu:=mu union {g11a}: fi: od: fi: fi: od: od: mu: end: #MaxMag1(n0,L,n,K): the ratio of the largest element #in the orbit of n0 starting at L to n0, or FAIL #allowing up to K iterations #For example, try: #MaxMag1(10,[(3*n+1)/2,n/2],n,1000); MaxMag1:=proc(n0,L,n,K) local gu: gu:=Orb(n0,L,n,K): if gu=FAIL then RETURN(FAIL): else max(op(gu[1]))/n0: fi: end: #MaxMag(n0,L,n,K): the max. ratio of the largest element #in the orbit of n0 starting at L to n0, together with the champion #or FAIL #allowing up to K iterations #For example, try: #MaxMag(100,[(3*n+1)/2,n/2],n,1000); MaxMag:=proc(n0,L,n,K) local aluf, shi,n1,lu: aluf:=1: shi:=1: for n1 from 1 to n0 do lu:=MaxMag1(n1,L,n,K) : if lu>shi then aluf:=n1: shi:=lu: fi: od: aluf,shi: end: #Merkhak(n0,L,n): the distance of integer n0 from 1 #using Collatz-rule L phrased in terms of n. #Warning (it tacitly assumed that all orbits go to 1) #For example, try #Merkhak(17,[(3*n+1)/2,n/2],n); Merkhak:=proc(n0,L,n) option remember: if n0=1 then 0: else Merkhak(Hafel(n0,L,n),L,n)+1: fi: end: #MerkhakT(n0,L,n): the number of steps to reduce integer n0 from 1 #using Collatz-rule L phrased in terms of n. #Warning (it tacitly assumed that all orbits go to 1) #For example, try #MerkhakT(17,[(3*n+1)/2,n/2],n); MerkhakT:=proc(n0,L,n) local lu: option remember: if n0=1 then RETURN(0): fi: lu:=Hafel(n0,L,n): if lu=0 then m1:=subs(n=n0,m): if not m1=[0$nops(m1)] then gu:=gu union {[op(1..k,m1)]}: fi: fi: od: gu1:={seq(Can1(g), g in gu)}: gu:={}: for g1 in gu1 do if IsPrim(g1) then gu:=gu union {g1}: fi: od: gu: end: #OrbsF(N0,L,n,K): Given a Collatz-like rule L phrased in terms of #n, and a pos. integer N0, and a large pos. integer K #finds the set of ultimate orbits obtained by running all initial #values up to N0 (using brute force), including FAIL if it #sems to run off to infinity. For example, try: #Orbs(1000,[3*n+1,n/2],n,100); OrbsN:=proc(N0,L,n,K) local n0,StillToDo,gu,lu: StillToDo:={seq(n0,n0=1..N0)}: gu:={}: while StillToDo<>{} do n0:=op(1,StillToDo): lu:=Orb(n0,L,n,K): StillToDo:=StillToDo minus {op(lu[1]),op(lu[2])}: gu:=gu union {Can1(lu[2])}: od: gu: end: #MostStuborn1(L,n,N0): Given a Collatz-style rule L, #phrased in terms of n, and a pos. integer N0 #finds the integer between 1 and N0 that takes #longest to reach 1, followed by the record #(it tacitly assumes that all trajectories go to 1) #For example, try: #MostStuborn1([(3*n+1)/2,n/2],n,100); MostStuborn1:=proc(L,n,N0) local n0,aluf,shi,hale: aluf:={}: shi:=0: for n0 from 1 to N0 do hale:=Merkhak(n0,L,n) : if hale>shi then shi:=hale: aluf:={n0}: elif hale=shi then aluf:=aluf union {n0}: fi: od: aluf,shi: end: #MostStuborn2(L,n,N0): Given a Collatz-style rule L, #phrased in terms of n, and a pos. integer N0 #finds the integer between 1 and N0 that takes #longest to reach a smaller integer, followed by the record #(it tacitly assumes that all trajectories go to 1) #For example, try: #MostStuborn2([(3*n+1)/2,n/2],n,100); MostStuborn2:=proc(L,n,N0) local n0,aluf,shi,hale: aluf:={}: shi:=0: for n0 from 1 to N0 do hale:=MerkhakT(n0,L,n) : if hale>shi then shi:=hale: aluf:={n0}: elif hale=shi then aluf:=aluf union {n0}: fi: od: aluf,shi: end: