「くどー問題」再び

くどー問題とその数値計算

Mathlogでくどーさんがとある数学の問題を記事にされています：

私はこれを「くどー問題」と呼んでいます。問題の内容は上記記事を見ていただくのがよいと思いますが、一応引用させていただきます：

$1$から$n(\ge 3)$までの整数が書かれたカードが1枚ずつある。この時、以下の条件を満たしながらカードを$1$枚ずつ取り除いていくことを考える。
$$
$\underset{―}{\mathrm{条}\mathrm{件}}$

1. どの時点においても「カードに書かれた数の$3$乗和」が「カードに書かれた数の和」で割り切れる
2. $1$と$n$のカードは最後の$2$枚になるまで残す
   $$

この時、条件を満たしながらカードがなくなるまで取り除くことができるような$n$は無数に存在するか？

くどーさんは条件を満たす$n$として$3,6$を見つけています。また$n\le 10$では条件を満たす$n$は存在しないことを数値計算で確かめてもらった記憶がある、と述べています。

これに関し私は以前数値計算を行いました： 「くどー問題」の数値計算

この記事では$3\le n\le 500$の範囲で数値計算を行い条件を満たす数を探索しましたが、$3,6$以外見つけることはできませんでした。

今回はJuliaというコンピュータ言語を用いて$3\le n\le 5000$までの範囲の殆どの数（※20個を除く）に関して検証を行いました。結果は以下です：

本記事ではこれに関して説明します。前回同様、今回も整数論的考察はないです。また前回の記事は読まなくても大丈夫です。

用語の定義

まず議論に必要な用語を定義します。

以上からくどーさんの提示した問題は

数値計算の概要

計算アルゴリズムは単純です。集合$\{1,2,\dots ,n\}$から要素を1つづつ取り除いて条件を満たすか確認し、$\mathrm{Remain}(n;m)$を各$m$で逐一構成することで$\mathrm{mmax}(n)$を計算します。

前回の記事ではすべての$\mathrm{Delete}(n;m)$を計算しました。しかしそれだと計算量が膨大になります（前章の$\mathrm{Delete}(19;15)$と$\mathrm{Remain}(19;15)$を比較してみてください）。そこで今回は基本的に$\mathrm{Remain}(n;m)$のみ計算し、そこに至るまでの経路（$=\mathrm{Delete}(n;m)$）は必要な場合だけ計算しました。　

結果

$\mathrm{mmax}(n)$の$n$依存性

まず$\mathrm{mmax}(n)$を$n$に対してプロットしてみます。

$\mathrm{mmax}$と$n$の関係

図1は各$n$に対する$\mathrm{mmax}$のグラフです。青の点が$\mathrm{mmax}(n)$、緑の線は$\mathrm{mmax}(n)=n-2$のグラフです。この線に青点が乗れば、その$n$は最後までカードを削除できる数字です。オレンジの上向き三角形は探索が終了していない$n$に対する$\mathrm{mmax}$の下限を示しています。つまりこれらの$n$では$\mathrm{mmax}$はこれより大きいことしか判明していません。

図1の詳細は下のリンクのグラフからアクセスすることができます：

このグラフは拡大・縮小することができ、かつ各データポイントにカーソルを重ねればそのデータの$(n,\mathrm{mmax}(n))$が表示されます（計算の終わってないデータでは$(n,inf(\mathrm{mmax}(n)))$が表示される）。

これらより以下がわかります：

これに関してコメントです：

• 1.に関して

  $7\le n\le 5000$の間では以下の20個の$n$（図1のオレンジ三角）:

  $1322,1456,1567,1586,1979,2203,2261,2522,2639,2882,3135,3266,3503,3535,3551,3699,3799,3815,4562,4928$

  を除き、条件を満たす数字は存在しないことがわかっています。上記の$n$では多くの削除のパターンが存在するため計算量が非常に大きくなり探索が完了していません。

