BB5 (状態数が5,文字数1 の Busy Beaver )

注意：この notebook はやや時間のかかる (5分程度の) プログラムを含んでいます．
notebook 全体でなく，1個ずつ評価することをお勧めします．なお評価は「初期条件を変えなければ不要」です．

遷移関数は Cat on Mat さんからお借りしました．　後で気づきましたが Wiki にもまとまって載っていました．

状態数がnで1文字のBBの書ける最大の文字数をB[n]とします. Cat on Mat さんによると
B[3]=6, B[4]=13, B[5]=4098, B[6]=46×(10^1439)  (WikiによるとB[6]=3.5×(10^18267)) だそう．B[5]，B[6]については正確な数字はまだ分かっていないそうです．

サイトのリンクを張っておきます.

Mathematica組み込みのTuringMachineの仕様は TuringMachine[rule,初期状態,総step数] となります．
ruleは{q,s}→{q’,s’,dir} (q:現在の状態,s:ヘッドの下の文字,q’:次の状態,s’:ヘッドが書き込む文字,dir:ヘッドの進む方向.+1,-1,0)
初期状態は{{(Subscript[q, 0]),pos},tape} ((Subscript[q, 0]):最初の状態,pos:最初のヘッドの位置,tape:最初のテープの状態) となります．

例えば{1,a}→{2,b,+1} は「状態1でheadの文字がa」なら「ヘッドの下に文字bを書いて，ヘッドは右に1つ進み，内部状態は状態2へ移る」という事です．

同じ例は遷移図の方では「《1》→(ab右)→《2》」 と表されます．

tapeの初期状態は変えることができますが，両端に空白がそれぞれ1個以上「最後まで残る」様にして下さい．

補助 Program

ClearAll["Global`*"]

turing[out_,comment_]:=(*TuringMachine の結果の表示に使う*)
Manipulate[
Block[{now,len,tape,control,pos,state,contents,boxes},
len=Length[out[[1,2]]];
now=out[[step]];
tape=now[[2]];
control=now[[1]];
pos=control[[2]];
state=control[[1]];
boxes={Line[{{1,0},{1,2}}],Table[Line[{{i,0},{i+1,0},{i+1,2},{i,2}}],{i,1,len}]};
contents={Table[Text[Style[tape[[i]],Tiny,Bold],{i+.5,1}],{i,1,len}],(*Tinyに変更*)
Green,Polygon[{{pos,-1},{pos+1,-1},{pos+1,-.5},{pos+.5,0},{pos,-.5}}],Black,Text[Style[state,Medium],{pos+0.5,-0.7}]};
Graphics[{boxes,contents},PlotRange->{{0,len+2},{-1.2,2}}]],Style[comment,12,Italic],{step,1,Length[out],1}](*PlotRange->Allでは何故か駄目*)

turing[out_]:=turing[out,""]

Main Program

Import["https://mixedmoss.com/mathematica/turing/jpg/BB5.jpg"]

rule={{s," "}->{2,X,-1},{s,X}->{s,X,-1},{2," "}->{3,X,+1},{2,X}->{2,X,+1},{3," "}->{s,X,-1},{3,X}->{4,X,+1},{4," "}->{s,X,-1},{4,X}->{5,X,+1},{5," "}->{F,X,+1},{5,X}->{3," ",+1}};

これはCat on Mat さんから借りてきた遷移関数ですが，(Wikiにも同様のものがあります) Cat on Mat さん自身は始めはPythonで組んで5分ぐらい掛かり，次にCでプログラムして14秒で4098個を数えて終了したそうです．しかしB[5] はまだ正確には分かっていないそうなので，B[5]はもっと大きい可能性があります．Cat　on Mat さんの所には Python と C のプログラムも置いてあります．

