バブルソートの遅さを体感してみる

バブルソートには思い出がある。大学生のときにPC-9801落ちゲーを作った。上から落ちてくるトランプのカードを5×5のマトリックスに積み上げて、縦横斜めでポーカーの役を作る。やりこむと「縦列で数字を揃えようとしても最後で失敗する確率が高くてスコアの期待値は低い。むしろ横のフラッシュを多く目指せ。フラッシュ作りは一気にやらないと失敗のペナルティが高すぎる」というような経験則が蓄積していく楽しいゲームだった。いまでもやりたいぐらい。たぶんオリジナルは「琉球」という名前のアーケードゲームだった。

そのとき、ポーカーの役判定のロジックを書き下ろしたときにやったのが、バブルソートだということを、かなり後になってから知った。当時はソートという言葉すら知らなかった。

カードは最大で5枚だし、チェックすべき列も1度に最大4列だから、どんなにひどいことをやっても処理時間はかからないということは分かっていた。かといって、人間が並べ替えるように、大きいのとか小さいのを順にピックアップして並べ直す、という処理を書くのは面倒そうに感じた。それで直感的に「隣り合うやつらを比べて置き換えて、それを要素数の回数繰り返せば綺麗に並ぶんちゃうか?」と思ってfor文を書いたらうまくいった。ソートという言葉も知らないままだったけど、それで問題なく役判定はできたから、そこは深く考えずに終わらせた。

カードの種類をチェックしてフラッシュ判定関数を作ってフラグを立てたり、ジョーカーを含む場合に一番特点が高くなるように解釈するやり方を考えたりして、それがうまく行ったときに、ああ、プログラミングってパズルに近い感覚があって楽しいなと思った。でも、ソートだけで1冊本が書けるほどの議論があるとはまったく想像がつかなかった。

いまでもまだそういう感覚がある。そんなに大きなデータばかり並べ替えるわけでもなし、O(n^2)だからって別にどうでもええんちゃうか、と。「最悪」という感覚が分からない。こういう感覚は、きっとO(n^2)を体感したことがないからじゃないかと思った。

グラフは見たことがあるし、データ量の2乗に比例して時間がかかるという理屈も分かる。でも、それって「実際オレのどう人生に関係するのよ?」みたいな。386-20MHzのCPUですら、ぼくの役判定のルーチンは判定に0.1秒もかかってなかったと思う。

もしかすると、なんらかのアプリのユーザーとして、O(n^2)の腐った実装にイライラさせられている可能性はある。でも、ソースコードと対比して見ているわけではないので、O(n^2)のソートを最悪と感じる感覚は分からない。

それで、バブルソートの遅さを体験してみることにした。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int ShowList(int *list, int size) {
  int i;
  for (i = 0; i < size; i++) {
    printf("%d\n",list[i]);
  }
}

int Sort(int *list, int size) {
  int i, j, tmp;
  for (i = 0; i < size; i++) {
    for (j = 0; j < size - i; j++) {
      if (list[j + 1] < list[j]) {
	tmp = list[j + 1];
	list[j + 1] = list[j];
	list[j] = tmp;
      }
    }
  }
}

int main(int argc, char *argv[]) {
  int i, n, *list, size;

  if (argc < 1) exit(1);
  size = atoi(argv[1]);

  list = malloc(sizeof(int) * size);

  srand(time(NULL));
  for (i = 0; i < size; i++) {
    n = rand() % 1000;
    list[i] = n;
  }
  
  //  ShowList(list, size);
  Sort(list, size);
  //  ShowList(list, size);
  return 0;
}
require 'benchmark'

n = 50000

Benchmark.bm do |x|
  1000.step(n, 1000) do |i|
    x.report("sort #{i}") do 
      system("./bubble " + i.to_s)
    end
  end
end
$ ruby bench.rb
      user     system      total        real