• 2.に関して

  以下に$\mathrm{mmax}=0,0,2,2,\dots$の繰り返しが起こる典型的な図を示します：

$\mathrm{mmax}=0,0,2,2,\dots$の構造を持つ部分の拡大図

このような周期的な変化は、$n$がある程度大きくなれば普遍的に現れます。よろしければインタラクティブな図を拡大して確認してみてください。

• 3.に関して
• $\mathrm{mmax}=n/2$付近に存在する$n$は、$|\mathrm{Remain}(n;m)|$および$\mathrm{Delete}(n;m)$も特徴的な性質を持ちます。これに関しては次章で述べます。
  ギャップに関しては、探索が終了していない$n$の$\mathrm{mmax}$が0から$n/2$の中間付近の値を持つ可能性はあり、そうであればギャップはそれほど明確ではなくなります。しかしながら、探索の終了していない$n$は$|\mathrm{Remain}(n;m)|$がある$m$から指数的に大きくなり、経験的にはその場合$n/2$より有意に小さな$\mathrm{mmax}$をもつことはないです。なのでギャップの存在はロバストであるように思います。あくまで予想ではありますが。

2.3.は何らかの理論的背景を持つと思うのですが、いまのところよくわかりません。

「カードの削除数」と「残りのカードの可能な組み合わせの数」の関係

次に$m$と$|\mathrm{Remain}(n;m)|$の関係を見てみます。すなわち各$n$における「カードの削除数」と「残りのカードの可能な組み合わせの数」との関係を調べてみます。

$n$がある程度大きい時、$m$ vs $|\mathrm{Remain}(n;m)|$は次の4つのパターンに分けられます：

1.は図1の$\mathrm{mmax}\simeq 0$に存在する点です。殆どの$n$はこれに分類されます。
2.の典型的な$n$は例えば$n=4407$（$\mathrm{mmax}=2146$）であり、$m$ vs $\mathrm{Remain}(4407;m)$をプロットすると次の図のようになります：
$n=4407$における$m$ vs $|\mathrm{Remain}(n;m)|$
図3のように殆どの$m$に対して$\mathrm{Remain}(4407,m)$は$1$、時に$2$になることがあり、そしてそのような小さい値を保ったまま最終的にゼロになります（※ゼロの点は示されていません。一番右側の点が$\mathrm{mmax}=2146$）。

$\mathrm{mmax}=n/2$より有意に大きな$\mathrm{mmax}$をもつ$n$の場合、2.と同様に$\mathrm{Remain}$が小さい値を保ったままゼロになる場合と、途中で指数的に大きくなる場合の2つに分けられます。3.は前者であり、例えば$n=1035,1961$はその例です。$n=1961$の図をプロットすると以下のようになります：
$n=1961$における$m$ vs $|\mathrm{Remain}(n;m)|$

これは図3と同様の振る舞いです。$n=4407$よりは$|\mathrm{Remain}|$が2になる回数は多いですが、結局小さい値を保ったままゼロになります。

一方4.の例は$n=703,735,799,839,935,1203,1379,3959,4255$などです。この場合$\mathrm{mmax}$が指数関数的に増加したのちにゼロになるまでの振る舞いに違いがあります。そのうち典型的な$n=703,1203,4255$を図示します（$y$軸logプロットです）：

$n=703$における$m$ vs $|\mathrm{Remain}(n;m)|$ ($y$軸 log scale)
$n=1203$における$m$ vs $|\mathrm{Remain}(n;m)|$ ($y$軸 log scale)
$n=4255$における$m$ vs $|\mathrm{Remain}(n;m)|$ ($y$軸 log scale)

図5,6,7に対応するインタラクティブなグラフへのリンクも貼っておきます：

• log plot
  $n=703$, log plot
  $n=1203$, log plot
  $n=4255$, log plot
• normal plot
  $n=703$, normal plot
  $n=1203$, normal plot
  $n=4255$, normal plot

