2011 年度「計算数学」 2011-12-02

N の値に応じて 関数 init() における配列の初期化を変えるのを忘れそうな人は、 「N が 100 未満だったら」「そうでなかったら」 と分けて初期化するように、関数 init() を書き直してください。

§37 Shell ソート

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

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

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

普通の日本語ではこれを「m-1 おきに取り出した」と呼ぶかもしれない。 このそれぞれの部分配列の中でソートを行ない、 それぞれの部分配列の中では小さい順にならんでいるようにすることを、 ここでは「m おきにソートする」と呼ぼう。

(37.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

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

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

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

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

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

(37.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 おきのソートが完了するのだから、 主張は正しい。)

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

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

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

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

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

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

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

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

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

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

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

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

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

(37.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 である部分に格子点があることがわかる。)

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

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

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

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

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

§38 課題2

(38.1) 課題1と同様のことを、Shell ソートについて行なえ。 なお、そのほうが見やすいと思うならば、 プログラムの出力に改行を加えても構わない。 課題1で書いた、 挿入ソートのソースファイルのコピーを元にする。 ソースファイル名は shell.c がよいだろう。

(38.2) ※ Shell ソートのプログラムを書く際には、 最初から配列を画面いっぱいにとること。 N が 10 ぐらいでは意味がない。 また、初めのうちは、 途中で m が変わるごとに 「40 おきにソートします」 などと出力させるとよいかもしれない。 (これらの間にソートされてゆく過程が出力されることになる。)

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

    for (m = 前述の大きな数; m ???; m = m/3) {
        挿入ソートの二重ループをここに入れ、一部手直しをする
    }

(38.4) N を大きくしての実験は、 N を 10000, 100000, 1000000 の三通りに変えて行なえ。

(38.5) 件名は「kadai2」(←全て半角文字、アルファベットは小文字、 途中にスペースをいれない)とせよ。 これ以外の注意点は課題1と同じである。

§39 発展問題

(39.1) ここは、時間などに余裕がある人のためのオプション項目である。

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

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

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

(39.5) 上で書いたプログラム(= m が 1743392200 から始まるもの)は、 int 型が 32 ビットの場合のものである。 それ以外の場合でも適切に動作する Shell ソートのプログラムを書け。

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

(39.7) (ビット数と、int 型で扱える数の範囲について簡単に説明する。 1 ビットには 0 か 1 かのどちらかが格納できる。 よって、 センターのパソコンのように int 型が 32 ビットであれば、 232 = 4294967296 通りの数が格納できる。 普通、このうち半分を非負の数に、残りの半分を負の数にあてるので、 センターの Linux では、 int 型であらわせる最小の数は -2147483648, 最大の数は 2147483647 となっている。)


岩瀬順一