BBS-BBS/47

トップ 差分 一覧 Farm ソース 検索 ヘルプ RSS ログイン

キミならどう書く 2.0 - ROUND 2 - まとめ - さいとう (2006年07月14日 23時33分15秒)

キミならどう書く 2.0 - ROUND 1 - はどうだったでしょうか?個人的には LL Ring のページに「なぜ解説がないのか?」と思ってしまうのですが、単なるアルゴリズムで終わらないように前回同様に今回も詳しく見ていきます。

今回も、

ライブラリなど花拳繍腿、自作こそ王者の技よ

というわけでスクラッチばかりです。

キミならどう書く 2.0 - ROUND 2 - のお題は「Collatz 予想」または「角谷予想」や「3x+1 問題」と呼ばれているものです。

日本語の Wikipedia に記載がないので、英語版で検索してみると以下の URL に詳細が記載されています。

つまり、Collatz は任意の自然数 n に対して、

  • n が偶数の時は 2 で割る
  • n が奇数の時は 3 をかけて、1 を足す

という操作を行っていくと、必ず 1 に収束する、という予想をしたというもので、未解決問題集のひとつとされています。

下記の URL を見れば、膨大な計算が現在も続けられていることも分かります。

さて、上記の n に対して 100 を実行した際の最大ステップとその際の n を求めよという解釈で進めます。

普通に作ってみる

自己再帰で組むのが楽そうなので、作ってみました。

#! /usr/bin/gawk -f
BEGIN {
  if ( ARGV[1] == "" ) {
    num = 100;
  } else {
    num = ARGV[1];
  }
  for ( i = 1; i <= num; i++ ) {
    count = 0;
    collatz( i );
    if ( count > max ) {
      max = count;
      max_num = i;
    }
  }
  printf( "%d max=%d\n", max_num, max );
}
function collatz( num ) {
  if ( num == 1 ) {
    result = 1;
    count++;
  } else {
    if ( num % 2 == 0 ) {
      result = collatz( num / 2 );
    } else {
      result = collatz( 3 * num + 1 );
    }
    count++;
  }
  return result;
}

少し荒削りな部分もありますが、簡単に作れます。

  • ARGV でも値 n を読み込めるようにします。
    • 100 だけじゃつまらないので、100 以外も簡単に指定できるようにしておきます。
  • 1〜100 は for ループで回します。
    • この際に最大ステップも求めています。
  • 関数として Collatz の式を定義しています。
    • ここでステップ数もカウントしています。

実行してみましょう

答えは、97 の時に 119 ステップです。バックグラウンドでいろいろ動作している中で計算させましたが、非常に高速に結果を出しています。

$ time gawk -f collatz.awk 100
97 max=119

real    0m0.045s
user    0m0.016s
sys     0m0.028s

計算速度に関しては面白いことが分かります。次に n を 10, 100, 1000, 10000, 100000 で計算させた時の時間を求めてみましょう。

Time n = 10 n = 100 n = 1000 n = 10000 n = 100000
real 0m0.031s 0m0.055s 0m0.428s 0m6.308s 1m6.532s
user 0m0.000s 0m0.016s 0m0.336s 0m4.896s 1m2.028s
sys 0m0.020s 0m0.036s 0m0.036s 0m0.060s 0m0.120s

100000 を超えたあたりから莫大な計算時間が発生していることが分かります。

つまり、この「Collatz の予想」の面白いところでもあるのですが、n が 1 に収束する場合でも、途中で非常に大きな値をとって、収束に時間のかかることがあるという可能性があります。

ここからが面白いところ?

では、実際にどうなっているのでしょうか?少しの改良で途中の結果をログとして出力させることができますので、以下のように変更して実行してみることにします。

#! /usr/bin/gawk -f
BEGIN {
  if ( ARGV[1] == "" ) {
    num = 100;
  } else {
    num = ARGV[1];
  }
  for ( i = 1; i <= num; i++ ) {
    count = 0;
    collatz( i );
    if ( count > max ) {
      max = count;
      max_num = i;
    }
  }
  printf( "%d max=%d\n", max_num, max );
}
function collatz( num ) {
  printf( "%d ", num );  # ←
  if ( num == 1 ) {
    result = 1;
    count++;
    printf( "\n" );      # ←
  } else {
    if ( num % 2 == 0 ) {
      result = collatz( num / 2 );
    } else {
      result = collatz( 3 * num + 1 );
    }
    count++;
  }
  return result;
}

矢印の部分の printf を追加してみました。

さて、ここで得られた結果で n = 97 の部分をグラフにしてみます。縦軸には最大値、横軸にはステップ数を取っています。

アッサリと収束するどころか、10000 に近い値になっていることが分かります。この最大値は n^a (a=2.1...) を超えないと予想されています。つまり 97^2.1 = 14866.91 を超えてないことになります。

さて、本当なのでしょうか?

先ほどのログから 9232 が最大であることから、確かに超えません。

実は、これが証明できれば、無限に発散することがないため、無限ループにならない限り、この Collatz の予想は正しいことになるわけです。

32 bit を超えて

この話を書いていると、小飼さんの Blog で以下のようなものがありました。

木村さんが言うには、gawkdobule で数値を扱っていて、double の mantissa は 52bit (+けちビット) あるので、通常の gawk でも 113383 は計算できます。

"One True AWK" で調べると確かに 113383 で引っかかるので、以前は 32 bit 整数でリミットしていたものと思われます。

ずーっと 113383 まで計算させるのが面倒なので、以下のように改造してみます。