$n=703$では指数的増加ののちすぐに$\mathrm{mmax}$がゼロになります。$n=1203$では一度小さくなってからもしばらくある程度の値を保ち、そして最後は急激に小さくなってゼロになります。$n=4255$では指数増加する領域が2つあり、小さくなってからもある程度の値を保つ領域を経てゼロになります。

exponentの計算

ここで指数的増加のexponentを求めてみます。指数的増加を示す領域において2点をとりexponent $b$（$\sim e^{bm}$で振る舞う）を計算します。そのような領域が複数存在する場合はそれぞれの領域で別々に$b$を求めます。さらに$\alpha :=e^b$とすると、$\alpha$は$m$に対する$|\mathrm{Remain}(n;m)|$の平均的な増加率：$|\mathrm{Remain}(n;m+1)|/|\mathrm{Remain}(n;m)|$の指標となります。 以下の表が$b,\alpha$の計算結果です：

$$\begin{array}{cccccc}n& \mathrm{region}& (m_1,|\mathrm{Remain}(n;m_1)|)& (m_2,|\mathrm{Remain}(n;m_2)|)& b& \alpha \\ 703& & (100,53)& (318,83938)& 0.03380& 1.034\\ 935& 1& (421,34262)& (442,63172)& 0.02913& 1.030\\ & 2& (445,35519)& (469,57558)& 0.02011& 1.020\\ & 3& (469,57558)& (498,98108)& 0.01839& 1.019\\ 1203& & (498,48)& (560,393)& 0.03391& 1.034\\ 3959& 1& (1810,27)& (1868,89)& 0.02057& 1.021\\ & 2& (2017,160)& (2120,689)& 0.01418& 1.014\\ 4255& 1& (1900,100)& (2070,1373)& 0.01541& 1.016\\ & 2& (2204,3108)& (2382,33128)& 0.01329& 1.013\end{array}$$

$b,\alpha$を求めるのに2点$(m_1,|\mathrm{Remain}(n;m_1))$と$(m_2,|\mathrm{Remain}(n;m_2))$を用いています。"region"は指数的増加領域が複数ある場合の各領域を示しています。$\alpha$の値より、指数増加領域では$m$が1増えるとき$|\mathrm{Remain}(n;m)|$は平均してだいたい1から3%程度増加することがわかります。

注意してほしいのは、この増加率は$|\mathrm{Delete}(n;m)|$の増加率とは違うことです。前回の記事で調べた$n$の範囲（$n\le 500$）では$|\mathrm{Delete}(n;m+1)|/|\mathrm{Delete}(n;m)|$は1.2から1.4程度でした。今回はちゃんと計算していませんが、$n$が$1000$から$5000$くらいでもおそらく同じようなものだと思います。これが示しているのは、違う$\mathrm{Remain}(n;m)$の要素の違う数字の削除が同じ$\mathrm{Remain}(n;m+1)$の要素をもたらすことがけっこうある、ということです。実際今回の数値計算において$\mathrm{Remain}(n;m)$の要素に対して削除作業をすると同じ$\mathrm{Remain}(n;m+1)$の要素に行き着くものが多くあるので、そのredundancyを取り除く作業をしています（「定義1」の4.と5.に示した$\mathrm{Remain}(19;15)$と$\mathrm{Delete}(19;15)$を比較してみてください）。

$\mathrm{Delete}(n;m)$の$m$依存性

これは前回の計算でもわかっていたのですが、$\mathrm{mmax}$がある程度大きくなるような$n$では$\mathrm{Delete}(n;m)$は非常に特徴的な振る舞いをします。言葉で言うより$\mathrm{Delete}(n;m)$を見たほうが早いので、幾つか具体的に例を示してみます。まず$n=19,59,119$の$m=\mathrm{mmax}$における$\mathrm{Delete}(n;\mathrm{mmax})$（の一部）を示します：

      - n=19, m=14(=mmax):
