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

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

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

#C言語での配列シャッフル実装
ご存知の通りC言語には、配列シャッフルに使える標準ライブラリー関数が存在しない。（C++には [std::shuffle](https://msdn.microsoft.com/library/dn688973.aspx) がある）

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

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

```c
/* 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` をはるかに上回る上限整数の中からひとつ値を選択して返すものとする。

##定義
次は実装。至極単純で、例えば以下のようなものでよい[^std]。

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

```c
#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" をご一読の上、考えていただきたい。

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

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

```c
#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` は事後条件で、これがプログラムの形式的な仕様である[^program]。その他の細かな記法は "The Science of Programming" に従っている。

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

```text
{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])`[^subst_array] を定義とする `swap` を使うことで先の `Q` は以下のように書き換えられる。

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

```text
{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` の場合も綺麗に書ける。

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

###`foldr`で、さらに整理
`swap`を使った定義を眺めると、事後条件の構造がだんだん見えてくる。

```text
{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](https://www.haskell.org/hoogle/?hoogle=foldr) 関数を使うと、事後条件 `Q` は以下のような畳み込みで定義できる。

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

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

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

```text
{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)` であるとする。

```math
\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` について展開すると、

```math
\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` を展開すると、

```math
\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 通りが得られることがわかる。

```math
\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（プログラミングの科学）