#! /usr/bin/gawk -f
BEGIN {
  if ( ARGV[1] == "" ) {
    max = 100;
  } else {
    max = ARGV[1];
  }
  count = 0;
  max_num = 0;
  collatz( max );
  printf( "%d max=%d(%d)\n", max, count, max_num );
}
function collatz( num ) {
  if ( max_num < num ) {
    max_num = num;
  }
  if ( num == 1 ) {
    result = 1;
    count++;
  } else {
    if ( num % 2 == 0 ) {
      result = collatz( num / 2 );
    } else {
      result = collatz( 3 * num + 1 );
    }
    count++;
  }
  return result;
}

計算させると、

$ gawk -f collatz.awk 113383
113383 max=248(2482111348)
               ~~~~~~~~~~

ちなみに、同じく小飼さんの Blog に出ている以下のものも計算できます。

ここに出ているものくらいであれば、gawk で計算可能範囲です。

篩を考える

先のキミならどう書く 2.0 - ROUND 1 - では、エラトステネスの篩という便利な篩があり、頭の中で素数を求める場合にも目の粗い篩を使っているわけですが、同じようにできないでしょうか?ここでは n が大きな値の探索に限ってみます。

まず、偶数の場合ですが、次が n/2 になるので、既に探索済みのものになり、調べる必要がなくなります。

と思ったら、3n+1 問題の中で、

これは自身なし

n=10000000000

g(9780657630) = 1133

4222.760u 4223.620s 2:23:12.06 98.3% 0+0k 0+0io 95pf+0w

というものがあり、9780657630 が偶数です。

しかし、9780657631 も 1133 ステップであることから、同じステップ数のものは複数存在する可能性があることを示唆しています。

次に、例えば、4n+1 の形の奇数の場合には、次が 3*(4n+1)+1 = 4*(3n+1) となるため、2 で 2 回割れます。つまり 3*(4n+1)+1 = 4*(3n+1) → 2*(3n+1) → 3n+1 となり、探索の必要がないでしょう。

if 文で組み込むと以下のようになります。

#! /usr/bin/gawk -f
BEGIN {
  if ( ARGV[1] == "" ) {
    num = 100;
  } else {
    num = ARGV[1];
  }
  for ( i = 1; i <= num; i++ ) {
    if ( i % 2 != 0 || ( n - 1 ) % 4 != 0 ) {
      count = 0;
      collatz( i );
      if ( count > max ) {
        max = count;
        max_num = i;
      }
    }
  }
  printf( "%d max=%d\n", max_num, max );
}
function collatz( num ) {
  if ( num == 1 ) {
    result = 1;
    count++;
  } else {
    if ( num % 2 == 0 ) {
      result = collatz( num / 2 );
    } else {
      result = collatz( 3 * num + 1 );
    }
    count++;
  }
  return result;
}

これを元に少し振るいにかけて 100000 の速度を見てみます。

Time ノーマル 篩あり
real 1m8.703s 1m9.386s
user 1m2.380s 1m3.044s
sys 0m0.364s 0m0.380s

あまり変わりませんね。むしろ遅いくらい。(w

実はこれ、if 文の "||" (OR) が影響しているようです。

if 文のところを

    if ( i % 2 != 0 ) {

にすると、

Time ノーマル 篩あり 偶数除去
real 1m8.703s 1m9.386s 0m44.331s
user 1m2.380s 1m3.044s 0m33.282s
sys 0m0.364s 0m0.380s 0m0.424s

と劇的 (?) に高速になりました。

こうしたあらかじめ人間が数学的な根拠から行う高速化としては、3n+1 問題:まだまだ高速化あたりに出ていますので、参考にしてみてください。

確かに Haskell や高速 CPU を用いると高速に解が求まるのかもしれませんが、その前に考えるということをしてみませんか?

最後に

私は数学者でも数学科の卒業でもなければ、大学時代に数学の単位を落としかけたことがある程度なので、この「Collatz の予想」がどのくらい凄いことなのかは分かりません。ただし、Collatz の予想を解くよりも、その背景の方が気になるので、いろいろ書いてみました。

LL だからこそ手軽に調べられる利点を生かして、いろいろと調べてみました。個人的には、先のキミならどう書く 2.0 - ROUND 1 - においても、かなり可読性が高い (C 言語ライクである) 部類に属していると思っており、特定の言語しか習得していなくても、ある程度読める言語であることが分かります。逆に、他の言語 (関数型言語は別ですが) を習得していれば、記述することも容易な言語ではないでしょうか。

さて、今回の最後はそうしたことを日本の awker 用に書いてくださった gawk のメンテナー Arnold Robbins の言葉で閉めたいと思います。お付き合いありがとうございます。

Awk is often called a "little" language. And in terms of number of features, that's true. But the features it does have, both singly and in combination, provide both power and elegance. You can do a lot of things with small, easy to read (and easy to write!) programs. If you invest the time to learn awk, you may find that you don't really *need* to invest a lot more time learning some of the other languages out there.

awk はしばしば「小さい」言語と呼ばれます。 特徴の点では、それは本当です。しかし、それが単独と組み合わせで持つ特徴は、パワーとエレガントの両方を提供してくれます。 小さく、読みやすい (そして、書くのも簡単です !) プログラムで多くのことができます。awk を学ぶ時間を費やすなら、ずっと多くの時間を他の言語のいくつかを学びながら費やす必要はないと本当にわかるでしょう。

…決まった。(コンプリート! (w