配列のランダムシャッフルと形式的証明

本記事は以下を記述する。興味のない人はここで回れ右。

• C言語での配列シャッフル実装
• その形式的証明

C言語での配列シャッフル実装

ご存知の通りC言語には、配列シャッフルに使える標準ライブラリー関数が存在しない。（C++には std::shuffle がある）

よってC言語向けの配列シャッフル関数を定義する。（よくある再発明に過ぎない）

宣言

C言語でのシャッフルAPIは、以下のインターフェースを持つだろう。

/* shuffle.h */
extern void* shuffle(void* destination,
                     void* source,
                     unsigned int size, unsigned int count,
                     unsigned int (*random)(void* context),
                     void* random_context);

source ポインターは、各要素が size の大きさを持ち count 要素を格納する配列へのポインター。
同サイズのメモリー領域をポイントする destination を渡す。シャッフル結果は destination に格納され、戻り値としても返される。 source と部分的に重複した場合の動作は未定義とする。
シャッフルは、乱数生成関数 random に基づいて実行する。 random 関数は random_context を伴って呼び出すたびに、要素数 count をはるかに上回る上限整数の中からひとつ値を選択して返すものとする。

定義

次は実装。至極単純で、例えば以下のようなものでよい¹。

# include <stdint.h>
# include <string.h>
# include "shuffle.h"

static void swap(void* array, unsigned size, unsigned i, unsigned j);

void* shuffle(void *dst, void* src, unsigned size, unsigned n,
              unsigned (*random)(void*), void* cookie) {
  memmove(dst, src, size * n);
  if (n > 1) {
    while (n != 1) {
      swap(dst, size, n - 1, random(cookie) % n);
      n -= 1;
    }
  }
  return dst;
}

void swap(void* array, unsigned size, unsigned i, unsigned j) {
  uint8_t* const p = array;
  uint8_t t[size];

  memcpy(t,            p + size * i, size);
  memcpy(p + size * i, p + size * j, size);
  memcpy(p + size * j, t,            size);
}

shuffle 本体のループは if からまとめて以下でよかろうというご意見はあろうとおもう。が、ここではあえて n != 1 をループのガードとした。この是非は参考文献に挙げた "The Science of Programming" をご一読の上、考えていただきたい。

for ( ; n > 1; --n) {
  swap(dst, size, n - 1, random(cookie) % n);
}

使用例

先の shuffle 実装は、たとえば以下のように使用する。

# include <stdio.h>
# include <stdlib.h>
# include "shuffle.h"

static unsigned my_random(void* _) {
  return random();
}

# define LENGTHOF(a) sizeof(a)/sizeof(a[0])

int main() {
  int arr[] = {1, 2, 3, 4, 5, 6, 7, 8};
  shuffle(arr, arr, sizeof(arr[0]), LENGTHOF(arr), my_random, NULL);

  { /* dump */
    int i;
    for (int i = 0; i < LENGTHOF(arr); ++i) {
      printf("%d ", arr[i]);
    }
    puts("");
  }

  return 0;
}

形式的証明

先の shuffle を、形式的に定義すると以下のようになる。
P は事前条件、 Q は事後条件で、これがプログラムの形式的な仕様である²。その他の細かな記法は "The Science of Programming" に従っている。

{P: b=B ∧ |r|=|b|-1 ∧ |b|>0}
n:= |b|
do n≠1 →
  p:= r[n-2] mod n
  b[p], b[n-1]:= b[n-1], b[p]
  n:= n-1
od
{Q: (n=1 ∧ b=B)
  ∨ (n=2 ∧ b=(B; 1:B[r[0] mod 2]; r[0] mod 2:B[1]))
  ∨ (n=3 ∧ b=((B; 2:B[r[1] mod 3]; r[1] mod 3: B[2]); 1:(B; 2:B[r[1] mod 3]; r[1] mod 3: B[2])[r[0] mod 2]; r[0] mod 2:(B; 2:B[r[1] mod 3]; r[1] mod 3: B[2])[1]))
  ...}

