2023 年度「計算数学b」 2023-06-16

§8.1 Shell ソート

(8.1.1) ここから、素朴でないソートについて学び始める。 自明でないアルゴリズムなので、まず、 「どんな手順でやるか」がわかりにくいことが多い。 手順を理解したうえで、 「なぜこれでソートできるのか?」 「ちょっと見たところ手間がかかるのに、なぜこのほうが速いのか?」 を、数学を用いて考察する(のだが、 時間がない。興味のある者だけ、プリントを見て自分で考えよ)。

(8.1.2) 初めに、Shell ソートを扱う。 Shell は人名なので大文字で始まる。

(8.1.3) (正の)自然数 m を一つ選んで固定する。 配列 a[1], ..., a[N] を、 その中から m おきに取り出した、 次の m 本の部分配列に分割する。

このそれぞれの部分配列の中でソートをおこない、 それぞれの部分配列の中では小さい順にならんでいるようにすることを、 ここでは「m おきにソートする」と呼ぼう。

(8.1.4) 具体例で説明する。 次を「3 おきにソートする」としよう。

    18 08 13 93 78 32 81 65 02 40

まずは a[1] から始まる部分配列を考える。 “いまは見ない”部分は「..」とする。

    18 .. .. 93 .. .. 81 .. .. 40

これをこの中で小さい順に並べ替えるとこうなる。 (どうやって並べ替えるかは後述。)

    18 .. .. 40 .. .. 81 .. .. 93

こんどは a[2] から始まる部分配列を考え、それを小さい順に並べ替える。

    .. 08 .. .. 78 .. .. 65 .. ..
                  ↓
    .. 08 .. .. 65 .. .. 78 .. ..

最後に a[3] から始まる部分配列を考え、それを小さい順に並べ替える。

    .. .. 13 .. .. 32 .. .. 02 ..
                  ↓
    .. .. 02 .. .. 13 .. .. 32 ..

よって、最終的にはこうなる。

    18 08 02 40 65 13 81 78 32 93

(8.1.5) たとえば配列の大きさ N が 1000 のとき、

のが Shell ソートである。 各ステップは挿入ソートでおこなう。 この数列 364, 121, 40, 13, 4, 1 は漸化式 m1 = 1, mn+1 = 3 * mn + 1 で定義される数列を逆順に見たものである。 この数列に特別な意味があるわけではない。 1 で終わる別の単調減少数列を選んでもよい。 ただし、あとでわかるように、 隣り合った項が互いに素であるほうが望ましいようである。

(8.1.6) 最後におこなう 1 おきのソートとは普通のソートのことだから、 上のアルゴリズムで正しくソートされることは間違いない。 よって、これから考えるべき問題は、 「なぜこれが速いのか?」である。

(8.1.7) 自然数 m, M が 0 < m < M を満たすとき、 次の意外な事実が成り立つ。

M おきにソートしてから m おきにソートした結果は、 M おきにソートされている。

