哲学対話イベントで使える Udon ギミックを作った話 #6 - キャリー付き乗算
こちらに記事を移行しました。
力学モデルの UdonSharp 実装において、マルチスレッドを実現するために Audio Filter を使用しましたが、そのために既存のライブラリが使用できず、乱数を自ら実装する必要がありました。
キャリー付き乗算を使用しました。これは厳密さが求められる分野には適しませんが、簡単な実装で素早く動くため、今回の用途には役立ちます。以下にコードをそのまま記載します。
ただし、ビット演算が含まれているため、初心者には分かりにくい部分もあるかと思います。できる限り詳しく説明しますので、ご覧ください。
キャリー付き乗算の疑似乱数 - MWC
private void Mwc(int seed)
{
_x = 0xffff0000 | (uint)(seed & 0xffff);
}まずは初期設定を行います。`seed` を使って `_x` に値をセットします。
これは、32 ビットある `seed` のうち下位 16 ビットだけを使用し、上位 16 ビットを 1 で埋めて `_x` に収める関数です。MWC は状態変数が 0 になるとずっと 0 しか出力しなくなるので、これを防ぐ意図があります。
`(seed & 0xffff)` を見てみましょう。
`0xffff` は二進数にすると `1111 1111 1111 1111` で、下位 16 ビットがすべて 1 の整数です。
`&` はビット演算の論理積（AND）という操作で、同じ位置にあるビットが互いに 1 のとき、1 を返します。
つまり、このコードは `seed` の下位 16 ビット部分を取り出しています。
行全体を見ると、`0xffff0000 | (uint)(seed & 0xffff)` とあります。
`0xffff0000` は二進数にすると `1111 1111 1111 1111 0000 0000 0000 0000` で、上位 16 ビットがすべて 1、下位 16 ビットがすべて 0 の、32 ビットの値です。
`|` は論理和（OR）という操作で、同じ位置にあるビットがどちらかでも 1 であれば、1 を返します。
つまり、これは取り出した 16 ビットの値と `0xffff0000` を組み合わせています。`seed` は下位 16 ビットだけを使い、上位 16 ビットを 1 で埋めています。
private int Next()
{
_x = ((_x & 0xffff) * 62904) + (_x >> 16);
return (int)(_x & 0xffff);
}この関数は、実際に新しい疑似乱数を計算します。
1 行目を見てみましょう。
`(_x & 0xffff) * 62904)`： `_x` の下位 16 ビットに 62904 を掛けています。
`(_x >> 16)`：`>>` は右シフトというビット演算で、右側にビットをずらす操作です。これを使って、`_x` の上位 16 ビット部分を取り出しています。
全体を見ると、`_x` の下位 16 ビットに 62904 を掛けて、`_x` の上位 16 ビットを足しています。
2 行目を見てみましょう。
先ほど計算した値の下位 16 ビットを取り出しています。
出力は 16 ビットの整数です。
一見、謎の計算をしています。どういうことか説明します。
実は、`_x & 0xffff` は割り算の余りを求めています。`_x mod (2 ^ 16)` です。
`_x >> 16` も実際のところ割り算で、 `_x / (2 ^ 16)` を求めています。
このように、ビット演算は掛け算や割り算を効率よく求めるために使われています。
なぜこうなるのでしょうか？実は、ここに基数の原理というべきものが隠されています。
右に 1 桁ずらす操作を 10 進法で考えてみましょう。
考える値は 10 です。
10 を右に 1 桁ずらしてみましょう。1 になります。
これで、`10 / (10 ^ 1) = 1` が求まりました。
左に 1 桁ずらす操作を 10 進法で考えてみましょう。
考える値は 1 です。
1 を左に 1 桁ずらしてみましょう。10 になります。
これで、`1 * (10 ^ 1) = 10` が求まりました。
下位 1 桁を取り出す操作を考えてみましょう。
考える値は 12 です。
下位一桁を取り出してみましょう。2 です。
これで、`12 mod (10 ^ 1) = 2` が求まりました。
つまり、やろうとしていることはこうです。
`carry = _x / (2 ^ 16)`
`x = _x mod (2 ^ 16)`
`(x * 62904 + carry) mod (2 ^ 16)`
簡単のために一部を変更して、式をよく見てみましょう。
x が 4、carry が 4、簡単のため 62904 を 6 とし、10 進法で考えます。
`x * 6 + carry` なので、`4 * 6 + 4` で、28 です。
`28 mod 10` は 8 です。得られた乱数は 8 です。
次の乱数を求めたいとき、carry は `28 / 10` で求めます。2 です。
`x * 6 + carry` なので、`8 * 6 + 2` で、50 です。
`50 mod 10` は 0 です。得られた乱数は 0 です。
次の乱数を求めたいとき、carry は `50 / 10` で、5 です。
...
というふうに計算できます。
実際のコードを見て、どんな計算をしているのか考えてみましょう。
`_x = 100000` とします。
`carry = _x / (2 ^ 16)`
`x = _x mod (2 ^ 16)`
`(x * 62904 + carry) mod (2 ^ 16)` で、乱数 58113 が得られます。
以上のようにして、キャリー付き乗算を用いた疑似乱数生成を実装しました。
次回予告
次回は、軽量化で試みて失敗したグリッド処理と近傍探索（kd-Tree探索）について扱います。