続・ UnityEngine.Random の実装と性質～ rotation(Uniform) の謎

はじめに

前回の記事から調査が進んだので、 UnityEngine.Random の実装と性質についての追加情報をまとめました。

調査環境は Unity 2022.2.0b14 (Windows) です。

なお、調査は入出力の観測にて行っているため、実際の計算式や手順と異なる場合があります。数学的には等価だと思いますが、浮動小数点数の計算誤差や環境などによって結果がわずかに異なる場合があります。ご了承ください。

前回までのおさらい

前回の記事では、 InitState() ・ value ・ Range(int, int) ・ Range(float, float) について調査しました。

public class UnityRandom  
{  
    private uint x, y, z, w;  
  
    private uint Next()  
    {  
        uint t = x ^ x << 11;  
        x = y;  
        y = z;  
        z = w;  
        return w = w ^ w >> 19 ^ t ^ t >> 8;  
    }  
  
  
    public void InitState(int seed)  
    {  
        x = (uint)seed;  
        y = x * 1812433253 + 1;  
        z = y * 1812433253 + 1;  
        w = z * 1812433253 + 1;  
    }  
  
    public float value  
        => (Next() & 0x7fffffu) / 8388607.0f;  
  
    public int Range(int min, int max)   
        => min + (int)(Next() % (uint)(max - min));  
  
    public float Range(float min, float max)  
    {  
        float t = value;  
        return t * min + (1 - t) * max;  
    }  
}

`InitState()` についての補足

前回の記事で InitState() の調査を行った際、どうやって線形合同法のパラメータを求めたかを紹介します。

最初の内部状態を $x_0$ 、次に遷移したときの内部状態を $x_1$ 、さらにその次を $x_2$ とします。
このとき、遷移関数から以下の関係式が得られます。

$\begin{aligned} x_1 &= m x_0 + a \\ x_2 &= m x_1 + a \end{aligned}$

ここで、 $m$ と $a$ は線形合同法のパラメータであり、 $m$ は乗算の定数、 $a$ は加算の定数です。

この 2 式を組み合わせて、

$\begin{aligned} x_2 - x_1 &= m(x_1 - x_0) \\ m &= (x_2 - x_1) (x_1 - x_0)^{-1} \end{aligned}$

とすることで $m$ を得ることができます。
なお、 $(x_1 - x_0)^{-1}$ は通常の割り算ではなく $2^{n}$ を法とする逆数を掛けることで計算します。今回は 32 bit なので $n = 32$ です。
プログラム的には以下のコードで求めることができます。

public static uint MultiplicativeInverse(uint x)  
{  
    if ((x & 1) == 0)  
        throw new ArgumentException("x must be odd");  
  
    uint y = x;  
    y = y * (2 - y * x);  
    y = y * (2 - y * x);  
    y = y * (2 - y * x);  
    y = y * (2 - y * x);  
  
    return y;  
}

$m$ がわかればあとは簡単で、最初の式から、

$a = x_1 - m x_0$

とすれば $a$ が求まります。

というわけで、連続する 3 つの値が観測できれば、線形合同法のパラメータそのものを復元できるといえます。

ところで、ビルドしたアプリやエディタの初期状態についてはよく分からない、としましたが、現在のスクリプトリファレンスを確認すると、

It is statically initialized with a high-entropy seed from the operating system, and stored in native memory where it will survive domain reloads.

とあるため、 /dev/random 的なもの (何らかの CSPRNG) で初期化されているようです。

以前の記事を書いたときには何も記載がなかったような気がします… こういった情報が増えていくのはありがたいことです。

今回の進捗

今回はそれ以外の応用的なメソッドについて調べてみました。

`insideUnitCircle`

単位円内のランダムな点を取る insideUnitCircle は、以下のように計算できます。

public Vector2 insideUnitCircle  
{  
    get  
    {  
        var angle = Range(0f, MathF.PI * 2);  
        var radius = MathF.Sqrt(Range(0f, 1f));  
        return new Vector2(  
            radius * MathF.Cos(angle),   
            radius * MathF.Sin(angle));  
    }  
}

念のためですが、 MathF は System 名前空間に生えている標準搭載のほうで、 Unity の Mathf とは異なります。
計算自体はほとんど変わりませんが、非 Unity 環境との互換性を考えて今回は MathF にしています。

`insideUnitSphere`

単位球内のランダムな点を取る insideUnitSphere は、以下のように計算できます。

