【gnuplot】コラッツ予想のデータをグラフで可視化 Part 1

こんにちは，Hiro（@hiroloquy）です．

今回は，以前gnuplotで作ったコラッツ予想（Collatz Conjecture）のプログラムを掘り下げてみます．

以前書いた記事では，gnuplotでプログラムを作成し，ターミナル上で計算途中の数字を表示することができました．

【未解決問題】gnuplotでコラッツ予想を作ってみた

この記事はこんな人にオススメ算数，数学が好きな人プログラミングが好きな人『予備校のノリで学ぶ「大学の数学・物理」』 (略...

hiroloquy.com

2020.08.03

しかし，それだけではgnuplotである必要はなく，gnuplotの強みであるグラフ描画を使っていない…
そこで今回，グラフを使って2つのデータ（再帰回数，計算中の最大値）を可視化してみます！

注意 本記事はコラッツ予想自体の解説記事ではありません，gnuplotのグラフ作成の具体例としてコラッツ予想を選んだ程度なので，gnuplotのグラフ描画に関する解説記事だと思ってください．

【復習】コラッツ予想とプログラム

コラッツ予想の中身は次のとおりで，四則演算しか使わないシンプルな未解決問題です．
（以降，2つの計算手順をまとめてコラッツ操作と呼ぶことにします．）

詳しく知りたい方はこちらの動画をご覧ください．

そしてgnuplotで作成したプログラムの動作はこちらのツイートで確認できます．
動画のように，コラッツ操作で求めた数を1秒毎に表示させられます．

pause(1)で表示間隔を1秒にして途中経過を眺められるように

偶数が連続すると気持ちいい pic.twitter.com/cl6TgL3XLH

— Hiro (@Sm_pgmf) July 22, 2020

しかし，これだけではC言語やpythonでやっても変わりません…。そこで，gnuplotに実装されているグラフ描画を活用するために，コラッツ予想に関するデータを可視化してみました．まあpythonには描画用ライブラリがありますが。

今回はコラッツ操作の再帰回数と計算中の$N$の最大値を，グラフで描画するプログラムを作りました．作ったプログラムの解説とともに，描画したグラフを見ていきます．

今回作ったプログラム

早速，作ったプログラムのコードを紹介します．

プログラムの構成・順序は以下のとおり．

1. 変数定義
2. 1からrangeNまでのすべての自然数でコラッツ予想を調べる
3. 調べて得たデータをfileNameで指定したファイル（date.txt）に保存
4. date.txt内の情報を使って3つのグラフ（iterN.png, iterN_marked.png, maxN.png）を描画

CollatzGraph.pltで使っている関数等の簡単な説明と注意事項はこちら．

• l.8 : system "mkdir フォルダ名"
フォルダ名で指定した名前のフォルダを作成可能．既にフォルダがあるならFile existsと表示される．もしフォルダ名に変数を利用したいなら，l.8のようにsprintf("mkdir フォルダ名", 変数)を使用．
  注意 MacOSはこれで大丈夫ですが，Windowsの方はstackoverflowの記事 (make an output directory in Gnuplot) を参照
• l.13 – 36 : do for [i=初期値:最終値:刻み幅] {...}
  for文と一緒．刻み幅は省略すると1．
  補足 C言語ならfor文に終了の条件式を与えるが，gnuplotでは具体的な値を与える．
  注意 C言語ならi++，i--とカウンタiを増減できるが，gnuplotではiを増やすことしかできない．i--にしたいなら{...}内部でj=（初期値+最終値)-iと別のカウンタを作る必要あり．
• l.39-40, 44 : set print "txtファイル" / print 変数 or "文字列" / unset print
printはターミナル上に変数・文字列を表示させる（=printf）ことができるが，printを使う前に，txtファイルをset printするとそのtxtファイルに変数・文字列を記録する（=fprintf）ことも可能．
  注意 set printの有無でprintの役割が異なるため，もしprintfとしてprintを使いたいなら，unset printでfprintfモードを解除する必要あり．
• l.60, 65-66, 74-76 : plot "ファイル名" using (X軸データ):(Y軸データ) ...
  データファイルのグラフを描画する．usingでデータ列を指定．
  このコマンドにはいろんな使い方・オプションがあり，ここでは説明しきれないため，詳しく知りたい方は新潟工科大学 竹野研究室のページを参照

このプログラムでは3つのグラフ（iterN.png, iterN_marked.png, maxN.png）が出力される．以下ではこれらの詳細についてまとめます．また，実行環境は以下のとおり．

• macOS Catalina 10.15.6
• gnuplot version 5.4 patchlevel 0

コラッツ操作の再帰回数iterNに関するグラフ

ここでは1に到達するまでのコラッツ操作の再帰回数iterNを調べます．
各Nでの再帰回数はiterN[ ]およびdata.txtの2列目に保存されており，plotコマンドでdata.txtを呼び出して描画します．

とりあえずiterNを描画する（iterN.png）

l.58 – 60のset output，plotコマンドによりiterN.pngが出力されます．
ここでは，x軸の値を10の累乗（rangeN=1e5の場合，rangeN/10=1e4）で割ることで，目盛りを10を基数とする指数表記にすることができます．

plot filePath using ($1/(rangeN/10)):2 w p pt 7 lt 6 ps pointSize

$\times 10^4$はl.55でlabelを使ってx軸近くに表示させています．

set label 1 at screen 0.99, 0.04 right sprintf("× 10^{%.0f}", log10(rangeN/10)) font ", 16"