[15, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 2, 18] 
[15, 4, 5, 6, 7, 8, 9, 16, 18, 12, 3, 17, 10, 11, 13] 
[15, 4, 5, 6, 7, 8, 9, 16, 18, 12, 3, 17, 10, 11, 14]
[15, 4, 5, 6, 7, 8, 9, 16, 18, 12, 3, 17, 10, 14, 11]
[15, 4, 5, 6, 7, 8, 9, 16, 18, 12, 10, 17, 3, 11, 13]
[15, 4, 5, 6, 7, 8, 9, 16, 18, 12, 10, 17, 3, 11, 14]
[15, 4, 5, 6, 7, 8, 9, 16, 18, 12, 10, 17, 3, 14, 11]
[15, 4, 5, 6, 7, 9, 8, 10, 11, 12, 13, 14, 16, 2, 18]
[15, 4, 5, 6, 7, 9, 8, 16, 18, 12, 3, 17, 10, 11, 13]
[15, 4, 5, 6, 7, 9, 8, 16, 18, 12, 3, 17, 10, 11, 14]
[15, 4, 5, 6, 7, 9, 8, 16, 18, 12, 3, 17, 10, 14, 11]
[15, 4, 5, 6, 7, 9, 8, 16, 18, 12, 10, 17, 3, 11, 13]
[15, 4, 5, 6, 7, 9, 8, 16, 18, 12, 10, 17, 3, 11, 14]
[15, 4, 5, 6, 7, 9, 8, 16, 18, 12, 10, 17, 3, 14, 11]

 - n=59 m=21(=mmax)
[30, 29, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] 
[30, 29, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 53] 
[30, 29, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 27, 26, 28] 
[30, 29, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 27, 26, 53] 
[30, 29, 10, 11, 52, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28] 
[30, 29, 10, 11, 52, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 38] 
[30, 29, 10, 11, 52, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 32, 58, 57, 56, 55] 
[30, 29, 10, 11, 52, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 54, 42, 57, 55, 44] 
[30, 29, 10, 11, 52, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 54, 42, 57, 55, 56] 
[30, 29, 10, 11, 52, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 54, 42, 57, 56, 43] 
...

- n=119, m=51(=mmax)
[60, 10, 11, 12, 13, ..., 56, 57, 58, 59] 
[60, 10, 11, 12, 13, ..., 56, 57, 58, 94] 
[60, 10, 11, 12, 13, ..., 37, 38, 39, 94, ..., 56, 57, 58, 40]
...
    

$n=19$ではすべての$\mathrm{Delete}(n;\mathrm{mmax})$、$n=59,119$では一部の$\mathrm{Delete}(n;\mathrm{mmax})$を示しました。これをみると

次に$\mathrm{Delete}(n;m)$の$m$依存性を$m$が小さな領域で見てみます。例えば$\mathrm{mmax}\simeq n/2$上に存在する$n$である$n=1239$に関してこれを見てみると以下のようになります：
$$\begin{array}{rl}& m=1:\mathrm{Delete}(1239;1)=\{[620]\}\\ & m=2:\mathrm{Delete}(1239;2)=\{[620,31]\}\\ & m=3:\mathrm{Delete}(1239;3)=\{[620,31,32]\}\\ & m=4:\mathrm{Delete}(1239;4)=\{[620,31,32,33]\}\\ & ⋮\\ & m=275:\mathrm{Delete}(1239;275)=\{[620,31,32,33,\dots ,303,304]\}\\ & m=276:\mathrm{Delete}(1239;276)=\{[620,31,32,33,\dots ,303,304,305],[620,31,32,33,\dots ,303,304,1015]\}\\ & ⋮\end{array}$$
これも最初の削除可能な数字が$620\simeq 1239/2$であり、次に$31,32,...$と続き、実に$m=275$まで連続した数字が並びます。しかも$m=275$までは可能な削除列は1つだけであり、$m=276$になって初めて2つ現れます。

