Math Kook

This notebook contains the commands used in the creation of Math Kook (K. Knight, Threeven Press, 2022).

This code is provided as-is, with a hearty “good luck” from the author, but without guarantees or support.

Chapter 1. What’s the problem?

The 3n+1 problem

C1[x_]:=If[EvenQ[x],x/2,(3x+1)]

NestWhileList[C1,7,#!=1&]

NestWhileList[C1,84,#!=1&]

NestWhileList[C1,1438,#!=1&]

Counter-examples

Length[NestWhileList[C1,28828795202195990691,#!= 1&]]

569936821221962380720^3+(-569936821113563493509)^3 +(-472715493453327032)^3

Chapter 2. Let’s prove it!

An equivalent formulation

T[x_]:=If[EvenQ[x],x/2,(3x+1)/2]

NestList[C1,18,20]

NestList[C1,33,20]

NestList[C1,55,20]

NestWhileList[C1,17,#!=1&]

NestWhileList[T,17,#!=1&]

NestList[T,27,70]

The 3n+1 tree

CollatzTree3nPlusOne[n_,cutoff_,height_,depth_]:=
If[(depth==1)||(n>cutoff),n,
If[2n>cutoff,
If[(n>4)&&(Mod[n,3]==2),
n[CollatzTree3nPlusOne[(2n-1)/3,cutoff,If[OddQ[n],height+1,height],depth-1]],n],
If[(n>4)&&(Mod[n,3]==2),
n[CollatzTree3nPlusOne[(2n-1)/3,cutoff,If[OddQ[n],height+1,height],depth-1],
CollatzTree3nPlusOne[2n,cutoff,If[OddQ[n],height+1,height],depth-1]],
n[CollatzTree3nPlusOne[2n,cutoff,If[OddQ[n],height+1,height],depth-1]]]]]

TreeForm[CollatzTree3nPlusOne[1,1000,0,14],VertexRenderingFunction->({White,Disk[#,.2],Red,Text[#2,#1]}&)];

Network`GraphPlot`ExprTreePlot[CollatzTree3nPlusOne[1,1000000,0,11],Bottom]

How long to reach 1?

BarChart[Table[Length[NestWhileList[T,n,#1≠1&]]-1,{n,2,50}],PlotTheme->"GrayColor",ChartLabels->{Placed[Style[#,FontSize->Scaled[.012]]&/@Range[2,50],Below]}]

BarChart[Table[Length[NestWhileList[T,n,#1≠1&]]-1,{n,2,50}],PlotTheme->"GrayColor",ChartLabels->{Placed[Style[#,FontSize->Scaled[.015]]&/@Range[2,50],Above]}]

BarChart[Table[n^1.1,{n,1,50}],PlotTheme->"GrayColor"]

|

BarChart[Table[Log[1.1,n],{n,1,50}],PlotTheme->"GrayColor"]

BarChart[Table[n^2/50,{n,1,50}],PlotTheme->"GrayColor"]

Our first pattern

Table[{n,Length[NestWhileList[T,n,#≠1&]]-1},{n,1,51}]

NestWhileList[T,49,#≠1&]

NestWhileList[T,50,#≠1&]

NestWhileList[T,51,#≠1&]

NestWhileList[T,28,#≠1&]

NestWhileList[T,29,#≠1&]

NestWhileList[T,12,#≠1&]

NestWhileList[T,13,#1≠1&]

NestWhileList[T,20,#≠1&]

NestWhileList[T,21,#≠1&]

Larger numbers take longer

BarChart[Table[Length[NestWhileList[T,n,#1≠1&]]-1,{n,1,500}],PlotTheme->"GrayColor"]

Table[{n, Length[NestWhileList[T,n,#1≠1&]]-1},{n,1,500}];

Sequences trend downwards

ToBase10Big[x_]:=IntegerDigits[x,10,225]

ArrayPlot[Map[ToBase10Big,b=NestWhileList[T,a=826474739491445506384651754684,#>1&]],ColorRules->{0->White},ColorFunction->"GrayLevel"]

a

Length[Map[If[EvenQ[#],1,Nothing]&,b]]

Length[Map[If[OddQ[#],1,Nothing]&,b]]

A more organized census

Table[{n,NestList[T,n,4]},{n,0,15}]

With[{m=4},Table[If[Or[n<2,Last[NestList[T,n,m]]<n],n,Nothing],{n,0,2^m-1}]]

Table[{2^m,Length[Table[If[Or[n<2,Last[NestList[T,n,m]]<n],True,Nothing],{n,0,2^m-1}]]/2^m*1.},{m,4,16}]

Structure

Log[3,2.]

Coin flipping

Total[Table[Binomial[10,n],{n,6,10}]]/2^10.

Total[Table[Binomial[100,n],{n,60,100}]]/2^100.

SetPrecision[Total[Table[Binomial[1000,n],{n,600,1000}]]/2^1000.,30]

Histogram[{Table[With[{b=Total[Table[RandomInteger[1],{n,1,10}]]},If[b< 6,b,Nothing]],{x,1,10000}],Table[With[{b=Total[Table[RandomInteger[1],{n,1,10}]]},If[b>= 6,b,Nothing]],{x,1,10000}]},{"Raw",800},ChartStyle->{LightGray,Black},Ticks->{{#+.5,#}&/@Range[0,10],Automatic}]

Histogram[{Table[With[{b=Total[Table[RandomInteger[1],{n,1,100}]]},If[b< 60,b,Nothing]],{x,1,10000}],Table[With[{b=Total[Table[RandomInteger[1],{n,1,100}]]},If[b>= 60,b,Nothing]],{x,1,10000}]},{"Raw",800},ChartStyle->{LightGray,Black},ChartLabels->{Placed[Style[#,FontSize->Scaled[.005]]&/@Range[0,100],Below]}]

Histogram[{Table[With[{b=Total[Table[RandomInteger[1],{n,1,1000}]]},If[b< 600,b,Nothing]],{x,1,10000}],Table[With[{b=Total[Table[RandomInteger[1],{n,1,1000}]]},If[b>= 600,b,Nothing]],{x,1,10000}]},{"Raw",800},ChartStyle->{LightGray,Black}]

Every operation sequence seen once and only once

Table[{b,(16b-1)/3},{b,1,10}]

Chapter 3. Let’s disprove it!

A heuristic argument

2^20-1

Table[1./2^(y-2),{y,{11,12,13,14,15,20,100}}]

Sum[1./2^(y-2),{y,11,Infinity}]

3^7

Where to look for 3n+1 cycles

Denom[k_,x_]:=2^k-3^x

Duval1[n_,prev_]:=Take[Flatten[NestList[Join,prev,Ceiling[n/Length[prev]]]],n]

Duval2[prev_]:=If[Last[prev]==1,Duval2[Take[prev,Length[prev]-1]],prev]

Duval3[prev_]:=ReplacePart[prev,Length[prev]->1]

Duval[n_,prev_]:=Duval3[Duval2[Duval1[n,prev]]]

Necklaces[k_,x_]:=Map[Reverse,Select[NestWhileList[Duval[k,#]&,{0},#!={1}&],Total[#]==x&&Length[#]==k&]]

Necklaces[5,3]

Necklaces[8,5]

BetaAnd[l_,i_,x_]:=If[Last[Take[l]]==1,2^(i-1)*3^(x-Count[l,1]),0]

BetaValue[necklace_]:=Total[Table[BetaAnd[Take[necklace,i],i,Count[necklace,1]],{i,1,Length[necklace]}]]

Map[BetaValue,Necklaces[5,3]]

LoopsByMember[k_,x_]:=Map[BetaValue[#]/Denom[k,x]&,Necklaces[k,x]]

LoopsByMember[5,3]

Map[{BetaValue[#],#}&,Necklaces[5,3]]

LoopForOperationSequence[n_,seq_,res_]:=If[seq=={},Reverse[res],If[First[seq]==0,LoopForOperationSequence[n/2,Rest[seq],Insert[res,n/2,1]],LoopForOperationSequence[(3n+1)/2,Rest[seq],Insert[res,(3n+1)/2,1]]]]

LoopForOperationSequence[19/5,{1,1,1,0,0},{19/5}]

LoopForOperationSequence[23/5,{1,1,0,1,0},{23/5}]

How many shots do we get?

Necklaces[8,5]

CountNecklaces[k_,x_]:=(1/k) Sum[MoebiusMu[j] * Binomial[k/j,x/j],{j,Divisors[GCD[k,x]]}]

Table[With[{k=kx[[1]],x=kx[[2]]},{k,x,CountNecklaces[k,x],1.000000000000000000*Binomial[k,x]/k}],{kx,{{6,4},{10,7},{19,12},{42,27}}}]

Expected number of cycles

1.-4/5 * 4/5

ExactTotalNecklaces[k_,x_]:=If[k==1,1,(1/k)Total[Map[MoebiusMu[#]*Binomial[k/#,x/#]&,Divisors[GCD[k,x]]]]]

ExpectedNumberOfCycles[q_,k_]:=ExactTotalNecklaces[k,Floor[k/Log[2,q]]]/Denom[q,k]

ExpectedNumberOfCycles[3,8]

ApproxTotalNecklaces[k_,x_]:=If[k==1,1,1.0*Binomial[k,x]/k]

ApproxTotalNecklaces[5,3]

GenDenom[q_,k_]:=If[k==1,100000,2^k-q^Floor[k/Log[2,q]]]

ProbOfCycle[q_,k_]:=If[ApproxTotalNecklaces[k,Floor[k/Log[2,q]]]==0,0.0,1.-(1-1/GenDenom[q,k])^ApproxTotalNecklaces[k,Floor[k/Log[2,q]]]]

ProbOfCycle[3,5]

ExpectedNumberOfCycles[q_,k_]:=ApproxTotalNecklaces[k,Floor[k/Log[2,q]]]/GenDenom[q,k]

ExpectedNumberOfCycles[3,5]

ExpectedNumberOfCycles[3,1000]*1.

ArrayPlot[Table[ProbOfCycle[2a+1,k],{a,1,1},{k,1,20}],Mesh->All]

3n+1

Table[{k,Floor[k/Log[2,3]]},{k,5,10}]

Table[{k,ExpectedNumberOfCycles[2a+1,k]*1., ProbOfCycle[2a+1,k]},{a,1,1},{k,5,10}]

Table[{k,ExpectedNumberOfCycles[2a+1,k]*1., ProbOfCycle[2a+1,k]*1.},{a,1,1},{k,100,100}]

Table[{k,ExpectedNumberOfCycles[2a+1,k]*1., ProbOfCycle[2a+1,k]*1.},{a,1,1},{k,1000,1000}]

Total expected cycles

Map[Total[#]&, Table[ExpectedNumberOfCycles[2a+1,k]*1.,{a,1,1},{k,5,20}]]

Map[Total,Table[ExpectedNumberOfCycles[2a+1,k]*1.,{a,1,1},{k,5,100}]]

Map[Total,Table[ExpectedNumberOfCycles[2a+1,k]*1.,{a,1,1},{k,5,1000}]]

Map[Total,Table[ExpectedNumberOfCycles[2a+1,k]*1.,{a,1,1},{k,5,10000}]]

Sum[(Sqrt[2 Pi k]/(Sqrt[2 Pi(k/Log[2,3])] Sqrt[2 Pi (k - (k/Log[2,3]))]))*((k/E)^k /(((k/Log[2,3])/E)^(k/Log[2,3]) * ((k-k/Log[2,3])/E)^(k-k/Log[2,3])))/(k 2^k),{k,1,100}]*1.

Simplify[(((k/2.7182818284)^k  Sqrt[2 3.1415926535k])/
(((k/1.5849625007211563)/2.7182818284)^(k/1.5849625007211563) * 
Sqrt[2 Pi(k/1.5849625007211563)]*
((k-k/1.5849625007211563)/2.7182818284)^(k-k/1.5849625007211563) * 
Sqrt[2 Pi (k-(k/1.5849625007211563))]))/(k 2^k)]

E^-0.034688

Sum[(0.82673 * 0.96591^k)/(k Sqrt[k]),{k,1,Infinity}]

Chapter 4. The qn+1 problem

Where to look for 3n-1 cycles

NegDenom[q_,k_]:=If[(k==1)&&(q>3),100000,2^k-q^Ceiling[k/Log[2,q]]]

NegExpectedNumberOfCycles[q_,k_]:=-ExactTotalNecklaces[k,Ceiling[k/Log[2,q]]]/NegDenom[q,k]

NegProbOfCycle[q_,k_]:=If[ExactTotalNecklaces[k,Ceiling[k/Log[2,q]]]==0,0.0,1.-(1-1/-NegDenom[q,k])^ExactTotalNecklaces[k,Ceiling[k/Log[2,q]]]]

Table[{k,Ceiling[k/Log[2,3]]},{k,5,10}]

Table[{k, NegExpectedNumberOfCycles[2a+1,k]*1.,NegProbOfCycle[2a+1,k]},{a,1,1},{k,5,20}]

The 3n-1 tree

CollatzTree3nMinusOne[n_,cutoff_,height_,depth_]:=
If[(depth==1)||(n>cutoff),n,
If[2n>cutoff,
If[(n>2)&&(Mod[n,3]==1),
n[CollatzTree3nMinusOne[(2n+1)/3,cutoff,If[OddQ[n],height+1,height],depth-1]],n],
If[(n>2)&&(n!=7)&&(Mod[n,3]==1),
n[CollatzTree3nMinusOne[(2n+1)/3,cutoff,If[OddQ[n],height+1,height],depth-1],
CollatzTree3nMinusOne[2n,cutoff,If[OddQ[n],height+1,height],depth-1]],
n[CollatzTree3nMinusOne[2n,cutoff,If[OddQ[n],height+1,height],depth-1]]]]]

Network`GraphPlot`ExprTreePlot[CollatzTree3nMinusOne[1,100000,0,12],Bottom]

Network`GraphPlot`ExprTreePlot[CollatzTree3nMinusOne[5,100000,0,10],Bottom]

The 5n+1 problem

Table[{k,Floor[k/Log[2,5]]},{k,5,10}]

Table[{k,ExpectedNumberOfCycles[2a+1,k]*1., ProbOfCycle[2a+1,k]},{a,2,2},{k,5,14}]

T5[n_]:=If[EvenQ[n],n/2,(5n+1)/2]

NestList[T5,13,7]

NestList[T5,17,7]

CollatzTree5nPlusOne[n_,cutoff_,height_,depth_]:=
If[(depth==1)||(n>cutoff),n,
If[2n>cutoff,
If[(n>2)&&(Mod[n,5]==3),
n[CollatzTree5nPlusOne[(2n-1)/5,cutoff,If[OddQ[n],height+1,height],depth-1]],n],
If[(n>2)&&(n!=3)&&(n!=33)&&(n!=43)&&(Mod[n,5]==3),
n[CollatzTree5nPlusOne[(2n-1)/5,cutoff,If[OddQ[n],height+1,height],depth-1],
CollatzTree5nPlusOne[2n,cutoff,If[OddQ[n],height+1,height],depth-1]],
n[CollatzTree5nPlusOne[2n,cutoff,If[OddQ[n],height+1,height],depth-1]]]]]

Network`GraphPlot`ExprTreePlot[CollatzTree5nPlusOne[1,1000000,0,14],Bottom]

Network`GraphPlot`ExprTreePlot[CollatzTree5nPlusOne[13,1000000,0,12],Bottom]

Network`GraphPlot`ExprTreePlot[CollatzTree5nPlusOne[17,1000000,0,12],Bottom]

Cycles for the qn+1 problem

ArrayPlot[Table[ExpectedNumberOfCycles[2a+1,k],{a,1,12},{k,1,20}],Mesh->All]

Map[{#[[1]],#[[2]],ExpectedNumberOfCycles[#[[1]],#[[2]]]}&,{{3,2},{3,5},{3,8},{5,7},{7,3},{15,4},{5,14},{11,7}}]

The start number 1

NestList[If[EvenQ[#],#/2,(9 # + 1)/2]&,1, 20]

NestList[If[EvenQ[#],#/2,(90 # + 1)/2]&,1, 20]

NestList[If[EvenQ[#],#/2,(2 # + 1)/2]&,1, 20]

NestList[If[EvenQ[#],#/2,(16# + 1)/2]&,1, 20]

The 181n+1 problem

ArrayPlot[Table[ExpectedNumberOfCycles[2a+1,k],{a,89,91},{k,1,20}],Mesh->All]

ExpectedNumberOfCycles[181,15]

ProbOfCycle[181,15]

{2^15,181^2}

NestList[If[EvenQ[#],#/2,(181 # + 1)/2]&,27, 15]

The 1093+1 problem

NestList[If[EvenQ[#],#/2,(7 # + 1)/2]&,3, 33]

C1093[n_]:=If[EvenQ[n],n/2,(1093n+1)/2]

NestList[C1093,3,10]

Height

Uaux[n_,q_,a_]:=If[EvenQ[n],Uaux[n/2,q,1],If[a==1,n,Uaux[q*n+1,q,1]]]

U[n_,q_]:=Uaux[n,q,0]

NestWhileList[U[#,3]&,27,#!=1&]

height[n_,q_]:=Length[NestWhileList[U[#,q]&,n,#!=1&]]-1

height[27,3]

Solve[Mod[2^(m+1)-1,3]==0,{m},Integers]

Table[(4^m-1)/3,{m,0,10}]

Solve[Mod[2^m-5,9]==0,{m},Integers]

Table[(2^(5+6m)-5)/9,{m,0,10}]

Limit[(3x)/(4^x-1),{x->Infinity}]

No height-2 numbers for 1093n+1

Solve[Mod[2^(m+1)-1,1093]==0,{m},Integers]

Solve[1093n==2^364-1,{n},Integers]

height[34379397369058857864804839477245555752983459087404499479704167927278448675053737507245365026675479641821555,1093]

Table[Mod[2^(364a)-1,1093],{a,2,10}]

Table[Mod[2^(364a)-1,1093^2],{a,2,10}]

MultiplicativeOrder[2,1093]

MultiplicativeOrder[2,1093^2]

Map[If[MultiplicativeOrder[2,#]==MultiplicativeOrder[2,#^2],#,Nothing]&,Table[Prime[n], {n,2,500}]]

MultiplicativeOrder[2,3511]

MultiplicativeOrder[2,3511^2]

MultiplicativeOrder[3,11]

Table[Mod[3^x,11],{x,1,20}]

MultiplicativeOrder[2,11]

Chapter 5. Circuit cycles

Enumerating all possible circuits

Grid[Table[2^a 2^x-3^x,{x,1,6},{a,0,6}]]

Format[primeFactorForm[n_Integer]]:=CenterDot@@Superscript@@@FactorInteger[n]/. _[x_]:>x

Grid[Table[primeFactorForm[2^a 2^x-3^x],{x,1,6},{a,0,6}]]

ArrayPlot[Table[If[(2^x 2^b-3^x>0)&&(3^x-2^x)>=(2^x 2^b - 3^x),0,1],{x,1,70},{b,1,70}]]

A systematic approach

With[{p=5},ArrayPlot[Table[With[{k=x+y},If[Divisible[2^k-3^x,p]&&Not[Divisible[3^x-2^x,p]],1,0]],{x,1,20},{y,1,20}],{Mesh->Automatic}]]

With[{p=7},ArrayPlot[Table[With[{k=x+y},If[Divisible[2^k-3^x,p]&&Not[Divisible[3^x-2^x,p]],1,0]],{x,1,20},{y,1,20}],{Mesh->Automatic}]]

With[{p=11},ArrayPlot[Table[With[{k=x+y},If[Divisible[2^k-3^x,p]&&Not[Divisible[3^x-2^x,p]],1,0]],{x,1,20},{y,1,20}],{Mesh->Automatic}]]

With[{p=13},ArrayPlot[Table[With[{k=x+y},If[Divisible[2^k-3^x,p]&&Not[Divisible[3^x-2^x,p]],1,0]],{x,1,20},{y,1,20}],{Mesh->Automatic}]]

ArrayPlot[Table[With[{k=x+y},If[Or[Divisible[2^k-3^x,5]&&Not[Divisible[3^x-2^x,5]],Divisible[2^k-3^x,7]&&Not[Divisible[3^x-2^x,7]],Divisible[2^k-3^x,11]&&Not[Divisible[3^x-2^x,11]],Divisible[2^k-3^x,13]&&Not[Divisible[3^x-2^x,13]]],1,0]],{x,1,20},{y,1,20}],{Mesh->Automatic}]

Length[Flatten[Table[With[{k=x+y},If[Or[Divisible[2^k-3^x,5]&&Not[Divisible[3^x-2^x,5]],Divisible[2^k-3^x,7]&&Not[Divisible[3^x-2^x,7]],Divisible[2^k-3^x,11]&&Not[Divisible[3^x-2^x,11]],Divisible[2^k-3^x,13]&&Not[Divisible[3^x-2^x,13]]],1,0]],{x,1,100},{y,1,100}]]]

Length[Flatten[Table[With[{k=x+y},If[Or[Divisible[2^k-3^x,5]&&Not[Divisible[3^x-2^x,5]],Divisible[2^k-3^x,7]&&Not[Divisible[3^x-2^x,7]],Divisible[2^k-3^x,11]&&Not[Divisible[3^x-2^x,11]],Divisible[2^k-3^x,13]&&Not[Divisible[3^x-2^x,13]]],1,Nothing]],{x,1,100},{y,1,100}]]]

3709./10000

Chapter 6. Ready for prime time

Prime factors

FactorInteger[46189]

FactorInteger[3553]

46189/3553

Map[FactorInteger,{20,21,22,23,88,89,90,91,131,132,133,134}]

Prime chaos

ArrayPlot[Table[If[PrimeQ[50*x+y+1],1,0],{x,0,0},{y,0,49}],{Mesh->True}]

ArrayPlot[Table[If[PrimeQ[50*x+y+1],1,0],{x,0,19},{y,0,49}],{Mesh->True}]

Prime patterns

ArrayPlot[Table[If[PrimeQ[50*x+y+1],1,0],{x,10000,10019},{y,0,49}],{Mesh->True}]

ArrayPlot[Table[If[PrimeQ[50*x+y+1],1,0],{x,1000000000,1000000019},{y,0,49}],{Mesh->True}]

Primes are never ending

{2*3*7*43+1,FactorInteger[2*3*7*43+1]}

Counting the primes between 1 and n

Primorial[n_]:=Product[Prime[i],{i,1,n}]

Primorial[3]

ActualPrimesUpToPrimorial[m_]:=Map[If[PrimeQ[#],#,Nothing]&,Range[1,Primorial[m]]]

ActualPrimesUpToPrimorial[4]

Length[ActualPrimesUpToPrimorial[4]]

PossiblePrimesUpToPrimorial[m_]:= pputp[Table[Prime[n],{n,m}],Range[2,Primorial[m]]]

pputp[primes_,result_]:=If[primes=={},result,pputp[Rest[primes],Map[If[(#!=First[primes])&&Divisible[#,First[primes]],Nothing,#]&,result]]]

PossiblePrimesUpToPrimorial[4]

Length[PossiblePrimesUpToPrimorial[4]]

Complement[PossiblePrimesUpToPrimorial[4],ActualPrimesUpToPrimorial[4]]

FactorInteger[143]

MaxPrimesUpToPrimorialPart1[m_]:=(m-1)/Primorial[m]

MaxPrimesUpToPrimorialPart2[m_] := Product[(Prime[i]-1)/Prime[i],{i,1,m}]

MaxPrimesUpToPrimorial[m_]:=MaxPrimesUpToPrimorialPart1[m]+MaxPrimesUpToPrimorialPart2[m]

Table[{i,1.*MaxPrimesUpToPrimorialPart1[i]},{i,1,10}]

Table[{i,1.*MaxPrimesUpToPrimorialPart2[10^i]},{i,1,5}]

Table[{i,1.*MaxPrimesUpToPrimorial[10^i]},{i,1,5}]

Table[{i,1.*MaxPrimesUpToPrimorial[i]*Primorial[i]},{i,1,20}]

Table[Product[(10^i-1)/10.^i,{i,1,n}],{n,1,10}]

How fast do primes thin out?

Counting powers of 2 between 1 and 10^y

Table[With[{x=Range[10^y]},{y,Count[x,n_/;IntegerQ[Log[2,n]]]}],{y,1,6}]

Counting perfect squares

Table[With[{x=Range[10^y]},{y,Count[x,n_/;IntegerQ[Sqrt[n]]]}],{y,1,6}]

Counting primes

Table[With[{x=Range[10^y]},{y,Count[x,n_/;PrimeQ[n]]}],{y,1,6}]

Table[{digits,Length[Table[If[PrimeQ[x],x,Nothing],{x,10^(digits-1),10^digits-1}]]/(10^digits-10^(digits-1)+1.)},{digits,2,6}]

Table[{x,(1/x)*0.4545},{x,2,8}]

Logarithms

SetPrecision[Log[10,500.],10]

SetPrecision[Log[10,50.],10]

jaggedplot= ListStepPlot[{{1,1},{10,2},{100,3},{120,3}},PlotStyle->{Black}];

smoothedplot=Plot[Log[10.,x]+1,{x,1,140},PlotTheme->"GrayColor"];

Show[jaggedplot,smoothedplot]

Table[{x,Log[10.,x]+1,Floor[Log[10.,x]]+1},{x,1,120}]

Log[10.,e]

Chapter 7. Back to the conjecture

If it were a whole number

Log[10,2^1000-1.]

NestList[#/2&,2^8,8]

NestList[(3#+1)/2&, 2^8-1,8]

Factors of Powers

Table[{x,3^x-2^x,FactorInteger[3^x-2^x]},{x, 1, 10}]

Modular arithmetic

Table[Mod[2^x,Prime[p]],{p,3,8},{x,1,14}]

Modular chaos

p=7933

PrimitiveRoot[7933]

Mod[2^1239,7933]

Table[{2^x,Mod[2^(2^x),7933]},{x,1,10}]

Mod[7596*1455*1047*2072*16*4*2,7933]

Position[Table[Mod[2^x,p],{x,1,p}],4621]

Table[Mod[2^x,7933],{x,1,20}]

Internet security

public = 7933

PrimitiveRoot[public]

mySecret = 56

iTextYou = Mod[2^mySecret,public]

yourSecret = 327

youTextMe = Mod[2^yourSecret,public]

shared = Mod[youTextMe^mySecret,public]

shared = Mod[iTextYou^yourSecret,public]

Mod[(2^327)^56,7933]

Solve[Mod[2^x,7933]==2886,Integers]

Solve[Mod[2^x,7933]==5375,Integers]

Where were we?

Table[Mod[3^x,Prime[p]],{p,3,8},{x,1,14}]

How about the denominator?

FactorInteger[2^72-3^24]

Column by column

Table[{x,2^(x+3)-3^x},{x, 1, 10}]

FactorInteger[-217]

Fourth column

Table[{x,Mod[16*2^x,7],Mod[3^x+5,7]},{x,1,10}]

Chapter 8. Is it provable?

Chapter 9. Divisibility by size

Reduce[1.6^x > 2(1.5)^x ,{x},Reals]

Reduce[3^(x/2)>1.6^x,{x},Reals]

A smaller numerator

Table[With[{x=Floor[k/Log[2,3]]},{k, x,3^x-2^x,2^(k-x)-1,2^k-3^x}],{k,5,19}]

Co-primes

a=With[{q=68},Table[If[CoprimeQ[x,y],1,0],{x,1,q},{y,1,q}]]; ArrayPlot[a,ColorRules->{0->White,1->Black}]

a=Total[With[{q=68},Flatten[Table[If[CoprimeQ[x,y],1,0],{x,1,q},{y,1,q}]]/(q*q)]]

a*1.

FactorInteger[3787238]

What a bunch of squares have to do with a circle

Table[{n,Sum[1./i^2,{i,1,n}]},{n,1,10}]

Table[{n,Product[Prime[i]^2./(Prime[i]^2-1),{i,1,n}]},{n,2,10}]

Squarefree numbers

HarmonicQ[n_]:=1==Last[NestWhileList[#/3&,With[{h=Last[NestWhileList[#/2&,n,EvenQ]]},h],Mod[#,3]==0&]]

HarmonicQ[54]

Map[If[HarmonicQ[#],#,Nothing]&,Range[1000000]]

Length[Map[If[HarmonicQ[#],#,Nothing]&,Range[1000000]]]

2^30

A fast co-prime algorithm

CoprimeQ[72893629935748932717391403489324200172323161876213872631,1812372931216312973629812377]

With[{x=72893629935748932717391403489324200172323161876213872631,y=1812372931216312973629812377},NestWhileList[If[#[[1]]>#[[2]],{Mod[#[[1]],#[[2]]],#[[2]]},{#[[1]],Mod[#[[2]],#[[1]]]}]&,{x,y},And[#[[1]]!=0,#[[2]]!=0,#[[1]]!=#[[2]]]&]]

With[{x=72893629935748932717391403489324200172323161876213872631,y=1812372931216312973629812381},NestWhileList[If[#[[1]]>#[[2]],{Mod[#[[1]],#[[2]]],#[[2]]},{#[[1]],Mod[#[[2]],#[[1]]]}]&,{x,y},And[#[[1]]!=0,#[[2]]!=0,#[[1]]!=#[[2]]]&]]

39% of a proof

ArrayPlot[Table[With[{k=b+x},If[CoprimeQ[x,b]&&!CoprimeQ[3^x-2^x,2^k-3^x],1,0]],{x,1,100},{b,1,100}]]

Total[Flatten[Table[If[CoprimeQ[x,b],1,0],{x,1,100},{b,1,100}]]]

FactorInteger[3^11-2^11]

FactorInteger[2^41-3^11]

Notes and references

1./1^2+1/2^2+1/3^2+1/4^2+1/5^2

Sum[1./n^2,{n,1,100}]

Chapter 10. Separation of powers

Plot[2^x 2^3 - 3^x,{x,0, 6},PlotStyle->Black]

A compact Gersonides proof

Grid[{Range[20],Table[Mod[3^k,8],{k,1,20}]}]

Herschfeld rules out 7s

Grid[{Range[20],Table[Mod[2^k,3],{k,1,20}]}]

Grid[{Range[20],Table[Mod[3^k,4],{k,1,20}]}]

Herschfeld tackles 5

Grid[{Range[20],Table[Mod[3^k,64],{k,1,20}]}]

Grid[{Range[20],Table[Mod[2^k,17],{k,1,20}]}]

Contours

Floor[1000000./Log[2,3]]

Diagonal of a square

41/29 > Sqrt[2]

SetPrecision[Sqrt[2.],20]

Circumference of a circle

Sqrt[10.]

Convergents[Sqrt[2],10]

Logarithms are irrational

2^Log[2,3.]

SetPrecision[Log[2,3.],20]

SetPrecision[1/Log[2,3.],20]

Rhin and Ellison

With[{k=1000000,x=630029},(0.00220823)*3^x / x^13.3]

With[{k=1000000,x=630029},2.^k/E^(x/10)]

With[{k=1000000,x=630029},2.^k-3^x]

With[{k=1000000,x=630029},(2.^k-3^x)/((0.00220823)*3^x / x^13.3)]

A penny is 26.5mm across. If 2^k-3^x were the width of a penny, Rhin’s estimate would be 10^491 light years. Our galaxy is about 105,000 light years across

26.5 *10^509 / 1000000 / (9.46 * 10^12) / 3.26 / 10^9

E^0.1

2/1.1

Reality[p_]:=2^Ceiling[p Log[2,3]]-3^p

Baker[p_]:=With[{q=p Log[2,3.]},3^p 1.11 / q^2]

Trivial[p_]:=With[{q=p Log[2,3.]},1.]

Rhin[p_]:= 3^p 0.00220823 / p^13.3

Ellison[p_]:=With[{q=p Log[2,3.]},2.^q E^(-q/10.)]

Dirichlet[p_]:=With[{q=p Log[2,3.]},3^p 1.11 / q]

Baker:     2^q - 3^p > 3^p 1.11a / q^(b-1)               ~ O(3^p / p^c) with c unknown
Rhin:     2^q - 3^p > 3^p 0.00220823 / p^13.3      ~ O(3^p / p^13.3)
Ellison:  2^q - 3^p > 2^q e^(-q/10)                           ~ O(1.8^q), or O(2.56^p)
Dirichlet:    2^q - 3^p < 3^p 1.11 / q                         ~ O(3^p / p) infinitely often

Ellison is the worst asymptotically, but it is tighter than Rhin until p > 571.

Reduce[Ellison[p]>Rhin[p],{p},Reals]

minp=2; maxp=700;

RealityPlot=ListLinePlot[Table[{p,Log[10,Reality[p]]},{p,minp,maxp}],ScalingFunctions->{"Log10","Log10"},PlotStyle->Gray];

BakerPlot=Plot[Log10[Baker[p]],{p,minp,maxp},ScalingFunctions->{"Log10","Log10"},PlotStyle->Black];

RhinPlot=Plot[Log[10,Rhin[p]],{p,minp,maxp},ScalingFunctions->{"Log10","Log10"},PlotStyle->Black];

EllisonPlot=Plot[Log[10.,Ellison[p]],{p,minp,maxp},ScalingFunctions->{"Log10","Log10"},PlotStyle->Black];

DirichletPlot=Plot[Log[10,Dirichlet[p]],{p,minp,maxp},ScalingFunctions->{"Log10","Log10"},PlotStyle->Red];

Show[RealityPlot,EllisonPlot,RhinPlot,PlotRange->All]

5n+1 circuits

Format[primeFactorForm[n_Integer]]:=CenterDot@@Superscript@@@FactorInteger[n]/. _[x_]:>x

Grid[Table[primeFactorForm[2^a 2^x-5^x],{x,1,4},{a,0,6}]]

Germain primes

Table[Mod[x^3,13],{x,1,12}]

Table[Mod[x^3,7],{x,1,12}]

Chapter 11. Let’s take a census

Empirical counts of height-h numbers

Height of any odd number:

T[n_]:=If[EvenQ[n],n/2,3n+1]

v[n_]:=If[EvenQ[n],v[n/2],n]

u[n_]:=If[EvenQ[n],v[n],v[3n+1]]

h[n_]:=If[n==1,0,1+h[u[n]]]

h[27]

Compile the list of heights for first 2^19 odd numbers. This takes a while.

heights=Prepend[Table[{2n+1,h[2n+1]},{n,1,2^18}],{1,0}];

Extract numbers of height h less than or equal to x:

pihx[h_,x_]:=Map[If[#[[2]]==h && #[[1]]<=x,#[[1]],Nothing]&,heights]

pihx[1,2^18]

Less than x = 1000

pihx12=Table[{h,pihx[h,1000]},{h,1,100}];

Length[Sort[Flatten[Map[#[[2]]&,pihx12]]]]

BarChart[Map[Length[#[[2]]]&,Take[pihx12,65]],ChartLabels->Range[65],PlotTheme->"GrayColor",BaseStyle->Directive[FontSize->8]]

Less than x = 10^6

pihx18=Table[{h,pihx[h,10^6]},{h,1,150}];

Missing one! Must have a huge height.

BarChart[Map[Length[#[[2]]]&,Take[pihx18,150]],PlotTheme->"GrayColor"]

Height-2 numbers

Let’s calculate how many height-2 numbers there are, by drawing a chart (separate).

There are two types.
	(a) Ones that satisfy: 3((3x+1)/(2^(2m-1)))+1 = 2^(6n-2)   		odd -> even* -> odd -> 2^n
	(b) And ones that satisfy: 3((3x+1)/(2^(2m)))+1 = 2^(6n+2)		odd -> even* -> odd -> 2^n

For (a):
	((2^(6n-2)-1)/3) * 2^(2m-1) - 1
	---------------------------------------  = x
          		                 3

Change the "=" to "<=", and simplify:
	 (2^(6n-2)-1) * 2^(2m-1) <= 9x + 3

If we ignore the first "-1" and the last "3", we can take the log of both sides and look for solutions in positive integer (n,m) combinations.

Less than x = 2^12

Length[Solve[6n-2+2m-1<=Log[2,2^12*9],{n,m},PositiveIntegers]]

Length[Solve[6n+2+2m<=Log[2,2^12*9],{n,m},PositiveIntegers]]

Length[pihx[2,2^12]]

Less than x = 2^18

Length[Solve[6n-2+2m-1<=Log[2,2^18*9],{n,m},PositiveIntegers]]

Length[Solve[6n+2+2m<=Log[2,2^18*9],{n,m},PositiveIntegers]]

Length[pihx[2,2^18]]

methodSolveH2[x_]:=Length[Solve[2^(6n-2)2^(2m-1)<= 9x + 3,{n,m},PositiveIntegers]]+Length[Solve[2^(6n+2) 2^(2m)<= 9x + 3,{n,m},PositiveIntegers]]

Grid[Transpose[Prepend[Table[With[{x=2^(6y)},{StringJoin["2^",ToString[6y]],If[y>3,"unk",Length[pihx[2,x]]],methodSolveH2[x]}],{y,1,8}],{"x","Actual","Solve"}]]]

Height-3 numbers

18n-14 over, 1 up, 6m-3 over, 1 up, 2p over
18n-14 over, 1 up, 6m-1 over, 1 up, 2p-1 over
18n-8 over, 1 up, 6m-5 over, 1 up, 2p-1 over
18n-8 over, 1 up, 6m-1 over, 1 up, 2p over
18n-2 over, 1 up, 6m-5 over, 1 up, 2p over
18n-2 over, 1 up, 6m-3 over, 1 up, 2p-1 over
18n-10 over, 1 up, 6m-4 over, 1 up, 2p-1 over
18n-10 over, 1 up, 6m over, 1 up, 2p over
18n-4 over, 1 up, 6m-2 over, 1 up, 2p over
18n-4 over, 1 up, 6m over, 1 up, 2p-1 over
18n+2 over, 1 up, 6m-4 over, 1 up, 2p over
18n+2 over, 1 up, 6m-2 over, 1 up, 2p-1 over

Random testing

height3sample[path_,n_,m_,p_]:=If[path==1,(((2^(18n-14)-1)/3 2^(6m-3)-1)/3 *2^(2p)-1)/3,
If[path==9,(((2^(18n-4)-1)/3 2^(6m-2)-1)/3 *2^(2p)-1)/3,
If[path==10,(((2^(18n-4)-1)/3 2^(6m)-1)/3 *2^(2p-1)-1)/3,
If[path==11,(((2^(18n+2)-1)/3 2^(6m-4)-1)/3 *2^(2p)-1)/3,
If[path==12,(((2^(18n+2)-1)/3 2^(6m-2)-1)/3 *2^(2p-1)-1)/3, Nothing]]]]]

height3sample[1,2,2,2]

h[1272582597]

methodSolveH3[x_]:=Length[Solve[2^(18n-14)2^(6m-3)2^(2p)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-14)2^(6m-1)2^(2p-1)<=27*x+12 ,{n,m,p},PositiveIntegers]]+ Length[Solve[2^(18n-8)2^(6m-5)2^(2p-1)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-8)2^(6m-1)2^(2p)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-2)2^(6m-5)2^(2p)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-2)2^(6m-3)2^(2p-1)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-10)2^(6m-4)2^(2p-1)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-10)2^(6m)2^(2p)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-4)2^(6m-2)2^(2p)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n-4)2^(6m)2^(2p-1)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n+2)2^(6m-4)2^(2p)<=27*x+12 ,{n,m,p},PositiveIntegers]]+Length[Solve[2^(18n+2)2^(6m-2)2^(2p-1)<=27*x+12 ,{n,m,p},PositiveIntegers]]

methodSolveH3[64]

Grid[Transpose[Prepend[Table[With[{x=2^(6y)},{StringJoin["2^",ToString[6y]],If[y>3,"unk",Length[pihx[3,x]]],methodSolveH3[x]}],{y,1,8}],{"x","Actual","Solve"}]]]

Counting solutions of Diophantine inequalities

Example:  2x + 3y <= 20.

Bogan-Dov’s 1972 bound on positive integer solutions to linear Diophantine inequality, using the “pyramid” method, also re-invented by myself.

boganDov[coeff_,n_]:=With[{h=Length[coeff]},(n - Total[coeff])^h / (h! Apply[Times,coeff])]

Bound gives lower bound of 18.75 positive integer solutions for this inequality, when there are actually 27.

{boganDov[{2,3},20.],Length[Solve[2n+3m<=20,{n,m},PositiveIntegers]]}

The bound is bad when the number of solutions is small, but good when larger. Note that if n < Total[coeff], then it predicts a negative number of solutions, but more usefully, it identifies that there are no solutions in reality.

Here are undercounted solutions for height-3 numbers, compared with accurate counts.

Table[{k,boganDov[{18,6,2.},k],Length[Solve[18p+6n+2m<=k,{n,m,p},PositiveIntegers]]},{k,1,50}]

Undercounting lattice points

Using the bound of Bogan-Dov.

methodAreaH2[x_]:=boganDov[{6,2},Log[2,x]+2 Log[2,3]+3] + boganDov[{6,2},Log[2,x]+2 Log[2,3]-2]

methodAreaH3[x_]:=boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+17]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+16]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+14]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+9]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+7]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+6]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+15]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+12]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+6]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+5]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+2]+boganDov[{18,6,2},Log[2,x]+3 Log[2,3]+-1]

Closed form for infinite sums

n^k / k!

Sum[n^k/k!,{k,0,Infinity}]

Sum[n^k/k!,{k,0,z}]

With[{n=2000},{n,E^n/2.,Sum[1.*n^k/k!,{k, 0, n}]}]

(n-k)^k / k!

Sum[(n-k)^k/k!,{k, 0, Infinity}]

Sum[(n-k)^k/k!,{k,0, n}]

With[{n=10},Table[1.*(n-k)^k/k!,{k, 0, 10}]]

Total[With[{n=10},Table[1.*(n-k)^k/k!,{k, 0, 10}]]]

@Gary from mathstackexchange worked out the closed form: (1 / (1+omega)) e^(omega n).

omega=ProductLog[1.]

ClosedForm1[n_]:=(1./(1+omega)) E^(omega n)

(1./(1+omega))

Undercounting census teams

Instructions:
	Move to the right 6n
	Go up
	Move to the right 6m
	Go up
	etc.

methodOneteam1[x_,h_]:=(Log[2,x]/6-h)^h / h!

methodOneteam2[x_,h_]:=(Log[2,x]+h Log[2,3]-6h)^h / (6^h * h!)

methodOneteam3[x_,h_]:=2(Log[2,x]+h Log[2,3]-6h)^h / (6^h * h!)

teamsInGroup[h_,i_]:=Binomial[h,i]-Binomial[h-2,i-1]

methodOneteam4[x_,h_]:=Sum[teamsInGroup[h,i]*(Log[2,x]+h Log[2,3]-6h)^h / (6^h * h!),{i,0,h+1}]

methodOneteam5[x_,h_]:=Sum[teamsInGroup[h,i]*(Log[2,x]+h Log[2,3]+2i-6h)^h / (6^h * h!),{i,0,h+1}]

Applying the infinite sum to methodOneteam1 gives O(x^0.13).

exponent=omega *Log[2,E]/6.

Getting rid of the multiplicative constant.

Reduce[(1/(1+omega))x^exponent > x^0.13,{x},Reals]

Reduce[(1/(1+omega))x^exponent > x^0.1,{x},Reals]

(10^1000)^0.13

x^0.13 is approximately sqrt(sqrt(sqrt(x))).

Table[With[{x=2^(6y)},{6y,(1./(1+omega))x^exponent}],{y,1,30}]

Height 2 with one census team: sample numbers found

Table[{n,m,((4*2^(6n)-1)/3*2^(6m)-1)/3},{m,1,5},{n,1,5}]

Height 3 with one census team: sample numbers found (not always whole numbers)

Table[{n,m,p,(((4*2^(6n)-1)/3*2^(6m)-1)/3*2^(6p)-1)/3},{p,1,2},{m,1,2},{n,1,2}]

Undercounting results

h=2

Grid[Transpose[Prepend[Table[With[{x=2^(6y)},{StringJoin["2^",ToString[6y]],If[y>3,"unk",Length[pihx[2,x]]],methodSolveH2[x*1.],methodAreaH2[x*1.],methodOneteam1[x*1.,2],methodOneteam2[x*1.,2],methodOneteam3[x*1.,2],methodOneteam4[x*1.,2],methodOneteam5[x*1.,2]}],{y,1,15}],{"x","Actual","Solver","Area","Oneteam1","Oneteam2","Oneteam3","Oneteam4","Oneteam5"}]]]

h=3

Grid[Transpose[Prepend[Table[With[{x=2^(6y)},{StringJoin["2^",ToString[6y]],If[y>3,"unk",Length[pihx[3,x]]],methodSolveH3[x*1.],methodAreaH3[x*1.],methodOneteam1[x*1.,3],methodOneteam2[x*1.,3],methodOneteam3[x*1.,3],methodOneteam4[x*1.,3],methodOneteam5[x*1.,3]}],{y,1,15}],{"x","Actual","Solver","Area","Oneteam1","Oneteam2","Oneteam3","Oneteam4","Oneteam5"}]]]

Histogram of predicted versus actual counts

It seems impossible to make a single plot, because x^0.13 is too small to see.
Here the true heights are for x = 2^18, but the predicted heights for x = 2^90.

x = 2^18

real18=Table[Length[pihx[h,2^18]],{h,1,150}];

guessed18=Table[If[h>3,0,methodOneteam1[2^18.,h]],{h, 1, 150}];

guessed90=Table[If[h>15,0,methodOneteam1[2^90.,h]],{h, 1, 150}];

BarChart[MapThread[{#1,#2}&,{real18,guessed90}]]

guessed90=Table[If[h>5,0,methodOneteam5[2^30.,h]],{h, 1, 150}];

BarChart[MapThread[{#1,#2}&,{real18,guessed90}]]

Chaos radiates from the number 1

CollatzTree1[n_,cutoff_,height_,depth_]:=
If[(depth==1)||(n>cutoff),n,
If[2n>cutoff,
If[(n>4)&&EvenQ[n]&&(Mod[n,3]==1),
n[CollatzTree1[(n-1)/3,cutoff,If[OddQ[n],height+1,height],depth-1]],n],
If[(n>4)&&EvenQ[n]&&(Mod[n,3]==1),
n[CollatzTree1[(n-1)/3,cutoff,If[OddQ[n],height+1,height],depth-1],
CollatzTree1[2n,cutoff,If[OddQ[n],height+1,height],depth-1]],
n[CollatzTree1[2n,cutoff,If[OddQ[n],height+1,height],depth-1]]]]]

TreeForm[CollatzTree1[1,10000,0,1000],VertexLabeling->Automatic,VertexRenderingFunction->None,EdgeRenderingFunction->({Black,Line[#]}&)]

Notes and references

state1=19000; state2=17000;state3=6000;

apportion[seats_]:=With[{totalpop=state1+state2+state3},{1.*seats*state1/totalpop,1.*seats*state2/totalpop,1.*seats*state3/totalpop}]

Grid[Table[{i,apportion[i]},{i,45,46}]]

Chapter 12. If there’s a cycle, it’s long

Relating cycle lengths with member sizes

(3^8-2^8)/(2^13-3^8)*1.

A[n_]:=(3n+1)/2

B[n_]:=n/2

B[B[B[B[B[A[A[A[A[A[A[A[A[(3^8-2^8)/(2^13-3^8)*1.]]]]]]]]]]]]]

A[A[A[A[A[A[A[A[(3^8-2^8)/(2^13-3^8)*1.]]]]]]]]

B[A[B[A[A[B[A[B[A[A[B[A[A[8.8909]]]]]]]]]]]]]

A formula for the bottom member of the high cycle

BetaAnd[l_,i_,x_]:=If[Last[Take[l]]==1,2^(i-1)*3^(x-Count[l,1]),0]

BetaValue[necklace_]:=1.*Total[Table[BetaAnd[Take[necklace,i],i,Count[necklace,1]],{i,1,Length[necklace]}]]

Rotations[l_]:=Table[RotateLeft[l,i],{i,1,Length[l]}]

InsertBetas[a_]:=Map[{BetaValue[#],#}&,a]

SelectMin[a_]:=First[MinimalBy[a,First]]

Stats[l_]:={Length[l],Length[Select[l,#==1&]]}

Rep[x_,n_]:=Flatten[Apply[Join,Table[{x},{m,1,n}]]]

CircuitLoop[x_,k_]:=Join[Rep[1,x],Rep[0,k-x]]

CircuitLoop[5,8]

AddOnes[x_,l_]:=With[{k=Length[l]},Table[Floor[(k/x)(i-1)+1.],{i,1,x}]]

HighLoop[x_,k_]:=With[{indices=AddOnes[x,Rep[0,k]]},Table[If[MemberQ[indices,i],1,0],{i,1,k}]]

HighLoop[5,8]

HighLoop[12,19]

HighLoop[41,65]

HighLoop[8,13]

BetaValue[HighLoop[8,13]]/(2^13-3^8)

Example:  x=971, k=1539

BetaValue[CircuitLoop[971,1539]]/(2^1539-3^971.)

BetaValue[HighLoop[971,1539]]/(2^1539-3^971.)

Table[(2.^(Ceiling[x*Log[2,3]]+1)-3^x)/(2.^Ceiling[x*Log[2,3]]-3^x),{x,80,100}]

While the numerator increases reliably by about 1.5 times:

Table[BetaValue[HighLoop[x, Ceiling[x*Log[2,3]]+1]]/BetaValue[HighLoop[x, Ceiling[x*Log[2,3]]]],{x,80,100}]

Increasing k from “x log 3” to “x log 3 + 1” causes the altitude of the highest (x)-loop to decrease.

How does the altitude change as we increase k beyond x log 3?

Here's an example with k=32+, x= 20.

loop=HighLoop[20,32]

loop=HighLoop[20,100]

Here is a log10 plot of altitude versus “extra k”.

ListPlot[With[{k=32,x=20},Table[{i,Log[10,BetaValue[Join[HighLoop[x,k+i],Table[0,{j,1,i}]]]/(2^(k+i)-3^x)]},{i,0,40}]]]

When k=40, x=20, we just have a repetition of the trivial 1-2-1 loop, with altitude = 1.

HighLoop[20,40]

BetaValue[HighLoop[20,40]]/(2^40-3^20)    (* same as 1-2-1 loop *)

All loops with x = 1 to 1000

MaximalBy[Table[With[{k=Ceiling[x* Log[2,3]]},{k,x,BetaValue[HighLoop[x,k]]/(2.^k-3.^x)}],{x,1,1000}],#[[3]]&]

Because beta(k,x) for the constructed loop involves floors and ceilings, it is hard to work with.

So we compute an overestimate of beta(k,x), by getting rid of the floors when we distribute the (3n+1)/2 operations among the n/2 operations, and getting rid of the ceiling in selecting k.

Sum[2^(((5 * Log[2,3]+1)/5)*i) 3^(5-i-1),{i,0,4}]

Sum[2^(((5 * Log[2,3]+1)/5)*i) 3^(5.-i-1),{i,0,4}]

BetaOverestimate[x_]:=Sum[2^((x * Log[2,3] + 1.)/x*i) 3^(x-i-1),{i,0,x-1}]

BetaOverestimate[5]

How good is the upper bound for beta(k,x)?

Table[With[{ups=Prime[x]},With[{k=Ceiling[ups Log[2,3]]},{k,ups,BetaExact[ups], BetaOverestimate[ups]}]],{x,1,20}]

This overestimate looks pretty tight.  As x increases, the estimate is somewhere between 1.4x and 2x too high:

Table[With[{ups=Prime[x]},With[{k=Ceiling[ups Log[2,3]]},{k,ups,BetaOverestimate[ups]/BetaExact[ups]}]],{x,80,100}]

Here’s a plot of actual log(beta), in blue dots, versus a solid-line upper bound.

maxc=200;betadotsplot=ListPlot[Table[With[{k=Ceiling[x Log[2,3]]},{x,Log[10,BetaValue[HighLoop[x,k]]]}],{x,180,maxc}],ScalingFunctions->{"Log10","Log10"}];

betaupperboundplot=Plot[Log[10,3^x*x/2. ],{x,1,maxc},ScalingFunctions->{"Log10","Log10"}];

betalowerboundplot=Plot[Log[10,3^x * x/6.],{x,1,maxc},ScalingFunctions->{"Log10","Log10"}];

Show[betadotsplot,betaupperboundplot]

The looseness is the upper bound comes from “getting rid of the floors” being pessimistic.  For every big jump, there’s going to be a corresponding small jump somewhere else in the sum.

Upper bounding a particular loop’s altitude

We can use our overestimate for beta(k,x) in order to calculate an upper-bound the altitude of the (1539,971)-loop.

Our upper bound is ~476,000, compared to the loop’s actual altitude of ~238,000:

BetaOverestimate[971.]/(2^1539-3^971)

BetaValue[HighLoop[971.,1539]]/(2^1539-3^971)

aixtable=Table[Table[If[(x==0)&&(a==0),"x      a->",If[x==0,a,If[a==0,x,Sum[a^(i/x)*1.,{i,0,x-1}]/x]]],{a,{0,1.035,1.05,1.1,1.5,2.,2.0706,3.,3.864,4.}}],{x,{0,10.,100.,1000.,10000.,100000.,1000000.}}];

Grid[aixtable,Alignment->Left,Frame->All,ItemStyle->"Text"]

Note that this sum can be re-written in closed form as  (a-1) / (xth root of a - 1).  I’m not sure if the closed-form is easier to compute, though.

Sum[a^(i/x),{i,0,x-1}]

Fast, incremental computation of altitude bounds

In checking altitudes of highest-loops for x = 1 to 1,000,000, we have recomputed the highest loop (and its altitude) for each x:
   ...
   (3^12345 (12345/2)) / (2^ceiling(12345 log 3) - 3^12345
   (3^12346 (12346/2)) / (2^ceiling(12346 log 3) - 3^12346
   ...

However, given 3^x, we can quickly, incrementally, compute 3^(x+1) = 3^x * 3.
Given x/2, we can quickly compute (x+1)/2 = x/2 + 1/2.
Given x log 3, we can quickly compute (x+1) log 3 = x log 3 + log 3.

We can also keep track of the best altitude so far.  Here’s a 6-tuple state that we can incrementally update:

{ x , 3^x , x/2 , xlog3 , (3^x * x/2) / (2^k-3^x) , highest-altitudes-so-far }

Updater[state_]:=With[{s1=state[[1]]+1,s2=state[[2]]*3,s3=state[[3]]+1/2,s4=state[[4]]+Log[2,3]},With[{s5=(s2*s3)/(2.^Ceiling[s4]-s2)},{s1,s2,s3,s4,s5,If[s5>state[[6]][[1]][[2]],Insert[state[[6]],{s1,s5},1],state[[6]]]}]]

Confirmation

Let’s confirm that the incremental version gives the same results as the non-incremental version, at least for x=1 to 100,000.

BetaOverestimateClosed[x_]:=(x/2) 3.^x

MaximalBy[Table[With[{k=Ceiling[x* Log[2,3]]},{k,x,BetaOverestimateClosed[x]/(2^k-3^x)}],{x,1,100000}],#[[3]]&];

First[Nest[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},100000-1][[6]]]

Empirical detection of high-altitude loops

Let’s start with x = 1 to 100.

res100=Map[{#[[1]],#[[5]]}&,NestList[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},99]]

For each x, we plot the altitude of the highest (x)-loop.

We can see that for loops up to length (x=100, k=159), every loop has a member lower than 5000.  So there are no integer loops shorter than x=5000.

bar100plot=BarChart[Map[#[[2]]&,res100],ChartLabels->Table[If[Mod[i,5]==0,i, ],{i,1,100}],PlotTheme->"GrayColor"];

step100=Nest[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},100-1][[6]]

step100plot=ListStepPlot[step100,Left,PlotRange->{{0,100},{0,5200}},PlotTheme->"GrayColor"];b=1

Show[bar100plot,step100plot,PlotRange->{{0,100},{-200,5200}}]

Show[bar100plot,PlotRange->{{0,100},{-200,5200}}]

(* Nest[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},10000000-1][[6]] *)

(* ListStepPlot[Insert[res10m,{10000000,res10m[[1]][[2]]},1],Left,ScalingFunctions->{"Log10","Log10"}] *)

(* res100m=Nest[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},100000000-1][[6]] *)

(* ListStepPlot[Insert[res100m,{100000000,res100m[[1]][[2]]},1],ScalingFunctions->{"Log","Log"}] *)

The plot continues

maxb=100000

As we increase x, let’s write down each loop that establishes a new “record” for altitude.

res=Nest[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},maxb-1][[6]]

Conclusion so far: no loop with x < 100,000 has a bottom member whose value exceeds 1.1 * 10^10.

So any integer loop has to have at least 100,000 odd terms.

resplot=ListStepPlot[res,Left,ScalingFunctions->{"Log10","Log10"}];

resdots=Map[{#[[1]],#[[5]]}&,NestList[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},maxb-1]];
resdotsplot=ListPlot[resdots,ScalingFunctions->{"Log10","Log10"}];

diophantineplot=Plot[0.71 x^2 ,{x,1,maxb},ScalingFunctions->{"Log10","Log10"}];

bakerplot=Plot[(x/2) 1.11 (x Log[2,3])^2,{x,1,maxb},ScalingFunctions->{"Log10","Log10"}];

rhinplot=Plot[(x/2) x^13.3 / 0.00220823 ,{x,1,maxb},ScalingFunctions->{"Log10","Log10"}];

lowerboundplot=Plot[x/3.,{x,1,maxb},ScalingFunctions->{"Log10","Log10"}];

circuit = Table[(3^x-2^x)/(2.^Ceiling[x Log[2,3]]-3^x),{x,1,maxb}];
circuitplot=ListPlot[circuit,ScalingFunctions->{"Log10","Log10"},PlotStyle->Orange];

ellisonplot=Plot[(x/2) 1.17^x,{x,1,1000},ScalingFunctions->{"Log10","Log10"}];

x-axis:  number of up-moves x
y-axis:  altitude of highest/lowest loop with up-moves <= x

First, the altitudes of highest loops (blue) and circuit loops (in orange).  The solid line shows the cumulative “record-holders” for loops with highest altitude.

Show[resdotsplot,resplot,circuitplot]

resplotbw=ListPlot[res,ScalingFunctions->{"Log10","Log10"},PlotTheme->"GrayColor"];

Show[resplotbw]

Show[resdotsplot,lowerboundplot,ellisonplot,rhinplot]

Show[diophantineplot,resplot]

Pushing it to x=11,000,000. (NB: This takes hours!)

(* maxb=11000000

As we increase x from 1 to 10,000,000, let’s record each loop that establishes a new “record” for altitude.

(* res11m=Nest[Updater,{1,3,1/2,Log[2,3],1,{{1,1}}},maxb-1][[6]]

(* resplot=ListPlot[res11m,ScalingFunctions->{"Log10","Log10"},PlotTheme->"GrayColor"];

(* diophantineplot=Plot[(x Log[2,3.])^2 /2.,{x,1,maxb},ScalingFunctions->{"Log10","Log10"},PlotTheme->"GrayColor"];

(* lowerboundplot=Plot[x/12.,{x,1,maxb},ScalingFunctions->{"Log10","Log10"},PlotTheme->"GrayColor"];

(* rhinplot=Plot[(x/2) x^13.3 / 0.00220823 ,{x,1,20},ScalingFunctions->{"Log10","Log10"},PlotTheme->"GrayColor"];

(* ellisonplot=Plot[(x/2) 1.17^x,{x,1,1000},ScalingFunctions->{"Log10","Log10"},PlotTheme->"GrayColor"];

x-axis:  number of up-moves x
y-axis:  altitude of highest/lowest loop with up-moves <= x

The solid jagged line shows the cumulative “record-holders” for loops with highest altitude.

(* Show[resplot,diophantineplot,lowerboundplot,rhinplot]

Notes and References

12*20*20*18*20+9*20*18+ 13*20+1

13*20*20*18*20

13*20*20*18*20/365.366

5125.6-(3113+(222/365))

Chapter 13. What can compute 3n+1?

De Mol’s tag system

DeMol[pattern_]:=With[{a=Rest[Rest[pattern]],b=First[pattern]},Switch[b,0,Join[a,{1,2}],1,Join[a,{0}],2,Join[a,{0,0,0}]]]

DeMol[{0,0,0,0,0,0}]

ArrayPlot[NestList[DeMol,{0,0,0,0,0,0},10],ColorRules->{0->Black,1->Gray,2->LightGray}]

An amazing assault on the 3n+1 problem

f[n_]:=2n

t[n_]:=2n+1

o[n_]:=3n

i[n_]:= 3n+1

z[n_]:=3n+2

l[n_]:= 1

r[n_]:=n

fm[{x1_,x2_}]:={{1,1},{1,0}}.{x1,x2}

tm[{x1_,x2_}]:={{1,3},{3,4}}.{x1,x2}+{1,1}

om[{x1_,x2_}]:={{7,2},{2,5}}.{x1,x2}+{1,0}

im[{x1_,x2_}]:={{2,1},{1,1}}.{x1,x2}+{1,0}

zm[{x1_,x2_}]:={{2,2},{2,4}}.{x1,x2}+{0,2}

lm[{x1_,x2_}]:={{1,5},{0,0}}.{x1,x2}

rm[{x1_,x2_}]:={2,2}

lm[tm[rm[{0,0}]]]

{{"<t>", r[t[l[0]]], lm[tm[rm[{0,0}]]]},{"<2>",r[z[l[0]]],lm[zm[rm[{0,0}]]]},{"<ft>",r[t[f[l[0]]]],lm[fm[tm[rm[{0,0}]]]]},{"<f2>",r[z[f[l[0]]]],lm[fm[zm[rm[{0,0}]]]]},{"<1f>",r[f[i[l[0]]]],lm[im[fm[rm[{0,0}]]]]},{"<1>",r[i[l[0]]],lm[im[rm[{0,0}]]]},{"<ff>",r[f[f[l[0]]]],lm[fm[fm[rm[{0,0}]]]]},{"<f>",r[f[l[0]]],lm[fm[rm[{0,0}]]]},{"<>",r[l[0]],lm[rm[{0,0}]]}}

A cellular automaton

BaseForm[FromDigits["341",6],10]

KorecTable={{6,6,6,1,1,2,2},{6,0,0,1,1,2,2},{4,3,3,4,4,5,5},{6,0,0,1,1,2,2},{4,4,4,4,4,5,5},{6,0,0,1,1,2,2},{4,3,3,4,4,5,5}}

KorecNext[pair_]:=With[{a=pair[[1]],b=pair[[2]]},Part[Part[KorecTable,If[a==6,1,a+2]],If[b==6,1,b+2]]]

KorecInitial={6,6,6,6,6,1,1,6,6,6,6,6,6,6,6,6,6,6,6,6,6}

listify[a_,b_]:={a,b}

KorecNextState[state_]:=With[{a=AppendTo[Rest[state],First[state]]},Drop[Insert[Map[KorecNext,MapThread[listify,{state,a}]],6,1],-1]]

NestList[KorecNextState,KorecInitial, 20]

The Farkas rule

T[n_]:=If[EvenQ[n],n/2,(3n+1)/2]

Farkas[n_]:=Switch[Mod[n,6],0,n/2,1,(n-1)/3,2,n/2,3,(3n+1)/2,4,(n-1)/3,5,(3n+1)/2]

NestWhileList[Farkas,3,#!=1&]

NestWhileList[T,3,#!=1&]

NestWhileList[Farkas,7,#!=1&]

NestWhileList[T,7,#!=1&]

The original Collatz rule

UnsameQorCutoff[z__]:=And[UnsameQ[z],Length[List[z]]<400]

FindLoops[rule_]:=Table[With[{u=NestWhileList[rule,m,UnsameQorCutoff,All]},If[And[First[u]==Last[u],First[u]==First[MinimalBy[u,Abs]]],u,Nothing]],{m,1,200}]

Original[n_]:=Switch[Mod[n,3],0,2n/3,1,(4n-1)/3,2,(4n+1)/3]

OriginalInv[n_]:=Switch[Mod[n,4],0,3n/2,1,(3n+1)/4,2,3n/2,3,(3n-1)/4]

NestList[Original, 4,10]

FindLoops[Original]

NestList[Original, 8,100]

NestList[OriginalInv, 8,100]

A barely expanding system

NearFlatUniverse[n_]:=Switch[Mod[n,5],0,n/5,1,(2n+3)/5,2,(3n+4)/5,3,(521n+2)/5,4,(n+1)/5]

3126/3125.

NestList[NearFlatUniverse, 3,70]

NestList[NearFlatUniverse, 4,20]

NestList[NearFlatUniverse, 5,20]

NestList[NearFlatUniverse, 7,20]

NestList[NearFlatUniverse, 8,20]

NestList[NearFlatUniverse, 138324323434,100]

LoopsQ[n_,f_]:=Not[DuplicateFreeQ[NestList[f,n,1000]]]

NumbersThatReachALoop2[f_,upto_]:=Table[{n,LoopsQ[n,f]},{n,1,upto}]

TotalNumbersThatReachALoop2[f_,upto_]:=Length[Select[NumbersThatReachALoop2[f,upto],First[Rest[#]]==True&]]

T5[n_]:=If[EvenQ[n],n/2,(5n+1)/2]

Table[{i,TotalNumbersThatReachALoop2[NearFlatUniverse,i]/i*1.},{i,{100,1000}}]

Table[{i,TotalNumbersThatReachALoop2[T5,i]/i*1.},{i,{100,1000}}]

3n+1 without the +1

This is a hand-designed 30-way generalized 3n+1 rule that simulates 3n+1, but does it without any +1’s. All conditions replace n with (Subscript[a, i]) n.
  n = 375/2 n if n = {2, 8} mod 30
  n = 125/2 n if n = {0} mod 30
  n = 5/3 n if n = {3,9,15,21,27} mod 30
  n = 2/5 n if n = {5, 10, 20, 25} mod 30
  n = 3/4 n if n = {6,6,12,16,18,24} mod 30
The other 12 remainders are not used because they are never reached by any start number (2^n).

simulator30way[n_]:=With[{m=Mod[n,30]},If[MemberQ[{2,8},m],(375/2)n,If[MemberQ[{0},m],(125/2)n,If[MemberQ[{15,3,9,27,21},m],(5/3)n,If[MemberQ[{10,20,5,25},m],(2/5)n,If[MemberQ[{4,16,6,12,18,24},m],(3/4)n,0]]]]]]

NestList[simulator30way,2^3,96]

Map[If[IntegerQ[Log[2,#]],Log[2,#],Nothing]&,NestList[simulator30way,2^3,96]]

Map[If[IntegerQ[Log[2,#]],Log[2,#],Nothing]&,NestList[simulator30way,2^31,202215]]

Properties of the 30-way rule

(375/2)^2 (125/2) (5/3)^5(2/5)^4(3/4)^6.

Cycles revisited

OrderInvariantTuples[vals_,n_]:=With[{i={#2,#1,Length@vals}&@@@Partition[Prepend[Array[i,n],1],2,1]},Partition[vals[[Flatten[Table@@{i[[All,1]],Sequence@@i}]]],n]];

OrderInvariantTuples[{375/2,125/2,5/3,2/5,3/4},2]

Flatten[Table[Map[With[{m=Apply[Times,#]},If[m==1,{k,#},Nothing]]&,OrderInvariantTuples[{375/2,125/2,5/3,2/5,3/4},k]],{k,1,30}],1]

Monks’ simulator

monks={1/11,136/15,5/17,4/5,26/21,7/13,1/7,33/4,5/2,7/1}

generalAction[fractions_,n_]:=With[{f=SelectFirst[fractions,IntegerQ[#*n]&]},If[MissingQ[f],0,f*n]]

run[rule_,n_,iter_]:=Map[primeFactorForm,NestList[generalAction[rule,#]&,n,iter]]

monksAction[n_]:=generalAction[monks,n]

NestList[monksAction,2, 10]

NestList[monksAction,1021019, 10]

Map[If[IntegerQ[Log[2,#]],Log[2,#],Nothing]&,NestList[monksAction,2^7,500]]

Length[NestWhileList[monksAction,2^7,#!=2^1&]]-1

Format[primeFactorForm[n_Integer]]:=CenterDot@@Superscript@@@FactorInteger[n]/. _[x_]:>x

Map[primeFactorForm,NestList[monksAction,2^7,100]]

Fractran

universal= {583/559,629/551,437/527,82/517,615/329,371/129,1/115,53/86,43/53,23/47,341/46,41/43,47/41,29/37,37/31,299/29,47/23,161/15,527/19,159/7,1/17,1/13,1/3}

37485 * 2^(2^4)

rununiversal[c_,n_,iter_]:=Map[primeFactorForm,NestList[generalAction[universal,#]&,c*2^(2^n),iter]]

rununiversal2[c_,n_,iter_]:=NestList[generalAction[universal,#]&,c*2^(2^n),iter]

rununiversal2[37485,3,187]

rununiversal[37485,3,3]

rununiversal2[2268945,1,400]

Chapter 14. What can 3n+1 compute?

Rule 110

update[state_]:=MovingMap[If[((Part[#,1]==0)&&(#[[2]]==0)&&(#[[3]]==0))||((#[[1]]==1)&&(#[[2]]==0)&&(#[[3]]==0))||((#[[1]]==1)&&(#[[2]]==1)&&(#[[3]]==1)),0,1]&,
state,{3,Center},0]

ArrayPlot[NestList[update,IntegerDigits[RandomInteger[2^1000],2,1000],500],ColorRules->{0->White,1->Black}, ColorFunction->"Rainbow"]

Engineering a 3n+1 behavior

T[n_,q_]:=If[EvenQ[n],n/2,(q*n+1)/2]

ArrayPlot[Map[IntegerDigits[#,2,50]&,NestList[T[#,3]&,398902141,100]],ColorRules->{0->White,1->Black},ColorFunction->"Rainbow"]

ArrayPlot[Map[IntegerDigits[#,2,50]&,NestList[T[#,3]&,308728191,100]],ColorRules->{0->White,1->Black},ColorFunction->"Rainbow"]

((128 *1)+1)/2187 * 128 - 1

Reduce[2187c - 128d== 1,{c},Integers]

Table[{f, 35 + 128f,598+2187f},{f,0,5}]

163*2^7-1

ArrayPlot[Map[IntegerDigits[#,2,50]&,NestList[T[#,3]&,20863,20]],ColorRules->{0->White,1->Black},ColorFunction->"Rainbow"]

Reduce[5b-3c==1,PositiveIntegers]

Reduce[9b-6c==1,PositiveIntegers]

With[{x=2187,y=128},NestWhileList[If[#[[1]]>#[[2]],{Mod[#[[1]],#[[2]]],#[[2]]},{#[[1]],Mod[#[[2]],#[[1]]]}]&,{x,y},And[#[[1]]!=0,#[[2]]!=0,#[[1]]!=#[[2]]]&]]

Reduce[11c-128d==1,PositiveIntegers]

Reduce[2187c-128d==1,PositiveIntegers]

Reduce[7c-11d==1,PositiveIntegers]

Reduce[11c-7d==1,PositiveIntegers]

Other behaviors

T[n_,q_]:=If[EvenQ[n],n/2,(3n+1)/2]

ArrayPlot[Map[IntegerDigits[#,2,600]&,NestList[T[#,3]&,2^328-1,328]],ColorRules->{0->White,1->Black},ColorFunction->"Rainbow"]

ArrayPlot[Map[IntegerDigits[#,2,600]&,NestList[T[#,3]&,2^328-42,328]],ColorRules->{0->White,1->Black},ColorFunction->"Rainbow"]

ArrayPlot[Map[IntegerDigits[#,2,600]&,NestList[T[#,3]&,2^328-42,2000]],ColorRules->{0->White,1->Black},ColorFunction->"Rainbow"]

A universal computer can simulate itself

Log[10,2^(2^16.)]

Log[2,Log[2,2^(2^16)]]

Practice gliders for 3n+1

T[n_]:=If[EvenQ[n],n/2,(3n+1)/2]

(* Flatten[Table[Table[With[{a=NestWhileList[T,n,(#!=1)&&(#!=m*n)&]},If[Last[a]==1,Nothing,{n,n*m,Length[a]-1,BaseForm[n,m],BaseForm[n*m,m]}]],{n,1,100000}],{m,3,13}],1]

(* Flatten[Table[Table[With[{a=NestWhileList[T,n,(#!=1)&&(#!=m*n)&]},If[Last[a]==1,Nothing,{n,n*m,Length[a]-1,BaseForm[n,m],BaseForm[n*m,m]}]],{n,1,100000}],{m,3,100}],1]

(* Flatten[Table[Table[With[{a=NestWhileList[T,n,(#!=1)&&(#!=m*n)&]},If[Last[a]==1,Nothing,{n,n*m,Length[a]-1,BaseForm[n,m],BaseForm[n*m,m]}]],{n,1,1000000}],{m,3,23}],1]

Map[BaseForm[#,10]&,NestList[T,671,32]]

Map[BaseForm[#,22]&,NestList[T,671,32]]

Map[BaseForm[#,10]&,NestList[T,1342,33]]

Map[BaseForm[#,11]&,NestList[T,1342,33]]

Pay attention to nature

T[n_]:=If[EvenQ[n],n/2,(3n+1)/2]

ToBase10Big[x_]:=IntegerDigits[x,10,55]

myplot=Table[ArrayPlot[Map[ToBase10Big,b=NestWhileList[T,a=RandomInteger[10^50],#>1&]],ColorRules->{0->White},ColorFunction->"GrayLevel"],{x,1,20}];

ResourceFunction["PlotGrid"][{myplot}]

Notes and references

Solve[Mod[n,3]==2 && Mod[n,5]==3 && Mod[n,7]==2,{n},PositiveIntegers]

Epilogue

Solve[18x-7y==35,{x,y},PositiveIntegers]

Table[{k, Floor[k / Log[2,3]],2^k-3^(Floor[k / Log[2,3]]),k^1.8},{k,2,20}]

{Sum[2^i 3^(x-i-1),{i,0,x-1}],Sum[2^i 5^(x-i-1),{i,0,x-1}]}

{Gamma[4],Gamma[4.5],Gamma[5]}

smoothed[x_]:=x+ 1/4-((2x+1)/4 )* Cos[Pi x]

{smoothed[7],smoothed[7.5],smoothed[8]}

NestList[smoothed,7,13]

NestList[smoothed,7.5,20]

Plot[smoothed[x],{x,0,31}]