public Vector3 insideUnitSphere  
{  
    get  
    {  
        var p = MathF.Asin(Range(-1f, 1f));  
        var t = Range(0, MathF.PI * 2);  
        var r = MathF.Cbrt(value);  
  
        return new Vector3(  
            r * MathF.Cos(p) * MathF.Cos(t),   
            r * MathF.Cos(p) * MathF.Sin(t),   
            r * MathF.Sin(p));  
    }  
}

p が xz 平面上の回転、 t が xy 平面上の回転、 r が中心からの距離です。
MathF.Cbrt は立方根を求める関数です。

`onUnitSphere`

単位球上のランダムな点を取る onUnitSphere は、以下のように計算できます。

public Vector3 onUnitSphere  
{  
    get  
    {  
        var p = MathF.Asin(Range(-1f, 1f));  
        var t = Range(0f, MathF.PI * 2f);  
  
        return new Vector3(  
            MathF.Cos(p) * MathF.Cos(t),   
            MathF.Cos(p) * MathF.Sin(t),   
            MathF.Sin(p));  
    }  
}

`rotation`

ランダムな Quaternion を取得する rotation は、以下のように計算できます。

public Quaternion rotation  
{  
    get  
    {  
        float x = Range(-1f, 1f);  
        float y = Range(-1f, 1f);  
        float z = Range(-1f, 1f);  
        float w = Range(-1f, 1f);  
  
        float reciprocalLength = 1 / MathF.Sqrt(x * x + y * y + z * z + w * w);  
        reciprocalLength *= w < 0 ? -1 : 1;  
  
        return new Quaternion(  
            x * reciprocalLength,  
            y * reciprocalLength,  
            z * reciprocalLength,  
            w * reciprocalLength);  
    }  
}

これの計算方法はスクリプトリファレンスに書いてあったので調査が簡単でした。
なお、 w は常に正の値をとるようです。

`rotationUniform`

よりランダムな (?) Quaternion を取得する rotationUniform は、以下のように計算できます。

public Quaternion rotationUniform  
{  
    get  
    {  
        float a = Range(0f, 1f);  
        float b = Range(0f, MathF.PI * 2);  
        float c = Range(0f, MathF.PI * 2);  
  
        float x = MathF.Sin(b) * MathF.Sqrt(1 - a);  
        float y = MathF.Cos(b) * MathF.Sqrt(1 - a);  
        float z = MathF.Sin(c) * MathF.Sqrt(a);  
        float w = MathF.Cos(c) * MathF.Sqrt(a));  
  
        float signOfW = w < 0 ? -1f : 1f;  
  
        return new Quaternion(  
            x * signOfW,  
            y * signOfW,  
            z * signOfW,  
            w * signOfW);  
    }  
}

これも計算手法がスクリプトリファレンスに書いてあったのですが、

Employs Hopf fibration to return a random Quaternion within a uniformly distributed selection space.

Hopf fibration とは……? こういう数学はさっぱりで Wikipedia を見ても全く分からなかったので、観察と直感で調べました。理論について詳しい方がいらっしゃればぜひご教授ください…

`rotation` と `rotationUniform` の違いについて

スクリプトリファレンス曰く、

(rotationUniform) Gives higher quality results compared to the more naive approach employed by rotation, though at a 40% performance cost.

とのことで、 rotationUniform は品質が良く、 rotation は高速な代わりに相対的に品質が良くない、とされています。
これは実際のところ、どの程度違うのでしょうか?

実際に点を描画してみる

rotation はぱっと見た感じ、一辺の長さが 2 の 4 次元超立方体内のベクトルを生成 → 正規化(4 次元超球上のベクトルに変換) → クォータニオンに、としている気がします。
となると、超立方体の角に相当する部分が濃くなりそうな雰囲気があります (これは 3 次元立方体 → 球の変換を考えれば理解しやすそうです) が、これが実際どのように影響するかというと……?

まずは、普通に見るとどうなるか見てみましょう。

これは、 Vector3.forward のベクトルをランダムに回転させた位置に点を打つ作業を行ったものです。左の赤いのが rotation, 右の緑が rotationUniform です。
まったく違いが判りません……。

それでは、 Quaternion.w < 0.01 のものだけ抽出するとどうなるでしょうか。

w < 0.01 のやつだけ抽出したのですが、こうすると rotation 側は X っぽい模様が見える…ような… (この処理が正しいか微妙) pic.twitter.com/XJabzhbLe6
— Andante (@andanteyk) 2023年4月1日