Q はプログラムの定義を書き下したものであるため、自明である（サボりすぎか……）。

補助関数で事後条件Qを整理

先の Q は配列の代入が恐ろしいことになってしまい読みづらい。補助関数をつかって整理していく。

swap関数を使って整理

swap(b, i, j) = (b; i:b[j]; j:b[i])³ を定義とする swap を使うことで先の Q は以下のように書き換えられる。

{Q: (n=1 ∧ b=B)
  ∨ (n=2 ∧ b=swap(B, 1, r[0] mod 2))
  ∨ (n=3 ∧ b=swap(swap(B, 2, r[1] mod 3), 1, r[0] mod 2))
  ...}

swap により見通しが良くなった。 n=4 の場合も綺麗に書ける。

(n=4 ∧ b=swap(swap(swap(B, 3, r[2] mod 4)
                           2, r[1] mod 3)
                           1, r[0] mod 2))

foldrで、さらに整理

swapを使った定義を眺めると、事後条件の構造がだんだん見えてくる。

{Q: (n=1 ∧ b=B)
  ∨ (n=2 ∧ b=swap(B, 1, r[0] mod 2))
  ∨ (n=3 ∧ b=swap(swap(B, 2, r[1] mod 3),
                          1, r[0] mod 2))
  ∨ (n=4 ∧ b=swap(swap(swap(B, 3, r[2] mod 4)
                               2, r[1] mod 3)
                               1, r[0] mod 2))
  ...}

つまり、 n が増えるたび swap の呼び出しが一つ増える。そして、その最も内側にある呼び出しが一番右の要素どこに移動するかを定め、以降、順に左に向かって入れ替え対象の要素をずらして行く。

ここで Haskell にもある foldr 関数を使うと、事後条件 Q は以下のような畳み込みで定義できる。

{Q: b=foldr (\(r,n) b → swap(b, |b|-1, r)) B (zip r [2..])}

Haskell に馴染みのない方に粗く説明しておく。
foldr は二引数関数、初期値、リストを受け取る関数。リストの右端の値から順に、初期値とに、渡した二引数関数を適用し、順次たたみ込んでいく。

結局、初出のコードは以下のように清書できる。

{P: b=B ∧ |r|=|b|-1 ∧ |b|>0}
do n≠1 →
  p:= r[n-2] mod n
  b[p], b[n-1]:= b[n-1], b[p]
  n:= n-1
od
{Q: b=foldr (\(r,n) b → swap(b, n-1, r)) B (zip r [2..])}

全ての順列が shuffle により、出る

乱数生成器である r により、 shuffle は入力列の順列からなる集合要素のうち、いずれかひとつを返せる性能があることを示したい。

入力列が (1, 2, 3, 4) であるとする。