sort 1000  0.000000   0.000000   0.010000 (  0.015362)
sort 2000  0.000000   0.000000   0.050000 (  0.057531)
sort 3000  0.000000   0.000000   0.110000 (  0.108889)
sort 4000  0.000000   0.000000   0.150000 (  0.156021)
sort 5000  0.000000   0.000000   0.220000 (  0.219979)
sort 6000  0.000000   0.000000   0.260000 (  0.262262)
sort 7000  0.000000   0.000000   0.390000 (  0.392773)
sort 8000  0.000000   0.000000   0.460000 (  0.462353)
sort 9000  0.000000   0.000000   0.590000 (  0.587144)
sort 10000  0.000000   0.000000   0.760000 (  0.764877)
sort 11000  0.000000   0.010000   0.890000 (  0.877772)
sort 12000  0.000000   0.000000   1.040000 (  1.045481)
sort 13000  0.000000   0.000000   1.230000 (  1.231655)
sort 14000  0.000000   0.000000   1.460000 (  1.458790)
sort 15000  0.000000   0.000000   1.630000 (  1.644733)
sort 16000  0.000000   0.000000   1.860000 (  1.857299)
sort 17000  0.010000   0.000000   2.140000 (  2.139211)
sort 18000  0.000000   0.000000   2.350000 (  2.351165)
sort 19000  0.000000   0.000000   2.670000 (  2.668415)
sort 20000  0.000000   0.000000   2.900000 (  2.899825)
sort 21000  0.000000   0.000000   3.230000 (  3.232370)
sort 22000  0.000000   0.000000   3.640000 (  3.637609)
sort 23000  0.000000   0.000000   4.070000 (  4.096720)
sort 24000  0.010000   0.000000   4.570000 (  4.624907)
sort 25000  0.000000   0.000000   5.080000 (  5.126958)
sort 26000  0.010000   0.010000   5.650000 (  5.669311)
sort 27000  0.000000   0.000000   6.460000 (  6.711002)
sort 28000  0.000000   0.000000   6.920000 (  6.979870)
sort 29000  0.000000   0.020000   7.310000 (  7.411398)
sort 30000  0.010000   0.000000   8.170000 (  8.197811)
sort 31000  0.000000   0.010000   8.600000 (  8.698379)
sort 32000  0.000000   0.010000   9.430000 (  9.533037)
sort 33000  0.000000   0.020000  10.380000 ( 10.395539)
sort 34000  0.000000   0.010000  10.830000 ( 10.993836)
sort 35000

最初のうちサクサク進む処理は目に見えて遅くなり、要素数が1万を越えた辺りからは、ゴキブリホイホイに足をとられたゴキブリのような感じで、グググググと処理がスローダウンするのが体感できる。

もっとも速いはずのintのソートですら、1、2万個で重くなる。と、考えてみて、やっと思い当たる節が出てきた。メールボックスの並べ替えだと、ふつうに数万のオーダーになる。

もしかして、20年前とか30年前にパソコンで相手にしていた量のデータでは、ソートなんて何でも良かったんじゃないだろうか。モデムで電話をかけてくれるアドレス帳アプリを作ったことがあるけど、そのときもバブルソートと知らずに書いた並べ替え関数で何の問題もなかった。

ところが、いまは数万個のオブジェクトなんてたいした量でもなくて、日常的に処理するデータ量なのじゃないか。

ぼくがゲームを作った当時にいちばん気にしていたのはダウンロード時間だった。当時、プログラマの友達に「ソートを自分で書くバカがいるか、ライブラリ使えよ」と言われたけど、パソコン通信時代にバイナリを小さくすることは結構重要だったので、ぼくは並べ替えごときでライブラリをリンクするのは納得がいかなかった。qsortなんて意味が分からない。文字出力も音楽もグラフィックも98のハードウェア依存のライブラリを自作していたので、ぼくのゲームは標準ライブラリを1つも使っていなかったし、そういうのは普通だったような気もする。実行バイナリは10KBしかなかった。これは、ダウンロードして遊んでくれる人に対して、30秒で済ませるか、5分待ってもらうかという違いだった。5分待たされてクソゲーなら怒るのがふつう。30秒なら「なんだ、小さい割に予想外にいい」と思ってるかもしれないと考えた。当時のゲームはダウンロード時間の長さでだいたいのデキが判断できた。

で、glibcのqsortを使ってみた。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int ShowList(int *list, int size) {
  int i;
  for (i = 0; i < size; i++) {
    printf("%d\n",list[i]);
  }
}

int cmp(const void *a, const void *b) {

  if (*((int*)a) > *((int*)b)) {
    return 1;
  } else {
    return -1;
  }
}

int main(int argc, char *argv[]) {
  int i, n, *list, size;

  if (argc < 1) exit(1);
  size = atoi(argv[1]);

  list = malloc(sizeof(int) * size);

  srand(time(NULL));
  for (i = 0; i < size; i++) {
    n = rand() % 1000;
    list[i] = n;
  }
  
  //  ShowList(list, size);
  qsort(list, size, sizeof(int), cmp);
  //  printf("\n");
  //  ShowList(list, size);
  return 0;
}

同じbench.rbを実行してみて、まあ驚いた。速いなんてもんじゃない。想定した範囲が一瞬で終わる。50万個のintでも0.02秒! バブルソートのほうは1万個でも1分以上かかるというのに。クイックソートは500万個でも2秒で終わる。ていうか1億個でも51秒で終わった。すごすぎる。