$1\leq N \leq 10^5$のときの N – iterN グラフは以下のとおり．
グラフは滅茶苦茶になると想像していたが，異なる自然数$N$で再帰回数が同じ（そうな）箇所がいくつもあった．部分部分で直線になっており，iterN[i]=iterN[i+1]=iterN[i+2]=...みたいな箇所があるのかも？

iterNが連続する箇所のみ赤色にする（iterN_marked.png）

iterN.pngで見られた部分的な直線について少し考えます．
見た目では直線に見えても，プロット点が重なっているだけの可能性もあるので，
iterNが異なる自然数i, i+1で同じかどうかを調べ，その箇所を赤くプロットするようにします．

l.30 – 35のif文で，iterNの変化がないかどうかを"i"と"o"で目印を付けてます．
iterNに変化がなくて直線になる箇所は"i"，変化した箇所は"o"にしています．
この目印はdata.txtの4列目に出力しています．

if(i>=2 && iterN[i] == iterN[i-1]){
  equalN[i-1] = "i"
  equalN[i] = "i"
} else {
  equalN[i] = "o"
}

そして"i"/"o"を使い，iterNが変化する箇所（"o"）では青色（lc 6），変化しない箇所（"i"）では赤色（lc 7）でプロットさせます．

plot fileName using ($1/(rangeN/10)):(stringcolumn(4) eq "o" ? $2 : 1/0) w p pt 7 lc 6 ps pointSize, \
     fileName using ($1/(rangeN/10)):(stringcolumn(4) eq "i" ? $2 : 1/0) w p pt 7 lc 7 ps pointSize

関数・データの一部分だけをプロットするには、三項演算子と1/0を使うのが便利です．1/0は定義されていない数であり、plotコマンド中では無視されます．これを利用すると，関数・データの中でも条件式を満たす箇所のみをプロットすることができます！
（参考：三項演算子を用いた条件分岐／米澤進吾 ホームページ）

lcの番号が何色を表すか（今回なら6→青色，7→赤色）がわからない場合はこちらの記事を参考に

【gnuplot】コマンド “test” を用いた出力サンプルの生成

この記事の内容gnuplotのtestコマンドを使って出力結果のサンプルを生成出力形式 (terminal) とカラーパ...

hiroloquy.com

2020.08.26

$1\leq N \leq 10^5$のときの N – iterN グラフ（色付け）は以下のとおり．

プロット点が大きく，重なって直線に見えているだけかと思っていたが，if文で調べるとiterNがちゃんと連続していることがわかった．

部分的な直線が想像よりも多く，グラフが真っ赤．
ちなみにdata.txtを調べるとiterNの最大値はN=77,031のときに350だった．

赤色の部分が何個あるかについても調べられますが，次回の記事でグラフにしてまとめます．

コラッツ操作中の最大値maxNに関するグラフ (maxN.png)

再帰回数iterNを調べるのはここで終わり，次に最大値maxNについて考える．

コラッツ操作は$1\leq N \leq10^5$の範囲で最大350回も計算されています．
それほど計算される中で$N$がどれだけ大きくなるのか？

そこで，コラッツ操作中に最大値をl.22 – 24で調べ，それをmaxNおよびdata.txtの3列目に保存します．l.69 – 76 でグラフの描画を行い，グラフをmaxN.pngとして出力させます．

ここでのプログラムの個人的ポイントはl.75で，usingで読み取ったデータをwith labelsで表示させる箇所です．

今回出力するmaxN.pngではいくつかピークが立ちます．
そこで，ピークのデータラベルをまとめて描画すれば，ある$N$がかなり大きな数まで増加する様子が確認できます．

ただし，with labelsを使うと各Nのデータラベルが表示されて真っ黒になるので，三項演算子を使って表示するデータラベルを絞ることができます．
今回はrangeN=1e5のときに5（実際は$5\times10^8$）を閾値にして条件式$3/scale>5を書いています．

plot filePath using ($1/(rangeN/10)):($3/scale) w l lt 6 lw lineWidth, \
     filePath using ($1/(rangeN/10)):($3/scale>5 ? $3/scale : 1/0):(sprintf("(%d, %d)", $1, $3)) w labels offset char -2, 1

plot using with labelsはデータラベルのみをプロットし，グラフの線・ 点はプロットされません．グラフを表示するためのplot with line (l.74) も必要．

$1\leq N \leq 10^5$のときの N – maxN グラフは以下のとおり．

グラフを見ると，ピークはいくつか立っているが，大体の$N$ではそこまで大きくなっていない．とはいえ，小さいピークはちらほらあるため，$N$の範囲・ スケールを変えても同様のグラフが得られそう．フラクタル的性質があるのかも？

$1\leq N \leq 10^5$の範囲でのmaxNの最大値はN=77,671のときに1,570,824,736．これは元のNの20,224.1倍！　このときの再帰回数は231回と少ないが，ここまで大きくなるなら奇数時の操作$3N+1$が一時的に連続したのだろう．

iterN_marked.pngではN=77,031のときに計算回数350回と多かったが，コラッツ操作$N/2$と$3N+1$の分量・ 出現タイミングの加減で$N$はそこまで大きくならずに済んだっぽい．

また，iterNの最大値がN=77,031のときに350，maxNの最大値がN=77,671のときに1,570,824,736という結果から，「再帰回数が多いほど到達する最大値は大きい」ということではないらしい．

さいごに

今回は以前gnuplotで作ったコラッツ予想のプログラムを改良し，コラッツ操作の再帰回数iterNと計算中の最大値maxNをグラフで描画するプログラムを作りました．コラッツ予想について少し知れて満足….

ここにまとめきれなかったグラフもあるので，後日まとめます．個人的には，using中で$10^n$で割って指数表記にするのと，with labelsを使うのが初めてだったので，工夫したグラフを描くいい勉強になった．