不思議な計算ループ

あまりに駄文なので、もう、どうしよう。でも書いちゃったしなー。ということで。

なにを計算するプログラムでしょう？　というお題があって、そりゃあもう気になるじゃあないですか、あなた。

/* precondition: n >= 0 */
int mystery(int n, int *p)
{
    int q, r, x, y;

    r = n + 1;
    q = 0;
    while (r >= 4) {
        x = r / 4;
        y = r % 4;
        q = q + x;
        r = x + y;
    }
    *p = r - 1;
    return q;
}

不思議な計算ループ - Think Stitch - PRINCIPIA

C での複数値返し

初見で感じたのが、 C は複数の値をまとめて返せなくて不便だなあってこと。ただ C でも構造体を使えば値はまとめて返せるんで、それで見通しがよくなるんじゃないかなー…　とおもっていた。ので書いてみた。結果、注目すべき点から目がそらされてしまう感じの残念なコードになってしまった。

typedef struct pair { unsigned q, r; } pair_t;

pair_t mystery(unsigned n) {
  pair_t v = { 0, n + 1 };
  while (v.r >= 4) {
    unsigned x, y;
    x = v.r / 4;
    y = v.r % 4;
    v.q = v.q + x;
    v.r  = x + y;
  }
  return v;
}

本来注目すべき点

ところでもとのコードだけれど、本来注目すべきは while の条件部分と x と y の計算部分。

unsigned mystery(unsigned n, unsigned *r) {
  unsigned p, x, y;
  p = 0;
  *r = n + 1;
  while (*r & ~3) {
    x = *r >> 2;
    y = *r & 3;
    p = p + x;
    *r = x + y;
  }
  *r = *r - 1;
  return p;
}

すごい。下位二ビットと残り上位ビットのシフト演算の組み合わせで、ああいう計算ができるプログラムになっているだなんて、ねえ。

それ、 Haskell で

ま、見通しが残念な書き換えになっちゃったとか、本来注目すべきは〜だとかは横において。みんな大好き Haskell で書き直してみよう！　こうなる。

mystery n = aux (0, n+1)  where
  aux (q, r) | r >= 4    = aux (q+x, x+y)  where (x, y) = divMod r 4
  aux (q, r) | otherwise = (q, r-1)

while ループの部分は、補助関数 aux を再帰的に呼び出す形でおきかえられる。
末尾再帰になっているのが美しい。美しいんだけれど、関数型ならではの別の方法、そう、みんな大好きリストがある。

無限リストを使った定義

mystery' n = (q, r-1)  where
  (q, r) = head (dropWhile (\(_, r) -> r >= 4) (aux (0, n+1)))
  aux (q, r) = (q, r) : aux (q+x, x+y) where (x, y) = divMod r 4

む。リストが見えない。見えづらいけれど aux (q, r) がリストだ。
こちらの定義での関数 aux はタプル (q, r) を取る一引数の関数 *1 で、 (q, r) を先頭にして再帰的にリスト aux (q+x, x+y) をつないだ無限リストを返す。ここで x と y の値は divMod r 4、つまり x が div r 4、 y が mod r 4 になっている。
最初に戻ると、こちらの定義 mystery' n はタプル (q, r-1) という値だと言っている。
ここで q と r の値は二行めの定義で決めていて、無限リスト aux (0, n+1) の要素のいくつかを dropWhile で落としたあと、その先頭要素を head で取った値（タプルだ）、となっている。
dropWhile の条件はラムダ式 \(_, r) -> r >= 4。タプルの第二要素の値 r が 4 以上のもの、と言っている。
蛇足ながら、このラムダ式はポイントフリースタイルで (>=4).snd と書き直せる。 snd はタプルの第二要素を返す関数。 (>=4) はセクションで Integral な型の値を一つとり、それが 4 以上なら真となる関数。それを . で関数合成する。 f.g という関数合成なあれ。ここで g が snd。タプルを取ってその第二要素を返す。その結果を、関数 (>=4) で評価する、という形になっている。わかりにくい！　たのしい！

結局のところ aux (0, n+1) は、(0, n+1) からはじまって、先の while ループで計算する途中の q, r の列を無限につないだリストになっている。
無限リストの計算は終わらないんだけれど、それを dropWhile で r が 4 以上のタプルをばんばん落としていって、 4 以下になった最初のタプルを返す、という仕組み。

iterate

mystery' での補助関数 aux は、関数 (:) を使って値とリストを自分でのりづけすることで無限リストを得ていたけれど、 Haskell には値 v を与えると、その値 v に関数 f を適用した f v をつないで、さらにその f v に f を適用した f (f v) をつないで、以下繰り返して無限リストを返す関数 iterate がある。
再帰呼び出しや無限リストの生成を自分で書き下すのはダサいそうだ。
iterate を使うと、 mystery' の aux (q, r) はこう書き換えられる。

aux (q, r) = iterate f (q, r)  where
  f (q, r) = let (x, y) = divMod r 4 in (q+x, x+y)

ループの回数

mystery' で定義した aux は、結果がでるまでの途中計算をぜんぶ記憶していると見ることができる。実際、次のようにして答えと、答えを求めるまでの中間計算をタプルにして出すこともできる。

mystery'' n = let
  vs = iterate f (0,n+1)
  f = \(q,r) -> let (x,y) = divMod r 4 in (q+x,x+y)
  (as,bs) = span ((>=4).snd) vs
  (q,r) = head bs
  in ((q,r-1),as++[(q,r)])

65535 を mystery に食わせたら 8 回の計算で結果がでるんだって。へえ。

*1:型は (Integral t) => (t, t) -> (t, t)、つまりタプルを取ってタプルを返す関数