グラフで比較すると、こんな感じ。

クイックソートのほうで、もっと上の方を見てみると、こんなグラフ。100万個から1800万個辺りまで100万個ステップ。

このグラフを見て、そもそもnlognのグラフのイメージがぼくには分かってないことに気がづいた。クイックソートがO(nlogn)という知識はあったけど、実感としてもグラフとしての理解していなかった。青がn、黄色がlogn、オレンジがnlognのグラフは以下の通り。

まだまだもやもやするけど、とりあえず、バブルソートが最悪というニュアンスは分かった気がする。

ところで、よくよくCPUモニターを見てみると、コアが2つとも活用されてる気がするのだけど、1億個の並べ替えをやると、どうもあるタイミングで片方が休んだりしている。マルチパス処理? 何やってるんだ?

とりあえず、glibcのqsort.cをメモがわりに。

/* Copyright (C) 1991,1992,1996,1997,1999,2004 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Written by Douglas C. Schmidt (schmidt@ics.uci.edu).

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
   02111-1307 USA.  */

/* If you consider tuning this algorithm, you should consult first:
   Engineering a sort function; Jon Bentley and M. Douglas McIlroy;
   Software - Practice and Experience; Vol. 23 (11), 1249-1265, 1993.  */

#include <alloca.h>
#include <limits.h>
#include <stdlib.h>
#include <string.h>

/* Byte-wise swap two items of size SIZE. */
#define SWAP(a, b, size)						      \
  do									      \
    {									      \
      register size_t __size = (size);					      \
      register char *__a = (a), *__b = (b);				      \
      do								      \
	{								      \
	  char __tmp = *__a;						      \
	  *__a++ = *__b;						      \
	  *__b++ = __tmp;						      \
	} while (--__size > 0);						      \
    } while (0)

/* Discontinue quicksort algorithm when partition gets below this size.
   This particular magic number was chosen to work best on a Sun 4/260. */
#define MAX_THRESH 4

/* Stack node declarations used to store unfulfilled partition obligations. */
typedef struct
  {
    char *lo;
    char *hi;
  } stack_node;

/* The next 4 #defines implement a very fast in-line stack abstraction. */
/* The stack needs log (total_elements) entries (we could even subtract
   log(MAX_THRESH)).  Since total_elements has type size_t, we get as
   upper bound for log (total_elements):
   bits per byte (CHAR_BIT) * sizeof(size_t).  */
#define STACK_SIZE	(CHAR_BIT * sizeof(size_t))
#define PUSH(low, high)	((void) ((top->lo = (low)), (top->hi = (high)), ++top))
#define	POP(low, high)	((void) (--top, (low = top->lo), (high = top->hi)))
#define	STACK_NOT_EMPTY	(stack < top)


/* Order size using quicksort.  This implementation incorporates
   four optimizations discussed in Sedgewick:

   1. Non-recursive, using an explicit stack of pointer that store the
      next array partition to sort.  To save time, this maximum amount
      of space required to store an array of SIZE_MAX is allocated on the
      stack.  Assuming a 32-bit (64 bit) integer for size_t, this needs
      only 32 * sizeof(stack_node) == 256 bytes (for 64 bit: 1024 bytes).
      Pretty cheap, actually.

   2. Chose the pivot element using a median-of-three decision tree.
      This reduces the probability of selecting a bad pivot value and
      eliminates certain extraneous comparisons.

   3. Only quicksorts TOTAL_ELEMS / MAX_THRESH partitions, leaving
      insertion sort to order the MAX_THRESH items within each partition.
      This is a big win, since insertion sort is faster for small, mostly
      sorted array segments.

   4. The larger of the two sub-partitions is always pushed onto the
      stack first, with the algorithm then concentrating on the
      smaller partition.  This *guarantees* no more than log (total_elems)
      stack size is needed (actually O(1) in this case)!  */