これらの傾向は$\mathrm{mmax}$が大きくなる$n$に関してかなり普遍的に当てはまります。またこれはなんとなくなのですが、そのような$n$は一桁目が$9$や$7$であることが多いように見えます（インタラクティブなグラフで確認してみてください）。

まとめ

くどーさんのご提案された問題: そろそろ誰かに解いてほしい自作問題(未解決) をJuliaによる数値計算により調べました。その結果以下の結果を得ました：

これらの現象の理論的背景はいまのところわかりません。

くどーさんのオリジナルの問題「条件を満たしながらカードがなくなるまで取り除くことができるような𝑛は無数に存在するか？」ですが、これに関しても今のところ何も言えないです。数値計算でこれにアプローチするには、条件を満たす$n$を数多く求め、ある$n$までで条件を満たす数の出現頻度を求めることが必要かと思います。また条件を満たす数の何らかのルールを見出すことが重要です。しかし今のところ3,6以外で条件を満たす数が求まっていないです。まずは6より大きな数で条件を満たすものを見つけるのが目標かと思います。ただ数値計算では、アルゴリズムに何らかの根本的な工夫をするか、または大規模計算機を使わないとこれ以上はなかなか厳しいかもしれません。

あと逆方向に探索することもできることを付記しておきます。これは$\{1,n\}$のカードの集合にカードを付け加えながら3乗和が1乗和で割り切れるか調べ、最終的に集合が$1$から$n$のカードすべてを持つかどうかを確認する探索法です。前回の記事ではこれを行ったのですが、非常に多くのパターンが現れるため数値計算が重くなりあまり有効ではありませんでした。ただ、前回は違う「付け加え」で結局同じ状態に行き着く場合の重複を除かずに計算していたので、これを取り除くと計算が軽くなるかもしれません。計算が非常に重くなる$n$に対しては、両方向から計算し、$m\simeq n/2$のところで両者を比較することで完全なカードの取り除きが可能か調べるのもよいかもしれません。

個人的な感想ですが、この問題が何らかの物理系と対応がつくと面白いように思います。テキトーな話ですが、「割り算の余り」と関係するといえば磁場中の量子系のエネルギー準位が思い浮かびます。Aharanov=Bohm効果、Diracの量子化条件のように、磁場中の半古典的な量子化条件では「割り切れる」ことが重要です。また量子ホール効果のエネルギー準位にも「余り」が登場します。しかしながら、$1$から$n$の3乗和が現れ、それを1乗和で割った余りがエネルギー準位や量子化条件等に現れるような物理系は存在するでしょうか...。

おしまい。${}_{◼}$

Appendix

Juliaのコード

Juliaは2012年に最初のバージョンがローンチされた比較的新しいコンピュータ言語ですjulia。Pythonと似通った言語でありながらCやfortranに匹敵するほど計算が高速であるという特徴を持ちます。Pythonを使える人ならJuliaもすぐに使えるようになると思います。

参考までに今回使用したコードをGistにアップロードしましたjulia_code。使い方を以下に説明します

      julia> file_open = h5open("KudoProb_n-27_mmin-1_mmax-100_finished.dat","r")
🗂️ HDF5.File: (read-only) KudoProb_n-27_mmin-1_mmax-100_finished.dat
├─ 🔢 list_length_pos_list
└─ 🔢 m_final
    

      julia> list_length_pos_list = read(file_open,"list_length_pos_list")
2×9 Matrix{Int64}:
1  2  3  4  5  6  7  8  9
1  1  1  2  1  1  1  1  0
    

      julia> m_final = read(file_open,"m_final")
9
    

      julia> file_open2=h5open("KudoProb_n-143_mmin-1_mmax-100.dat","r")
🗂️ HDF5.File: (read-only) KudoProb_n-143_mmin-1_mmax-100.dat
├─ 🔢 list_length_pos_list
├─ 🔢 m_final
└─ 🔢 pos_list
    

      julia calc_KudoProb_parallel.jl (n_threadの値) (nの値) (mmin_fileの値) (mmax_fileの値) (mmaxの値)