斜めから見た画像はこちらで、

上から見た画像はこちらです。

rotation のほうには赤い X のような偏りが見えています。これがもしかすると 4 次元超立方体の辺や角の部分かもしれません。
対して rotationUniform のほうにも濃い点ができていますが、これはおそらく正しいはずです。
理由を説明しましょう。こちらの文献曰く、 $\lambda$ 軸周りに $\theta$ だけ回転させるクォータニオンはこのような値を持っています。

$\begin{split} x &= \lambda_x \mathrm{sin}(\theta / 2) \\ y &= \lambda_y \mathrm{sin}(\theta / 2) \\ z &= \lambda_z \mathrm{sin}(\theta / 2) \\ w &= \mathrm{cos}(\theta / 2) \\ \end{split}$

ということは $w \approx 0$ の場合には $\theta \approx 180 °$ であるはずなので、 Vector3.forward を 180° 回転させた位置に (上から見た画像の画面下方向に) 集中するはず…です。

ヒストグラムにしてみる

ところで、直感的な映像による確認もよいですが、簡単な統計データとしても見てみましょう。

クォータニオンは一辺 -1 <= x <= 1 の範囲の 4 次元超立方体内に収まるので、一辺を 16 等分して 65536 個の領域に分割し、どの領域に何個クォータニオンが入ったかのヒストグラムをつくりました。
簡単のため Unity の AnimationCurve を悪用してグラフを作りました。

まず、以下が rotation のほうのヒストグラムです。 Unity が途中で線の描画を諦めてしまっている点にご留意ください…

次に、以下が rotationUniform のほうのヒストグラムです。

これを見ると、特に両端のエリア (Quaternion.x が -1 や +1 に近いエリア) において rotationUniform のほうが均等に分布していそうな印象を受けます。
また、 rotation では部分的に突出したデータが見られます。グラフの縦軸にも注目してみてください。

標準偏差を取ると、 rotation では 62.2241 ・ rotationUniform では 54.0865 となりました。 (なお、平均値は 16 です。)
数値的にも rotationUniform のほうがばらつきが少ないことが示唆されます。

次に、 eulerAngles を使って 3 次元にした場合のヒストグラムを見てみましょう。
同様に一辺 16 等分して 4096 個の領域に分割してヒストグラムを作成しました。

まず、 rotation のほうのヒストグラムです。

次に、 rotationUniform のほうのヒストグラムです。

こちらも、やはり rotationUniform のほうが突出した部分が少なく、円形の理想的な形に近い分布になっていそうです。

標準偏差を取ると、 rotation では 21.9124 ・ rotationUniform では 19.6578 となりました。 (同じく、平均値は 16 です。)

パフォーマンス

最後に、パフォーマンスの差について調べてみましょう。
私の環境 (12th Gen Intel(R) Core(TM) i7-12700F) で $2^{20}$ 個の Quaternion を生成するのにかかった時間は、 rotation は 61.167 ms ・ rotationUniform は 83.261 ms でした。
rotation のほうが 36 % ほど速く、スクリプトリファレンスに載っていた数値と概ね合致しました。

まとめ

以上を私見を交えてまとめると、

正味の違いはほぼ分からないが、詳しく見ると rotation で偏りがわずかに見られる
基本的には rotationUniform を使うべき
- 何万オーダーで生成するような、パフォーマンス第一の場合は rotation を検討する

といった感じでしょうか。

`ColorHSV()`

ランダムな色を取得する ColorHSV() は、以下のように計算できます。

public Color ColorHSV(  
    float hueMin, float hueMax,  
    float saturationMin, float saturationMax,   
    float valueMin, float valueMax,   
    float alphaMin, float alphaMax)  
{  
    // note: (max, min) の順。実際は Lerp(value, min, max) のように実装されています  
    float h = Range(hueMax, hueMin);  
    float s = Range(saturationMax, saturationMin);  
    float v = Range(valueMax, valueMin);  
    float a = Range(alphaMax, alphaMin);  
  
    Color HsvToRgb(float h, float s, float v, float a)  
    {  
        float hStar = h * 6;  
  
        float c = v * s;  
        float x = c * (1 - Math.Abs(hStar % 2f - 1f));  
  
        return hStar switch  
        {  
            >= 0 and  
            < 1 => new Color((v - c) + c, (v - c) + x, (v - c) + 0, a),  
            < 2 => new Color((v - c) + x, (v - c) + c, (v - c) + 0, a),  
            < 3 => new Color((v - c) + 0, (v - c) + c, (v - c) + x, a),  
            < 4 => new Color((v - c) + 0, (v - c) + x, (v - c) + c, a),  
            < 5 => new Color((v - c) + x, (v - c) + 0, (v - c) + c, a),  
            < 6 => new Color((v - c) + c, (v - c) + 0, (v - c) + x, a),  
            _ => throw new ArgumentException()  
        };  
    }  
  
    return HsvToRgb(h, s, v, a);  
}