(8.1.8) (証明の概略: 「m おきの転倒数」を # {(i, j) | i < j かつ j-im の倍数、かつ a[i] > a[j]} で定義する。 いま、配列は M おきにソートされているとする。 m おきのソートが完了していなければ、 m おきの転倒数は 0 でなく、 転倒している組 (a[i], a[i+m]) が存在する。 このような i を一つ選んで固定し、 組 (a[i+k*M], a[i+m+k*M]) (k は整数)のうち、転倒している組すべてに互換をほどこす。 この操作のあと、m おきの転倒数が減少していること、 および M おきにソートされたままであることを観察せよ。 この操作を m おきの転倒数が 0 になるまでくり返してゆけば m おきのソートが完了するのだから、 主張は正しい。)

(8.1.9) 仮に M が 5, m が 3 だったとしよう。 5 おきにソートし、それから 3 おきにソートしたあと、

であることはもちろんだが、a[i-5]a[i] にも注意すると

が言える。 a[i-10]a[i] であることからは

が出てくる。 ここに出てくる項をよく見ると、 「c が 8 以上ならば a[i-c]a[i]」が成り立つことがわかる。

(8.1.10) このあとに 1 おきのソートを挿入ソートでおこなえばソートが完了するが、 上の考察から、任意の i に対し #{j | j < i かつ a[j] > a[i]} < 8 であることがわかるから、 どの要素を挿入する際も高々 8 回の比較、高々 7 回の交換で十分である。 計算量をおおざっぱに調べてみよう。

Ο(N2) であることに変わりはないが、 挿入ソートの計算量の N2/4 と比べると、 N2 の係数が小さくなった分だけ、 速いはずである。

(8.1.11) ※ 実験してみると、 3 おきにソートする際の計算量はこれよりはるかに小さい。 それは、 その前に 5 おきにソートしたことで、小さいものが左に、 大きいものが右に寄っているためであろう。

(8.1.12) 上の例では「5 おき」「3 おき」にソートしたのがよい“下準備”となった。 そこから 8 という数が出てきて、 最後の「1 おき」のソートの計算量が 8N となったのだった。 一般の場合を考えよう。 自然数 m, M が 0 < m < M を満たすとき、 M おきにソートしてから m おきにソートしたあとでは 0 以上の整数 x, y に対し a[i-(m*x+M*y)]a[i] が成り立つ。

(8.1.13) 次の定理はすでに習ったであろう。

a と b を互いに素な自然数とし、 d を整数とするとき、 方程式 ax+by = d は整数解 (x, y) をもつ。

(8.1.14) ことばを変えて言えば、 a と b とが互いに素なら直線 ax+by = d の上には格子点 (= x 座標も y 座標も整数である点)が存在する、 ということである。 (x, y) が解であるとき、任意の整数 k に対し (x+kb, y-ka) も解である。 これ以外の解がないことはすぐわかる。 よって、格子点は直線 ax+by = d の上に x 座標でみて b おき、y 座標でみて a おきに並んでいる。

(8.1.15) この定理から次が得られる。

a と b を互いに素な自然数とするとき、 整数 d が (a-1)(b-1) 以上ならば方程式 ax+by = d は 0 以上の整数の解 (x, y) をもつ。

(8.1.16) (証明の概略: 上の定理から、直線 ax+by = d の上に格子点があることはわかる。 この直線が二直線 x = -1, y = -1 と交わる点の、x 座標の差を考えよう。 d = (a-1)(b-1) のとき、 直線 y = -1 とは x = b-1+1/a で交わる。 これと x = -1 との差は b+1/a となり、b よりも大きい。 直線 ax+by = d は、d が大きくなるほど上方へ移動するので、 d がより大きければ、この差はさらに大きくなる。 よって、d が (a-1)(b-1) 以上であるならば、 直線 ax+by = d のうち x > -1 かつ y > -1 である部分に格子点があることがわかる。)

(8.1.17) ※ 方程式 ax+by = (a-1)(b-1)-1 は 0 以上の整数の解 (x, y) を持たない。 これは (b-1, -1) と (-1, a-1) とが隣り合った解であることからわかる。

(8.1.18) ※ さっきの例で出てきた「8」は (M-1)(m-1) = 4 * 2 にあたる。

(8.1.19) 上の例では「5 おき」「3 おき」「1 おき」にソートしたが、 この授業の実習では (8.1.5) で述べた単調増加数列を逆にたどりつつ、 もっと何度も「m おきのソート」をおこなう。 mN なら 「m おきのソート」とは何もしないことだから、 m が大きすぎても構わない。 よって、m の値は、(8.1.5) で述べた単調増加数列の項のうち、 いま使っているコンパイラの int 型で表せる最大の数から始めればよい。 その値は、 簡単なプログラムで求められるように 1743392200 である。 このようにした場合の Shell ソートの理論的な計算量は私には全くわからない。

(8.1.20) 【重要】 上では m 本の部分配列をそれぞれ m おきにソートすると説明したが、 実際には

とするほうが楽だろう。これなら三重ループで済むからである。

§8.2 課題2

(8.2.1) 課題1と同様のことを、Shell ソートについておこなえ。 課題1で書いた、 挿入ソートのソースファイルのコピーを元にせよ。 ソースファイル名は shell.c がよいだろう。 ファイル名には小文字が便利なので、頭文字を小文字にした。

(8.2.2) ※ Shell ソートのプログラムを書く際には、 最初から配列を画面いっぱいにとること。 N は 78 がよいだろう。 挿入ソートの完成後にコメントにした関数 print() の呼び出しは復活させる。 また、初めのうちは、 途中で m が変わるごとに 「40 おきにソートします.」 などと出力させるとよい。 (注意:「ソート」と出力させるには、 ソースファイルには「ソ\ート」と書くのだった。 これはコンパイラの都合である。(1.11.12) 参照。)

(8.2.3) 挿入ソートのプログラムは for が二重になっていたが、 それをさらに for ループで囲み、 そのループで m を減らしてゆく。その際、(8.1.5) で述べた単調増加数列を逆にたどる際の漸化式は m = (m-1)/3 でもよいが、 C言語では int 型の正の数どうしの割り算では余りは切り捨てられるから、 m = m/3 と書けばよい。 ソートをおこなう部分の全体は次のようになる。

  for (m = 1743392200; m ???; m = m/3) {
    printf("%d おきにソ\ートします.\n", m);
    /* 挿入ソートの二重ループがここにくる。一部手直しが必要 */
  }

(8.2.4) N を大きくしての実験は、 N を 10000, 100000, 1000000 の三通りに変えて五回ずつ。

(8.2.5) 件名は「??? kadai2」(←全て半角文字、 ??? は自分の履修者番号、その後ろに半角スペース一つ、 kadai は小文字、 kadai2 の間にはスペースを入れない)とせよ。 これ以外の注意点は課題1と同じである。

(8.2.6) このプリントの説明に従っていないプログラムについては、採点できない場合があります。

§8.3 時間などに余裕がある人のための発展問題

(8.3.1) 5 おき、3 おき、1 おきにソートする Shell ソートのプログラムを書き、(8.1.7), (8.1.9) の主張が正しいことを確認せよ。 さらに、それぞれの段階における、 そこまでの比較回数・交換回数を出力するようにせよ。 それを見て、(8.1.10) の計算が正しいか確かめよ。

(8.3.2) Shell ソートの m の取り方を別のものに変えたらどうなるか。 効率を調べてみよ。 N/2 から始めて 2 で割ってゆくものは、 N が 2 のベキ乗のとき効率がよくないと言われている。

(8.3.3) 前に (8.1.19) に書いた定数 1743392200 を求めるプログラムを書け。

(8.3.4) 挿入ソートを、N を 100000 にして動かし、 かかった時間を時計で測定せよ。それを元に、 1000000 のときの所要時間を推測せよ。 それと今回の Shell ソートの場合とを比べてみよ。

(8.3.5) 関数 print() の改変版 print2(int i, int j) を書け。 この関数は、 print() と同様に配列の要素を出力するが、 a[i]a[j] だけはその前あるいは後ろに何か文字を出力する、 というものである。 関数 swap() から呼び出すときは print2() のほうを使えば、 いま交換した要素がどれだかわかり、見やすくなる。 ((8.5.6) も参照。)

§8.4 ソートプログラムの効率の限界

(8.4.1) 互いに異なる三つの数 a, b, c を比較して小さい順に並べるとしよう。 下の図は、最初に a < b かどうかを調べ、「はい」なら上へ、 「いいえ」なら下へ進む、と読む。一番右が結論である。

                          +--------------- 「a < b < c」
                          |
            +- (b < c ?) -+
            |             |             +- 「a < c < b」
            |             +- (a < c ?) -+
            |                           +- 「c < a < b」
-(a < b ?) -+
            |                           +- 「b < a < c」
            |             +- (a < c ?) -+
            |             |             +- 「b < c < a」
            +- (b < c ?) -+
                          |
                          +--------------- 「c < b < a」

高々3回の比較で結論に至るが、これを2回に減らすことは不可能である。 それは、2回の分岐では4通りの結論しかありえないのに対し、 結論は 3! = 6 通りあるからである。

(8.4.2) 同様の議論により、 高々 k 回の比較で N 個の元をソートするアルゴリズムがあれば、 2kN! が成り立つ。 自然対数をとると k log(2) ≥ log(N!) = log(1) + log(2) + log(3) + ... + log(N-1) + log(N) となる。 この右辺は、関数 y = log(x) の積分と比較することにより、 N log(N) - N + 1 よりも大きく、 (N+1) log(N+1) - N よりも小さいことがわかる。 N が十分大きければほぼ N log(N) とみなせる。 よって、 N が十分大きいとき、 “ほぼ k ≥ N log(N) / log(2) = N log2(N)”となる。 情報科学では対数の底として 2 を選ぶことが多いらしい。 その約束を採用するなら、“ほぼ k ≥ N log(N)”と書ける。 オーダーで言えば、Ο(N log(N)) が限界となる。 これより真に小さいオーダーの計算量で済むアルゴリズムは存在しない。 (現行のコンピュータを使う限り。)

(8.4.3) 次回から、Ο(N log(N)) のソートアルゴリズムを学ぶ。

修正:関数 myrand() を次のように変えてください

この授業の課題は学術メディア創成センター演習室のパソコンでこなしてもらうつもりで準備してきました。 RAND_MAX の値が大きいコンパイラの場合は修正が必要です。 関数 myrand() を次のように修正してください。 こうすれば、RAND_MAX が 32677 でも、もっと大きくても、うまくゆきます。

int myrand(void) {      /* 1 から 1073741824 までに値をとる乱数を返す */
  int r;

  r = rand() % 32678;
  return rand() % 32678 * 32768 + r + 1;
}

重要なお断り

以下のプログラム例は、 ここの演習室のコマンドプロンプトでは期待通りには動きません。 画面制御エスケープシーケンスに対応していないからです。 最近の Windows のコマンドプロンプトや Mac のターミナルは対応しているようです。

§8.5 画面制御エスケープシーケンス

(8.5.1) いままで、printf() で画面に文字を出力してきたが、 文字の色は白のみで、上から下へと順に書いてゆくことしかできなかった。 改行は別として。

(8.5.2) ここのコマンドプロンプトは、 printf("\033[2J"); のようにして \033[ で始まる文字列を画面に出力することで、画面を“制御”できる。 これらの文字列を「画面制御エスケープシーケンス」と(私は)呼んでいる。 これは C 言語の規則ではなく、コマンドプロンプトの仕様である。

(8.5.3) (一時期の Windows では、この画面制御エスケープシーケンスが使えなかった。)

(8.5.4) (\n のように \ で始まるものを C 言語のエスケープシーケンスと呼ぶことがある。 これは C 言語の規則なので、混同しないように。)

(8.5.5) 出力する文字に色などをつけるもの。 例えば \033[33m で、その次から出力する文字が黄色になる。 元に戻すには \033[m である。 次のプログラムは、この 33 をいろいろな数に変えて、その効果を見るものである。

#include <stdio.h>

int main() {
  int n;

  for (n = 0; n < 256; n++) {
    printf("\033[%dm %03d \033[m", n, n);
    if (n % 16 == 15) {
      printf("\n");
    }
  }
  printf("\n");
  printf("\033[05;45;92m 05;45;92 の重ねがけ \033[m"); 
}

最後の行で示したように、\033[05;45;92m として、 05, 45, 92 の三つの効果を重ねがけすることもできる。

(8.5.6) ソートプログラムを書くにあたってのこれらの使いみちの一つを示す。 関数 print() の改変版 print2(int i, int j) を作る。 この関数は、print() と同様に配列の要素を出力するが、 a[i]a[j] だけは別の色などで出力する、というものである。 関数 swap() から呼び出すときは print2() のほうを使えば、 いま交換した要素だけ、別の色などで出力させることができ、見やすくなる。

(8.5.7) 関数 print() を次で置き換えると、 配列の値を棒グラフのように表示する。 N の値はコマンドプロンプトの行数以下とすること。

void print(void) {
  int i, k;

  for (i = 1; i <= N; i++) {
    printf("\033[%d;1H", i);      /* カーソル移動。\033[2;4H なら 2 行 4 列へ */
    printf("%02d:", a[i]);
    for (k = 0; k < a[i]; k++) {  /* a[i] の値と同じ数だけ # を出力 */
      printf("#");
    }
    printf("\033[0K");            /* 行末まで削除 */
  }
}

(8.5.8) ほかには、画面全体を消去する printf("\033[2J"); に使い道がある。 次の節のプログラムに使用例を載せておいた。 これ以外は、ネットを検索して探してみられたい。

§8.6 エコーなし文字入力

(8.6.1) C 言語に元々ついているキーボードからの入力関数は、 scanf() も含めて、すべて、入力された文字を画面に書く。 これを「エコー」と呼ぶ。 ここのコマンドプロンプトでは「エコーなし」の入力も可能だが、 どのコンピュータでも動くという保証はない。

(8.6.2) エコーなし入力関数 getch() を使ったサンプルプログラム。 (「'」はシフトしながら 7 のキーを打つ。)

#include <stdio.h>
#include <conio.h>    /* getch() を使うのに必要 */

int main() {
  int i, number;
  char c;             /* 入力された文字は char 型に代入する */

  printf("\033[2J");              /* 全画面消去 */
  printf("\033[1;1H");            /* カーソルを (1, 1) へ */
  printf("u で増加、d で減少、e で終了\n");

  number = 0;
  c = 0;
  for (     ; c != 'e';   ) {     /* 「e という文字」は 'e' と書く */
    printf("\033[2;1H");          /* カーソルを (2, 1) へ */
    printf("%10d\n", number);

    c = getch();                  /* エコーなしの入力 */
    if (c == 'u') {
      number++;
    } else if (c == 'd') {
      number--;
    } else {
      printf("\a");               /* 警告(音が鳴ったりします)*/
    }
  }

  printf("あなたは %d を入力しました.\n", number);
}

岩瀬順一