\begin{eqnarray}
(1, 2, 3, 4) => \left\{
\begin{array}{ll}
(1, 2, 3, 4), & if r[2] = 3 + 4k \\
(1, 2, 4, 3), & if r[2] = 2 + 4k \\
(1, 4, 3, 2), & if r[2] = 1 + 4k \\
(4, 2, 3, 1), & if r[2] = 0 + 4k.
\end{array} \right.
\end{eqnarray}

さらに r[2] = 3 + 4k について展開すると、

\begin{eqnarray}
(1, 2, 3, 4) => \left\{
\begin{array}{ll}
(1, 2, 3, 4), & if r[1] = 2 + 3l \\
(1, 3, 2, 4), & if r[1] = 1 + 3l \\
(3, 2, 1, 4), & if r[1] = 0 + 3l.
\end{array} \right.
\end{eqnarray}

同様に r[1] = 2 + 3l を展開すると、

\begin{eqnarray}
(1, 2, 3, 4) => \left\{
\begin{array}{ll}
(1, 2, 3, 4), & if r[0] = 1 + 2m \\
(2, 1, 3, 4), & if r[0] = 0 + 2m.
\end{array} \right.
\end{eqnarray}

このように展開を考えることで以下の全順列 24 通りが得られることがわかる。

\begin{eqnarray}
(1, 2, 3, 4) => \left\{
\begin{array}{ll}
(1, 2, 3, 4), & if r = (1 + 2m, 2 + 3l, 3 + 4k) \\
(2, 1, 3, 4), & if r = (0 + 2m, 2 + 3l, 3 + 4k) \\
(1, 3, 2, 4), & if r = (1 + 2m, 1 + 3l, 3 + 4k) \\
(3, 1, 2, 4), & if r = (0 + 2m, 1 + 3l, 3 + 4k) \\
(3, 2, 1, 4), & if r = (1 + 2m, 0 + 3l, 3 + 4k) \\
(2, 3, 1, 4), & if r = (0 + 2m, 0 + 3l, 3 + 4k) \\
(1, 2, 4, 3), & if r = (1 + 2m, 2 + 3l, 2 + 4k) \\
(2, 1, 4, 3), & if r = (0 + 2m, 2 + 3l, 2 + 4k) \\
(1, 4, 2, 3), & if r = (1 + 2m, 1 + 3l, 2 + 4k) \\
(4, 1, 2, 3), & if r = (0 + 2m, 1 + 3l, 2 + 4k) \\
(4, 2, 1, 3), & if r = (1 + 2m, 0 + 3l, 2 + 4k) \\
(2, 4, 1, 3), & if r = (0 + 2m, 0 + 3l, 2 + 4k) \\
(1, 4, 3, 2), & if r = (1 + 2m, 2 + 3l, 1 + 4k) \\
(4, 1, 3, 2), & if r = (0 + 2m, 2 + 3l, 1 + 4k) \\
(1, 3, 4, 2), & if r = (1 + 2m, 1 + 3l, 1 + 4k) \\
(3, 1, 4, 2), & if r = (0 + 2m, 1 + 3l, 1 + 4k) \\
(3, 4, 1, 2), & if r = (1 + 2m, 0 + 3l, 1 + 4k) \\
(4, 3, 1, 2), & if r = (0 + 2m, 0 + 3l, 1 + 4k) \\
(4, 2, 3, 1), & if r = (1 + 2m, 2 + 3l, 0 + 4k) \\
(2, 4, 3, 1), & if r = (0 + 2m, 2 + 3l, 0 + 4k) \\
(4, 3, 2, 1), & if r = (1 + 2m, 1 + 3l, 0 + 4k) \\
(3, 4, 2, 1), & if r = (0 + 2m, 1 + 3l, 0 + 4k) \\
(3, 2, 4, 1), & if r = (1 + 2m, 0 + 3l, 0 + 4k) \\
(2, 3, 4, 1), & if r = (0 + 2m, 0 + 3l, 0 + 4k).
\end{array} \right.
\end{eqnarray}

参考書籍

• The Science of Programming（プログラミングの科学）

1. C89 では uint8_t vs[var]; のような宣言は不正と思っていたが gcc -Wall --std c89 -c shuffle.c -o shuffle.o が無警告でコンパイルできたので swap 関数を定義して利用した。 ↩

2. ある状態 P から状態 Q に至る差分を埋めるのが、ふだん、わたしたちが書いているコード、プログラムである。そして数式による表明 P と Q をつなぐ手続き型のコード S。この (P, S, Q) の三つ組をホーア・トリプルという。 ↩

3. tSoP での配列の代入は (b; i:v) として与えられる。すなわち元の配列は元のまま残り、ただしそのインデックス i を参照したときは値 v を返す。この観点で (b; i:b[j]; j:b[i]) は、 i と j を参照したときに b[j] と b[i] を返しほかはそのまま、つまり swap を意味している。 ↩