ほぼ Wikipedia に載っているものと同じです。

ちなみに、これは UnityCsReference で実際のコードを確認することもできます。

crazy ...?https://t.co/RAvrolDoxw
— Andante (@andanteyk) 2023年4月1日

補足・雑記など

その他、気付いたことや雑記などです。

浮動小数点数の計算について

今回実験をしていて気が付いたのですが、実行環境が Unity 上 (mono) か Visual Studio 上 (dotnet) かで出力される浮動小数点数の値が微妙に異なる場合があるようです。
Unity 上では絶対に合わなかったテストケースが dotnet に移植したプロジェクトでは合った場合があったためです。

全く根拠のない推測ですが、 mono のほうはいわゆる -ffast-math のようなアグレッシブな最適化をしているのではないかと思います。

`UnityEngine.Random` に値を仕込む方法

調査にあたっては、まず UnityEngine.Random から任意の値を出力できるようにする必要があります。いわゆる「乱数調整」みたいなものです。これをどうやって行ったのかを記載しておきます。

Random の内部構造は Xorshift であることは本記事の冒頭や前回の記事でお伝えした通りです。
Xorshift は内部状態をそのまま吐き出すタイプなので、かなり楽に調整ができます。

具体的には、

4 個の内部状態に欲しい値をそれぞれセットする
乱数を 4 回巻き戻す

という手順を通すだけです。

Xorshift を巻き戻す処理は具体的には以下で行えます。どうやって求めたかはこちらの記事を参照してください。

public void Prev()  
{  
    w ^= z ^ z >> 19;       // t ^ t >> 8  
    w ^= w >> 8; w ^= w >> 16;  
    w ^= w << 11; w ^= w << 22;     // w = x  
    uint tx = w;  
    w = z;  
    z = y;  
    y = x;  
    x = tx;  
}

そして、こういう感じで内部状態をセット + 巻き戻し処理を行います。
HackState() は欲しい値を 4 つセットすると、それに対応する内部状態を返してくれます。

public static uint[] HackState(Span<uint> values)  
{  
    var emu = new UnityRandom() { 
        x = values[0], 
        y = values[1], 
        z = values[2], 
        w = values[3] };  
    emu.Prev();  
    emu.Prev();  
    emu.Prev();  
    emu.Prev();  
    return new uint[] { emu.x, emu.y, emu.z, emu.w };  
}

あとは、リフレクションで内部状態を書き込むメソッドを作成して Random.state = state; で書き込めば、 UnityEngine.Random の出力値を 4 つぶんまではコントロールできます。

public static UnityEngine.Random.State SetState(uint[] values)  
{  
    // ボックス化が必要なので object で受ける  
    object state = new UnityEngine.Random.State();  
  
    var flags = BindingFlags.NonPublic | BindingFlags.Instance | BindingFlags.SetField;  
  
    state.GetType().GetField("s0", flags).SetValue(state, (int)values[0]);  
    state.GetType().GetField("s1", flags).SetValue(state, (int)values[1]);  
    state.GetType().GetField("s2", flags).SetValue(state, (int)values[2]);  
    state.GetType().GetField("s3", flags).SetValue(state, (int)values[3]);  
  
    return (UnityEngine.Random.State)state;  
}

おわりに

というわけで、現状存在するすべてのメソッド・プロパティの挙動を調べることができました。満足です。
(浮動小数点数の下位桁が微妙にずれる問題はありますが、めちゃくちゃ面倒なので見なかったことに……)

付録として、上記のソースをまとめたものを添付しておきます。

UnityRandom.cs

はじめに

前回までのおさらい

InitState() についての補足

今回の進捗

insideUnitCircle

insideUnitSphere

onUnitSphere

rotation

rotationUniform

rotation と rotationUniform の違いについて