void
_quicksort (void *const pbase, size_t total_elems, size_t size,
	    __compar_d_fn_t cmp, void *arg)
{
  register char *base_ptr = (char *) pbase;

  const size_t max_thresh = MAX_THRESH * size;

  if (total_elems == 0)
    /* Avoid lossage with unsigned arithmetic below.  */
    return;

  if (total_elems > MAX_THRESH)
    {
      char *lo = base_ptr;
      char *hi = &lo[size * (total_elems - 1)];
      stack_node stack[STACK_SIZE];
      stack_node *top = stack;

      PUSH (NULL, NULL);

      while (STACK_NOT_EMPTY)
        {
          char *left_ptr;
          char *right_ptr;

	  /* Select median value from among LO, MID, and HI. Rearrange
	     LO and HI so the three values are sorted. This lowers the
	     probability of picking a pathological pivot value and
	     skips a comparison for both the LEFT_PTR and RIGHT_PTR in
	     the while loops. */

	  char *mid = lo + size * ((hi - lo) / size >> 1);

	  if ((*cmp) ((void *) mid, (void *) lo, arg) < 0)
	    SWAP (mid, lo, size);
	  if ((*cmp) ((void *) hi, (void *) mid, arg) < 0)
	    SWAP (mid, hi, size);
	  else
	    goto jump_over;
	  if ((*cmp) ((void *) mid, (void *) lo, arg) < 0)
	    SWAP (mid, lo, size);
	jump_over:;

	  left_ptr  = lo + size;
	  right_ptr = hi - size;

	  /* Here's the famous ``collapse the walls'' section of quicksort.
	     Gotta like those tight inner loops!  They are the main reason
	     that this algorithm runs much faster than others. */
	  do
	    {
	      while ((*cmp) ((void *) left_ptr, (void *) mid, arg) < 0)
		left_ptr += size;

	      while ((*cmp) ((void *) mid, (void *) right_ptr, arg) < 0)
		right_ptr -= size;

	      if (left_ptr < right_ptr)
		{
		  SWAP (left_ptr, right_ptr, size);
		  if (mid == left_ptr)
		    mid = right_ptr;
		  else if (mid == right_ptr)
		    mid = left_ptr;
		  left_ptr += size;
		  right_ptr -= size;
		}
	      else if (left_ptr == right_ptr)
		{
		  left_ptr += size;
		  right_ptr -= size;
		  break;
		}
	    }
	  while (left_ptr <= right_ptr);

          /* Set up pointers for next iteration.  First determine whether
             left and right partitions are below the threshold size.  If so,
             ignore one or both.  Otherwise, push the larger partition's
             bounds on the stack and continue sorting the smaller one. */

          if ((size_t) (right_ptr - lo) <= max_thresh)
            {
              if ((size_t) (hi - left_ptr) <= max_thresh)
		/* Ignore both small partitions. */
                POP (lo, hi);
              else
		/* Ignore small left partition. */
                lo = left_ptr;
            }
          else if ((size_t) (hi - left_ptr) <= max_thresh)
	    /* Ignore small right partition. */
            hi = right_ptr;
          else if ((right_ptr - lo) > (hi - left_ptr))
            {
	      /* Push larger left partition indices. */
              PUSH (lo, right_ptr);
              lo = left_ptr;
            }
          else
            {
	      /* Push larger right partition indices. */
              PUSH (left_ptr, hi);
              hi = right_ptr;
            }
        }
    }

  /* Once the BASE_PTR array is partially sorted by quicksort the rest
     is completely sorted using insertion sort, since this is efficient
     for partitions below MAX_THRESH size. BASE_PTR points to the beginning
     of the array to sort, and END_PTR points at the very last element in
     the array (*not* one beyond it!). */

#define min(x, y) ((x) < (y) ? (x) : (y))

  {
    char *const end_ptr = &base_ptr[size * (total_elems - 1)];
    char *tmp_ptr = base_ptr;
    char *thresh = min(end_ptr, base_ptr + max_thresh);
    register char *run_ptr;

    /* Find smallest element in first threshold and place it at the
       array's beginning.  This is the smallest array element,
       and the operation speeds up insertion sort's inner loop. */

    for (run_ptr = tmp_ptr + size; run_ptr <= thresh; run_ptr += size)
      if ((*cmp) ((void *) run_ptr, (void *) tmp_ptr, arg) < 0)
        tmp_ptr = run_ptr;

    if (tmp_ptr != base_ptr)
      SWAP (tmp_ptr, base_ptr, size);

    /* Insertion sort, running from left-hand-side up to right-hand-side.  */

    run_ptr = base_ptr + size;
    while ((run_ptr += size) <= end_ptr)
      {
	tmp_ptr = run_ptr - size;
	while ((*cmp) ((void *) run_ptr, (void *) tmp_ptr, arg) < 0)
	  tmp_ptr -= size;

	tmp_ptr += size;
        if (tmp_ptr != run_ptr)
          {
            char *trav;

	    trav = run_ptr + size;
	    while (--trav >= run_ptr)
              {
                char c = *trav;
                char *hi, *lo;

                for (hi = lo = trav; (lo -= size) >= tmp_ptr; hi = lo)
                  *hi = *lo;
                *hi = c;
              }
          }
      }
  }
}