B[5]をMathematicaで計算しようとしたのですが，4098個ということは，テープの長さも少なくとも4098区画は必要です．すると重くなり，MathematicaのTuringMachine コマンドでは１日以上掛けても計算できませんでした．経験上，1日で駄目なものは一週間でも駄目なので，自分で少し改良することにしました．

[1]　default の TuringMachine コマンドをそのまま使用.　しかし終わらない…

まずはテープ全体を全ての区画が表示できる様に500回ぐらいでやってみます.

short=Table[" ",45];
bb5=TuringMachine[rule,{{s,11},short},500];

turing[bb5](*テープの部分をクリックしてdragして下さい.大きくなります.*)

次は終了状態まで達するように実行したい所ですが，危険なので，コメントアウトしてあります.　ちなみに下の回数mは改良型のマシーンを動かして得た回数です.　速いマシーンをお持ちの方は自己責任で評価してみて下さい．私は10万回が一晩かけても終わりませんでした．

(*longtape=Table[" ",8000];
m=11798826;
bb5long=TuringMachine[rule,{{s,1000},longtape},m]*)

[2] TuringMachine コマンドを改良して マシンが終了状態になるまで動かす．

上記のように一晩かけても終わらなかったので自分で少し改良することにしました．基本的なアイデアは現在のヘッドの前後でテープを短く切ります．(下では名前はcutとしてます．幅は 2w+1 です．) このcutに対してTuringMachineコマンドでn(n<=w)ステップ進めます．このときの出力の第2要素( newcut ) を長いテープの元の場所にはめ込んで newtape を作ります．(ヘッドは一度に１つしか進めないのでヘッドはcut からはみ出しません．ですから長いテープでTuringMachineを働かせたのと同じ結果になります．） これをNestWhileListで繰り返します．するとnステップ毎(以下1ブロックと言う)の結果が得られます．これで終了状態Fまで到達できました．(^o^)　個数はやはり4098個でした．回数は11798826回 (約1200万回), 時間は最短で4分30秒掛かりました．(Pythonと同じですね.) しかし私は７年ぐらい前に買った CPUがCore i5-4670　という骨董品を使っているので，今のPC/Mac ならず～と速いと思います．

改良Program

(*{状態q,ヘッド位置h,全tape}={{q,h},{" ",...,X,...,X,...," "}}から,nステップ後のセット(newset)を作る*)
n=100;w=100;(*n≤wでないと駄目.[3]参照*)
turingMachine[{{q_,h_Integer},tape_List}]:=Block[{cut,newset,newcut,newh,newq,newtape},
cut=Take[tape,{h-w,h+w}];
newset=TuringMachine[rule,{{q,w+1},cut},n]//Last;
newq=newset[[1,1]];
newh=h+newset[[1,2]]-(w+1);
newcut=newset[[2]];
newtape=Join[Take[tape,h-w-1],newcut,Take[tape,{h+w+1,Length[tape]}]];
{{newq,newh},newtape}
]

これを終了状態になるまで NestWhileList で繰り返します．私の骨董品のPCでは270秒掛かりました．(ちなみにn=m=50なら306秒,n=m=200なら280秒でした.)
                 !historyの評価は時間とメモリーがかかります!

longtape=Table[" ",8000];
notyetF[x_]:=UnsameQ[x[[1,1]],F];
history=NestWhileList[turingMachine,{{s,1000},longtape},notyetF];

終了状態まで行きました．!(^_^)!   そして X は4098個．総ステップ数 m=11798826 でした．

Count[Last[history][[2]],X]

Length[history]

FirstPosition[TuringMachine[rule,history[[%-1]],n][[All,1,1]],F]//First

m=(%%-2)*n+(%-1)

[3] 結果の分析

各ステップ (n,2n,3n...) 時に於いて，左端と右端のXの位置と個数を求めます．

posX1=Map[FirstPosition[#[[2]],X]&,Drop[history,1]]//Flatten;

posX2=Map[Length[longtape]+1-FirstPosition[Reverse[#[[2]]],X]&,Drop[history,1]]//Flatten;

kosuu=Map[Count[#[[2]],X]&,Drop[history,1]]

Max[kosuu]
Position[kosuu,%]//Flatten

面白いことに個数は単調増加してません．最大個数は6142個もあります．しかも何回も達しています. 全て終了の直前です．でもここでは終了できなかったようです．

さらにXの密度 (密度が1の時は隙間なく埋まっている状態) を求めます. 後のグラフを見ると分かりやすいです．

mitsudo=kosuu/(posX2-posX1+1)

Position[mitsudo,1]//Flatten

Xの範囲, 個数のグラフです．

ListPlot[{Labeled[posX1,"左端のXの位置"],Labeled[posX2,"右端のXの位置"],Labeled[kosuu,"Xの個数"]},AxesLabel->{"ブロック数"}]

Xは右の方にばかり広がっているようです．左には殆ど動いていません.

今度はXの密度のグラフです.　大きな波があります．単調減少ではありません．何度も急速に落下してます．

ListPlot[Labeled[mitsudo,"Xの密度"],AxesLabel->{"ブロック数"}]

次にヘッドの位置の変化を見てみます. ブロック単位の平均速度になります．下のグラフは面白いです．結節点の step数はほぼ公比3の等比数列になっています．

head=Map[#[[1,2]]&,history];

ListPlot[Labeled[head,"ヘッドの位置"],AxesLabel->{"ブロック数"}]

上の曲線と下の曲線の乖離が大きい点は速く動いていて，この２曲線の結節点 (目視では約5 箇所) では殆どヘッドは動いてないはずです．(平均速度の話)

さらにヘッドが「1step毎に」どれだけ動くか調べます．

先のグラフは「ブロックごと(1ブロック=n steps)」 のヘッドの位置の変化でしたが，今度は n×k (kは自然数) 経った時点から，「1回ごと」に 5000回までの変化を調べます．
全部は調べられないので開始直後と，密度=1 の点 （おそらくグラフの結節点と同じ) の前後, さらに 終了点の直前を調べます，
                   !評価にやや時間がかかります!

checklist={0,180,230,117810};ListPlot[TuringMachine[rule,history[[#+1]],5000][[All,1,2]], AxesLabel->{"step数"},PlotLabel->ToString[#]<>"ブロック終了時から"]&/@checklist//GraphicsRow

ブロック数 k=0 の時は，振動が徐々に大きくなり，突然すぼみまた徐々に振動が増えています．これほど分かり易くは無いですが，他の「密度=1 の点」の時も同様です． また終了前はず～と右に動いています．（これが[2]のプログラムで n≤w とした理由です)

各状態の出現頻度を見てみます．最初はブロック単位，次はステップ単位です．

history[[All,1,1]]

Tally[%]

圧倒的に状態sと状態2が占めます．しかもほぼ同数です． step単位で見ると次の様になります．開始から 500 steps(=5 blocks) の結果です.

bb5[[All,1,1]]

Tally[%]

終了直前の50ステップの様子です．”kiritori” を使って，ステップ単位の変化を見る事ができます．

(*kiritori[{開始ステップ，終了ステップ}］．*)

kiritori[{start_Integer,end_Integer}]:=Block[{grp,motodata,heads,hmin,hmax,modified},
grp=Quotient[start,n];
motodata=Drop[TuringMachine[rule,history[[grp+1]],end-grp*n],start-n*grp];
heads=motodata[[All,1,2]];
hmin=Min[heads];
hmax=Max[heads];
modified={#[[1]]-{0,hmin-2,0},Take[#[[2]],{hmin-1,hmax+1}]}&/@motodata;
turing[modified,"「step」は,"<>ToString[start]<>"steps からの増加量. テープは"<>ToString[hmin-1]<>"~"<>ToString[hmax+1]<>"区画."]
]

kiritori[{m-49,m}](*mは総ステップ数*)