2025-05-21

モンゴメリ乗算について～ a x b mod m の高速化

はじめに

$a \times b \bmod m$ を高速かつ正確にやりたいとき、ありますよね。
とくに「 $a \times b$ の結果がワードサイズ (ulong とか) を超えてオーバーフローする場合」にどうするのか、という問題があります。
本稿では、それを高速かつ正確に行うことができる「モンゴメリ乗算」と呼ばれる手法について紹介したいと思います。

モンゴメリ乗算を三行で説明すると、

変な係数がかかるかわりに、オーバーフローを気にせず $a \times b \bmod m$ が計算できる
ついでに自然と乗法逆数も使えるようになる
頑張れば任意の $m > 0$ に対して適用可能

という感じです。

私が調べている最中、数式ばっかで肝心のプログラムコードがない例を良く見かけたので、この記事ではなるべくコードも載せていこうと思います。
逆に、証明については省略している場合が多いです。興味のある方は各キーワードでググってください。

なお、以下特に断りがない限り、実装言語は C# 、数値型は ulong (C でいう uint64_t 、符号なし 64 ビット整数) であるものとします。とはいえ実装固有のものは少ないのでほかの言語や型にも応用可能です。

モンゴメリ乗算とは

モンゴメリ乗算 (Montgomery multiplication) は、先に述べたように $a \times b \bmod m$ を高速かつオーバーフローさせることなく実行する手法です。

モンゴメリリダクション

まずは基本となる操作である「モンゴメリリダクション」について説明します。

モンゴメリリダクション $MR(x)$ は、以下のような操作です。

$MR(x) = x R^{-1} \bmod m$

ここで、 $R$ は $R > m$ かつ $\gcd(R, m) = 1$ の整数、 $R^{-1}$ は乗法逆数 (モジュラ逆数、 $R R^{-1} \equiv 1 \bmod m$ を満たす値) です。乗法逆数の求め方については後述します。

さて、これを安直に実装しようとすると $x R^{-1} \bmod m$ で詰まると思います。それができれば苦労はしない。
ですが、うまく $R$ を設定することによって効率よく計算可能になります。

具体的には、 $R$ を 2 冪、特にワードサイズ (ulong) と同じ $R = 2^{64}$ にすることです。

$\begin{align*} result &= \lfloor x / R \rfloor - \lfloor (x m^{-1} \bmod R) m / R \rfloor \\ &= \lfloor x / 2^{64} \rfloor - \lfloor (x m^{-1} \bmod 2^{64}) m / 2^{64} \rfloor \end{align*}$

$MR(x) = \begin{cases} result + m & \text{if } result \lt 0 \\ result & \text{otherwise}\end{cases}$

$\begin{align*} x&: \text{unsigned 128 bit integer} \\ R&: 2^{64} \\ m&: \text{unsigned 64 bit odd integer} \\ m^{-1}&: \text{unsigned 64 bit odd integer, } m \times m^{-1} \equiv 1 \pmod R \end{align*}$

どのあたりが効率よくなっているのかというと、まずは $\lfloor x / R \rfloor$ 、より具体的には $\lfloor x / 2^{64} \rfloor$ です。これは $x$ が Uint128 としたときの上位 64 ビットにあたります。 (どうして 128 ビットなのかは後述します。)
後半も同様で、まず $x m^{-1} \bmod R$ は $x m^{-1}$ の下位 64 ビットなので、掛け算したらそのままオーバーフローさせれば OK です。 $\lfloor (...) / R \rfloor$ は同様に上位 64 bit を取り出せばよいです。
式をよく見ると、それ以外の除算・剰余算は出てこないことが分かります。つまり、 実質除算・剰余算をすることなく $\bmod m$ をとることができます 。

最後に $result$ の値域ですが、 $-m \lt result \lt m$ になります。負値になった時は $m$ を足して正の値に戻します。
なお、 $m$ は $R = 2^{64}$ と互いに素である以上、奇数である必要があります。

以上をプログラムに落とし込むとこんな感じになります。

public static ulong MontgomeryReduction  
    (ulong xlo, ulong xhi, ulong m, ulong mInv)  
{  
    ulong result = xhi - Math.BigMul(xlo * mInv, m, out _);  
    return result > xhi ? result + m : result;  
}

数式のいかつさに比べるとずいぶんシンプルに見えるのではないでしょうか。

Math.BigMul は 64bit x 64bit = 128bit の掛け算を行うメソッドで、第三引数に下位 64bit ・戻り値に上位 64bit をセットします。今回は下位ビットは不要なので捨てています。

Unity 環境 (.NET Standard 2.1) にはまだ実装されていません。つらい。
利用したい場合は Unity.Burst.Intrinsics.Common.umul128 (上位と下位が逆なので注意が必要です) もしくは自前の polyfill を使うことになります。

public static ulong BigMul(ulong a, ulong b, out ulong lo)  
{  
    ulong alo = a & 0xffffffff, ahi = a >> 32;  
    ulong blo = b & 0xffffffff, bhi = b >> 32;  
  
    ulong lolo = alo * blo;  
    ulong lohi = alo * bhi;  
    ulong hilo = ahi * blo;  
    ulong hihi = ahi * bhi;  
  
    lo = lolo + ((lohi + hilo) << 32);  
    ulong carry = ((lolo >> 32) + (lohi & 0xffffffff) + (hilo & 0xffffffff)) >> 32;  
    ulong hi = hihi + (lohi >> 32) + (hilo >> 32) + carry;  
  
    return hi;  
}

話を戻して、 result > xhi はオーバーフローを検出する構文です。
ひとつ前の行で result = xhi - (略) と計算していますが、この引き算が負のオーバーフローを起こした場合に result > xhi が成り立ちます。覚えておくと何かの役に立ちます、はい。

ところで、 Wikipedia と式が微妙に違うことに気づかれた方は慧眼です。
Wikipedia では以下のような式になっています：

$result = \lfloor x / R \rfloor + \lfloor (x m^{-1} \bmod R) m / R \rfloor \\ m \times m^{-1} \equiv -1 \pmod{2^{64}}$

1 つめの式の引き算が足し算になっていて、 $m$ の逆数の定義が $-1$ と合同になっています。
これでも間違いではないのですが、 $result$ の値域が $0 \le result \lt 2m$ になりオーバーフロー検出が難しくなってしまうため、本稿では引き算の式を採用します。結果自体は変わりません。
加えて、 $m$ の逆数に関しても $1$ と合同のほうが高速に計算できる場合があるので、このようにしています。

乗法逆元

前の項で出てきた $m^{-1}$ や $R^{-1}$ の求め方について説明します。

まず、乗法逆元 (モジュラ逆元、 multiplicative inverse とも) というのは、乗法における逆元であり、もとの数に掛けると単位元 $1$ になるような数です。それはそう。
例えば、 $x$ の加法 (足し算) の逆元は皆さんおなじみ $-x$ で、 $x + (-x) = 0$ が成り立ちます。同様に、乗法 (掛け算) の逆元は「実数の範囲では」 $1/x$ (または $x^{-1}$ )になり、 $x \times (1/x) = 1$ となります。
整数の範囲では $1/x$ をとることはできませんが、 $\bmod m$ の世界では似たような操作を考えることができます。

$x$ と $m$ が互いに素な、つまり最大公約数が 1 ( $\gcd(x, m) = 1$ ) なとき、必ず乗法逆元 $x^{-1}$ が存在し、 $x \times x^{-1} \equiv 1 \pmod m$ を満たします。

具体例を考えると、 $3 \pmod 7$ の乗法逆数は $5$ です。 $3 \times 5 = 15 \equiv 1 \bmod 7$ となることから確かめられます。

なお、 $x$ と $m$ が互いに素でない場合は乗法逆元が存在しません。例えば、 $2 \pmod 6$ は何を掛けても $0, 2, 4$ のどれかにしかならないので乗法逆元が存在しません。

さて、加法逆元は - 単項演算子ですぐに求められますが、乗法逆元はどうやって求められるのでしょうか？
いくつか方法があるので紹介します。

拡張ユークリッドの互除法

一般に適用可能な手法が「拡張ユークリッドの互除法」 (Extended Euclidean algorithm) です。
拡張ユークリッドの互除法を使うと、 $x$ と $m$ の最大公約数 $\gcd(x, m)$ と、存在していれば乗法逆数 $x^{-1} \bmod m$ を求められます。一度に求められてお得です。

// TODO: a, mod must be smaller than 2^63  
public static (ulong gcd, ulong inverse)  
    Egcd(ulong a, ulong mod)  
{  
    ulong x0 = 0, x1 = 1, y0 = mod, y1 = a;  
  
    while (y1 != 0)  
    {  
        ulong q = y0 / y1;  
        (x0, x1) = (x1, x0 - q * x1);  
        (y0, y1) = (y1, y0 - q * y1);  
    }  
  
    ulong inverse = x0 >= mod ? x0 + mod : x0;  
    return (y0, inverse);  
}

前述したように、 gcd が 1 のときに限り inverse の値が有効になります。
場合によっては例外を投げたり、 0 を返すような実装にしてもいいかもしれません。

なお、 TODO に書いた通り q * x1 の計算でオーバーフローすると正しい値が求められない場合があります。
このあたりは後で説明する別の手法を使うと解決できます。

Jeffrey Hurchalla 法 ( $m$ が 2 冪の場合)

$m$ が 2 冪の場合 (ここでは $m = 2^{64}$ のとき) 、 Jeffrey Hurchalla 氏が提案した手法を使うとより高速に求められます。 *1

public static ulong MultiplicativeInverse(ulong a)  
{  
    Debug.Assert((a & 1) != 0);  
  
    ulong x0 = (3 * a) ^ 2;  
    ulong y = 1 - a * x0;  
  
    ulong x1 = x0 * (1 + y);  
    y *= y;  
  
    ulong x2 = x1 * (1 + y);  
    y *= y;  
  
    ulong x3 = x2 * (1 + y);  
    y *= y;  
  
    ulong x4 = x3 * (1 + y);  
  
    return x4;  
}

Debug.Assert にある通り、 $m$ が 2 冪の場合、偶数は乗法逆元を持ちません ( $\gcd(x, m) \ge 2$ になるため) 。
あとはまるでビット黒魔術ですね。さっぱりわかりません。

最初の (3 * a) ^ 2 では下位 5 bit の乗法逆元が求められます。
あとは x1 で 10 bit, x2 で 20 bit, x3 で 40 bit, x4 で 80 bit ぶん求められるようです。

これは従来のニュートン法を用いた手法と比べると、命令レベルの並列実行性が高いために高速になるそうです。
比較については下記サイトが詳しいです。

Integer multiplicative inverse via Newton's method

以上から、 $m^{-1} \bmod R$ は Jeffrey Hurchalla 法で求めるのがよさそうです。
$R^{-1}$ は手動で求める必要はありませんが、もし欲しければ $MR(1)$ で求めることができます。

モンゴメリ乗算

さて、話をモンゴメリ乗算に戻しましょう。

モンゴメリ乗算で $a \times b \bmod m$ をどうやるのかというと、以下のようにします。

$\begin{align*} A :&= MR(a R^{2}) \\ B :&= MR(b R^{2}) \\ C :&= MR(AB) \\ result :&= MR(C) = ab \bmod m \end{align*}$

どういうことか説明しましょう。
まず $A$ では $a \times R^{2}$ を計算し、それをモンゴメリリダクションしています。
$MR(x) = x R^{-1} \bmod m$ だったことを思い出すと、 $MR(a R^{2}) = a R^{2} R^{-1} = aR \pmod m$ となることが分かります。普通の数に $R^{2}$ を掛けてモンゴメリリダクションし、 $aR$ の形式にすることを一般に「モンゴメリ表現への変換」と言うそうです。
$R^{2}$ は $m$ に依存する定数なので、事前に計算しておきます。 (求め方は後述します。)
すると、実際の処理としては $a \times R^{2}$ を 64bit x 64bit = 128 bit の掛け算で求めて、これをモンゴメリリダクションすることになります。だから、モンゴメリリダクションの引数が ulong 2 個分 (xlo, xhi) だったのですね。

乗算してからモンゴメリリダクションするところまでをまとめてメソッドにしておきましょう。

/// a x b mod m  
public static ulong MontgomeryMultiplication(ulong a, ulong b, ulong m, ulong mInv)  
{  
    ulong xhi = Math.BigMul(a, b, out ulong xlo);  
    return MontgomeryReduction(xlo, xhi, m, mInv);  
}

次に、 $B$ も同様にモンゴメリ表現への変換を行います。

続いて、 $C$ では $A \times B$ を計算してモンゴメリリダクションします。
これでどうなるのかというと、

$\begin{align*} MR(AB) &= MR(aR \times bR) \\ &= MR(ab R^2) \\ &= ab R^2 R^{-1} \\ &= ab R \pmod m \end{align*}$

となります。 $R$ が 1 つ残っているので、モンゴメリ表現を維持しています。

最後に、 $C$ に対してモンゴメリリダクションを行うことで $R$ を打ち消して $ab$ だけにします。
これで $a \times b \bmod m$ の計算が完了しました。

なお、乗算を一回だけ行う場合は、以下のようにするとモンゴメリ乗算 2 回で求めることができて効率的です。

$result := MR(MR(ab)R^{2})$

モンゴメリ加算・減算

モンゴメリ乗算について説明しましたが、実は加算・減算も行うことができます。
といっても、単に足したり引いたりするだけです。 (オーバーフローへの対処は必要です。)

$A + B = aR + bR = (a+b)R \\ A - B = aR - bR = (a-b)R$

public static ulong MontgomeryAddition(ulong a, ulong b, ulong m)  
{  
    ulong add = a + b;  
    return add < a || add >= m ? add - m : add;  
}  
  
public static ulong MontgomerySubtraction(ulong a, ulong b, ulong m)  
{  
    ulong sub = a - b;  
    return a < b ? sub + m : sub;  
}

加算のほうの補足説明として、 add >= m でオーバーフローしなかったけど m 以上になった場合、 add < a でオーバーフローした ( $2^{64}$ 以上になった) 場合を処理しています。

$R^{2}$ の計算

さて、モンゴメリ表現への変換のために $R^{2} = 2^{128} \bmod m$ を使いましたが、これはどのように計算すればよいでしょうか？
そもそも多倍長乗算をしたくないからモンゴメリ乗算をしようとしているのに、ここで BigInteger に頼っては本末転倒です。どうにか求める方法を考えてみましょう。

とりあえず手始めに、 $R = 2^{64} \bmod m$ を求めてみましょう。
普通に考えると ulong で扱える領域を超えているので難しそうですが、実は $m$ を引くことで計算できるようになります。

$R \equiv R - m \pmod m$

で、 $2^{64}$ は ulong 的にはオーバーフローして $0$ なので、プログラムとしては以下のようにして求められます。

ulong rmod = (0ul - mod) % mod;

次に、 $R^{2}$ を考えましょう。
しかし、ここですぐにモンゴメリ乗算は使えません。なぜなら、

$MR(R \times R) = R \times R \times R^{-1} = R$

に戻ってしまうためです。

したがって、私は以下のようにいくつか加算してから乗算するようにしています。

$\begin{align*} r :&= R + R = 2R \\ r :&= r + r = 4R \\ r :&= MR(r \times r) = 2^{4}R \\ r :&= MR(r \times r) = 2^{8}R \\ r :&= MR(r \times r) = 2^{16}R \\ r :&= MR(r \times r) = 2^{32}R \\ r :&= MR(r \times r) = 2^{64}R = 2^{128} \pmod m\\ \end{align*}$

この定数計算はよくやるので、 $m^{-1} \bmod R$ の計算とまとめてメソッドにしておくと便利です。

/// returns (m^-1 mod R, R mod m, R^2 mod m)  
public static (ulong modinv, ulong rmod, ulong r2mod)   
    MontgomeryConstant(ulong mod)  
{  
    ulong modinv = MultiplicativeInverse(mod);  
    ulong rmod = (0ul - mod) % mod;  
  
    ulong r2mod = rmod;  
    r2mod = MontgomeryAddition(r2mod, r2mod, mod);  
    r2mod = MontgomeryAddition(r2mod, r2mod, mod);  
    r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);  
    r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);  
    r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);  
    r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);  
    r2mod = MontgomeryMultiplication(r2mod, r2mod, mod, modinv);  
  
    return (modinv, rmod, r2mod);  
}

冪剰余 $a^{b} \bmod m$

モンゴメリ乗算は単発の計算にも便利ですが、真価を発揮するのは $m$ が固定値で何度も計算するときです。
一番わかりやすいのは冪剰余 $a^{b} \bmod m$ の計算でしょう。

繰り返し二乗法を使います。
一応解説しておくと、 $b$ を二進数に分解して (例えば $13 = 2^{0} + 2^{2} + 2^{3}$ )、 $a$ を自乗して $a^{2^{0}}, a^{2^{1}}, a^{2^{2}}, a^{2^{3}}, ...$ としながら、 $b$ のビットが立っているところだけ掛けて答えにする ( $a^{0} \times a^{2^{2}} \times a^{2^{3}} = a^{13}$ ) 、といったアルゴリズムです。

コードにするとこういう感じですね。

public static ulong ModPow(ulong value, ulong exponent, ulong mod)  
{  
    var (modinv, rmod, r2mod) = MontgomeryConstant(mod);  
  
    ulong power = MontgomeryMultiplication(value, r2mod, mod, modinv);  
    ulong result = 1;  
  
    while (exponent > 0)  
    {  
        if ((exponent & 1) != 0)  
        {  
            result = MontgomeryMultiplication(result, power, mod, modinv);  
        }  
  
        power = MontgomeryMultiplication(power, power, mod, modinv);  
        exponent >>= 1;  
    }  
  
    return result;     
}

ポイントとしては、 power はモンゴメリ表現 $aR$ で持っていますが、 result はそうではない $r$ ( $rR$ ではない) というところです。
power のほうは $MR(aR \times aR) = a^{2} R$ となってモンゴメリ表現を維持し、 result のほうは $MR(r \times aR) = ra$ となって通常の表現を維持します。
result をモンゴメリ表現で持っておくと、あとでリダクションしないといけなくなるため、それを節約しています。

個人的には、「モンゴメリ表現」と考えてしまうとどうも固定観念を持ってしまいがちなので、「 $R^{-1}$ がついてくるかわりにオーバーフローしない乗算」と考えると便利だと思います。

発展問題

$m$ が偶数のときのモンゴメリ乗算

$R$ と $m$ は互いに素でなければならない以上、偶数の $m$ を使ってモンゴメリ乗算することはできません。
でも $a \times b \bmod m$ の計算はしたいです。これだけのために多倍長乗算はしたくないです。どうすればいいでしょうか？

実は方法があります。

中国の剰余定理

中国の剰余定理 (Chinese Remainder Theorem; CRT) を応用すると、異なる法を持つ数値群からある法の数値を復元することができます。

分かりづらいので式にすると、

$x \equiv a \pmod{m_1} \\ x \equiv b \pmod{m_2}$

があって $a, b$ が既知であり、かつ $m_1$ と $m_2$ が互いに素なら、

$x \equiv c \pmod{m_1 m_2}$

を満たす $c$ (と $x$ ) を求めることができます。

これが何の役に立つのか、というと、

$m = M \times 2^{d}$ に分解する ( $M$ は奇数)
とのそれぞれを計算する
1. $A$ は法が奇数なのでモンゴメリ乗算
2. $B$ は法が 2 冪なのでそのまま計算 (最後にマスクをかけるだけ)
中国の剰余定理より $A \times B \bmod (M \times 2^{d})$ を復元する

という流れで、偶数でもモンゴメリ乗算 (のようなこと) ができるようになります。

具体的な復元手順を示しましょう。
まず、 $m = M \times 2^{d}$ に分解します。 $M$ は奇数です。
このとき、 BitOperations.TrailingZeroCount が使えます。

例によって Unity では使えません。つらい。
以下の polyfill を使うか、ループ + シフトで頑張る必要があります。

private static readonly byte[] TrailingZeroLookup = new byte[64] {  
    0, 1, 59, 2, 60, 40, 54, 3, 61, 32, 49, 41, 55, 19, 35, 4, 62, 52, 30, 33, 50, 12, 14, 42, 56, 16, 27, 20, 36, 23, 44, 5, 63, 58, 39, 53, 31, 48, 18, 34, 51, 29, 11, 13, 15, 26, 22, 43, 57, 38, 47, 17, 28, 10, 25, 21, 37, 46, 9, 24, 45, 8, 7, 6  
};  
  
public static int TrailingZeroCount(ulong value)  
{  
    if (value == 0)  
    {  
        return 64;  
    }  
  
    return TrailingZeroLookup[((value & (0 - value)) * 0x03F566ED27179461) >> 58];  
}

謎のテーブルルックアップについては、 "de bruijn trailing zero count" とググってください。

次に、 $A = a \times b \bmod M$ を計算します。
$M$ は奇数なので、モンゴメリ乗算で計算できます。

続いて、 $B = a \times b \bmod 2^{d}$ を計算します。
2 冪なので、普通に $a \times b$ して最後に & ((1ul << d) - 1) のマスクをかけるだけです。

最後に、中国の剰余定理より復元を行います。
ここで、中国の剰余定理を発展させた Garner のアルゴリズムを使うと効率よく復元できます。

Garner のアルゴリズム

Garner のアルゴリズムは、

$\begin{cases} x \equiv a_1 \pmod{m_1} \\ x \equiv a_2 \pmod{m_2} \\ ... \\ x \equiv a_k \pmod{m_k} \end{cases}$

の $a_1, a_2, ..., a_k$ が既知のときに、

$x = x_1 + x_2 m_1 + x_3 m_1 m_2 + ... + m_k p_1 ... p_{k-1} \pmod {m_1 m_2 ... m_k}$

を満たす $x$ を計算するアルゴリズムです。なんのこっちゃ、と思われたかもしれませんが、要するに上記の条件を満たす最小の $x$ を計算できるということです。
なお簡単のため、 $m_1, m_2, ..., m_k$ は互いに素であるとします。

まず初項 $x_1$ は、 $a_1 \equiv x_1 \pmod {m_1}$ です。そのまま。

次に $x_2$ は、 $a_2 \equiv x_1 + x_2 m_1 \pmod {m_2}$ から求めます。
式を整理すると、

$\begin{align*} a_2 &\equiv x_1 + x_2 m_1 \pmod {m_2} \\ a_2 - x_1 &\equiv x_2 m_1 \\ (a_2 - x_1) (m_1^{-1} \bmod m_2) &\equiv x_2 \end{align*}$

となって $x_2$ が求まります。

$x_3$ も同様に、

$\begin{align*} a_3 &\equiv x_1 + x_2 m_1 + x_3 m_1 m_2 \pmod {m_3} \\ a_3 - x_1 &\equiv x_2 m_1 + x_3 m_1 m_2 \\ a_3 - x_1 &\equiv m_1 (x_2 + x_3 m_2) \\ (a_3 - x_1)(m_1^{-1} \bmod m_3) &\equiv x_2 + x_3 m_2 \\ (a_3 - x_1)(m_1^{-1} \bmod m_3) - x_2 &\equiv x_3 m_2 \\ ((a_3 - x_1)(m_1^{-1} \bmod m_3) - x_2)(m_2^{-1} \bmod m_3) &\equiv x_3 \end{align*}$

といった感じで計算していきます。

最後に、当初の式 $x = x_1 + x_2 m_1 + x_3 m_1 m_2 + ... + m_k p_1 ... p_{k-1}$ に $x_i$ を代入すれば、 $x$ を求めることができます。

それでは、 Garner のアルゴリズムを今回の問題に適用してみましょう。

まず、 $x \equiv A \pmod M$ です。したがって、 $x_1 = A \pmod M$ です。はい。
次に、 $x \equiv B \pmod{2^{d}}$ より、

$\begin{align*} x_2 &= (B - x_1) (M^{-1} \bmod 2^d) \pmod{2^{d}} \\ &= (B - A) (M^{-1} \bmod 2^d) \pmod{2^{d}} \end{align*}$

が成り立ちます。

最終的に、

$\begin{align*} x &= x_1 + x_2 M \pmod{M \times 2^d}\\ &= A + ((B - A) (M^{-1} \bmod 2^d) \bmod 2^d) M \end{align*}$

となります。

プログラムに落とし込むと以下のようになります。

public static ulong ModMul(ulong a, ulong b, ulong mod)  
{  
    // if mod is even then  
    if ((mod & 1) == 0)  
    {  
        int evenBits = TrailingZeroCount(mod);  
        ulong oddMod = mod >> evenBits;  
        var (modinv, rmod, r2mod) = MontgomeryConstant(oddMod);  
        ulong mask = (1ul << evenBits) - 1;  
  
  
        // MR(a * b)  
        ulong thi = BigMul(a, b, out ulong tlo);  
        ulong tOdd = thi - BigMul(tlo * modinv, oddMod, out _);  
        if (tOdd > thi)  
        {  
            tOdd += oddMod;  
        }  
  
        // == (a * b) & mask  
        ulong tEven = tlo & mask;  
  
        // MR(ab * R^2)  
        thi = BigMul(tOdd, r2mod, out tlo);  
        tOdd = thi - BigMul(tlo * modinv, oddMod, out _);  
        if (tOdd > thi)  
        {  
            tOdd += oddMod;  
        }  
  
        // Garner's algorithm  
        ulong t = tOdd + (((tEven - tOdd) * modinv) & mask) * oddMod;  
        return t;  
    }  
    else  
    {  
        // if mod is odd then  
        // do normal Montgomery Multiplication  
        var (modinv, rmod, r2mod) = MontgomeryConstant(mod);  
  
        var mont = MontgomeryMultiplication(MontgomeryMultiplication(  
            a, b, mod, modinv), r2mod, mod, modinv);  
        return mont;  
    }  
}

ここでポイントとなるのは t = tOdd + (((tEven - tOdd) * modinv) & mask) * oddMod の部分です。
一般に、乗法逆元 $m^{-1} \bmod 2^{d}$ の下位 $d - k$ ビットもまた乗法逆元 $m^{-1} \bmod 2^{d-k}$ となり、 $(m \times m^{-1}) \bmod 2^{d-k} \equiv 1$ を満たします。したがって、乗法逆元を掛けてからまとめてマスクをとっても問題ありません。

モンゴメリ乗算中に 2 で割る

モンゴメリ乗算において、 2 で割りたくなることはしばしばあります。つまり $a / 2 = a \times 2^{-1} \pmod m$ です。
正式には $2^{-1}R \bmod m$ を計算してモンゴメリ乗算 $MR(a \times 2^{-1}R)$ しないといけないのですが、これを簡便に求める方法があります。

public static ulong MontgomeryDivision2(ulong a, ulong mod)  
{  
    return (a & 1) != 0 ?  
        (a >> 1) + (mod >> 1) + 1 :  
        (a >> 1);  
}

$a$ が偶数のときはそのまま $a / 2$ を、奇数のときは $(a + m) / 2$ を計算しています。
$a$ が奇数のとき、 $m$ も奇数ですので最下位ビットで必ず繰り上がりが発生します。そのため +1 しているというわけです。

乗法逆元の応用

Extended Binary GCD

拡張ユークリッドの互除法では、 q の計算で除算を利用していました。
除算は一般に重い演算です。どうにかして回避する方法はないでしょうか？

実はあります。それが Extended Binary GCD アルゴリズムです。
Extended Binary GCD は "Binary" とついている通り、 2 で割ることで収束を目指します。 2 で割る、ということは高速なシフト演算に置換できるということで、高速化が望めます。

// Extended Binary GCD.  
// returns x^-1 if available, 0 otherwise  
public static ulong Ebgcd(ulong x, ulong mod)  
{  
    ulong a = x, u = 1, b = mod, v = 0;  
  
    while (a != 0)  
    {  
        if ((a & 1) == 0)  
        {  
            a >>= 1;  
            // (u / 2) mod m  
            u = MontgomeryDivision2(u, mod);  
        }  
        else  
        {  
            if (a < b)  
            {  
                (a, b) = (b, a);  
                (u, v) = (v, u);  
            }  
  
            a = (a - b) >> 1;  
            // ((u - v) / 2) mod m  
            u = MontgomeryDivision2(MontgomerySubtraction(  
                u, v, mod), mod);  
        }  
    }  
  
    if (b != 1)  
    {  
        return 0;       // x is not invertible  
    }  
  
    return v;  
}

もう一つの利点として、 u, v に算術オーバーフローが発生しない計算だけで構成されているので、拡張ユークリッドの互除法のように結果が壊れる心配をしなくてよいことが挙げられます。

なんでこれを最初に紹介しなかったのかというと、 Binary GCD 法は収束が遅くループ回数がかさむため、ナイーブに実装すると普通の拡張ユークリッドの互除法のほうが速い場合が多いためです。

ナイーブに、といったということは、カリカリにチューニングすれば高速になる可能性はあります。

Hybrid Extended GCD

これがチューニング例です。名前の通り、通常の拡張ユークリッドの互除法と Binary GCD 法のいいとこどりをした手法です。

Greatest common divisor, the extended Euclidean algorithm, and speed! – Daniel Lemire's blog

こちらで提案されていた GCD 用のアルゴリズムを拡張して逆元計算に使えるようにしたものがこちらになります。

public static ulong ExtendedHybridGcd(ulong x, ulong mod)  
{  
    ulong a = mod, b = x;  
    ulong u = 0, v = 1;  
  
    if (b == 0)  
    {  
        return 0;  
    }  
  
    // returns (a / b, a % b)  
    (ulong div, a) = Math.DivRem(a, b);  
    //u = u - v * div + mod;  
    u = mod - div;      // assumes u == 0, v == 1  
  
    if (a == 0)  
    {  
        return 0;  
    }  
  
    int ash = BitOperations.TrailingZeroCount(a) & 63;  
    int bsh = BitOperations.TrailingZeroCount(b) & 63;  
    if (bsh > 0 && ash > 0)  
    {  
        return 0;       // gcd >= 2  
    }  
  
    a >>= ash;  
    b >>= bsh;  
  
    for (int i = 0; i < ash; i++)  
    {  
        u = MontgomeryDivision2(u, mod);  
    }  
    for (int i = 0; i < bsh; i++)  
    {  
        v = MontgomeryDivision2(v, mod);  
    }  
  
    do  
    {  
        ulong aminusb = a - b;  
        ulong uminusv = MontgomerySubtraction(u, v, mod);  
  
        if (a > b)  
        {  
            (a, b) = (b, aminusb);  
            (u, v) = (v, uminusv);  
        }  
        else  
        {  
            b = b - a;  
            v = mod - uminusv;  
        }  
  
        bsh = BitOperations.TrailingZeroCount(aminusb) & 63;  
        b >>= bsh;  
  
  
        for (int i = 0; i < bsh; i++)  
        {  
            v = MontgomeryDivision2(v, mod);  
        }  
    } while (b != 0);  
  
    if (a != 1)  
    {  
        return 0;       // gcd(x, mod) == a > 1  
    }  
  
    return u;  
}

Method	Mean	Error	StdDev
ExtendedHybridGcd	50.71 ns	0.464 ns	0.387 ns
ExtendedGcd	104.09 ns	0.601 ns	0.562 ns
ExtendedBinaryGcd	129.84 ns	0.893 ns	0.836 ns

安直な実装に比べると 2 倍ほど速くなっています。
ただ、値を変えると結構性能がぶれることがあるので、お手元の環境で実際に測定してみてください。

フェルマーの小定理

乗法逆元を求めるにあたり、もう一種類やり方があります。
それはフェルマーの小定理を応用するものです。フェルマーの小定理というのは、 $m$ が素数のとき、

$a^{m} \equiv a \pmod{m}$

を満たす、という定理です。これを変形すると、

$a^{m-2} \equiv a^{-1} \pmod {m}$

となります。

ただ、 $m$ が素数という厳しめの条件がつくことと、繰り返し二乗法をもってしても多数の乗算が必要であることから、実用上はそんなに嬉しくないです。
しいて言えば BigInteger.ModPow(a, mod - 2, mod) とワンライナーで書けることぐらいでしょうか。

乗法逆元を使った倍数判定

$n$ が $d$ の倍数である、という判定は、ふつうは n % d == 0 でとると思います。
ですがこれは剰余算 (=除算) を含む、それなりに重い計算になります。もっと軽くしたいところです。

n % d == 0 と等価な計算を以下に示します。なお、 $d$ は奇数とします。

// n % d == 0?  
public static bool IsDiv(ulong n, ulong d)  
{  
    return n * MultiplicativeInverse(d) <= ~0ul / d;  
}

数式で書くと $n \times (d^{-1} \bmod 2^{64}) \le (2^{64} - 1) / d$ です。
ここで $d$ が定数なら、事前に乗法逆数と ~0ul / d を計算して埋め込んでおくことで高速化が見込めます。

実はこれ、コンパイラも似たようなことをやっています。
sharplab で覗いてみると、

public class C {  
    public bool M(ulong n) {  
        return n % 3 == 0;  
    }  
}

C.M(UInt64)  
    L0000: mov rcx, rdx  
    L0003: mov rdx, 0xaaaaaaaaaaaaaaab  
    L000d: mov rax, rcx  
    L0010: mul rdx  
    L0013: shr rdx, 1  
    L0016: lea rax, [rdx+rdx*2]  
    L001a: sub rcx, rax  
    L001d: sete al  
    L0020: movzx eax, al  
    L0023: ret

アセンブリを C# に直訳すると、

public bool M(ulong n)  
{  
    return (n - ((n * 0xaaaaaaaaaaaaaaabul) >> 1) * 3) == 0;  
}

0xaaaaaaaaaaaaaaabul は大体 $2^{64} \times 2 / 3$ なので、 $\lfloor (n \times 2/3) /2 \rfloor \times3 == n$ のような計算をしている、とみると理解しやすそうです。

当然見にくくなるので積極的におすすめはしませんが、一分一秒、一ナノ秒を争う状況では使える知識かもしれません。

実践: Miller-Rabin 素数判定法

冪剰余 $a^{b} \bmod m$ 以外にもモンゴメリ乗算が役立つ手法の一例として、今回は Miller-Rabin 素数判定法を紹介します。

Miller-Rabin 素数判定法は、確率的素数判定アルゴリズムです。名前の通り、素数かどうかを判定するために利用できます。
確率的、と言ったように、本来は乱数を生成して判定に利用します。一回の判定で最大 $1/4$ の確率で間違う (合成数を素数として判定する可能性がある) ので、例えば 32 回やって $1/4^{32}$ みたいに現実的に問題ないレベルまで持っていく、みたいなことをします。
しかし、 $n \lt 2^{64}$ までという制約をかけて、注意深く選ばれた数値をもとに計算すると 確実に素数かそうでないかを判定する ことができます。 *2
(私は最近知ってとても驚きました。愚直に乱数で処理していたのはいったいなんだったのか……)

まずはアルゴリズムの紹介がてら、ソースコードを示します。

// returns true if value is prime  
public static bool MillerRabin(ulong value)  
{  
    // 2 以下・偶数の場合の前処理  
    if (value <= 2)  
    {  
        return value == 2;  
    }  
    if ((value & 1) == 0)  
    {  
        return false;  
    }  
  
    ulong n1 = value - 1;  
    int s = TrailingZeroCount(n1);  
    ulong d = n1 >> s;  
  
    // 「注意深く選ばれた」定数  
    // 2^64 までならこの 7 つで対応可能  
    ReadOnlySpan<ulong> MillerRabinConstants =   
        [2, 325, 9375, 28178, 450775, 9780504, 1795265022];  
  
    foreach (var a in MillerRabinConstants)  
    {  
        // 割り切れた場合は続行  
        if (a % value == 0)  
        {  
            continue;  
        }  
  
        // a^d mod value == 1 or value-1 なら続行  
        ulong t = ModPow(a, d, value);  
        if (t == 1 || t == n1)  
        {  
            continue;  
        }  
  
        // t^(2^i) mod value について、  
        // == value-1 の要素があれば続行、なければ合成数  
        int i;  
        for (i = 0; i < s; i++)  
        {  
            t = ModMul(t, t, value);  
            if (t == n1)  
            {  
                break;  
            }  
        }  
        if (i == s)  
        {  
            return false;  
        }  
    }  
  
    // 全てのテストを通過したら素数  
    return true;  
}

さて、ここでは冪剰余 $a^{d} \bmod m$ と $t^{2} \bmod m$ が登場します。これをモンゴメリ乗算を使って高速化しよう、というわけです。

ModPow と ModMul を安直に BigInteger で計算した場合と、ちゃんとモンゴメリ乗算した場合とで比較してみましょう。

BenchmarkDotNet v0.14.0, Windows 11 (10.0.22631.5189/23H2/2023Update/SunValley3)  
12th Gen Intel Core i7-12700F, 1 CPU, 20 logical and 12 physical cores  
.NET SDK 10.0.100-alpha.1.24623.5  
  [Host]     : .NET 10.0.0 (10.0.24.62010), X64 RyuJIT AVX2  
  DefaultJob : .NET 10.0.0 (10.0.24.62010), X64 RyuJIT AVX2

Method	Mean	Error	StdDev
MillerRabin	13.265 us	0.1002 us	0.0888 us
MillerRabin_Montgomery	2.732 us	0.0289 us	0.0271 us

実に 5 倍弱の改善となったことが分かります。

余談

$n \lt 2^{64}$ の素数判定においては、 Baillie-PSW 素数判定法というさらに高速な手法があります。
ちょっと実装は大変ですが、 MillerRabin_Montgomery の 2 倍以上高速に (私の環境では約 1μs で) 判定することができます。
以下の記事に詳しいので、気になる方は調べてみてください。

64bit数の素数判定

その他の応用

例えば、楕円曲線法 (Elliptic Curve Method; ECM) による素因数分解では、そもそもが $\bmod m$ の世界での理論のため、モンゴメリ乗算を呼吸するかのように使いこなしていたりします。

おわりに

モンゴメリ乗算、理解できるとめちゃくちゃに便利なので、今回知ることができて良かったです。
本記事を通じて知ったという方がいらっしゃれば、理解の助けになっていれば幸いです。

*1:Hurchalla, J. (2022). An Improved Integer Modular Multiplicative Inverse (modulo $2^ w$ ). arXiv preprint arXiv:2204.04342.

*2:https://miller-rabin.appspot.com/

2025-05-02

TextMeshPro で自作画像をフォントにする

はじめに

Unity の TextMeshPro で、自分が作った文字画像を表示したいと思ったことはありませんか？私はあります。
その際、どういった手順で作業すればいいかを書き残しておきます。

※絵文字とかの Sprite 埋め込みとは異なります。

以下の画像を、

こうしたい、ということです。

手順

テクスチャを描く

まずはお好みのツールで文字テクスチャを描きます。
今回は ASCII を想定しているので、 256x256 の画像を 16x16 ずつに区切り、その中に該当する文字コードを一文字ずつ描きます。

注意すべき点としては、必ず文字の周囲に 2px 以上のマージンを設けることです。
この例では 16x16 から上下左右マージン 4px を設けて 8x8 の範囲に描くようにしています。
ぎちぎちに詰めて描きたくなるものですが、そうすると描画時に問題になる可能性が高いです。

Unity へインポート

テクスチャができたら .png で出力し、 Unity でインポートします。
この時、以下のようにインポート設定を変更しておきます。

Texture Type: Default
Alpha is Transparency: ✅
Generate Mipmap: 🟩

FontAsset 生成

次に、適当なフォント (なんでもよい、Liberation Sans でも可) を用意して、通常通り FontAsset を生成します。
具体的には、対象のフォントファイルを右クリック→ Create/TextMeshPro/FontAsset/Bitmap です。

そうしたら、生成された FontAsset の Inspector から Update Atlas Texture を押し、以下の設定にします。

Sampling Point Size: Custom Size, ＜作った文字の高さ＞px
Padding: 2px
Atlas Resolution: ＜作ったテクスチャのサイズ＞
Character Set: Custom Characters
Custom Character List: 20-7E
Render Mode: RASTER_HINTED

できたら Generate Font Atlas ボタンを押し、生成して Save します。

そうしたら、 FontAsset の Inspector から、

Generation Settings/Atlas Population Mode: Static

に変更しておきます。

この時点でアトラスがぐちゃぐちゃだったとしても、あとで置換するので問題ありません。

テクスチャを上書き

Inspector タブを右クリック→ Debug を選択し、デバッグメニューを開きます。

Atlas Textures: ＜作ったテクスチャ＞

に変更したら、Inspector タブを右クリック→ Normal に戻します。

次に、 FontAsset の子の Material を選択して、

Debug Settings/Font Atlas: ＜作ったテクスチャ＞

に変更します。

文字情報の置換

この時点では各テクスチャの座標と文字情報の紐づけが行われていないので、紐づけを行います。

以下のスクリプトを適当な Editor フォルダ以下に作成します。

// FontAssetHandCreator.cs  
using TMPro;  
using UnityEditor;  
using UnityEngine.TextCore;  
  
namespace OperatorOverload.Editor.Serialization  
{  
    public class FontAssetHandCreator  
    {  
        [MenuItem("CONTEXT/TMP_FontAsset/Create Hand")]  
        public static void CreateFontAsset(MenuCommand menuCommand)  
        {  
            int textureWidth = 256;     // テクスチャの幅(px)  
            int textureHeight = 256;    // テクスチャの高さ(px)  
            int characterAreaWidth = 16;    // ひとつの文字領域の幅(px)  
            int characterAreaHeight = 16;   // ひとつの文字領域の高さ(px)  
            int glyphX = 4; // 文字の x オフセット  
            int glyphY = 4; // 文字の y オフセット  
            int glyphWidth = 3; // 文字の幅  
            int glyphHeight = 5;    // 文字の高さ  
            float horizontalMargin = 1; // 文字間の間隔  
  
            var fontAsset = menuCommand.context as TMP_FontAsset;  
            fontAsset!.glyphTable.Clear();  
            fontAsset.characterTable.Clear();  
  
            for (int i = 0x20; i < 0x7f; i++)  
            {  
                var glyph = new Glyph((uint)i, new GlyphMetrics(glyphWidth, glyphHeight, 0, glyphHeight, glyphWidth + horizontalMargin),  
                    new GlyphRect(  
                        i % (textureWidth / characterAreaWidth) * characterAreaWidth + glyphX,  
                        textureHeight - (i / (textureWidth / characterAreaWidth) * characterAreaHeight + characterAreaHeight - glyphY),  
                        glyphWidth, glyphHeight));  
                var character = new TMP_Character((uint)i, glyph);  
  
                fontAsset.glyphTable.Add(glyph);  
                fontAsset.glyphLookupTable[glyph.index] = glyph;  
                fontAsset.characterTable.Add(character);  
                fontAsset.characterLookupTable[character.unicode] = character;  
            }  
  
            EditorUtility.SetDirty(fontAsset);  
        }  
    }  
}

そうしたら、 FontAsset の Inspector の上のほうを右クリック→ Create Hand を選択します。
すると文字情報が生成されるので、忘れずに保存します。

試してみる

うまくいっていれば、この時点でフォントが使えるようになっているはずです。

やったね！

おわりに

なんかうまくいかなくて n ヵ月放置→試してみる→ を繰り返すこと 3 回、ようやく実現できました。
また忘れそうなので残しておきます。

2025-02-21

.NET 環境と Unity 環境では ArrayPool の上限が異なる

はじめに

ArrayPool<T> は、 T 型の配列 T[] をいい感じにプールしてくれる機構です。 GC ゴミを削減したいときによくお世話になるクラスです。

ところで、 ArrayPool<T>.Shared.Rent() で得られるバッファは「指定した長さ以上」になることが明記されています。
以上、というのが具体的にいくつかというと、 基本的には 2 冪に切り上げた値 (e.g. 100 なら 128 、 200 なら 256 。最低 16) になります。 BitOperations.RoundUpToPowerOf2() と同じような感じです。

さて、 基本的には と書いた以上、例外条項が存在するということです。

いつ 2 冪以外の配列を返すのか

まず、 0 を指定した場合は空配列 ( Array.Empty<T>() ) への参照を返します。それはそう。もちろん本題はこれではないです。

.NET 9 環境

.NET 9 環境では、 $2^{30}$ (1073741824; 10 億ちょい) 以上の値を指定した場合はそれと同じ長さの配列を返します。
なぜかというと、内部プールで対応できる長さの上限が 1 << 30 のためです。
現行の実装の const int NumBuckets = 27; がそれに相当していて、 1 << (27 - 1 + 4) == 1 << 30 になります。 +4 は 16 のぶん。

それより大きな値を指定した場合は単に配列を確保するのと同じような挙動になります。 ( GC.AllocateUninitializedArray<T>() が使われるため、未初期化であることと引き換えに new T[] よりは多少速いです。)
そのような巨大な配列も合法的に Return() することはできますが、プールされずにそのまま GC 送りになります。

Unity 環境

一方で問題の Unity 環境 (Unity 6000.1.0b1) です。
Unity 環境では、 $2^{20}$ (1048576) 以上の値を指定した場合はそれと同じ長さの配列を返すようです。低すぎる！！！

Unity が採用している .NET Standard 2.1 が生まれたころの実装 (2018 年ごろ、多分このあたり ) を確認すると、確かに DefaultMaxArrayLength = 1024 * 1024 とかいう記述があります。

当然ながらそれより大きい値を指定すると新規確保されたうえで GC 送りの刑に処されます。
MiB 単位のバッファを借りることはしばしばあると思うのですが、 GC 削りのために ArrayPool<T> を使っていても一切無駄になります。ひどい！

余談

なんでこれに気づいたのかというと、確保したバッファが Unity 環境で 2 冪じゃなくて処理が崩壊したからです。

具体的に言うと、借りたバッファを拡大するときにこういう感じの処理を入れると思うのですが、

public void ResizeIfOver(int addSize)  
{  
    if (currentSize + addSize > currentArray.Length)  
    {  
        var oldArray = currentArray;  
        currentArray = ArrayPool<byte>.Shared.Rent(oldArray.Length + addSize);  
        oldArray.AsSpan().CopyTo(currentArray);  
        ArrayPool<byte>.Shared.Return(oldArray);  
    }  
}

このとき、 oldArray.Length + addSize が暗黙に 2 冪に切り上げられると仮定していました。
サイズが 2 冪で上がっていくので、拡大処理自体はあまり発生しない ( $O(\log n)$ みたいな感じ) になると思っていました。

ところがどっこい、 ArrayPool<T> の上限を超えてしまっていた場合、新規確保したバッファはジャスト oldArray.Length + addSize なので、拡大が毎回発生するうえ ( $O(\log n)$ に対して $O(n)$ みたいな) 、確保したバッファは全部 GC 送りされるので、実行が絶望的に遅くなります。

おわりに

ドキュメントに書いていない仕様を信じすぎてはいけない。それはそう。
というわけで気を付けましょう。生兵法は怪我のもと……

2025-02-14

MemoryPack をもっと効率よく使う～ IBufferWriter ・ ReadOnlySequence とは

はじめに

MemoryPack という爆速のシリアライゼーションライブラリがあります。 Cysharp さんには毎度お世話になっております。

ところで、ヘルプの Serialize API の項を見てみると、

byte[] Serialize<T>(in T? value, MemoryPackSerializerOptions? options = default)  
void Serialize<T, TBufferWriter>(in TBufferWriter bufferWriter, in T? value, MemoryPackSerializerOptions? options = default)

For performance, the recommended API uses BufferWriter. This serializes directly into the buffer. It can be applied to PipeWriter in System.IO.Pipelines, BodyWriter in ASP .NET Core, etc.

(訳) パフォーマンスのために、 BufferWriter を使用する API が推奨されます。これはバッファに直接シリアライズします。
System.IO.Pipelines.PipeWriter や、 ASP .NET Core の BodyWriter が該当します。

とあります。
要するに、パフォーマンスのためには byte[] ではなく BufferWriter を使うオーバーロードを使用するのが良い、とされています。
ですが、この BufferWriter とは何でしょうか？どうやって使えばいいのでしょうか？
このあたりはヘルプに載っていなかったので、調べてみました。
(実は常識なのかもしれませんが……)

`IBufferWriter<T>` とは

System.Buffers.IBufferWriter<T> は、 <T> 型のデータを書き込むためのインターフェイスです。それはそう。

インターフェイスの定義としてはこんな感じです。

namespace System.Buffers;  
  
public interface IBufferWriter<T>  
{  
    // `count` 個のデータが書き込まれたことを通知する  
    void Advance(int count);  
  
    // `sizeHint` 個以上のバッファを返す (0 のときは 1 以上)  
    Memory<T> GetMemory(int sizeHint = 0);  
  
    // `sizeHint` 個以上のバッファを返す (0 のときは 1 以上)  
    Span<T> GetSpan(int sizeHint = 0);  
}

これを使う側としては、例えばこういう感じになります。

// 何かのソースから writer を取得  
IBufferWriter<byte> writer = GetSomeWriter();  
  
// GetSpan で書き込みバッファを取得  
var span = writer.GetSpan(sizeHint: 16);  
  
// バッファに書き込む (`span.Length` は 16 "以上" なことに注意)  
for (int i = 0; i < 16; i++)  
{  
    span[i] = SomeData();  
}  
  
// 書き込んだデータ数の分だけ Advance する  
writer.Advance(16);  
  
// 同様にして、次の書き込みを行う...

そこまで難しくはありませんね。
Span<T> (あるいは Memory<T>) をもらってそこに直接書き込んでは Advance() する、を繰り返すだけです。

IBufferWriter<T> を実装しているクラスとしては、 System.Buffers.ArrayBufferWriter<T> と System.IO.Pipelines.PipeWriter があります。

`ArrayBufferWriter<T>`

ArrayBufferWriter<T> は、名前の通り配列をバックに持つ IBufferWriter<T> の実装です。
書き込んだデータは .WrittenMemory または .WrittenSpan プロパティから取得することができます。

例として、 MemoryPack を用いてディープコピーを行うコードはこんな感じになります。

var writer = new ArrayBufferWriter<byte>();  
  
// `Source` をシリアライズして `writer` に書き込む  
MemoryPackSerializer.Serialize(writer, Source);  
  
// 書き込んだデータを `WrittenSpan` から取得してデシリアライズ  
return MemoryPackSerializer.Deserialize<TestClass>(writer.WrittenSpan)!;

`PipeWriter`

PipeWriter は、データを書き込めるパイプラインを提供するクラスです。
これ自身は抽象クラスなので、実際には Pipe を経由して使うことになります。

例として、同じくディープコピーを行うコードはこんな感じになります。

// `pipe` を作成  
var pipe = new Pipe();  
  
// `pipe.Writer` で `PipeWriter` を取得、そこにシリアライズ  
MemoryPackSerializer.Serialize(pipe.Writer, Source);  
  
// Flush して読み込めるようにする  
await pipe.Writer.CompleteAsync();  
  
// `PipeReader` から読み込み結果を取得  
ReadResult readResult = await pipe.Reader.ReadAsync();  
  
// 結果のバッファ (`readResult.Buffer`) からデータを読み込んでデシリアライズ  
return MemoryPackSerializer.Deserialize<TestClass>(readResult.Buffer)!;

async/await が増えてちょっと重たい感じがしますね。
本来は通信とかファイル I/O などに使うもののようです。たぶん。

アロケーション問題

さて、 ArrayBufferWriter<T> は配列をバックに持っている、と書きました。
それでは、内部配列を超える量を書き込んだらどうなるのでしょうか？

ソースを見てみると、 Array.Resize() によってサイズを変更していることが分かります。
したがって、配列が拡大するたびに既存要素のコピーと新規配列のアロケーションが発生します。当たり前と言えば当たり前です。

Pipe については、今回やりたいことに対して結構大仰なことをやっている気がします。非同期操作も必要になりますし。

とはいえ体感とか予想の話なので、実際にベンチマークをとってみましょう。

[MemoryDiagnoser]  
public class BufferWriter  
{  
    private readonly TestClass Source = new TestClass(123, "Alice", "GDD", Enumerable.Range(0, 100).Select(i => i.ToString()).ToArray());  
  
    [Benchmark]  
    public TestClass ByteArray()  
    {  
        var bytes = MemoryPackSerializer.Serialize(Source);  
        return MemoryPackSerializer.Deserialize<TestClass>(bytes)!;  
    }  
  
    [Benchmark]  
    public async ValueTask<TestClass> Pipe()  
    {  
        var pipe = new Pipe();  
        MemoryPackSerializer.Serialize(pipe.Writer, Source);  
        await pipe.Writer.CompleteAsync();  
        var readResult = await pipe.Reader.ReadAsync();  
        return MemoryPackSerializer.Deserialize<TestClass>(readResult.Buffer)!;  
    }  
  
    [Benchmark]  
    public TestClass ArrayBufferWriter()  
    {  
        var writer = new ArrayBufferWriter<byte>();  
        MemoryPackSerializer.Serialize(writer, Source);  
        return MemoryPackSerializer.Deserialize<TestClass>(writer.WrittenSpan)!;  
    }  
  
    [Benchmark]  
    public TestClass With()  
    {  
        return Source with {   
            Name = new string(Source.Name),   
            Address = new string(Source.Address),   
            Status = Source.Status.Select(i => new string(i)).ToArray() };  
    }  
}  
  
[MemoryPackable]  
public partial record class TestClass(  
    int Id, string Name, string Address, string[] Status)  
{  
}

BenchmarkDotNet v0.14.1-nightly.20250107.205, Windows 11 (10.0.22631.4751/23H2/2023Update/SunValley3)  
12th Gen Intel Core i7-12700F 2.10GHz, 1 CPU, 20 logical and 12 physical cores  
.NET SDK 10.0.100-alpha.1.24623.5  
  [Host]     : .NET 9.0.2 (9.0.225.6610), X64 RyuJIT AVX2  
  Job-FLCZHK : .NET 9.0.2 (9.0.225.6610), X64 RyuJIT AVX2  
  
Affinity=00001111111111111111

Method	Mean	Error	StdDev	Gen0	Gen1	Allocated
ByteArray	2,550.8 ns	36.54 ns	34.18 ns	0.3891	0.0038	4.98 KB
Pipe	2,940.1 ns	43.75 ns	38.78 ns	0.6638	0.0114	8.5 KB
ArrayBufferWriter	2,773.8 ns	34.34 ns	30.44 ns	0.6104	0.0076	7.84 KB
With	963.5 ns	18.81 ns	30.91 ns	0.3138	0.0048	4.01 KB

ByteArray: byte[] を使ったシリアライズ・デシリアライズ
Pipe: Pipe を使ったシリアライズ・デシリアライズ
ArrayBufferWriter: ArrayBufferWriter<byte> を使ったシリアライズ・デシリアライズ
With: with でコピーしたやつ (理論値比較用)

ベンチマークより、 byte[] 経由のディープコピーよりも Pipe ・ ArrayBufferWriter を使ったほうがアロケーションが多くなっています。
推測ですが、やはり内部バッファの拡大に伴い無駄なアロケーションが増えているのだと思います。
逆に言えば、 byte[] のほうは MemoryPack 内部でかなり最適化されている (最終的な byte[] 以外には無駄がない) ことが推測されます。さすがです。

さらに言えるのが、全て理論値 (with によるコピー) よりもアロケーションが多いということです。
できることなら、無駄なアロケーションは減らしたいものですね。

`IBufferWriter<T>` を実装してみよう

さて、 GC 、ひいてはアロケーションを減らすといえば、 ArrayPool<T> ですね。
一応説明すると、配列をよしなにプールして再利用することでアロケーションやガベージを減らせる機構です。
単なる配列の代わりに ArrayPool<T> から貸し出された配列をバックに持つ IBufferWriter<T> を実装すれば、アロケーション問題を解決できそうですね。

`ArrayPoolBufferWriter<T>` の実装

というわけで実装してみました。

特に難しいところはありませんね。
GetSpan() (または GetMemory()) で内部バッファよりも大きい量をリクエストされたら、 ArrayPool<T>.Shared.Rent() で新規バッファを用意→旧バッファのデータをコピー→旧バッファを ArrayPool<T>.Shared.Return() で返却、として拡大します。

実装時のポイントとしては、 GetSpan() で返却するバッファのサイズを「sizeHint ジャスト」にするのではなく「現状返せるキャパシティの限界まで」とする点が挙げられます。
できるだけ大きなバッファを返すことで、呼び出し側が GetSpan() を呼ぶ回数を減らすことができます。
(sizeHint はあくまでヒントなので、 sizeHint 以上の要素にも合法的に書き込むことができます。)

それでは、ベンチマークをとってみましょう。

    [Benchmark]  
    public TestClass ArrayPoolBufferWriter()  
    {  
        using var writer = new ArrayPoolBufferWriter<byte>(1024);  
        MemoryPackSerializer.Serialize(writer, Source);  
        return MemoryPackSerializer.Deserialize<TestClass>(writer.WrittenSpan)!;  
    }

Method	Mean	Error	StdDev	Gen0	Gen1	Allocated
ByteArray	2,550.8 ns	36.54 ns	34.18 ns	0.3891	0.0038	4.98 KB
Pipe	2,940.1 ns	43.75 ns	38.78 ns	0.6638	0.0114	8.5 KB
ArrayBufferWriter	2,773.8 ns	34.34 ns	30.44 ns	0.6104	0.0076	7.84 KB
👉 ArrayPoolBufferWriter	2,478.5 ns	27.44 ns	24.32 ns	0.3090	0.0038	3.99 KB
With	963.5 ns	18.81 ns	30.91 ns	0.3138	0.0048	4.01 KB

アロケーションに関しては理論値 (With) よりも小さくすることができました。
(With のほうが大きいのは .Select() によるものだと思われます。)
速度も byte[] のものより速くなっています。やったね！

バッファ拡大時のコピーを回避する

ここでひとつ気になる点があるとすれば、バッファ拡大時の処理です。
現状、バッファの拡大を行う際には、

旧バッファの 2 倍以上のサイズを持った新バッファの貸出
旧バッファから新バッファへ、データのコピー
旧バッファの返却

といった手順で行っています。
ここで、「データのコピー」については改善の余地がありそうです。

単一の配列を持つデータ構造だとこれはどうしようもないのですが、複数の配列からなるデータ構造なら、このコピーを回避することができそうです。

block-beta  
columns 1  
block:id1  
    a1["Length: 128"]  
    a2[" "]  
end  
space  
block:id2  
    b1["Length: 256"]  
end  
a1 -- "copy" --> b1

このように、結果が単一の配列だと新しい配列へのコピーは避けられませんが、

block-beta  
columns 1  
block:id1  
    a1["Length: 128"]  
    a2[" "]  
end  
space  
block:id2  
    b1["Length: 128"]  
    b2["Length: 256"]  
end  
a1 -- "==" --> b1

結果が複数の配列でいいなら、以前の配列をそのまま使いまわすことができるのでコピーを削減できます。

ですが、複数の配列なんて受け付けてくれるのでしょうか？
ここで MemoryPack の Deserialize API を見てみましょう：

T? Deserialize<T>(ReadOnlySpan<byte> buffer)  
int Deserialize<T>(ReadOnlySpan<byte> buffer, ref T? value)  
T? Deserialize<T>(in ReadOnlySequence<byte> buffer)  
int Deserialize<T>(in ReadOnlySequence<byte> buffer, ref T? value)  
async ValueTask<T?> DeserializeAsync<T>(Stream stream)

ReadOnlySpan<byte> はいつものやつですが、 ReadOnlySequence<byte> を受け付けるオーバーロードもあります。

`ReadOnlySequence<T>` とは

ReadOnlySequence<T> は、文字通り読み取り専用シーケンスを表す構造体です。
単一の連続した領域を表す ReadOnlySpan<T> とは異なり、複数の領域からなるバッファを表すことができます。

内部実装はというと、概ね ReadOnlyMemory<T> の LinkedList のような実装になっています。
ノードは ReadOnlySequenceSegment<T> です。

図にすると、こういう感じになっています。

block-beta  
columns 9  
block  
    a["Segment\nLength: 128"] space  
    b["Segment\nLength: 256"] space  
    c["Segment\nLength: 512"] space  
    d["Segment\nLength: 1024"] space  
    e["Segment\nLength: 2048"]  
end  
a --"Next"--> b  
b --"Next"--> c  
c --"Next"--> d  
d --"Next"--> e

ここではサイズを 2 冪にしていますが、実際はなんでも大丈夫です。

なお、 ReadOnlySequenceSegment<T> は抽象クラスのため、自分で実装するときは具象クラスを作る必要があります。
とはいえ、手を入れる必要はほぼなく、単に継承すればそれだけで済みます。

public sealed class SequenceSegment<T> : ReadOnlySequenceSegment<T>  
{  
    internal SequenceSegment(  
        ReadOnlyMemory<T> memory,  
        ReadOnlySequenceSegment<T>? next,  
        int runningIndex)  
    {  
        Memory = memory;  
        Next = next;  
        RunningIndex = runningIndex;  
    }  
}

重ねて言いますが、 ReadOnlySequenceSegment<T> は抽象クラスです。
必然的にアロケーションが生じる設計となってしまっています。かなしい。

`ArrayPoolSegmentedBufferWriter<T>` の実装

というわけで実装してみました。

ポイントとしては、セグメントのサイズを 2 冪で増やしていく (最初が 64 なら次は 128, 256, 512, ...) ようにしたことです。
セグメントを細かく切りすぎても良くなさそうなのと、 ArrayPool<T> の性質上あまり同じサイズのバッファを同時かつ大量に借りるのは良くない説があるのとでこのようにしました。

ただ、あまり大きくすると Large Object Heap 送りされてそれはそれでよろしくない、みたいな話はあります。
比較的小さいバッファで済むとわかっているなら、均等に 64 KiB ぐらいで切っておく手もあるでしょう。
ただ、今回の目的 (コピーの回避) の真価は大きな配列に対して発揮されるものなので、このようにしました。

ベンチマークをとってみましょう。

    [Benchmark]  
    public TestClass ArrayPoolSegmentedBufferWriter()  
    {  
        using var writer = new ArrayPoolSegmentedBufferWriter<byte>(1024);  
        MemoryPackSerializer.Serialize(writer, Source);  
        return MemoryPackSerializer.Deserialize<TestClass>(writer.GetWrittenSequence())!;  
    }

Method	Mean	Error	StdDev	Gen0	Gen1	Allocated
ByteArray	2,550.8 ns	36.54 ns	34.18 ns	0.3891	0.0038	4.98 KB
Pipe	2,940.1 ns	43.75 ns	38.78 ns	0.6638	0.0114	8.5 KB
ArrayBufferWriter	2,773.8 ns	34.34 ns	30.44 ns	0.6104	0.0076	7.84 KB
ArrayPoolBufferWriter	2,478.5 ns	27.44 ns	24.32 ns	0.3090	0.0038	3.99 KB
👉 ArrayPoolSegmentedBufferWriter	2,621.5 ns	36.11 ns	33.77 ns	0.3204	-	4.09 KB
With	963.5 ns	18.81 ns	30.91 ns	0.3138	0.0048	4.01 KB

結果としては、それほど改善しませんでした。
コピーよりも非連続バッファを作ったり管理したりするオーバーヘッドが大きかった可能性があります。
あとアロケーションも若干ながら増えていますね。これは SequenceSegment がクラスなためかと思われます。

ただ、本領を発揮できるのは大きな配列を確保した場合 (= 大きなコピーが発生すると予想される場合) でしょう。試しに、 TestClass の string[100] を string[1000000] にして実験してみました。

Method	Mean	Error	StdDev	Gen0	Gen1	Gen2	Allocated
ByteArray	77.44 ms	1.513 ms	2.072 ms	3714.2857	3571.4286	714.2857	58.26 MB
Pipe	87.51 ms	1.739 ms	2.438 ms	4833.3333	4666.6667	833.3333	58.72 MB
ArrayBufferWriter	76.28 ms	1.478 ms	1.582 ms	4000.0000	3833.3333	1000.0000	77.01 MB
ArrayPoolBufferWriter	77.41 ms	1.490 ms	2.530 ms	3714.2857	3571.4286	714.2857	45.01 MB
ArrayPoolSegmentedBufferWriter	75.43 ms	1.480 ms	1.872 ms	3714.2857	3571.4286	714.2857	45.02 MB
With	61.49 ms	1.224 ms	1.941 ms	3777.7778	3666.6667	777.7778	45.01 MB

予想通り、 ArrayPoolSegmentedBufferWriter が最速という結果になりました。
ただし、 10 KiB ほどアロケーションが増えている点に注意が必要です。

余談

色々調べたところ、 dotNext ライブラリに ArrayPoolBufferWriter<T> とほぼ同じようなものがあるようです (PoolingArrayBufferWriter) 。
やはり皆考えることは同じですね。

なら公式に実装してくれ、というのはまた別のお話……

おわりに

IBufferWriter<T> のことがちょっとわかりました（原義）。
ついでに ReadOnlySequence<T> についても少し理解が深まりました。

本コードは GitHub で公開しているほか、本稿で作った ArrayPoolBufferWriter<T> ・ ArrayPoolSegmentedBufferWriter<T> は ArrayPoolCollection ライブラリに同梱されています。
このライブラリを使うと ArrayPool<T> を活用したコレクションやプール機構が使えるようになるので、良ければぜひ使ってあげてください。

2024-12-30

高速で均等なシャッフル手法～乱数を絞りつくす編

はじめに

配列の要素をランダムな順に並び替える「シャッフル」は、ゲーム分野 (ポーカーや麻雀など) のみならず、機械学習の前処理など、幅広く利用されています。
本稿では、このシャッフルを均等に、かつ高速に行う手法について紹介します。

ランダムソート

まずは悪い例から順を追って見ていきましょう。

こういったコードを見た、あるいは書いた経験のある方も多いかと思います。

static IEnumerable<T> LinqShuffle<T>(IEnumerable<T> source) {  
    return source.OrderBy(_ => Random.Shared.Next());  
}

乱数をキーとしてソートすることでシャッフルするという、とても簡単に実装できるコードです。

乱数の代わりに Guid を使うパターンもありますね。
本質的には乱数と同じです。

static IEnumerable<T> LinqShuffleGuid<T>(IEnumerable<T> source) {  
    return source.OrderBy(_ => Guid.NewGuid());  
}

ただ、この手法は結構非効率です。具体的に挙げていきましょう。

ソートなので遅い・均等ではない

OrderBy は .NET 9 時点では安定なクイックソートで実装されています。
詳しい方はクイックソートは安定ソートではないのでは？と思われたかもしれませんが、この実装では等値だった場合にはインデックスの差を返す比較関数を用いることで安定ソートにしてあります。
クイックソートの平均計算量は $O(n \log n)$ です。速いと言えば速いですが、 $O(n)$ で計算したいところです。

OrderBy は安定ソートのため、キー (ここでは Random.Shared.Next() の値) が重複した場合にもともとの順序が維持されます。
したがって、位置的に先頭に近い値は先頭付近に、末尾に近い値は末尾付近に出現しやすくなります。
直感的にはほぼありえないような確率に思えるかもしれませんが、それなりに現実的に起こりえます。
誕生日攻撃の計算式 $n(0.5; H) \approx 1.1774 \sqrt{H}$ より、 54562 個程度の要素にソートを行うと約 50% の確率で衝突が発生します。

OrderBy は、メモリ (ヒープ) を $O(n)$ で消費します。これは、列挙する際に内部的にもとの配列のコピーを作成するためです。また、固定サイズですが IOrderedEnumerable<TElement> 自体や比較関数などのアロケーションも行われます。大きな配列になってくると地味に重くなってきます。

以上に述べたように、ソートによるシャッフルは実装が楽なかわりに欠点が多く存在します。

ただ、基数ソートといった $O(n)$ のアルゴリズムを利用する、あるいは巨大な配列に対しては並列化可能なソートを利用することで速くなるかもしれません。検討の余地がある……かも……？

手元で試した限りでは、 10 万要素の IEnumerable<int> に対して source.AsParallel().OrderBy(...) で並列ソートした場合に、非並列のものより 2 倍程度早くなりました。ただ、小さいコレクションに対しては 100 倍程度遅くなる・ヒープも数倍～数十倍消費する・並列実行のため擬似乱数まわりの扱いが難しくなるなど難点が多いため、おすすめはできません。

危険なランダムソート

これもソートを使ったシャッフルの実装ですが、絶対にしてはいけません。

static IEnumerable<T> LinqShuffle_Danger<T>(IEnumerable<T> source) {  
    return source.OrderBy(_ => _,   
        Comparer<T>.Create((_, _) => Random.Shared.Next(int.MinValue, int.MaxValue)));  
}

これはキーではなく比較関数が返す結果をランダムにする手法なのですが、これを行ってしまうとソートの前提となる関係性が破綻してしまいます。
シャッフルの移動先が偏ったり、運が悪いとランダムに以下のような ArgumentException が発生したりします。
「ランダムに」というのが厄介で、テスト時に成功して本番でこける、といった事故も起こりかねません。

System.ArgumentException: 'Unable to sort because the IComparer.Compare()   
method returns inconsistent results. Either a value does not compare   
equal to itself, or one value repeatedly compared to another value   
yields different results. IComparer: 'System.Comparison`1[System.Int32]'.'

Fisher-Yates shuffle

詳しい方は Fisher-Yates shuffle をご存知かと思います。
このアルゴリズムは非常に効率的です。

static void FisherYates<T>(Span<T> source) {  
    for (int i = source.Length - 1; i >= 1; i--) {  
        int j = Random.Shared.Next(i + 1);  
        (source[i], source[j]) = (source[j], source[i]);  
    }  
}

Fisher-Yates shuffle の計算量は $O(n)$ で、ランダムソートより効率的です。
メモリを追加で消費することもありません。 (LINQ 版と同じように元の配列を変更しない (inside-out な) 実装にする場合はもちろん元の配列と同じぶん消費しますが、どちらにせよ理想的です。)
加えて、 擬似乱数生成器が理想的であれば 均等にシャッフルされます。ある要素がある位置に配置される確率が全て等しくなります。

このアルゴリズムは .NET 8 Preview 1 で追加された Random.Shuffle でも利用されています。

擬似乱数まわりの最適化

さて、 Fisher-Yates shuffle は十分に効率的なうえ、シンプルです。これ自体を改良するのはかなり難しいでしょう。
ここで改良の余地があるのは、 Random.Shared.Next() 、つまり擬似乱数生成の部分です。

擬似乱数生成器の選定

擬似乱数生成器そのものの選定も重要になってきます。

完全に均等なシャッフルを目指すなら CSPRNG (暗号論的擬似乱数生成器; RandomNumberGenerator ) を使う手もあるかもしれませんが、その分パフォーマンスは犠牲になります。
実用的には、十分に内部状態の大きい (より厳密には均等分布次元の大きい) 擬似乱数生成器を使用すべきでしょう。

なお、お金が関わるような場合 (ガチャとか) やゲームの流れを大幅に左右する場合 (麻雀とか) の場合は、 CSPRNG を使用すべきところだと思います。

なぜ内部状態の大きい擬似乱数生成器が必要なのか、について簡単に説明すると、 $n$ 個の要素のシャッフルの結果は $n!$ 通りある以上、擬似乱数生成器側でも $n!$ 通りの乱数が生成できる必要があるためです。
具体例を挙げると、トランプ（ジョーカー 2 枚を含む、 54 枚）では $54! \lt 2^{238}$ であるので、少なくとも 238 bit 以上の乱数を生成できる必要が出てきます。麻雀牌 (花牌は除く、同種の牌を区別するものとする) なら $136! \lt 2^{773}$ なので 773 bit 以上必要です。
しかもこれは理想的な実装を行った場合の話で、通常はそれ以上のビット数が必要になります。

Next() を何回も呼び出せば 238 bit ぐらい余裕で生成できるじゃないか、と思われたかもしれませんが、擬似乱数生成器の内部実装によっては「出現しない組み合わせ」が生じる可能性があります。
具体例を挙げると、 64 bit の線形合同法では、 64 bit までなら任意の bit 列を出力できますが、それより大きい bit 列の場合はほぼ確実に出現しない組み合わせが生じます。
より具体的に、以下の線形合同法 Lcg64 を用いて 2 個の連続した出力を観測するとき、最初が 0 なら次は必ず 1442695040888963407 になります。それ以外のペア、例えば [0, 0] などは絶対に出力されません。

static ulong Lcg64(ref ulong state)  
    => state = state * 6364136223846793005 + 1442695040888963407;

したがって、 64 bit の線形合同法を用いてトランプをシャッフルしようとした場合、絶対に生成されない組み合わせや、出やすい組み合わせが出てきてしまいます。この場合は、最低限 xoshiro256++ (256 bit) 、余裕をもって xoroshiro1024++ (1024 bit) などを使用すべきでしょう。

それでいて、もちろん高速性も重要です。
例えば、メルセンヌツイスタ mt19937 であれば 19937 bit まで生成できるので大抵の用途のシャッフルに耐えます ( $2081! \lt 2^{19937}$ ; 理論上は 2081 枚のカードを均等にシャッフル可能) 。ただ、速度はモダンな擬似乱数生成器に比べると遅いです。

主要な (?) アルゴリズムの内部状態の bit 数と、それによってシャッフルできるカードの枚数上限を示します。

Algorithm	bits	cards
`LCG` (線形合同法)	64	20
`xoroshiro128+`	128	34
`shioi128`	128	34
`seiran128`	128	34
`xoshiro256**`	256	57
`culumi`	256	57
`xoroshiro1024*`	1024	171
`mt19937` (メルセンヌツイスタ)	19937	2081
`SCP-1214-EX`	4749265984	182651279

また余談です。シャッフルに限った話ではありませんが、擬似乱数生成器の初期化にも気を配る必要があります。例えば、メルセンヌツイスタの初期化関数には 32 bit のシードを受け取るものがありますが (オリジナル実装の init_genrand(unsigned long s)) 、これを利用してしまうと高々 $2^{32}$ 通りの系列 (シャッフル結果) しか得られなくなってしまいます。初期化時には内部状態以上の情報量を持ったソース (CSPRNG など) を用いて、全域をまんべんなくランダムにする必要があります。

それなら最初から CSPRNG を使えばよくない？という話もあります。難しいですね。
まぁ現実的には実行速度や取り回しのしやすさ、再現性の担保のために普通の PRNG を使うことになるでしょう。

その際のポイントとしては、できる限り擬似乱数生成器インスタンスを使いまわすこと (都度初期化しないこと) が挙げられます。シード値が擬似乱数生成器の内部状態より小さい場合はなおさら。
擬似乱数生成器は使い続けることを前提に設計されており、初期化直後 (特に小さなシード値によるもの) はランダムでない (何らかの相関があったり、立っているビット数が少なかったりする) 値を出力する場合があります。

乱数を絞りつくす

話を戻して、 Fisher-Yates shuffle のコードをもういちど見てみましょう。

static void FisherYates<T>(Span<T> source) {  
    for (int i = source.Length - 1; i >= 1; i--) {  
        int j = Random.Shared.Next(i + 1);  
        (source[i], source[j]) = (source[j], source[i]);  
    }  
}

これを見ると、 $n$ 個の要素に対して $n - 1$ 回の Next() 呼び出しがあることがわかります。
Next() 、つまり乱数生成は相対的に重い処理であるため、この部分の最適化を図りたいです。

64 bit 環境向けの擬似乱数生成器は、大抵の場合一度に 64 bit の乱数を生成できますので、 $2^{64}$ 通りの乱数を得ることができます。
例えば 100 要素のシャッフルなら、一回あたり高々 100 通りぶんの乱数しか必要ないのですから、 1 回で $2^{64}$ ぶんの情報量を持つ乱数を消費してしまうにはもったいないです。
相対的に重い Next() 呼び出しの回数を減らすため、できる限り乱数を絞りつくす必要があります。

絞りつくすといいことがもう一つあります。乱数を絞りつくす実装では、必要な均等分布次元数 (≒ 内部状態の bit 数) を減らすことができます。例えば 20 要素のシャッフルに対して 19 回 Next を呼ぶ素朴な実装では、 Next が 64 bit の乱数を出力する場合 $64 \times 19 = 1216$ bit ぶん必要であるのに対し、絞りつくす実装では 1 回の Next 呼び出しでよい ( $20! \lt 2^{64}$ ) ので $64$ bit で済みます。

乱数を絞りつくす実際の工程はこんな感じです。

$n = \prod_{k=2}^{20}$ を求める
64 bit 擬似乱数 $r$ を生成する
$r$ を $0 \le s \lt n$ の範囲に均等に変換できるように調整する
$r$ を分割してインデックス 2 ～ 20 ぶんの乱数を得る
得た乱数で Fisher-Yates shuffle を行う
$n = \prod_{k=21}^{33}$ を求め、以下繰り返し

それぞれの工程について詳しく見ていきましょう。

$n$ を求める

$n = \prod_{k=2}^{20} = 2 \times 3 \times \ldots \times 20 = 2432902008176640000$ を求めます。
これは何からきているかというと、 ulong で表現可能な ( $2^{64}$ より小さい) 範囲で最大の階乗の数です。

64 bit 擬似乱数 $r$ を生成する

ulong 全域に一様分布する乱数 r を生成します。
System.Random には残念ながら NextUInt64() は生えていないので、自前でお好みのアルゴリズムを実装した擬似乱数生成器があると良いでしょう。高速なものをチョイスすればより高速に、均等分布次元の高いものをチョイスすればより大きな配列のシャッフルに使えます。

一応 Random.NextBytes() で頑張れば不可能ではありませんが、オーバーヘッドがあるかもなので素直に自作することをおすすめします。

実は内部的には NextUInt64() は実装されているのですが、 internal なので触れません。残念。

$r$ を $0 \le s \lt n$ の範囲に均等に変換できるように調整する

ここは一般的な擬似乱数生成における範囲変換と同様で、 Next(int max) などと同じイメージです。
ただ注意すべきなのは「均等に変換」というところです。
例えば、安直に r % n で変換した場合は n が 2 の冪乗でない限り最小値付近が最大値付近より出やすくなります。

具体的には、 r % n で範囲変換した場合、 $\lbrack 0, 2^{64} \bmod{n})$ の範囲の数は $\lfloor 2^{64} / n \rfloor + 1$ 個、 $\lbrack 2^{64} \bmod{n}, n)$ の範囲の数は $\lfloor 2^{64} / n \rfloor$ 個出現します。
したがって、確率に偏りが出ます。

そのため、r を再生成する・別の乱数と組み合わせて補正をかけるなどして、均等に出るようにします。

今回は、 Swift で提案されている手法をベースにして調整を行います。
この手法は、もともとは一様分布の乱数を特定の範囲に偏りなく変換するための手法です。
Math.BigMul を巧みに使うことによって重い除算や剰余算をする必要をなくし、また Lemire 氏が提案している方式に比べて連続再試行となる確率が低いという特徴があるため、高速に実行することができます。

具体的には、Swift 式は $i$ 試行目で再試行になる確率は $n / 2^{64i}$ と指数関数的に低くなっていきます。
対して Lemire 式の場合は $(2^{64} \bmod n) / 2^{64}$ と $i$ に依存しない定数となります。最初の 1 回の確率は Swift 式に比べて低くなりますが、試行回数を重ねても一定です。

加えて Lemire 式の場合、 2 回目の乱数を振る前に剰余算 % を実行する必要があります。これは場合によっては 1 回の乱数生成に匹敵するレベルの時間がかかります。
したがって、今回の用途では Swift 式のほうが有利と判断しました。

Swift 式の実装例を以下に示します。

// factorial は本文中の n に対応; 範囲の上限 (この値を含まない)  
  
// 64 bit 乱数 r を生成  
ulong r = rng.Next();  
  
ulong rlo = r * factorial;  
  
// r * factorial の下位 64 bit (rlo) を見て、繰り上がりの可能性があれば…  
// (後続の処理で足される最大値が factorial - 1 なので、  
//  rlo <= (2^64) - factorial なら繰り上がりは発生せず、処置不要)  
while (rlo > 0ul - factorial)  
{  
    // 追加で乱数 t を生成し、繰り上がるかを調べる  
    // 下記の筆算をやるイメージ  
    //   [rhi] . [rlo]        -> r * factorial の 上位 rhi / 下位 rlo  
    // +     0 . [thi] [tlo]  -> t * factorial の 上位 thi / 下位 tlo  
    // ---------------------  
    //   [carry   sum] [tlo]  -> rhi + carry が求めるべきもの  
    ulong thi = Math.BigMul(rng.Next(), factorial, out ulong tlo);  
    ulong sum = rlo + thi;  
    ulong carry = sum < thi ? 1ul : 0ul;  
  
    // sum == 0xffff...ffff であれば、今後繰り上がりの可能性があるのでもう一度  
    // そうでなければこれ以上繰り上がりは発生しないので、 carry を足して終了  
    if (sum != ~0ul)  
    {  
        // r に carry(1) を足す → rlo が factorial 増える →   
        // while の条件式から必ずオーバーフローするので rhi が 1 増える  
        r += carry;  
        break;  
    }  
  
    rlo = tlo;  
}  
  
// rhi は偏りなく 0 <= x < factorial の範囲に分布する一様乱数  
ulong rhi = Math.BigMul(r, factorial, out _);

お分かりいただけたでしょうか？
私は最初このアルゴリズムを見たとき感動しました。よく思いつきますね……

Lemire 式で実装する場合はこのようになります。

// 64 bit 乱数 r を生成  
ulong r = rng.Next();  
  
ulong rlo = r * factorial;  
  
// 事前チェック。常に下式は成立するので、  
// (0 - factorial) % factorial < factorial  
// この if で弾ければ時間のかかる % をスキップできる  
if (rlo < factorial)  
{  
    // 2^64 % factorial == (2^64 - factorial) % factorial  
    ulong mod = (0 - factorial) % factorial;  
  
    // 0 <= rlo < mod の場合、再抽選  
    while (rlo < mod)  
    {  
        r = rng.Next();  
        rlo = r * factorial;  
    }  
}  
  
// rhi は偏りなく 0 <= x < factorial の範囲に分布する一様乱数  
ulong rhi = Math.BigMul(r, factorial, out _);

体感ですが、通常時の範囲変換はこちらのほうが速い場合が多いです。
Lemire 式のほうが乱数を複数生成する確率が低いので、特に乱数生成が重い場合に有利になりがちです。
使い分け (とベンチマーク) が大切ということかもしれません。

$r$ を分割してインデックス 2 ～ 20 ぶんの乱数を得る

$r$ が計算できたら、各インデックスを取り出します。

int t2 = (int)Math.BigMul(r, 2ul, out r);   // [0, 2)  
int t3 = (int)Math.BigMul(r, 3ul, out r);   // [0, 3)  
// ...  
int t20 = (int)Math.BigMul(r, 20ul, out r);   // [0, 20)

64 bit . 64 bit の固定小数点数をイメージするとわかりやすいかもしれません。
最初の r が 0.r 、つまり 0 ～ 1 の乱数と見立てて、 2, 3, ... を掛けたときの上位 64 bit = 整数部分を得ることで 0 以上 2, 3, ... 未満の乱数を取得します。

論文 "Batched Ranged Random Integer Generation" *1 では、可変進数のような考え方をしていました。 $1!$ の位 (0 ~ 1) から始まり、 $2!$ の位 (0 ～ 2)、 $3!$ の位 (0 ～ 3)、…… といった感じです。

Fisher-Yates shuffle を行う

ここはベースのコードとほぼ同じです。
違う点があるとすれば、オリジナルのコードでは i <= x < source.Length の範囲でランダムなインデックスを生成していましたが、こちらでは 0 <= x < i の範囲で生成しています。
for 文を i++; で回すことによって、 $n$ を事前に計算してキャッシュしておけるようにするためです。
こういうことをしても大丈夫か、と不安になるかもしれないので、数学的帰納法っぽく証明？しておきます。

まず、長さ 1 の配列は、各要素 (といっても要素 [0] だけです) が均等な確率 (1) で各位置 ([0]) に存在するので、均等にシャッフルされていると言えます。
次に、長さ $k$ の配列があり、それは均等にシャッフルされているとします。この配列に $k+1$ 番目の要素を追加したうえでランダムに $i$ ( $1 \le i \le k + 1$ ) 番目の要素と交換したとき、均等にシャッフルされていると言えるでしょうか。
まず、追加した $k + 1$ 番目の要素は等しい確率 ( $\frac{1}{k+1}$ ) で全ての場所に移動するため、均等であると言えます。その他の要素は移動していないか、 $\frac{1}{k+1}$ の確率で末尾と交換されたかなので、各位置に均等な確率で存在している状態を維持します。
以上から、このシャッフル方式でも問題なくシャッフルできるといえます。ふわっとしていますがこんな感じでいかがでしょうか……

次の $n$ を求めて、必要なぶん繰り返す

$n = \prod_{k=21}^{33} = 21 \times 22 \times \ldots \times 33 = 3569119343741952000$ を求めて、同様に操作を繰り返します。
これを元の配列の長さと同じ分までやります。
途中まで必要であれば (長さが 25 だった場合など) 、そこまでで乗算を打ち切ってしまってよいです。

結果

ベンチマーク結果を示します。
BatchedSwift が上記の「乱数を絞りつくした」実装です。
DataClass は record DataClass(double x, double y, double z, double w) のクラスです。クラスと構造体で性能特性が違う可能性を考慮してテストしています。

Method	array	Mean	Error	StdDev
LinqSort	DataClass[1024]	52,252.68 ns	1,007.564 ns	1,237.379 ns
FisherYatesSwift	DataClass[1024]	8,163.16 ns	63.845 ns	59.721 ns
BatchedSwift	DataClass[1024]	6,221.22 ns	101.979 ns	90.401 ns
LinqSort	DataClass[32]	969.79 ns	11.481 ns	10.739 ns
FisherYatesSwift	DataClass[32]	261.47 ns	2.584 ns	2.417 ns
BatchedSwift	DataClass[32]	134.73 ns	2.616 ns	2.569 ns
LinqSort	Int32[1024]	49,661.06 ns	727.360 ns	680.373 ns
FisherYatesSwift	Int32[1024]	3,425.73 ns	46.750 ns	43.730 ns
BatchedSwift	Int32[1024]	1,532.79 ns	18.482 ns	16.384 ns
LinqSort	Int32[32]	878.98 ns	10.204 ns	7.966 ns
FisherYatesSwift	Int32[32]	94.11 ns	1.894 ns	2.528 ns
BatchedSwift	Int32[32]	32.54 ns	0.227 ns	0.212 ns

Linq とは比べ物にならないレベルで Fisher-Yates 群が速いです。それはそう。
また、 Batched は生の Fisher-Yates に比べて 1.5 ～ 2 倍程度早くなっていることが分かります。

小手先の高速化

アルゴリズムレベルではない、小手先の高速化手法について書きます。
うまくいったやつとそうではないやつがあるので注意してください。

`Next()` メソッドの (手動) インライン展開

乱数生成をインライン展開することで高速化を図ります。

例えば、

for (int i = 0; i < length; i++)  
{  
    ulong r = rng.Next();  
    // do something  
}

これを、こういう感じにします。

var state = rng.State;  
for (int i = 0; i < length; i++)  
{  
    ulong r = StaticNext(state);  
    // do something  
}  
rng.State = state;  
  
// ---  
  
[MethodImpl(MethodImplOptions.AggressiveInlining)]  
static ulong StaticNext(State state) { /* same as Next()*/ }

rng.State の更新を最後に移動させ、 Next() を静的関数に実装しなおしているのがポイントです。

関数呼び出しをスキップできるようになるほか、手動で工夫して展開すると都度メモリに書かずにレジスタ上で完結するようになるため、多少の高速化が見込めます。

もちろん、擬似乱数生成器とべったり癒着することになるので、一長一短です。

上限を削る

乱数の再生成が必要になるのは r > 0 - n の場合でしたね。つまり、 n が小さいほど乱数を再生成する確率が下がります。
今まで $n = \prod_{k=2}^{20}$ として計算していましたが、これを $n = \prod_{k=2}^{18}$ のようにしたらどうでしょうか？

均等分布次元が減るのと引き換えに、棄却率を下げて再生成 (=遅延) を減らそう、という試みです。
ちょっと試した感じでは $n \lt 2^{58}$ の制約をつけたときに速度のバランスが良い、ということがわかっていますが、均等分布次元を犠牲にするほどの劇的な加速は得られていませんので、微妙です。

タプルでの交換をやめる

現代の C# では、以下のコードで要素のスワップができます。

(span[a], span[b]) = (span[b], span[a]);

しかし、これはどうしてか、以下の従来のコードのほうとアセンブリの生成結果が異なる場合があります。

var t = span[a];  
span[a] = span[b];  
span[b] = t;

体感としては、タプルを使わないコードのほうが簡潔なアセンブリを生成する傾向があります。

具体例は Sharplab を確認してみてください。

SIMD 化

C# では Vector128 などを経由して SIMD 化することができます。

このコードで SIMD 化できそうなところとしては、

$n = \prod_{i}^{k}$ の計算
各インデックスの計算

が挙げられます。
ただ、ちょっと試した限りではオーバーヘッドのほうが大きく、高速化にはつながりませんでした。

配列アクセス時の範囲チェック削除

// values[m] = something;  
Unsafe.Add(ref MemoryMarshal.GetReference(values), m) = something;

こう書くと IndexOutOfRangeException を飛ばすコードがなくなります。
このコードは高速性と危険性が表裏一体なので、十分なデバッグをしてから最後に実装してください。

$n$ のキャッシュ

$n$ の値は実行ごとに変わらないので、キャッシュしておいたり事前計算して埋め込んでおいたりすることもできます。

静的キャッシュ (事前計算して switch) や動的キャッシュ (Dictionary に登録) など試してみましたが、手間の割に高速化しませんでした。なので初手の $n = 20!$ だけ埋め込むのがよさそうです。

Fisher-Yates shuffle の誤った実装例

誤った例ですので、真似しないでください！

例えば、乱数生成の範囲指定で +1 し忘れた場合 ( $0 \le r \lt i$ ) 、以下のようになります。

static void FisherYates_Wrong_OffByOne<T>(Span<T> source) {  
    for (int i = source.Length - 1; i >= 1; i--) {  
        int j = Random.Shared.Next(i); // instead of i + 1  
        (source[i], source[j]) = (source[j], source[i]);  
    }  
}

この場合、「サットロのアルゴリズム」という変種になり、円順列を生成するようになります。
また、ある要素がシャッフル後に同じ位置に配置される確率が 0 になります。

また、乱数生成の範囲指定で常に配列全部の範囲を指定した場合も、偏りが生じてしまいます。

static void FisherYates_Wrong_Entire<T>(Span<T> source) {  
    for (int i = 0; i < source.Length; i++) {  
        int j = Random.Shared.Next(source.Length); // instead of i + 1  
        (source[i], source[j]) = (source[j], source[i]);  
    }  
}

初心者が何も見ずに実装するとこうなる場合が多い気がします。

交換によって生じるパターン数が $n^{n}$ になる一方で、シャッフルによって生じるパターン数は $n!$ です。 $n^{n} \not\equiv 0 \mod{n!}$ ですので、必ず偏りが生じます。

実際にどのように偏るのかについては、 Wikipedia が詳しいです。

https://www.sega-mj.com/arcade/viewer/haiyama/viewer.html

ところで、シャッフルの実例としてコードを探していたところ、これを見つけました。
この「サンプルコード」の実装ではまさしく上記の間違ったシャッフル法が実装されています。
そのうえ範囲変換が剰余で実装されているのでそこでも偏っています。二重苦。

MergeShuffle

さて、高速化にあたって思いつく手法のひとつとして、並列化が挙げられます。
並列にシャッフルを実行するアルゴリズムとして、 MergeShuffle があります。 *2

実装例はこんな感じです。分割統治法のような感じですね。
$2^{n}$ 個の領域に分割してそれぞれ Fisher-Yates でシャッフルし、それらをマージしていく感じです。

public static void Shuffle<TRng, TSpan>(TRng rng, Span<TSpan> span)  
    where TRng : IRandom  
{  
    Divide(rng, dist, span);  
}  
  
private static void Divide<TRng, TSpan>(TRng rng, Span<TSpan> span)  
    where TRng : IRandom  
{  
    if (span.Length <= 16)  
    {  
        FisherYates(rng, dist, span);  
    }  
    else  
    {  
        Divide(rng, dist, span[..(span.Length / 2)]);  
        Divide(rng, dist, span[(span.Length / 2)..]);  
        Merge(rng, dist, span);  
    }  
}  
  
private static void Merge<TRng, TSpan>(TRng rng, Span<TSpan> span)  
    where TRng : IRandom  
{  
    int start = 0;  
    int mid = span.Length / 2;  
    int end = span.Length - 1;  
  
    ulong r = rng.Next();  
    int entropy = 64;  
  
    while (true)  
    {  
        // エントロピーがなくなったら補充  
        if (entropy == 0)  
        {  
            r = rng.Next();  
            entropy = 64;  
        }  
  
        // 1 bit 取り出す  
        ulong bit = r & 1ul;  
        r >>= 1;  
        entropy--;  
  
        // bit 1 なら [start] と [end] を交換  
        if (bit == 0)  
        {  
            if (start == mid)  
            {  
                break;  
            }  
        }  
        else  
        {  
            if (mid == end)  
            {  
                break;  
            }  
            (span[start], span[end]) = (span[end], span[start]);  
            mid++;  
        }  
  
        start++;  
    }  
  
    while (start < end)  
    {  
        // [0, start) の乱数を生成、それと [start] を交換  
        int index = (int)rng.Next((ulong)start);  
        (span[start], span[index]) = (span[index], span[start]);  
        start++;  
    }  
}

ただ、実際に実装してみると遅いです。 Fisher-Yates に処理を足している感じなのでそれはそう。実装が悪いだけかもしれませんが。
乱数生成やシャッフルをうまく並列化できれば、大きな配列に対して効果が見込めそう……ではあります。

Feistel 構造を利用したシャッフル

面白い性質を持ったシャッフル手法のひとつとして、 Feistel 構造を利用したシャッフルが挙げられます。

Feistel 構造は、ブロック暗号の構成法の一種です。 DES などで使われています。
簡単な実装例は以下のようになります。

// internal state  
uint left = ..., right = ...;  
  
// 4 rounds (any number of rounds)  
for (int round = 0; round < 4; round++)  
{  
    (left, right) = (right, left ^ Round(right));  
}  
  
// here, left and right are encrypted  
  
uint Round(uint x) => /* returns any value */;

ここでポイントとなるのは、 Round() には任意の関数を用いることができることです。
速度と品質を天秤にかけて、 (いい意味で) 適当な関数を設定できます。

さて、ブロック暗号とシャッフルに何の関係があるのか、と思った方もいるかと思います。

暗号化できるということは、復号もできます。それはそう。
そして復号ができるということは、ある種の全単射関数のように振る舞うということです。
どういうことかというと、例えば 4 bit の Feistel 構造を構成して連番 [0, 1, 2, ..., 15] を入力したとき、それを暗号化した後の値は [0, 12, 8, ..., 7] みたいになるのですが、これは連番と一対一対応する、すなわち連番の順序を「シャッフル」したものと同じになります。
ということはつまり、シャッフルに使える、というわけです。

具体的な流れとしては、

$n$ 要素の並べ替えをしたいとき、 $n \lt 2^{2b}$ を満たす $2b$ ビットの Feistel 構造をつくる
でループ
1. $i$ を暗号化して $f(i)$ を求める
2. $f(i) \lt n$ なら、 $f(i)$ 番目の要素を返す (yield return)

という感じです。
「 $2b$ ビットの Feistel 構造」は、単に uint のペア (64 bit) にビットマスクを掛ければよいです。具体的には、 18 bit が必要なら 0x1ff (9 ビット) のマスクを掛ければそれが 2 個なので 18 bit になります。

このシャッフルの利点は、 $i$ 番目の要素がどこに移動したかを $O(1)$ で取得できる点です。 Fisher-Yates の場合は全要素の処理が終わるまで座標は確定しませんが、 Feistel 構造なら $i$ を暗号化するだけなので長さに依存せずに座標を取得できます。なので、超巨大な配列からいくつかの要素をランダムに抽出したい、といった用途については効率的に行えるかもしれません。
また、面白い性質としては、復号することでシャッフルを「元に戻す」ことができます。

対して、欠点としては、要素数が 2 冪でない場合の処理が結構めんどくさいことが挙げられます。要素数が 2 冪でない場合、範囲外参照になる場合があるので、それを読み飛ばす必要が出てきます。 yield return するような実装ならこれは簡単なのですが、インプレースな (追加領域を確保せず、元の配列をいじるような) 実装は難しいです。
また、 1 つの要素を取得するのにかかる時間が比較的長くなってしまう問題があります。ちゃんとした暗号化 (シャッフル) をするためには最低でも 2 ラウンド必要ですし、きちんとしたハッシュ関数を使う必要があります。それに対して、 Fisher-Yates であれば 1 つあたり 1 回の乱数生成、より最適化すれば 1 回の乗算と数回に 1 回の乱数生成だけで済んでしまいます。

おわりに

いろいろなシャッフル手法と、高速で均等なシャッフルを行うにあたっての工夫についてまとめました。
バニラの Fisher-Yates より速い手法がある、というのを初めて知った時は驚きました。

最後に、高速で均等な実行ができる手法の実装例を挙げておきます。

Fast shuffle by batching

*1:Brackett‐Rozinsky, Nevin, and Daniel Lemire. "Batched ranged random integer generation." Software: Practice and Experience (2024).

*2:Bacher, Axel, et al. "Mergeshuffle: a very fast, parallel random permutation algorithm." arXiv preprint arXiv:1508.03167 (2015).

2024-12-19

高速で頑健なシェーダー乱数の比較と提案

はじめに

シェーダー (HLSL) で乱数を発生させたくなることがしばしばあります。
直接的にはホワイトノイズ、間接的にはパーリンノイズなどが挙げられます。

そういうシェーダーを書いたことがある方は、 frac(sin(...)) みたいな擬似乱数生成器を目にしていることでしょう。
ですが、乱数を発生させる手法は星の数ほど存在します。品質や速度もそれぞれ異なる擬似乱数生成器がたくさんあります。

本稿ではそれらのシェーダー乱数の性能比較を行っていきたいと思います。

シェーダー乱数の定義

さて、乱数とはいっても、シェーダーでよく使われる乱数は CPU ベースの (メルセンヌツイスタとかの) 擬似乱数生成器とは設計レベルで異なることが多いです。

まず、状態を持たない、ある種のハッシュ関数のような実装をすることがほとんどです。
大抵の場合、座標を引数にとることになるでしょう。毎フレーム変化する乱数が欲しい場合には、次元の一つに時間を入れればよいです。

そして、出力の値域です。
CPU ベースのものなら uint32_t や uint64_t 型全部をカバーするような場合がほとんどです。
しかし、今回の用途はシェーダ用なので float 型だったり、最悪 $\lbrack 0, 255 \rbrack$ ぐらい取れれば OK 、みたいな場合があります。
例えば、改良パーリンノイズなら 12 通りの乱数が得られれば動きます。ホワイトノイズも (HDR でなければ) 256 通りあれば十分でしょう。
今回は、返り値は float の $\lbrack 0, 1 )$ で、最低でも 256 通り (ただし、できれば $2^{16}$ 通り以上) の値が得られることを要件とします。

以上から、シグネチャはこんな感じになりそうですね。

// returns [0, 1)  
float rand(float4 pos)  
{  
    ...  
}

このシグネチャに合わない関数は、概ね以下の方針で改造するものとします。

入力が足りない場合
- float のみ：繰り返し呼ぶ f(f(f(f(x) + y) + z) + w)
- float2 など： f(v.xy + v.zw)
出力が多い場合 (float2 とかを返す場合)
- ひとつだけ使う： f(v).x
- 束ねる (総和をとる、 xor するなど)： dot(f(v), 1)

また、入力 pos は整数座標が入っていると仮定してもよいものとします。
つまり、 $(0.0, 0.0)$ と $(0.4, 0.9)$ が同じ値を返すことを許容します。
これはどうしてかというと、第一に uint ベースのアルゴリズムが多いため、第二にパーリンノイズなどの勾配ベクトルを生成する用途では整数格子なので小数以下はそもそも見ていないためです。
ついでに ±∞ や NaN 、非正規化数などへの対処は考えなくてよいものとします。

asint() 関数で直接 float から int へ変換することもできますが、基本的には普通に int (uint) へのキャストをすることにします。

乱数の品質については、擬似乱数生成器の品質を測れるテストスイートである PractRand でテストすることにします。
詳しくは後述します。

なお、本稿での float は CPU での float と同じ、 IEEE 754 規格の 32bit 単精度浮動小数点数とします。

half や fixed は精度が異なるためアルゴリズムが破綻しますので、考えないものとします。
特にモバイル向けのシェーダを書く場合、乱数関係の実装においては必ず float 精度またはそれに相当する指定を行ってください。

以上をまとめると、シェーダー乱数の要件は、

状態を持たない
シグネチャは float rand(float4 pos)
pos は整数としてよい
戻り値は $\lbrack 0, 1 )$ 、最低 256 通り以上

です。

PractRand テストについて

PractRand は、 2024 年現在主力となる擬似乱数生成器のテストスイートです。
diehard(er) や TestU01 といった他のテストスイートに比べると格段に検出力が高いほか、比較的早く問題を発見することができます。

具体的には、 TestU01 の BigCrush だとどんなに速くとも 3 時間程度かかるうえに偽陽性が出る (真の乱数であっても検定に失敗する) ことがあるのですが、 PractRand では最速で 1 秒程度で問題検出されるうえ、偽陽性も今のところ確認されていません。

PractRand のテストに落ちるということは、「少なくとも (バイト長) 出力したときに、なんらかの統計的・構造的問題が見られ、真の乱数ではないと見破れる」ことを示します。テストに落ちても乱数として「全く」役に立たないということではないですが、客観的な品質の指標として役立ちます。

今回は 64 KiB ( $2^{16}$ バイト) をテストの開始地点とし、 1 TiB ( $2^{40}$ バイト) までテストに通れば高品質だと定義することにします。
本稿では PractRand ver. 0.94 を使用しました。
ちなみに体感としては、 64 KiB で即座に失敗する擬似乱数生成器は人間でも目に見えるレベルで問題があります。

とても雑にまとめると、テストに落ちるまでのバイト長が長ければ乱数として強いということです。

本来の擬似乱数生成器のテストでは 32 TiB 程度まで見るほか、検定を expand モードにして実施したりしますが、時間がかかるので今回はそこまでしません。
参考までに、 1 TiB のテストには実測 1 ～ 2 時間かかりますので、 32 TiB のテストには単純計算で 32 ～ 64 時間 (1 日半～ 2 日半) かかります。

テスト用コードはこういう感じです。
0, -1, +1, -2, +2, -3, +3, ... という感じに進めていき、 65536 になったら 0 に戻して次の次元を加算に行きます。
また、コードに書いたように float の返り値の上位 16 bit だけをチェックするものとします。実用上、 float の下位ビットに何かあってもほぼ問題にならないと考えたためです。

class DummyRNG : public PractRand::RNGs::vRNG16 {  
public:  
    int x, y, z, w;  
  
    void increment() {  
        x = x >= 0 ? ~x : -x;  
        if (x >= 65536) {  
            x = 0;  
        }  
        else return;  
  
        y = y >= 0 ? ~y : -y;  
        if (y >= 65536) {  
            y = 0;  
        }  
        else return;  
  
        z = z >= 0 ? ~z : -z;  
        if (z >= 65536) {  
            z = 0;  
        }  
        else return;  
  
        w = w >= 0 ? ~w : -w;  
        if (w >= 65536) {  
            w = 0;  
        }  
    }  
  
    Uint16 raw16() {  
        increment();  
        // call some shader hashes  
        return (Uint16)(shader(x, y, z, w) * 65536);  
    }  
  
    void walk_state(PractRand::StateWalkingObject* walker) {  
        x = y = z = w = 0;  
    }  
  
    std::string get_name() const { return "shaderimpl"; }  
};

ちなみに、テストにかけるために全部 C++ に移植する羽目になり泣いていました。

シェーダー乱数一覧

シェーダー乱数を集めるにあたって、網羅的に調査されている素敵な論文 "Hash functions for gpu rendering" を見つけました。 *1
基本的にはこの論文からチョイスして、後は私が見つけたものについて調査します。

この論文の "Detailed hash results and code" に再現用のコードが載っているのですが、それぞれの出典元と比較したところ異なる実装が行われている場合がいくつか見受けられました。
本稿ではなるべく出典元の実装を尊重するため、一部で当該論文とは異なる実装を行っている場合があります。

なお、シグネチャが float rand(float4 v) ではないものに関しては、できるだけオリジナルに近い実装を rand_orig() 関数として実装した後、変換コードを rand() に書く、という運用をしました。

各乱数には、以下の表をつけています。

key	value
Instruction Mix	ピクセルシェーダーの命令数
GPU Duration	描画にかかった時間
FPS	実測 Frames per Second
PractRand Failed	PractRand で検定失敗するバイト長

Instruction Mix は、ピクセルシェーダーの命令数です。 NVIDIA Nsight Graphics で測定しました。

GPU Duration は、 1024x1024 のテクスチャに描画したときにかかった時間です。これも NSight Graphics で測定しました。

FPS は、 1024x1024 の quad (UI) を用意して描画したときの FPS です。しばらく動作させて最大の数値を記録しました。

PractRand Failed は、 PractRand でテストに失敗したときのバイト長 (2冪) です。
大きければ大きいほど品質が良いです。テストに失敗しなかった場合は > 2^41 (→ 2 TiB の検定にパスした) のように最低限確認できたところまで記載しています。

速度の参考までに、白塗りするだけのシェーダーはこんな感じでした。

key	value
Instruction Mix	0
GPU Duration	25.60 μs
FPS	2676
PractRand Failed	×

改良パーリンノイズの定数テーブル

static int perlin_permutation[512] = {  
    151,160,137,91,90,15,  
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,  
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,  
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,  
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,  
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,  
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,  
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,  
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,  
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,  
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,  
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,  
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,  
  
    151,160,137,91,90,15,  
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,  
    190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,  
    88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,  
    77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,  
    102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,  
    135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,  
    5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,  
    223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,  
    129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,  
    251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,  
    49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,  
    138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180  
};  
  
float perlinperm(float4 pos)  
{  
    return   
        perlin_permutation[  
            perlin_permutation[  
                perlin_permutation[  
                    perlin_permutation[int(pos.x) & 0xff]  
                    + (int(pos.y) & 0xff)]  
                + (int(pos.z) & 0xff)]  
            + (int(pos.w) & 0xff)] * (1.0 / 256.0);  
}

これは、改良パーリンノイズの勾配ベクトル生成に使われている擬似乱数生成器です。 *2
事前に定数テーブルを生成しておいて、それをうまく参照することで乱数を生成します。

定数テーブルは、 $\lbrack 0, 255 \rbrack$ の連番を適当にシャッフルしたものを 2 つつなげて長さ [512] にすることで生成します。
2 つつなげることで、都度 & をとる必要がなくなっています。

利点としては、定数テーブルをシャッフルしなおすことでシード値の設定が可能なこと、テーブルさえコピペしてしまえば実装が簡単なことが挙げられます。
欠点としては、座標 $\lbrack 0, 255 \rbrack$ までしか対応していないこと (それ以降はループします) 、それなりのサイズの定数テーブルが必要なため、定数ロードに時間がかかる可能性があることが挙げられます。

微妙にパターンが見える気がしますね。

key	value
Instruction Mix	21
GPU Duration	128.0 μs
FPS	2501
PractRand Failed	2¹⁶

案の定テストにも即座に落ちています。

mod289

float mod289_permute(float i)  
{  
    float im = mod(i, 289.0);  
    return mod(((im * 34.0) + 10.0) * im, 289.0);  
}  
  
float mod289(float4 i)  
{  
    return mod289_permute(  
        mod289_permute(  
            mod289_permute(  
                mod289_permute(i.w)  
                + i.z)  
            + i.y)  
        + i.x) * (1.0 / 289);  
}

これは、タイル化可能シンプレックスノイズ (psrdnoise) の論文で使用されている、勾配ベクトル生成用の擬似乱数生成器です。 *3

$x_{n+1} = 34 x_n^{2} + 10 x_n \mod{289}$ の遷移式を使って生成します。

これはさざ波のようなパターンが見える気がしますね。

key	value
Instruction Mix	39
GPU Duration	27.65 μs
FPS	2656
PractRand Failed	2¹⁶

案の定テストにも落ちています。

FNV1

uint fnv1_orig(float x, float y, float z, float w)  
{  
    const uint prime = 16777619;  
      
    uint ret = 2166136261;   
    uint key = x;  
  
    ret *= prime;  
    ret ^= ((key >> 0) & 0xff);               
    ret *= prime;  
    ret ^= ((key >> 8) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 16) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 24) & 0xff);  
  
    key = y;  
  
    ret *= prime;  
    ret ^= ((key >> 0) & 0xff);               
    ret *= prime;  
    ret ^= ((key >> 8) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 16) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 24) & 0xff);  
  
    key = z;  
  
    ret *= prime;  
    ret ^= ((key >> 0) & 0xff);               
    ret *= prime;  
    ret ^= ((key >> 8) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 16) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 24) & 0xff);  
  
    key = w;  
  
    ret *= prime;  
    ret ^= ((key >> 0) & 0xff);               
    ret *= prime;  
    ret ^= ((key >> 8) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 16) & 0xff);  
    ret *= prime;  
    ret ^= ((key >> 24) & 0xff);  
  
    return ret;  
}  
  
float fnv1(float4 v)  
{  
    return fnv1_orig(v.x, v.y, v.z, v.w) * 2.3283064365386962890625e-10;  
}

GPU でのハッシュ関数の実装について検討された論文に載っていました。 *4
FNV ハッシュについては Wikipedia を参照してください。

ちなみに、最後に掛けられている定数 2.3283064365386962890625e-10 は $2^{-32}$ (1.0 / 4294967296) です。

これも小さなさざ波のようなパターンが見える気がします。

key	value
Instruction Mix	50
GPU Duration	30.72
FPS	2656
PractRand Failed	2¹⁶

案の定テストに落ちています。

JenkinsHash

uint jenkins_orig(float4 v)  
{  
    uint ret = 0;  
  
    uint key = v.x;  
  
    ret += ((key >> 0) & 0xff);               
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 8) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 16) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 24) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
      
    key = v.y;  
  
    ret += ((key >> 0) & 0xff);               
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 8) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 16) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 24) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
  
    key = v.z;  
  
    ret += ((key >> 0) & 0xff);               
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 8) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 16) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 24) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
  
    key = v.w;  
  
    ret += ((key >> 0) & 0xff);               
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 8) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 16) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
    ret += ((key >> 24) & 0xff);  
    ret += ret << 10;  
    ret ^= ret >> 6;  
  
  
    ret += ret << 3;  
    ret ^= ret >> 11;  
    ret += ret << 15;  
  
    return ret;  
}  
  
float jenkins(float4 v)  
{  
    return jenkins_orig(v) * 2.3283064365386962890625e-10;  
}

1997 年に Bob Jenkins 氏が開発したハッシュ関数です。 *5
上記の文献での one_at_a_time 関数がそれです。

見た目上は綺麗なホワイトノイズになっている気がしますね。

key	value
Instruction Mix	93
GPU Duration	27.65 μs
FPS	2721
PractRand Failed	2²¹

しかし、テストには失敗しています。
こういう目視では確認できないような品質がチェックできるのが PractRand の強みです。

Blum Blum Shub

uint bbs65521_orig(uint value)  
{  
    value %= 65521;  
    value = (value * value) % 65521;  
    value = (value * value) % 65521;  
    return value;  
}  
  
float bbs65521(float4 v)  
{  
    return bbs65521_orig(bbs65521_orig(bbs65521_orig(bbs65521_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w)) * (1.0 / 65521);  
}

Blum-Blum-Shub は、 1986 年に発表された暗号論的擬似乱数生成器です。 *6

$x_{n+1} = x_n^{2} \mod{M}$ の漸化式で更新されます。
今回は $M = 65521$ として、品質向上のために 2 回処理してから返しています。

本来は $M$ に大きな素数 $p, q$ を用いて $M = p q$ となるように構成するのですが、今回は小型版のようなものなので「暗号論的」な強度はありません。

見た目上の違和感はありませんね。

key	value
Instruction Mix	53
GPU Duration	27.65 μs
FPS	2735
PractRand Failed	2¹⁶

でもテストには落ちています。

加えて、定数を 4093 に変えてやってみたバージョンを以下に示します。

なぜ 4093 なのかというと、 float 型の精度が 24 bit 分しかないためです。
2 つの数を乗算したときに 24 bit をオーバーしないためには、法 (mod の部分) がその平方根 ( $\sqrt{2^{24}} = 2^{12}$ ) 以下になるようにしなければなりません。
その条件下で最も大きい素数が $4093$ だからです。

float bbs4093_orig(float seed)  
{  
    float prime24 = 4093;  
    float s = frac(seed / prime24);  
    s = frac(s * s * prime24);  
    s = frac(s * s * prime24);  
    return s;  
}  
  
float bbs4093(float4 v)  
{  
    return bbs4093_orig(v.x + bbs4093_orig(v.y + bbs4093_orig(v.z + bbs4093_orig(v.w))));  
}

これも特にパターンは見受けられない気がします。
ただ、動画として見ると左側に線のようなものが見受けられました。

key	value
Instruction Mix	49
GPU Duration	27.65 μs
FPS	2667
PractRand Failed	2¹⁶

これもテストにはあっという間に落ちています。

CityHash

uint mur(uint a, uint h)  
{  
    a *= 0xcc9e2d51;  
    a = (a >> 17) | (a << 15);  
    a *= 0x1b873593;  
    h ^= a;  
    h = (h >> 19) | (h << 13);  
    return h * 5 + 0xe6546b64;  
}  
  
uint fmix(uint h)  
{  
    h ^= h >> 16;  
    h *= 0x85ebca6b;  
    h ^= h >> 13;  
    h *= 0xc2b2ae35;  
    h ^= h >> 16;  
    return h;  
}  
  
uint city_orig(uint4 s)  
{  
    uint len = 16;  
  
    uint a = s.y;  
    uint b = s.y;  
    uint c = s.z;  
    uint d = s.z;  
    uint e = s.x;  
    uint f = s.w;  
    uint h = len;  
  
    return fmix(mur(f, mur(e, mur(d, mur(c, mur(b, mur(a, h)))))));  
}  
  
float city(float4 s)  
{  
    return city_orig(uint4(s)) * 2.3283064365386962890625e-10;  
}

CityHash は、 2011 年に google が開発した非暗号論的ハッシュ関数です。 *7

今回は、その中でも 32 bit のハッシュ値を返す CityHash32 を使います。

パターンは特に見受けられませんね。

key	value
Instruction Mix	49
GPU Duration	27.65 μs
FPS	2695
PractRand Failed	2⁴¹

ここにきてようやく TiB の壁 ( $2^{40}$ ) を突破した擬似乱数生成器 (ハッシュ?) が現れました。
さすが本業のハッシュ関数、品質が違いますね。

ESGTSA

uint esgtsa_orig(uint s)  
{  
    s = (s ^ 2747636419) * 2654435769;  
    s = (s ^ s >> 16) * 2654435769;  
    s = (s ^ s >> 16) * 2654435769;  
    return s;  
}  
  
float esgtsa(float4 v)  
{  
    return (esgtsa_orig(esgtsa_orig(esgtsa_orig(esgtsa_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w))) * 2.3283064365386962890625e-10;  
}

ESGTSA は、 "Evolving sub-grid turbulence for smoke animation" という論文で提案されているアルゴリズムです。 *8
論文タイトルの頭文字をとっているわけですね。
もともとは、自然な煙のアニメーションを生成するための論文です。

keijiro 大先生によれば、 Unity の HDRP でもこの実装が使われているそうです。

これもパターンは見受けられませんね。

key	value
Instruction Mix	38
GPU Duration	26.62 μs
FPS	2721
PractRand Failed	2⁴⁰

1 TiB のテストでちょうど失敗しましたが、概ね高品質といえるでしょう。シンプルめの実装なので意外でした。

Fast

float fast_orig(float2 v)  
{  
    v = (1.0 / 4320.0) * v + float2(0.25, 0.0);  
    float state = frac(dot(v * v, float2(3571, 3571)));  
    return frac(state * state * (3571.0 * 2.0));  
}  
  
float fast(float4 v)  
{  
    return fast_orig(v.xy + v.zw);  
}

Unreal Engine が典拠らしいです。
私にはアクセス権がないので確認できませんでした……

float 型だけで (uint 型を経由せずに) 計算できるのがポイントでしょうか。

レンズの干渉模様みたいなものが見えますね。

key	value
Instruction Mix	16
GPU Duration	27.65 μs
FPS	2664
PractRand Failed	2¹⁶

テストには即座に落ちてしまいました。

Hash without Sine

float hashwosine(float4 p)  
{  
    p = frac(p * float4(0.1031, 0.1030, 0.0973, 0.1099));  
    p += dot(p, p.wzxy + 33.33);  
    return frac((p.x + p.y) * (p.z + p.w));  
}

この関数は、 2014 年に Dave_Hoskins 氏が Shadertoy 上で発表したものです。 *9

sin 関数が環境依存で値が変わる問題 (後述) に対処すべく、 sin なしで計算できることを念頭に設計されているようです。

微妙に縞模様が見えるような、そうでもないような……

key	value
Instruction Mix	31
GPU Duration	28.67 μs
FPS	2702
PractRand Failed	2¹⁶

テストには即座に落ちてしまいました。

license

// Hash without Sine
// MIT License...
/ Copyright (c)2014 David Hoskins.
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE./

HybridTaus

uint taus(uint z, int s1, int s2, int s3, uint m)  
{  
    uint b = ((z << s1) ^ z) >> s2;  
    return ((z & m) << s3) ^ b;  
}  
  
uint hybridtaus_orig(uint4 z)  
{  
    z.x = taus(z.x, 13, 19, 12, 0xfffffffe);  
    z.y = taus(z.y,  2, 25,  4, 0xfffffff8);  
    z.z = taus(z.z,  3, 11, 17, 0xfffffff0);  
    z.w = z.w * 1664525 + 1013904223;  
  
    return z.x ^ z.y ^ z.z ^ z.w;  
}  
  
float hybridtaus(float4 s)  
{  
    return hybridtaus_orig(uint4(s)) * 2.3283064365386962890625e-10;  
}

HybridTaus は、 2007 年に Howes らによって開発されました。 *10

Hybrid とあるように、 Tausworthe Generator (LFSR; メルセンヌツイスタなどの祖先) と LCG (線形合同法) のハイブリッドとなっています。

何も映っていません。
というのも、 x, y, z 成分が小さい場合 (といっても 65536 以下とかそういう場合でも) 結果への寄与がほぼ 0 になるためです。
(最終的に float に変換する際は上位ビットが使われるが、 x, y, z がほぼ下位ビットにしか伝播しないためです。)

key	value
Instruction Mix	25
GPU Duration	27.65 μs
FPS	2649
PractRand Failed	×

論外なのでテストしていません。
あくまでも「内部状態を更新していく擬似乱数生成器」としての運用を主眼としているものなので、ハッシュ用途には難しかったようです。

Interleaved Gradient Noise

float ign_orig(float2 v)  
{  
    float3 magic = float3(0.06711056, 0.00583715, 52.9829189);  
    return frac(magic.z * frac(dot(v, magic.xy)));  
}  
  
float ign(float4 v)  
{  
    return ign_orig(v.xy + v.zw);  
}

Interleaved Gradient Noise は、 2014 年に Jorge Jimenez 氏によって発表されたノイズです。 *11

ディザと乱数の中間のような性質を持っており、 Temporal Anti Aliasing (TAA) などに用いるとよい結果が得られるそうです。

乱数というよりディザっぽいですね。それはそう。

key	value
Instruction Mix	10
GPU Duration	26.62 μs
FPS	2632
PractRand Failed	2¹⁶

もちろんテストには落ちています。あくまでもディザ用であって乱数ではないということですかね。

IQ's Integer Hash 1

uint iqint1_orig(uint n)  
{  
    n ^= n << 13;  
    return n * (n * n * 15731 + 789221) + 1376312589;  
}  
  
float iqint1(float4 pos)  
{  
    return iqint1_orig(uint(pos.x) + iqint1_orig(uint(pos.y) + iqint1_orig(uint(pos.z) + iqint1_orig(uint(pos.w))))) * 2.3283064365386962890625e-10;  
}

iq 氏が 2017 年に Shadertoy 上で公開したハッシュ関数です。 *12

ちなみに、オリジナルのソースコードには以下のように記されています。

Do NOT use this hash as a random number generator.

特にパターンは見られず、概ね良好ですね。

key	value
Instruction Mix	30
GPU Duration	28.67 μs
FPS	2630
PractRand Failed	2¹⁷

しかしテストには落ちてしまっています。

license

// The MIT License
// Copyright © 2017 Inigo Quilez
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

IQ's Integer Hash 2

uint3 iqint2_orig(uint3 x)  
{  
    uint k = 1103515245;  
  
    x = ((x >> 8) ^ x.yzx) * k;  
    x = ((x >> 8) ^ x.yzx) * k;  
    x = ((x >> 8) ^ x.yzx) * k;  
      
    return x;  
}  
  
float iqint2(float4 pos)  
{  
    return (dot(iqint2_orig(pos.xyz), 1) + dot(iqint2_orig(pos.w), 1)) * 2.3283064365386962890625e-10;  
}

同じく、 iq 氏が 2017 年に Shadertoy 上で公開したハッシュ関数です。 *13

こちらも特にパターンは見られませんね。優秀そうです。

key	value
Instruction Mix	42
GPU Duration	26.62 μs
FPS	2622
PractRand Failed	2⁴²

実際にテストにも通っています。

license

IQ's Integer Hash 3

uint iqint3_orig(uint2 x)  
{  
    uint2 q = 1103515245 * ((x >> 1) ^ x.yx);  
    uint n = 1103515245 * (q.x ^ (q.y >> 3));  
  
    return n;  
}  
  
float iqint3(float4 pos)  
{  
    uint value = iqint3_orig(pos.xy) + iqint3_orig(pos.zw);  
    return value * 2.3283064365386962890625e-10;  
}

これも同じく、 iq 氏が 2017 年に Shadertoy 上で公開したハッシュ関数です。 *14

これもぱっと見はパターンがなく、よさそうですね。

key	value
Instruction Mix	24
GPU Duration	27.65 μs
FPS	2580
PractRand Failed	2¹⁶

しかし、テストには即座に落ちてしまいました。違いがわからない……

ただし、 2024 年に更新されたらしく、現在は以下のアルゴリズムに置き換わっています。
従来のアルゴリズムは https://www.shadertoy.com/view/XlGcRh から参照することができます。

uint iqint32_orig(uint2 p)  
{  
    p *= uint2(73333, 7777);  
    p ^= uint2(3333777777, 3333777777) >> (p >> 28);  
    uint n = p.x * p.y;  
    return n ^ n >> 15;  
}  
  
float iqint32(float4 pos)  
{  
    uint value = iqint32_orig(pos.xy) + iqint32_orig(pos.zw);  
    return value * 2.3283064365386962890625e-10;  
}

これもぱっと見は大丈夫そうですが……

key	value
Instruction Mix	34
GPU Duration	27.65 μs
FPS	2728
PractRand Failed	2¹⁸

これもテストには落ちてしまいました。改良前に比べれば長生きしているので、改良はされているようですね。

license

// The MIT License
// Copyright © 2017,2024 Inigo Quilez
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

JKISS32

uint jkiss32_orig(uint2 p)  
{  
    uint x = p.x, y = p.y;  
    uint z = 345678912, w = 456789123, c = 0;  
    int t;  
  
    y ^= y << 5;  
    y ^= y >> 7;  
    y ^= y << 22;  
  
    t = int(z + w + c);  
    z = w;  
    c = uint(t < 0);  
    w = uint(t & 2147483647);  
    x += 1411392427;  
  
    return x + y + w;  
}  
  
float jkiss32(float4 p)   
{  
    return (jkiss32_orig(p.xy) + jkiss32_orig(p.zw)) * 2.3283064365386962890625e-10;  
}

JKISS32 は、 2010 年に David Jones 氏が公開した乱数のベストプラクティス集に載っているアルゴリズムです。 *15

もともと George Marsaglia 大先生の KISS 擬似乱数生成器 *16 があって、それを改良したものだそうです。

乗算を使わずに計算できる利点があるとのことです。

なお、もともとは x, y, z, w, c が内部状態 (state) となっていたものをシェーダー用に改造した模様です。

見るからにダメです。
x がほぼ素通しになる設計上、横縞になってしまうようです。

key	value
Instruction Mix	15
GPU Duration	26.62 μs
FPS	2687
PractRand Failed	×

論外なのでテストはしていません。

LCG

uint lcg_orig(uint p)  
{  
    return p * 1664525 + 1013904223;  
}  
  
float lcg(float4 v)  
{  
    return (lcg_orig(lcg_orig(lcg_orig(lcg_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w))) * 2.3283064365386962890625e-10;  
}

言わずと知れた線形合同法です。この定数は Numerical Recipes によるものです。 *17

はっきりとパターンが見えますね。

key	value
Instruction Mix	14
GPU Duration	27.65 μs
FPS	2695
PractRand Failed	2¹⁶

当然テストにも即座に落ちています。

MD5

uint F(uint3 v)  
{  
    return (v.x & v.y) | (~v.x & v.z);  
}  
  
uint G(uint3 v)  
{  
    return (v.x & v.z) | (v.y & ~v.z);  
}  
  
uint H(uint3 v)  
{  
    return v.x ^ v.y ^ v.z;  
}  
  
uint I(uint3 v)  
{  
    return v.y ^ (v.x | ~v.z);  
}  
  
void FF(inout uint4 v, inout uint4 rotate, uint x, uint ac)  
{  
    v.x = v.y + rotl(v.x + F(v.yzw) + x + ac, rotate.x);  
    rotate = rotate.yzwx;  
    v = v.yzwx;  
}  
  
void GG(inout uint4 v, inout uint4 rotate, uint x, uint ac)  
{  
    v.x = v.y + rotl(v.x + G(v.yzw) + x + ac, rotate.x);  
    rotate = rotate.yzwx;  
    v = v.yzwx;  
}  
  
void HH(inout uint4 v, inout uint4 rotate, uint x, uint ac)  
{  
    v.x = v.y + rotl(v.x + H(v.yzw) + x + ac, rotate.x);  
    rotate = rotate.yzwx;  
    v = v.yzwx;  
}  
  
void II(inout uint4 v, inout uint4 rotate, uint x, uint ac)  
{  
    v.x = v.y + rotl(v.x + I(v.yzw) + x + ac, rotate.x);  
    rotate = rotate.yzwx;  
    v = v.yzwx;  
}  
  
uint K(uint i)  
{  
    return uint(abs(sin(float(i) + 1.0)) * float(0xffffffff));  
}  
  
uint4 md5_orig(uint4 u)  
{  
    uint4 digest = uint4(0x67452301, 0xefcdab89, 0x98badcfe, 0x10325476);  
    uint4 r, v = digest;  
    uint i = 0;  
  
    uint m[16];  
  
    m[0] = u.x;  
    m[1] = u.y;  
    m[2] = u.z;  
    m[3] = u.w;  
    m[4] = 0;  
    m[5] = 0;  
    m[6] = 0;  
    m[7] = 0;  
    m[8] = 0;  
    m[9] = 0;  
    m[10] = 0;  
    m[11] = 0;  
    m[12] = 0;  
    m[13] = 0;  
    m[14] = 0;  
    m[15] = 0;  
  
    r = uint4(7, 12, 17, 22);  
  
    FF(v, r, m[0], K(i++));  
    FF(v, r, m[1], K(i++));  
    FF(v, r, m[2], K(i++));  
    FF(v, r, m[3], K(i++));  
    FF(v, r, m[4], K(i++));  
    FF(v, r, m[5], K(i++));  
    FF(v, r, m[6], K(i++));  
    FF(v, r, m[7], K(i++));  
    FF(v, r, m[8], K(i++));  
    FF(v, r, m[9], K(i++));  
    FF(v, r, m[10], K(i++));  
    FF(v, r, m[11], K(i++));  
    FF(v, r, m[12], K(i++));  
    FF(v, r, m[13], K(i++));  
    FF(v, r, m[14], K(i++));  
    FF(v, r, m[15], K(i++));  
  
    r = uint4(5, 9, 14, 20);  
    GG(v, r, m[1], K(i++));  
    GG(v, r, m[6], K(i++));  
    GG(v, r, m[11], K(i++));  
    GG(v, r, m[0], K(i++));  
    GG(v, r, m[5], K(i++));  
    GG(v, r, m[10], K(i++));  
    GG(v, r, m[15], K(i++));  
    GG(v, r, m[4], K(i++));  
    GG(v, r, m[9], K(i++));  
    GG(v, r, m[14], K(i++));  
    GG(v, r, m[3], K(i++));  
    GG(v, r, m[8], K(i++));  
    GG(v, r, m[13], K(i++));  
    GG(v, r, m[2], K(i++));  
    GG(v, r, m[7], K(i++));  
    GG(v, r, m[12], K(i++));  
  
    r = uint4(4, 11, 16, 23);  
    HH(v, r, m[5], K(i++));  
    HH(v, r, m[8], K(i++));  
    HH(v, r, m[11], K(i++));  
    HH(v, r, m[14], K(i++));  
    HH(v, r, m[1], K(i++));  
    HH(v, r, m[4], K(i++));  
    HH(v, r, m[7], K(i++));  
    HH(v, r, m[10], K(i++));  
    HH(v, r, m[13], K(i++));  
    HH(v, r, m[0], K(i++));  
    HH(v, r, m[3], K(i++));  
    HH(v, r, m[6], K(i++));  
    HH(v, r, m[9], K(i++));  
    HH(v, r, m[12], K(i++));  
    HH(v, r, m[15], K(i++));  
    HH(v, r, m[2], K(i++));  
  
    r = uint4(6, 10, 15, 21);  
    II(v, r, m[0], K(i++));  
    II(v, r, m[7], K(i++));  
    II(v, r, m[14], K(i++));  
    II(v, r, m[5], K(i++));  
    II(v, r, m[12], K(i++));  
    II(v, r, m[3], K(i++));  
    II(v, r, m[10], K(i++));  
    II(v, r, m[1], K(i++));  
    II(v, r, m[8], K(i++));  
    II(v, r, m[15], K(i++));  
    II(v, r, m[6], K(i++));  
    II(v, r, m[13], K(i++));  
    II(v, r, m[4], K(i++));  
    II(v, r, m[11], K(i++));  
    II(v, r, m[2], K(i++));  
    II(v, r, m[9], K(i++));  
  
    return digest + v;  
}  
  
float md5(float4 v)  
{  
    uint4 result = md5_orig(v);  
    return dot(result, 1) * 2.3283064365386962890625e-10;  
}

暗号論的ハッシュ関数である (と現在言っていいかどうかは諸説ありますが) MD5 *18 をシェーダー用に改造したものです。 *19

さすがに大丈夫そうです。綺麗なホワイトノイズですね。

key	value
Instruction Mix	227
GPU Duration	78.85 μs
FPS	2748
PractRand Failed	> 2³⁸

もちろん、統計的な検定はパスできます。
ただしさすがに重く、めちゃくちゃ時間がかかりそうだったので途中で止めてしまいました。
腐っても暗号論的ハッシュ関数なので大丈夫だと思います。多分。

MurmurHash3

uint fmix32(uint h)  
{  
    h ^= h >> 16;  
    h *= 0x85ebca6b;  
    h ^= h >> 13;  
    h *= 0xc2b2ae35;  
    h ^= h >> 16;  
    return h;  
}  
  
uint murmur34_orig(uint4 seed)  
{  
    uint c1 = 0xcc9e2d51, c2 = 0x1b873593;  
    uint h = 0;  
    uint k = seed.x;  
  
    k *= c1;  
    k = rotl(k, 15);  
    k *= c2;  
  
    h ^= k;  
    h = rotl(h, 13);  
    h = h * 5 + 0xe6546b64;  
  
    k = seed.y;  
  
    k *= c1;  
    k = rotl(k, 15);  
    k *= c2;  
  
    h ^= k;  
    h = rotl(h, 13);  
    h = h * 5 + 0xe6546b64;  
  
    k = seed.z;  
  
    k *= c1;  
    k = rotl(k, 15);  
    k *= c2;  
  
    h ^= k;  
    h = rotl(h, 13);  
    h = h * 5 + 0xe6546b64;  
  
    k = seed.w;  
  
    k *= c1;  
    k = rotl(k, 15);  
    k *= c2;  
  
    h ^= k;  
    h = rotl(h, 13);  
    h = h * 5 + 0xe6546b64;  
  
    h ^= 16;  
    return fmix32(h);  
}  
  
float murmur34(float4 v)  
{  
    return murmur34_orig(v) * 2.3283064365386962890625e-10;  
}

MurmurHash3 は、ハッシュ関数のテストスイートである SMHasher 内に実装されているハッシュ関数です。
暗号論的ではないですが、そのぶん速度的にも品質的にも性能が良いことで有名ですね。

今回は、 32bit ベースの実装である MurmurHash3_x86_32 を使います。

大丈夫そうですね。綺麗なホワイトノイズです。

key	value
Instruction Mix	43
GPU Duration	39.94 μs
FPS	2683
PractRand Failed	2⁴¹

テストも問題なく、 TiB の壁を突破しました。

PCG

uint pcg_orig(uint v)  
{  
    uint state = v * 747796405 + 2891336453;  
    uint word = ((state >> ((state >> 28) + 4)) ^ state) * 277803737;  
    return word ^ word >> 22;  
}  
  
float pcg(float4 v)  
{  
    return pcg_orig(pcg_orig(pcg_orig(pcg_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w)) * 2.3283064365386962890625e-10;  
}

PCG は、 2014 年に O'neill 氏によって開発された擬似乱数生成器です。 *20

簡単に PCG について説明すると、かの有名な線形合同法に特殊な出力関数をかませることで品質を向上させつつ、内部理論が明らかな線形合同法のメリット (ジャンプが使えるなど) も得ることができるという一挙両得な擬似乱数生成器です。

PCG にはバリアントがたくさんありますが、中でも pcg_oneseq_32_rxs_m_xs_32_random_r を使います。
要するに、

内部状態 32 bit
単一シーケンス
出力関数が rxs_m_xs_32
- rightshift - xorshift - multiply - xorshift 32bit の意

という意味です。

key	value
Instruction Mix	38
GPU Duration	29.70 μs
FPS	2704
PractRand Failed	2³⁸

テストには失敗していますが、大健闘しています。
これも本業は擬似乱数生成器なのですごいです。
出力関数に重きを置いているのが功を奏している感じですね。

PCG2D

uint2 pcg2d_orig(uint2 v)  
{  
    v = v * 1664525 + 1013904223;  
  
    v.x += v.y * 1664525;  
    v.y += v.x * 1664525;  
  
    v = v ^ v >> 16;  
  
    v.x += v.y * 1664525;  
    v.y += v.x * 1664525;  
  
    v = v ^ v >> 16;  
  
    return v;  
}  
  
float pcg2d(float4 v)  
{  
    return (dot(pcg2d_orig(uint2(v.xy)), 1) + dot(pcg2d_orig(uint2(v.zw)), 1)) * 2.3283064365386962890625e-10;  
}

これは、 2020 年に Mark Jarzynski 氏らが設計したハッシュ関数です。 *21
名前の通り PCG を設計のコンセプトとしていますが、その実態は大きく異なり、ハッシュ関数となっています。

key	value
Instruction Mix	37
GPU Duration	27.65 μs
FPS	2664
PractRand Failed	2²⁷

見た目は大丈夫そうですが、テストには失敗しています。

PCG3D

uint3 pcg3d_orig(uint3 v)  
{  
    v = v * 1664525 + 1013904223;  
  
    v.x += v.y * v.z;  
    v.y += v.z * v.x;  
    v.z += v.x * v.y;  
  
    v = v ^ v >> 16;  
  
    v.x += v.y * v.z;  
    v.y += v.z * v.x;  
    v.z += v.x * v.y;  
  
    return v;  
}  
  
float pcg3d(float4 v)  
{  
    return (dot(pcg3d_orig(v.xyz), 1) + dot(pcg3d_orig(v.w), 1)) * 2.3283064365386962890625e-10;  
}

これも同じく、 PCG にインスパイアされて設計されたハッシュ関数です。
Unreal Engine が典拠らしいです。
3 入力 3 出力しかなかったので、無理やり 4 入力 1 出力に拡張しました。

key	value
Instruction Mix	38
GPU Duration	28.67 μs
FPS	2700
PractRand Failed	2⁴²

比較的軽量で、テストに通っています。

PCG3D16

uint3 pcg3d16_orig(uint3 v)  
{  
    v = v * 12829 + 47989;  
  
    v.x += v.y * v.z;  
    v.y += v.z * v.x;  
    v.z += v.x * v.y;  
  
    v.x += v.y * v.z;  
    v.y += v.z * v.x;  
    v.z += v.x * v.y;  
  
    v >>= 16;  
  
    return v;  
}  
  
float pcg3d16(float4 v)  
{  
    uint3 a = pcg3d16_orig(v.xyz);  
    uint3 b = pcg3d16_orig(float3(v.w, 0, 0));  
      
    return ((dot(a, 1) + dot(b, 1)) & 65535) * 1.52587890625e-5;  
}

PCG3D の線形合同法の定数を 16 ビットにして、出力も 16 ビットに絞ったバージョンです。
これも 3 入力 3 出力しかなかったので無理矢理拡張しました。

key	value
Instruction Mix	30
GPU Duration	27.65 μs
FPS	2706
PractRand Failed	2²⁵

見た目は問題なさそうですが、テスト結果によれば品質は落ちてしまっています。

PCG4D

uint4 pcg4d_orig(uint4 v)  
{  
    v = v * 1664525 + 1013904223;  
  
    v.x += v.y * v.w;  
    v.y += v.z * v.x;  
    v.z += v.x * v.y;  
    v.w += v.y * v.z;  
  
    v = v ^ v >> 16;  
  
    v.x += v.y * v.w;  
    v.y += v.z * v.x;  
    v.z += v.x * v.y;  
    v.w += v.y * v.z;  
  
    return v;  
}  
  
float pcg4d(float4 v)  
{  
    uint4 a = pcg4d_orig(v);  
    return dot(a, 1) * 2.3283064365386962890625e-10;  
}

これは、 PCG2D と同じく 2020 年に Mark Jarzynski 氏らが設計したハッシュ関数です。 *22

key	value
Instruction Mix	29
GPU Duration	27.65 μs
FPS	2652
PractRand Failed	2⁴²

それなりに速く、かつ TiB の壁を超える頑健性も示しています。

Pseudo

float pseudo_orig(float2 v)  
{  
    v = frac(v / 128.0) * 128.0 + float2(-64.340622, -72.465622);  
    return frac(dot(v.xyx * v.xyy, float3(20.390625, 60.703125, 2.4281209)));  
}  
  
float pseudo(float4 v)  
{  
    return pseudo_orig(v.xy + v.zw);  
}

これも Unreal Engine が典拠らしいです。

key	value
Instruction Mix	20
GPU Duration	27.65 μs
FPS	2605
PractRand Failed	2¹⁶

出力にもパターンが見える気がしますね。
案の定テストにも落ちています。

Ranlim32

uint ranlim32_orig(uint j)  
{  
    uint u, v, w1, w2, x, y;  
  
    v = 2244614371;  
    w1 = 521288629;  
    w2 = 362436069;  
  
    u = j ^ v;  
  
    u = u * 2891336453 + 1640531513;  
    v ^= v >> 13;  
    v ^= v << 17;  
    v ^= v >> 5;  
    w1 = 33378 * (w1 & 0xffff) + (w1 >> 16);  
    w2 = 57225 * (w2 & 0xffff) + (w2 >> 16);  
  
    v = u;  
  
    u = u * 2891336453 + 1640531513;  
    v ^= v >> 13;  
    v ^= v << 17;  
    v ^= v >> 5;  
    w1 = 33378 * (w1 & 0xffff) + (w1 >> 16);  
    w2 = 57225 * (w2 & 0xffff) + (w2 >> 16);  
  
    x = u ^ u << 9;  
    x ^= x >> 17;  
    x ^= x << 6;  
    y = w1 ^ w1 << 17;  
    y ^= y >> 15;  
    y ^= y << 5;  
  
    return (x + v) ^ (y + w2);  
}  
  
float ranlim32(float4 v)  
{  
    return ranlim32_orig(ranlim32_orig(ranlim32_orig(ranlim32_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w)) * 2.3283064365386962890625e-10;  
}

Ranlim32 は Numerical Recipes 3rd edition に記載されています。 *23

なんかとにかく全部乗せみたいな感じですね。線形合同法、 Xorshift などの要素が含まれています。

key	value
Instruction Mix	79
GPU Duration	27.65 μs
FPS	2595
PractRand Failed	2²⁸

見た目上は問題ありませんが、テストに落ちています。
他のハッシュ関数と比べると結構重いので厳しいところ。

Superfast

uint superfast4_orig(uint4 data)  
{  
    uint hash = 8, tmp;  
  
    hash += data.x & 0xffff;  
    tmp = (((data.x >> 16) & 0xffff) << 11) ^ hash;  
    hash = hash << 16 ^ tmp;  
    hash += hash >> 11;  
  
    hash += data.y & 0xffff;  
    tmp = (((data.y >> 16) & 0xffff) << 11) ^ hash;  
    hash = hash << 16 ^ tmp;  
    hash += hash >> 11;  
  
    hash += data.z & 0xffff;  
    tmp = (((data.z >> 16) & 0xffff) << 11) ^ hash;  
    hash = hash << 16 ^ tmp;  
    hash += hash >> 11;  
  
    hash += data.w & 0xffff;  
    tmp = (((data.w >> 16) & 0xffff) << 11) ^ hash;  
    hash = hash << 16 ^ tmp;  
    hash += hash >> 11;  
  
    hash ^= hash << 3;  
    hash += hash >> 5;  
    hash ^= hash << 4;  
    hash += hash >> 17;  
    hash ^= hash << 25;  
    hash += hash >> 6;  
  
    return hash;  
}  
  
float superfast4(float4 v)  
{  
    return superfast4_orig(v) * 2.3283064365386962890625e-10;   
}

Superfast は、 Paul Hsieh 氏によって開発されたハッシュ関数です。 *24

CPU においては CRC32 や FNV に比べて高速に実行できたらしいです。
現代においてはちょっと命令数が多めに見えますが果たして……？

key	value
Instruction Mix	43
GPU Duration	26.62 μs
FPS	2636
PractRand Failed	2¹⁹

見た目上は問題なさそうに見えますが、テストに落ちてしまっています。

TEA

uint2 scrambleTea(uint2 v, int rounds)  
{  
    uint y = v.x;  
    uint z = v.y;  
    uint sum = 0;  
  
    for (int i = 0; i < rounds; i++)  
    {  
        sum += 0x9e3779b9;  
        y += ((z << 4) + 0xa341316c) ^ (z + sum) ^ ((z >> 5) + 0xc8013ea4);  
        z += ((y << 4) + 0xad90777d) ^ (y + sum) ^ ((y >> 5) + 0x7e95761e);  
    }  
  
    return uint2(y, z);  
}  
  
float tea4(float4 v)  
{  
    return dot(scrambleTea(v.xy, 4) + scrambleTea(v.zw, 4), 1) * 2.3283064365386962890625e-10;   
}

TEA (Tiny Encription Algorithm) は、軽量なブロック暗号アルゴリズムです。 *25

ラウンド数を指定できます。オリジナルでは 32 だったようですが、重すぎるので今回は 4 とします。

key	value
Instruction Mix	87
GPU Duration	29.70 μs
FPS	2626
PractRand Failed	2²¹

織物のような模様が見えますね。
テストにも落ちています。

Trig

float trig_orig(float2 pos)  
{  
    return frac(sin(dot(pos, float2(12.9898, 78.233))) * 43758.5453123);  
}  
  
float trig(float4 pos)  
{  
    return frac(sin(dot(pos, float4(12.9898, 78.233, 42.234, 25.3589))) * 43758.5453123);  
}

有名なワンライナー乱数です。
名前はたぶん Trigonometric functions (三角関数) からきています。
あの The Book of Shaders にも載っています。

コードや定数でググれば多数の使用例が出てきます。例えば、 Unity の ShaderGraph の Random Range ノードでもこの実装が使われています。

利点としては、実装が簡単なことが挙げられます。
欠点としては、三角関数を利用しているため異なる環境下で再現性がないこと、案外わかりやすいパターンがあることが挙げられます。

key	value
Instruction Mix	11
GPU Duration	27.65 μs
FPS	2639
PractRand Failed	2¹⁶

特徴的な縞模様になっており、テストにも即座に落ちています。

Wang

uint wang_orig(uint v)  
{  
    v = (v ^ 61) ^ (v >> 16);  
    v *= 9;  
    v ^= v >> 4;  
    v *= 0x27d4eb2d;  
    v ^= v >> 15;  
    return v;  
}  
  
float wang(float4 v)  
{  
    return wang_orig(wang_orig(wang_orig(wang_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w))* 2.3283064365386962890625e-10;   
}

Wang は、 Thomas Wang 氏が 1997 年に発表したハッシュ関数です。 *26

key	value
Instruction Mix	41
GPU Duration	27.65 μs
FPS	2639
PractRand Failed	2³⁵

健闘しましたが、テスト結果は 1 TiB には届きませんでした。

Xorshift128

uint4 xorshift128_orig(uint4 v)  
{  
    v.w ^= v.w << 11;  
    v.w ^= v.w >> 8;  
    v = v.wxyz;  
    v.x ^= v.y;  
    v.x ^= v.y >> 19;  
  
    return v;  
}  
  
float xorshift128(float4 v)  
{  
    uint4 vv = xorshift128_orig(uint4(v));  
    return dot(vv, 1) * 2.3283064365386962890625e-10;  
}

Xorshift128 は、 George Marsaglia 氏が 2003 年に発表した擬似乱数生成器です。 *27

本業は擬似乱数生成器なので、 v.y, v.z は v.x, v.y そのままになっているなど、ハッシュ関数にするには若干怪しいところがあります。

案の定でした。

key	value
Instruction Mix	10
GPU Duration	26.62 μs
FPS	2601
PractRand Failed	×

論外なのでテストしていません。

Xorshift32

uint xorshift32_orig(uint v)  
{  
    v ^= v << 13;  
    v ^= v >> 17;  
    v ^= v << 5;  
    return v;  
}  
  
float xorshift32(float4 v)  
{  
    return xorshift32_orig(xorshift32_orig(xorshift32_orig(xorshift32_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w)) * 2.3283064365386962890625e-10;  
}

Xorshift32 は、同じく George Marsaglia 氏が 2003 年に発表した擬似乱数生成器です。 *28
これは半分ハッシュ関数のような設計になっているので、 xorshift128 よりはいい成績になりそうですが……

key	value
Instruction Mix	33
GPU Duration	30.72 μs
FPS	2601
PractRand Failed	2¹⁶

さざ波のような模様がありますね。
残念ながらテストにも即座に落ちてしまっています。

xxHash32

uint xxhash_orig(uint4 value)  
{  
    uint XXH_PRIME32_1 = 0x9E3779B1;  
    uint XXH_PRIME32_2 = 0x85EBCA77;  
    uint XXH_PRIME32_3 = 0xC2B2AE3D;  
  
    uint4 state = uint4(XXH_PRIME32_1 + XXH_PRIME32_2, XXH_PRIME32_2, 0, -XXH_PRIME32_1);  
  
    state += value * XXH_PRIME32_2;  
    state = rotl(state, 13);  
    state *= XXH_PRIME32_1;  
  
    uint h32 = rotl(state[0], 1) + rotl(state[1], 7) + rotl(state[2], 12) + rotl(state[3], 18);  
    h32 += 16;  
  
    h32 ^= h32 >> 15;  
    h32 *= XXH_PRIME32_2;  
    h32 ^= h32 >> 13;  
    h32 *= XXH_PRIME32_3;  
    h32 ^= h32 >> 16;  
  
    return h32;  
}  
  
float xxhash(float4 v)  
{  
    return xxhash_orig(v) * 2.3283064365386962890625e-10;   
}

xxHash は、軽量な非暗号論的ハッシュ関数です。 *29

今回はその中でも 32bit ベースの xxHash32 を用います。

key	value
Instruction Mix	42
GPU Duration	27.65 μs
FPS	2676
PractRand Failed	2²⁷

見た目上は問題なさそうですが、テストには落ちてしまいました。意外です。

license

BSD 2-Clause License (https://www.opensource.org/licenses/bsd-license.php)

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright notice, this
list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Wyhash

uint2 mul64(uint a, uint b)  
{  
    uint alo = a & 0xffff, ahi = a >> 16;  
    uint blo = b & 0xffff, bhi = b >> 16;  
  
    uint lo = alo * blo;  
    uint mid1 = ahi * blo;  
    uint mid2 = alo * bhi;  
    uint hi = ahi * bhi;  
  
    return uint2(  
        lo + (mid1 << 16) + (mid2 << 16),  
        hi + (mid1 >> 16) + (mid2 >> 16) + (((mid1 & 0xffff) + (mid2 & 0xffff) + (lo >> 16)) >> 16));  
}  
  
void wymix32(inout uint a, inout uint b)  
{  
    uint2 c = uint2(a ^ 0x53c5ca59u, 0);  
    c = mul64(c.x, b ^ 0x74743c1b);  
  
    a = c.x;  
    b = c.y;  
}  
  
uint wyhash_orig(float4 value, uint seed)  
{  
    uint see1 = 16;  
    wymix32(seed, see1);  
      
    seed ^= uint(value.x);  
    see1 ^= uint(value.y);  
    wymix32(seed, see1);  
      
    seed ^= uint(value.z);  
    see1 ^= uint(value.w);  
    wymix32(seed, see1);  
      
    wymix32(seed, see1);  
    wymix32(seed, see1);  
    return seed ^ see1;  
}  
  
float wyhash(float4 v)  
{  
    return wyhash_orig(v, 0xa0b428db) * 2.3283064365386962890625e-10;  
}

Wyhash は、 wangyi-fudan 氏によって 2019 年に開発された高速なハッシュ関数です。 *30
2 倍長乗算 (32bit x 32bit -> 64bit) をうまく使っているハッシュ関数です。
CPU ならこの乗算は専用命令があってすぐできるのですが、シェーダー上では intrinsics が存在しないので筆算法で解きました。

なお、同じ GPU でも OpenGL (GLSL) には umulExtended() 関数があり、 2 倍長乗算が一発で計算できるらしいです。うらやましい。
HLSL ではまだ検討段階らしいです。

2 倍長乗算をためしに Karatsuba 法で実装してみたらむしろ遅くなった、という悲しい事件がありました。ここに供養しておきます。

// csharp karatsuba multiplication  
public static (uint lo, uint hi) mulhi(uint a, uint b)  
{  
    uint alo = a & 0xffff, ahi = a >> 16;  
    uint blo = b & 0xffff, bhi = b >> 16;  
  
    uint hihi = ahi * bhi;  
    uint lolo = alo * blo;  
  
    uint aa = alo - ahi;  
    uint bb = bhi - blo;  
    uint mid = aa * bb;  
    uint midSign = (mid != 0 && ((aa ^ bb) >> 31) != 0) ? 1u : 0u;  
  
    uint lo = lolo + (mid << 16) + (lolo << 16) + (hihi << 16);  
    uint carry = (lolo >> 16) + (mid & 0xffff) + (lolo & 0xffff) + (hihi & 0xffff);  
    uint hi = hihi + (mid >> 16) + (lolo >> 16) + (hihi >> 16) + (carry >> 16) - (midSign << 16);  
  
    return (lo, hi);  
}

key	value
Instruction Mix	87
GPU Duration	28.67 μs
FPS	2636
PractRand Failed	2⁴²

2 倍長乗算のせいか命令数は増えてしまっていますが、テストに問題なくパスしています。

なお、 wyhash には rapidhash という後継があるようですが、 64 x 64 = 128 bit の 2 倍長乗算が必要なためシェーダーで高効率に実装するのは難しそうです。

lowbias32

uint lowbias32_orig(uint x)  
{  
    x ^= x >> 16;  
    x *= 0x7feb352d;  
    x ^= x >> 15;  
    x *= 0x846ca68b;  
    x ^= x >> 16;  
    return x;  
}  
  
float lowbias32(float4 v)  
{  
    return lowbias32_orig(lowbias32_orig(lowbias32_orig(lowbias32_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w)) * 2.3283064365386962890625e-10;  
}

これは、 Chris Wellons 氏によって 2018 年に発表されたハッシュ関数です。 *31
品質の良いハッシュを探索するツールによって発見されたそうです。

key	value
Instruction Mix	41
GPU Duration	28.67 μs
FPS	2638
PractRand Failed	2⁴²

比較的高速で、テストにもパスしています。

triple32

uint triple32_orig(uint x)  
{  
    x ^= x >> 17;  
    x *= 0xed5ad4bb;  
    x ^= x >> 11;  
    x *= 0xac4c1b51;  
    x ^= x >> 15;  
    x *= 0x31848bab;  
    x ^= x >> 14;  
    return x;  
}  
  
float triple32(float4 v)  
{  
    return triple32_orig(triple32_orig(triple32_orig(triple32_orig(uint(v.x)) + uint(v.y)) + uint(v.z)) + uint(v.w)) * 2.3283064365386962890625e-10;  
}

これも Chris Wellons 氏によって 2018 年に発表されたハッシュ関数です。 *32
lowbias32 に比べて命令数が増えているものの品質が向上しており、「ランダムな置換 (並べ替え) と統計的に区別できない」とまで言われています。

key	value
Instruction Mix	53
GPU Duration	27.65 μs
FPS	2603
PractRand Failed	2³⁹

なぜか lowbias32 よりもテスト結果が悪くなっています。どうして？

fihash

float fihash_orig(float2 v)  
{  
    uint2 u = asint(v * float2(141421356, 2718281828));  
    return float((u.x ^ u.y) * 3141592653) * 2.3283064365386962890625e-10;  
}  
  
float fihash(float4 v)  
{  
    return fihash_orig(v.xy + v.zw);  
}

2024 年に Lumi 氏によって作成されたハッシュ関数です。 *33

key	value
Instruction Mix	9
GPU Duration	28.67 μs
FPS	2642
PractRand Failed	2¹⁶

命令数がずば抜けて少ないです。それはそう。
見た目上は問題なさそうですが、テストには即座に落ちてしまっています。

ちなみに、 2018 年に James_Harnett 氏によってほぼ同じ構造のハッシュ関数が提案されていることを注記しておきます。
(シンプルな構造ゆえ、シンクロニシティ的に生まれたものかと思います。たぶん。)
https://www.shadertoy.com/view/MdcfDj

license

// The MIT License
// Copyright © 2024 Giorgi Azmaipharashvili
// Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

fast32hash

float4 fast32hash_orig(float2 v)  
{  
    const float2 offset = float2(26, 161);  
    const float domain = 71;  
    const float someLargeFloat = 951.135664;  
  
    float4 p = float4(v, v + 1);  
    p = p - floor(p * (1.0 / domain)) * domain;  
    p += offset.xyxy;  
    p *= p;  
    return frac(p.xzxz * p.yyww * (1.0 / someLargeFloat));  
}  
  
float fast32hash(float4 x)  
{  
    return fast32hash_orig(x.xy + x.zw).x;  
}

2011 年に Brian Sharpe 氏が発表したハッシュ関数です。 *34

uint を使わずに実装できるため、整数計算が遅い GPU でも高速に実行できることが期待されます。

key	value
Instruction Mix	17
GPU Duration	28.67 μs
FPS	2653
PractRand Failed	2¹⁶

明らかに繰り返しのパターンが見えますね。
テストにも即座に落ちてしまっています。

Philox

uint4 philox_round(uint2 key, uint4 ctr)  
{  
   uint2 lohi0 = mul64(ctr.x, 0xD2511F53);  
   uint2 lohi1 = mul64(ctr.z, 0xCD9E8D57);  
  
   return uint4(lohi1.y ^ ctr.y ^ key.x, lohi1.x, lohi0.y ^ ctr.w ^ key.y, lohi0.x);  
}  
  
uint2 philox_bumpkey(uint2 key)  
{  
    return key + uint2(0x9E3779B9, 0xBB67AE85);  
}  
  
uint4 philox_orig(uint2 key, uint4 ctr)  
{  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
    key = philox_bumpkey(key);  
    ctr = philox_round(key, ctr);  
     
    return ctr;  
}  
  
float philox(float4 v)  
{  
    uint4 u = v;  
    return philox_orig(uint2(0xf19cd101, 0x3d30), u).x * 2.3283064365386962890625e-10;  
}

Philox は、 2011 年に発表された擬似乱数生成器です。 *35
なかでも、 32 bit 環境でスタンダードな Philox-4x32-10 を使用しました。数字は 4 個 x 32 bit - 10 ラウンドを表しています。

key	value
Instruction Mix	294
GPU Duration	62.46 μs
FPS	2729
PractRand Failed	2⁴²

命令数が多い分、テストは問題なさそうです。
ちょっと重いのが難点かも。

AES-CTR

static const int AesSbox[256] = {  
    0x63, 0x7c, 0x77, 0x7b, 0xf2, 0x6b, 0x6f, 0xc5, 0x30, 0x01, 0x67, 0x2b, 0xfe, 0xd7, 0xab, 0x76,  
    0xca, 0x82, 0xc9, 0x7d, 0xfa, 0x59, 0x47, 0xf0, 0xad, 0xd4, 0xa2, 0xaf, 0x9c, 0xa4, 0x72, 0xc0,  
    0xb7, 0xfd, 0x93, 0x26, 0x36, 0x3f, 0xf7, 0xcc, 0x34, 0xa5, 0xe5, 0xf1, 0x71, 0xd8, 0x31, 0x15,  
    0x04, 0xc7, 0x23, 0xc3, 0x18, 0x96, 0x05, 0x9a, 0x07, 0x12, 0x80, 0xe2, 0xeb, 0x27, 0xb2, 0x75,  
    0x09, 0x83, 0x2c, 0x1a, 0x1b, 0x6e, 0x5a, 0xa0, 0x52, 0x3b, 0xd6, 0xb3, 0x29, 0xe3, 0x2f, 0x84,  
    0x53, 0xd1, 0x00, 0xed, 0x20, 0xfc, 0xb1, 0x5b, 0x6a, 0xcb, 0xbe, 0x39, 0x4a, 0x4c, 0x58, 0xcf,  
    0xd0, 0xef, 0xaa, 0xfb, 0x43, 0x4d, 0x33, 0x85, 0x45, 0xf9, 0x02, 0x7f, 0x50, 0x3c, 0x9f, 0xa8,  
    0x51, 0xa3, 0x40, 0x8f, 0x92, 0x9d, 0x38, 0xf5, 0xbc, 0xb6, 0xda, 0x21, 0x10, 0xff, 0xf3, 0xd2,  
    0xcd, 0x0c, 0x13, 0xec, 0x5f, 0x97, 0x44, 0x17, 0xc4, 0xa7, 0x7e, 0x3d, 0x64, 0x5d, 0x19, 0x73,  
    0x60, 0x81, 0x4f, 0xdc, 0x22, 0x2a, 0x90, 0x88, 0x46, 0xee, 0xb8, 0x14, 0xde, 0x5e, 0x0b, 0xdb,  
    0xe0, 0x32, 0x3a, 0x0a, 0x49, 0x06, 0x24, 0x5c, 0xc2, 0xd3, 0xac, 0x62, 0x91, 0x95, 0xe4, 0x79,  
    0xe7, 0xc8, 0x37, 0x6d, 0x8d, 0xd5, 0x4e, 0xa9, 0x6c, 0x56, 0xf4, 0xea, 0x65, 0x7a, 0xae, 0x08,  
    0xba, 0x78, 0x25, 0x2e, 0x1c, 0xa6, 0xb4, 0xc6, 0xe8, 0xdd, 0x74, 0x1f, 0x4b, 0xbd, 0x8b, 0x8a,  
    0x70, 0x3e, 0xb5, 0x66, 0x48, 0x03, 0xf6, 0x0e, 0x61, 0x35, 0x57, 0xb9, 0x86, 0xc1, 0x1d, 0x9e,  
    0xe1, 0xf8, 0x98, 0x11, 0x69, 0xd9, 0x8e, 0x94, 0x9b, 0x1e, 0x87, 0xe9, 0xce, 0x55, 0x28, 0xdf,  
    0x8c, 0xa1, 0x89, 0x0d, 0xbf, 0xe6, 0x42, 0x68, 0x41, 0x99, 0x2d, 0x0f, 0xb0, 0x54, 0xbb, 0x16,  
};  
  
uint aesctr_subword(uint x)  
{  
    return uint(AesSbox[x & 0xff]) ^  
        uint(AesSbox[(x >> 8) & 0xff]) << 8 ^  
        uint(AesSbox[(x >> 16) & 0xff]) << 16 ^  
        uint(AesSbox[x >> 24]) << 24;  
}  
  
uint4 aesctr_round(uint4 k, uint rcon)  
{  
    //uint keygenassist0 = aesctr_subword(k1);  
    //uint keygenassist1 = keygenassist0 >> 8 ^ keygenassist0 << 24 ^ rcon;  
    uint keygenassist2 = aesctr_subword(k.w);  
    uint keygenassist3 = keygenassist2 >> 8 ^ keygenassist2 << 24 ^ rcon;  
  
    uint4 t = uint4(  
        k.x,  
        k.y ^ k.x,  
        k.z ^ k.y ^ k.x,  
        k.w ^ k.z ^ k.y ^ k.x);  
  
    t ^= keygenassist3;  
      
    return t;  
}  
  
uint aesctr_multiply(uint x, uint y)  
{  
    uint result = 0;  
    for (uint mask = 1; mask < 256; mask <<= 1)  
    {  
        if (y & mask)  
        {  
            result ^= x;  
        }  
        x = x << 1 ^ ((x & 0x80) ? 0x1b : 0);  
    }  
    return result;  
}  
  
uint4 aesctr_mixcolumns(uint4 u)  
{  
    uint4 r;  
    r.x =   
        uint(aesctr_multiply(2, u.x & 0xff) ^ aesctr_multiply(3, ((u.x >> 8) & 0xff)) ^ ((u.x >> 16) & 0xff) ^ (u.x >> 24)) ^  
        uint((u.x & 0xff) ^ aesctr_multiply(2, (u.x >> 8) & 0xff) ^ aesctr_multiply(3, (u.x >> 16) & 0xff) ^ (u.x >> 24)) << 8 ^  
        uint((u.x & 0xff) ^ ((u.x >> 8) & 0xff) ^ aesctr_multiply(2, (u.x >> 16) & 0xff) ^ aesctr_multiply(3, (u.x >> 24))) << 16 ^  
        uint(aesctr_multiply(3, u.x & 0xff) ^ ((u.x >> 8) & 0xff) ^ ((u.x >> 16) & 0xff) ^ aesctr_multiply(2, (u.x >> 24))) << 24;  
    r.y =   
        uint(aesctr_multiply(2, u.y & 0xff) ^ aesctr_multiply(3, ((u.y >> 8) & 0xff)) ^ ((u.y >> 16) & 0xff) ^ (u.y >> 24)) ^  
        uint((u.y & 0xff) ^ aesctr_multiply(2, (u.y >> 8) & 0xff) ^ aesctr_multiply(3, (u.y >> 16) & 0xff) ^ (u.y >> 24)) << 8 ^  
        uint((u.y & 0xff) ^ ((u.y >> 8) & 0xff) ^ aesctr_multiply(2, (u.y >> 16) & 0xff) ^ aesctr_multiply(3, (u.y >> 24))) << 16 ^  
        uint(aesctr_multiply(3, u.y & 0xff) ^ ((u.y >> 8) & 0xff) ^ ((u.y >> 16) & 0xff) ^ aesctr_multiply(2, (u.y >> 24))) << 24;  
    r.z =   
        uint(aesctr_multiply(2, u.z & 0xff) ^ aesctr_multiply(3, ((u.z >> 8) & 0xff)) ^ ((u.z >> 16) & 0xff) ^ (u.z >> 24)) ^  
        uint((u.z & 0xff) ^ aesctr_multiply(2, (u.z >> 8) & 0xff) ^ aesctr_multiply(3, (u.z >> 16) & 0xff) ^ (u.z >> 24)) << 8 ^  
        uint((u.z & 0xff) ^ ((u.z >> 8) & 0xff) ^ aesctr_multiply(2, (u.z >> 16) & 0xff) ^ aesctr_multiply(3, (u.z >> 24))) << 16 ^  
        uint(aesctr_multiply(3, u.z & 0xff) ^ ((u.z >> 8) & 0xff) ^ ((u.z >> 16) & 0xff) ^ aesctr_multiply(2, (u.z >> 24))) << 24;  
    r.w =   
        uint(aesctr_multiply(2, u.w & 0xff) ^ aesctr_multiply(3, ((u.w >> 8) & 0xff)) ^ ((u.w >> 16) & 0xff) ^ (u.w >> 24)) ^  
        uint((u.w & 0xff) ^ aesctr_multiply(2, (u.w >> 8) & 0xff) ^ aesctr_multiply(3, (u.w >> 16) & 0xff) ^ (u.w >> 24)) << 8 ^  
        uint((u.w & 0xff) ^ ((u.w >> 8) & 0xff) ^ aesctr_multiply(2, (u.w >> 16) & 0xff) ^ aesctr_multiply(3, (u.w >> 24))) << 16 ^  
        uint(aesctr_multiply(3, u.w & 0xff) ^ ((u.w >> 8) & 0xff) ^ ((u.w >> 16) & 0xff) ^ aesctr_multiply(2, (u.w >> 24))) << 24;  
    return r;  
}  
  
uint aesctr_orig(uint4 x)  
{  
    uint4 seed[11];  
  
    seed[0] = x;  
  
    seed[ 1] = aesctr_round(seed[0], 0x01);  
    seed[ 2] = aesctr_round(seed[1], 0x02);  
    seed[ 3] = aesctr_round(seed[2], 0x04);  
    seed[ 4] = aesctr_round(seed[3], 0x08);  
    seed[ 5] = aesctr_round(seed[4], 0x10);  
    seed[ 6] = aesctr_round(seed[5], 0x20);  
    seed[ 7] = aesctr_round(seed[6], 0x40);  
    seed[ 8] = aesctr_round(seed[7], 0x80);  
    seed[ 9] = aesctr_round(seed[8], 0x1b);  
    seed[10] = aesctr_round(seed[9], 0x36);  
  
    uint4 ctr = uint4(1, 0, 0, 0);  
  
    uint4 t = ctr ^ seed[0];  
  
    for (int i = 1; i <= 9; i++)  
    {  
        uint4 u;  
  
        u.x = ((t.x >> 0) & 0xff) << 0 ^  
            ((t.y >> 8) & 0xff) << 8 ^  
            ((t.z >> 16) & 0xff) << 16 ^  
            ((t.w >> 24) & 0xff) << 24;  
        u.y = ((t.y >> 0) & 0xff) << 0 ^  
            ((t.z >> 8) & 0xff) << 8 ^  
            ((t.w >> 16) & 0xff) << 16 ^  
            ((t.x >> 24) & 0xff) << 24;  
        u.z = ((t.z >> 0) & 0xff) << 0 ^  
            ((t.w >> 8) & 0xff) << 8 ^  
            ((t.x >> 16) & 0xff) << 16 ^  
            ((t.y >> 24) & 0xff) << 24;  
        u.w = ((t.w >> 0) & 0xff) << 0 ^  
            ((t.x >> 8) & 0xff) << 8 ^  
            ((t.y >> 16) & 0xff) << 16 ^  
            ((t.z >> 24) & 0xff) << 24;  
  
        u.x = aesctr_subword(u.x);  
        u.y = aesctr_subword(u.y);  
        u.z = aesctr_subword(u.z);  
        u.w = aesctr_subword(u.w);  
  
        u = aesctr_mixcolumns(u);  
  
        t = u ^ seed[i];  
    }  
  
    if (++ctr.x == 0) {  
        if (++ctr.y == 0) {  
            if (++ctr.z == 0) {  
                ++ctr.w;  
            }  
        }  
    }  
  
    {  
        uint4 u;  
  
        u.x = ((t.x >> 0) & 0xff) << 0 ^  
            ((t.y >> 8) & 0xff) << 8 ^  
            ((t.z >> 16) & 0xff) << 16 ^  
            ((t.w >> 24) & 0xff) << 24;  
        u.y = ((t.y >> 0) & 0xff) << 0 ^  
            ((t.z >> 8) & 0xff) << 8 ^  
            ((t.w >> 16) & 0xff) << 16 ^  
            ((t.x >> 24) & 0xff) << 24;  
        u.z = ((t.z >> 0) & 0xff) << 0 ^  
            ((t.w >> 8) & 0xff) << 8 ^  
            ((t.x >> 16) & 0xff) << 16 ^  
            ((t.y >> 24) & 0xff) << 24;  
        u.w = ((t.w >> 0) & 0xff) << 0 ^  
            ((t.x >> 8) & 0xff) << 8 ^  
            ((t.y >> 16) & 0xff) << 16 ^  
            ((t.z >> 24) & 0xff) << 24;  
  
        u.x = aesctr_subword(u.x);  
        u.y = aesctr_subword(u.y);  
        u.z = aesctr_subword(u.z);  
        u.w = aesctr_subword(u.w);  
  
        t = u ^ seed[i];  
    }  
  
    return t;  
}  
  
float aesctr(float4 v)  
{  
    return aesctr_orig(uint4(asint(v.x), asint(v.y), asint(v.z), asint(v.w))) * 2.3283064365386962890625e-10f;  
}

AES (Advanced Encryption Standard) は、 2001 年にアメリカで標準暗号として定められた共通鍵暗号アルゴリズムです。 *36

今回は CTR モードを利用して実装しています。
CPU 上なら AES-NI 専用命令が使えてそれなりに速いのですが、シェーダー上となると完全にソフトウェア実装なのでめちゃくちゃ時間がかかります。

綺麗なんですがめちゃくちゃ重いです。

key	value
Instruction Mix	1021
GPU Duration	4970 μs
FPS	133
PractRand Failed	> 2³⁵

あまりにも重い (64 GiB のテストに半日かかった) のでテストはここまでで諦めました。

heptaplex-collapse noise

// @ENDESGA 2023  
  
uint heptaplex_orig(uint x, uint y, uint z)  
{  
    x = ~(~x - y - z) * ~(x - ~y - z) * ~(x - y - ~z);  
    y = ~(~x - y - z) * ~(x - ~y - z) * ~(x - y - ~z);  
    z = x ^ y ^ (~(~x - y - z) * ~(x - ~y - z) * ~(x - y - ~z));  
    return z ^ ~(~z >> 16);  
}  
  
float heptaplex(float4 v)   
{  
    return (heptaplex_orig(v.x, v.y, v.z) + heptaplex_orig(v.w, 0, 0)) * 2.3283064365386962890625e-10;  
}

heptaplex-collapse noise は、 2023 年に ENDESGA 氏が Shadertoy 上で発表したノイズです。 *37

ぱっと見は綺麗です。

key	value
Instruction Mix	46
GPU Duration	28.67 μs
FPS	2557
PractRand Failed	2¹⁹

テストは即死ではないものの、それなりの早さで落ちてしまいました。

IbukiHash

// IbukiHash by Andante  
// This work is marked with CC0 1.0. To view a copy of this license, visit https://creativecommons.org/publicdomain/zero/1.0/  
  
float ibuki(float4 v)  
{  
    const uint4 mult =   
        uint4(0xae3cc725, 0x9fe72885, 0xae36bfb5, 0x82c1fcad);  
  
    uint4 u = uint4(v);  
    u = u * mult;  
    u ^= u.wxyz ^ u >> 13;  
      
    uint r = dot(u, mult);  
  
    r ^= r >> 11;  
    r = (r * r) ^ r;  
          
    return r * 2.3283064365386962890625e-10;  
}

見慣れない関数かと思いますが、それはそうです。 私が今作りました。
命令数 (Instruction Mix) がなるべく少なく、かつテストには 1 TiB まで通るように設計しました。

key	value
Instruction Mix	26
GPU Duration	27.65 μs
FPS	2681
PractRand Failed	2⁴¹

設計思想

IbukiHash の全ての行に意味があります。
興味がない方は読み飛ばしていただいて大丈夫です。

const uint4 mult =   
    uint4(0xae3cc725, 0x9fe72885, 0xae36bfb5, 0x82c1fcad);

乗算の定数は、詳細は省きますがビットが分散しやすい値にする必要があります。
今回は分散の良い値を研究している論文から引用しました。 *38

絶対にこの定数である必要はないので変更することもできます。たとえば uint4 が返り値としてほしくなった場合に別の定数で計算するなど。
ただ、少なくとも最上位ビットが 1 でかつ最下位の 3 ビットが 0b101 で、それなりにビットが分散している (popcnt(mult) が 16 に近い奇数で、かつ 0x0000ffff みたいに固まっていない) 必要があります。

uint4 u = uint4(v);

asint() ではなく単純なキャストにしているのは、下位ビットに意味のある値を集中させるためです。
今回引数に取るのは座標なので、どうせ $\lbrack -2^{16}, 2^{16} \rbrack$ ぐらいの小さい範囲にしかなりません。そうなると、どこに意味のある情報 (ビット) を置くかが重要になります。
そして今後ビットを攪拌するにあたって、下位ビットに情報があったほうが好都合なためです。

asint() を使うと、 float 型のビットパターンの都合上上位ビットに情報が集まりがち、下位ビットが 0 になりがちであまりよろしくなかったのです。
測定結果を見た限りでは asint() のほうが多少速い感はありましたが、品質が落ちる（早くテストに落ちる）傾向がありました。なので致し方ありません。

u = u * mult;

u に定数を掛け、 $2^{32}$ で割った余りを得ます。
掛け算はハッシュ関数のなかでも非常に重要な要素で、下位ビットの情報を上位ビットに幅広く伝播させることができます。
めちゃくちゃ簡単に言えば、上位ビットのハッシュとしての「品質」を大きく高めることができる演算です。

$2^{n}$ を法としたとき、奇数定数との掛け算 x *= (a | 1) は全単射の操作です。
つまり、特定の値に偏ることがなく、まんべんなくビットを混ぜることができます。

ここで、 u の各要素にそれぞれ別の定数を掛けている理由は、同じ定数を掛けると u.x と u.y の値を交換しても同じ結果が得られるようになってしまい、 $y = x$ を軸に対称となってしまうためです。

HLSL では uint4 どうしの掛け算などをまとめて行えるのでアセンブリもそうやって最適化されるのかと思いきや、少なくとも中間言語 (DXIL) レベルでは別々に計算したとき (u.x *= mult.x; u.y *= mult.y; ...) と同じような命令が発行されているようです。

u ^= u.wxyz ^ u >> 13;

u に xorshift っぽいことをします。
まず u ^ u >> 13 ですが、乗算で品質が上がった上位ビットを下位ビットに右シフトで伝播させ、上も下も品質を向上させます。
加えて x ^= x >> a は乗算と同じく全単射の操作です。同じく特定の値に偏らずまんべんなくビットを混ぜられます。
シフト定数 13 は、 32 bit の半分 (16) より小さい最大の素数なので使っています。シフトを大きくすると下位ビットに伝播しやすくなるのですが、伝播させられる情報量自体は減ってしまうので、バランスよく選びました。

また、 u.wxyz とスウィズルした値を xor することによって、別のビットとの「絡み」を発生させます。今までは u.xyzw のそれぞれの要素の中でだけビットの情報伝播が行われていましたが、ここで他の要素と混ぜることでよりハッシュらしくします。

uint r = dot(u, mult);

次に、内積をとります。言い換えれば、

uint r = u.x * mult.x + u.y * mult.y + u.z * mult.z + u.w * mult.w;

です。乗算によって再度攪拌したのち、全要素を加算してひとつにまとめます。
ひとつにまとめることで今後の命令数を減らせます。
このまとめのタイミングが重要で、早すぎると xyzw の差別化が行えずにテストに落ち、遅すぎると命令数が増えて実行が遅くなります。

また、あえて dot() を使ったのは、命令の最適化が行えるからです。
シェーダーの中間言語である DXIL には、 32 bit 整数の乗算と加算を同時に行える mad 命令があります。 fma (fused multiply add) の整数版のようなイメージです。
dot() を使うと、この mad 命令を活用してこういう感じに展開してくれます。

; dot(u, mult)  
mul %1, u.x, mult.x        ; %1 = u.x * mult.x;  
mad %1, u.y, mult.y, %1    ; %1 = u.y * mult.y + %1;  
mad %1, u.z, mult.z, %1    ; %1 = u.z * mult.z + %1;  
mad %1, u.w, mult.w, %1    ; %1 = u.w * mult.w + %1;

ふつうに乗算と加算で計算するとなぜか mad 命令にしてくれないので、 dot() を使って命令の短縮を図ります。

 r ^= r >> 11;

また xorshift をして、 dot() で品質が向上した上位ビットの情報を下位ビットにも伝播させます。
ここのシフト定数 11 は、前回のシフト定数 13 と互いに素であることから選びました。
そのほうが品質が良くなる傾向があるようです。体感ですが……

r = (r * r) ^ r;

まず、 r = (r * r) ^ (r | 1) は全単射の操作です。
詳しい原理は私は分かっていませんが少なくとも uint ( $2^{32}$ を法とした演算) の範囲では総当たりで全単射になっていることを確認しています。
これが奇数定数の乗算 r *= (a | 1); よりもビットの攪拌性能が良いという噂があり、実際にテスト結果も向上したことから採用しました。
じゃあなんで全部これにしなかったのかというと、 mul だけで済む奇数定数乗算に比べて mul, xor と 1 命令増えてしまうためです。最後の一番大事なところに使いました。

ところで、 r | 1 が抜けているじゃないかと思った方は正しいです。
ですがここでは抜けていてもよいのです。 or がなくなることで全単射操作ではなくなります (最下位ビットが必ず 0 になります) 。ですが、 float の精度は 24 bit ですので、 uint でつくった 32 ビットのうち下位 8 ビットの情報は切り捨てられます。したがって大した問題にはなりません。
1 命令ぶん早くするために理論も犠牲にするという手法 (?) です。

return r * 2.3283064365386962890625e-10;

最後に、 $2^{-32}$ を掛けて $\lbrack 0, 1 )$ の float に戻します。
ここで、代わりに asfloat(0x3f800000 | r >> 9) - 1; とする手法もあります。
asfloat() を使うほうは、要するに $\lbrack 1.0, 2.0 )$ にビットパターン変換して 1 を引くことで $\lbrack 0, 1 )$ に変換する手法です。 asfloat() 法のほうが速くなる場合があるらしいのですが、変換後の最下位ビットが必ず 0 になるため 23 bit 精度になってしまう問題があります。
対して、掛け算のほうでは 24 bit 精度を維持できます。そのため、掛け算を選択しました。

以上が設計のお話です。
いつもは CPU 上で擬似乱数生成器を設計してみたりしているのですが、 GPU (シェーダー) 上となると速くて使える関数や命令が違うので、勝手が違って楽しかったです。

まとめ

みんな大好き比較グラフのお時間です。
FPS が微妙にあてにならない感があったので、命令数で比較することにします。

縦軸が Instruction Mix (命令数; 少ないほうが速いと期待される) 、横軸が PractRand Failed (40 以上は合格とみなしてよい) です。
つまり、右下に行けば行くほど軽くて品質が良いことになります。

IbukiHash は合格したシェーダー乱数のなかでは命令数が一番少なく、軽くて強いことが分かります。
PCG4D もいい線ですね。

また、品質を気にせず速さだけを求めるなら fihash でしょうか。
テストには落ちたものの見た目上は問題なさそうだったので、活用できるかもしれません。

元データの表も貼っておきます。
PractRand Failed → Instruction Mix の順にソートしてあります。

Algorithm	Instruction Mix	GPU Duration	FPS	PractRand Failed
PCG4D	29	27.65	2652	42
PCG3D	38	28.67	2700	42
lowbias32	41	28.67	2638	42
IQInt2	42	26.62	2622	42
Wyhash	87	28.67	2636	42
Philox	294	62.46	2729	42
ibuki	26	27.65	2681	41
MurmurHash3	43	39.94	2683	41
CityHash	49	27.65	2695	41
ESGTSA	38	26.62	2721	40
triple32	53	27.65	2603	39
PCG	38	29.7	2704	38
MD5	227	78.85	2748	> 38
Wang	41	27.65	2586	35
AESCTR	1021	4970	133	> 35
Ranlim32	79	27.65	2595	28
PCG2D	37	27.65	2664	27
xxHash32	42	27.65	2676	27
PCG3D16	30	27.65	2706	25
TEA	87	29.7	2626	21
JenkinsHash	93	27.65	2721	21
Superfast	43	26.62	2636	19
heptaplex-collapse	46	28.67	2557	19
IQInt32	34	27.65	2728	18
IQInt1	30	28.67	2630	17
fihash	9	28.67	2642	16
Interleaved Gradient Noise	10	26.62	2632	16
Trig	11	27.65	2639	16
LCG	14	27.65	2695	16
Fast	16	27.65	2664	16
fast32hash	17	28.67	2653	16
Pseudo	20	27.65	2605	16
PerlinPerm	21	128	2501	16
IQInt3	24	27.65	2580	16
Hash without Sine	31	28.67	2702	16
Xorshift32	33	30.72	2636	16
mod289	39	27.65	2656	16
BBS4093	49	27.65	2667	16
FNV1	50	30.72	2656	16
BBS65521	53	27.65	2735	16
Xorshift128	10	26.62	2601	0
JKISS32	15	26.62	2687	0
HybridTaus	25	27.65	2649	0

余談：GPU によって `sin()` の返り値が違う問題

シェーダーの sin() などの数学関数は 環境依存 であり、 GPU によって結果が異なる場合がある……と hashwosine の項で書きました。
これは本当なのでしょうか？

実際に試してみましょう。
一番お手軽なのが PC (NVIDIA GeForce RTX 3060 Ti) と Android (ASUS_I005DC; Adreno 660) 間での比較ですね。
Unity でコンピュートシェーダーを書いてビルドして試してみました。

ついでに、ググっていたら NVIDIA GPU における sin() はこういうコードで実装されている場合がある、とありました。

// see: https://developer.download.nvidia.com/cg/sin.html  
static float pseudoSin(float a)  
{  
    var c0 = new Vector4(0.0f, 0.5f, 1.0f, 0.0f);  
    var c1 = new Vector4(0.25f, -9.0f, 0.75f, 0.159154943091f);  
    var c2 = new Vector4(24.9808039603f, -24.9808039603f,  
                        -60.1458091736f, 60.1458091736f);  
    var c3 = new Vector4(85.4537887573f, -85.4537887573f,  
                        -64.9393539429f, 64.9393539429f);  
    var c4 = new Vector4(19.7392082214f, -19.7392082214f,  
                        -1.0f, 1.0f);  
  
    Vector3 r0, r1, r2;  
  
    r1.x = c1.w * a - c1.x;                
    r1.y = frac(r1.x);                    
    r2.x = (r1.y < c1.x) ? 1 : 0;          
    r2.y = (r1.y >= c1.y) ? 1 : 0;  
    r2.z = (r1.y >= c1.z) ? 1 : 0;  
    r2.y = Vector3.Dot(r2, new Vector3(c4.z, c4.w, c4.z));               
    r0 = new Vector3(c0.x - r1.y, c0.y - r1.y, c0.z - r1.y);  
    r0 = new Vector3(r0.x * r0.x, r0.y * r0.y, r0.z * r0.z);  
    r1 = new Vector3(c2.x * r0.x + c2.z, c2.y * r0.y + c2.w, c2.x * r0.z + c2.z);  
    r1 = new Vector3(r1.x * r0.x + c3.x, r1.y * r0.y + c3.y, r1.z * r0.z + c3.x);  
    r1 = new Vector3(r1.x * r0.x + c3.z, r1.y * r0.y + c3.w, r1.z * r0.z + c3.z);  
    r1 = new Vector3(r1.x * r0.x + c4.x, r1.y * r0.y + c4.y, r1.z * r0.z + c4.x);  
    r1 = new Vector3(r1.x * r0.x + c4.z, r1.y * r0.y + c4.w, r1.z * r0.z + c4.z);  
    return Vector3.Dot(r1, -r2);                  
}

これと値が一致するかも比べてみましょう。
0 ～ PI/2 (90°) まで、 256 等分した値をそれぞれの sin に与えて比較しました。

<GPU ネイティブの sin()> = <NVIDIA の sin()> ; <CPU の Mathf.Sin()> というフォーマットです。

Android



0 = 0 ; 0

0.006135901 = 0.006135881 ; 0.006135885

0.01227157 = 0.01227158 ; 0.01227154

0.01840678 = 0.01840681 ; 0.01840673

0.0245413 = 0.02454126 ; 0.02454123

0.03067489 = 0.03067487 ; 0.0306748

0.03680732 = 0.03680724 ; 0.03680722

0.04293814 = 0.04293823 ; 0.04293826

0.04906758 = 0.04906774 ; 0.04906768

0.05519516 = 0.05519527 ; 0.05519525

0.06132067 = 0.06132072 ; 0.06132074

0.06744387 = 0.06744397 ; 0.06744392

0.07356453 = 0.07356453 ; 0.07356457

0.07968242 = 0.07968247 ; 0.07968244

0.08579731 = 0.08579737 ; 0.08579732

0.09190897 = 0.09190899 ; 0.09190895

0.09801717 = 0.0980171 ; 0.09801714

0.1041217 = 0.1041216 ; 0.1041216

0.1102223 = 0.1102223 ; 0.1102222

0.1163187 = 0.1163187 ; 0.1163186

0.1224108 = 0.1224107 ; 0.1224107

0.128498 = 0.1284981 ; 0.1284981

0.1345806 = 0.1345807 ; 0.1345807

0.1406582 = 0.1406583 ; 0.1406582

0.1467304 = 0.1467305 ; 0.1467305

0.1527971 = 0.1527972 ; 0.1527972

0.1588581 = 0.1588582 ; 0.1588582

0.1649131 = 0.1649132 ; 0.1649131

0.1709619 = 0.1709619 ; 0.1709619

0.1770042 = 0.1770043 ; 0.1770042

0.1830399 = 0.1830398 ; 0.1830399

0.1890687 = 0.1890687 ; 0.1890687

0.1950904 = 0.1950904 ; 0.1950903

0.2011047 = 0.2011046 ; 0.2011046

0.2071115 = 0.2071114 ; 0.2071114

0.2131102 = 0.2131104 ; 0.2131103

0.2191012 = 0.2191013 ; 0.2191012

0.2250838 = 0.2250839 ; 0.2250839

0.2310581 = 0.2310581 ; 0.2310581

0.2370236 = 0.2370237 ; 0.2370236

0.2429802 = 0.2429802 ; 0.2429802

0.2489276 = 0.2489277 ; 0.2489276

0.2548656 = 0.2548657 ; 0.2548657

0.2607941 = 0.2607941 ; 0.2607941

0.2667128 = 0.2667128 ; 0.2667128

0.2726214 = 0.2726214 ; 0.2726214

0.2785197 = 0.2785197 ; 0.2785197

0.2844076 = 0.2844076 ; 0.2844076

0.2902848 = 0.2902847 ; 0.2902847

0.2961508 = 0.2961509 ; 0.2961509

0.3020059 = 0.302006 ; 0.3020059

0.3078496 = 0.3078496 ; 0.3078497

0.3136817 = 0.3136817 ; 0.3136818

0.319502 = 0.3195021 ; 0.319502

0.3253103 = 0.3253103 ; 0.3253103

0.3311063 = 0.3311063 ; 0.3311063

0.3368898 = 0.3368899 ; 0.3368899

0.3426607 = 0.3426607 ; 0.3426607

0.3484187 = 0.3484187 ; 0.3484187

0.3541636 = 0.3541635 ; 0.3541635

0.3598951 = 0.3598951 ; 0.3598951

0.3656131 = 0.365613 ; 0.365613

0.3713173 = 0.3713173 ; 0.3713172

0.3770073 = 0.3770074 ; 0.3770074

0.3826834 = 0.3826834 ; 0.3826835

0.388345 = 0.3883451 ; 0.388345

0.393992 = 0.3939921 ; 0.3939921

0.3996242 = 0.3996242 ; 0.3996242

0.4052413 = 0.4052413 ; 0.4052413

0.4108432 = 0.4108432 ; 0.4108432

0.4164295 = 0.4164296 ; 0.4164296

0.4220003 = 0.4220003 ; 0.4220003

0.4275551 = 0.4275551 ; 0.4275551

0.4330939 = 0.4330938 ; 0.4330938

0.4386163 = 0.4386162 ; 0.4386162

0.4441222 = 0.4441222 ; 0.4441222

0.4496114 = 0.4496114 ; 0.4496113

0.4550835 = 0.4550836 ; 0.4550836

0.4605387 = 0.4605387 ; 0.4605387

0.4659764 = 0.4659765 ; 0.4659765

0.4713967 = 0.4713967 ; 0.4713967

0.4767992 = 0.4767992 ; 0.4767992

0.4821838 = 0.4821838 ; 0.4821838

0.4875502 = 0.4875501 ; 0.4875502

0.4928982 = 0.4928982 ; 0.4928982

0.4982277 = 0.4982277 ; 0.4982277

0.5035384 = 0.5035384 ; 0.5035384

0.5088302 = 0.5088301 ; 0.5088302

0.5141028 = 0.5141028 ; 0.5141028

0.5193561 = 0.519356 ; 0.519356

0.5245896 = 0.5245897 ; 0.5245897

0.5298036 = 0.5298036 ; 0.5298036

0.5349975 = 0.5349976 ; 0.5349976

0.5401714 = 0.5401715 ; 0.5401715

0.545325 = 0.545325 ; 0.545325

0.550458 = 0.550458 ; 0.550458

0.5555702 = 0.5555702 ; 0.5555702

0.5606616 = 0.5606616 ; 0.5606616

0.5657318 = 0.5657318 ; 0.5657318

0.5707808 = 0.5707808 ; 0.5707808

0.5758082 = 0.5758082 ; 0.5758082

0.580814 = 0.580814 ; 0.580814

0.5857979 = 0.5857979 ; 0.5857979

0.5907598 = 0.5907597 ; 0.5907597

0.5956992 = 0.5956993 ; 0.5956993

0.6006164 = 0.6006165 ; 0.6006165

0.605511 = 0.6055111 ; 0.605511

0.6103828 = 0.6103828 ; 0.6103828

0.6152316 = 0.6152316 ; 0.6152316

0.6200572 = 0.6200572 ; 0.6200572

0.6248595 = 0.6248595 ; 0.6248595

0.6296383 = 0.6296383 ; 0.6296383

0.6343933 = 0.6343933 ; 0.6343933

0.6391245 = 0.6391245 ; 0.6391245

0.6438316 = 0.6438316 ; 0.6438316

0.6485144 = 0.6485144 ; 0.6485144

0.6531729 = 0.6531729 ; 0.6531729

0.6578068 = 0.6578067 ; 0.6578067

0.6624157 = 0.6624157 ; 0.6624158

0.6669999 = 0.6669999 ; 0.6669999

0.6715589 = 0.671559 ; 0.671559

0.6760927 = 0.6760927 ; 0.6760927

0.680601 = 0.680601 ; 0.680601

0.6850836 = 0.6850836 ; 0.6850837

0.6895405 = 0.6895406 ; 0.6895406

0.6939715 = 0.6939715 ; 0.6939715

0.6983762 = 0.6983763 ; 0.6983763

0.7027548 = 0.7027547 ; 0.7027547

0.7071067 = 0.7071068 ; 0.7071068

0.7114323 = 0.7114322 ; 0.7114322

0.7157309 = 0.7157308 ; 0.7157308

0.7200025 = 0.7200025 ; 0.7200025

0.7242471 = 0.7242471 ; 0.7242471

0.7284644 = 0.7284644 ; 0.7284644

0.7326542 = 0.7326543 ; 0.7326543

0.7368165 = 0.7368165 ; 0.7368166

0.7409511 = 0.7409511 ; 0.7409512

0.7450578 = 0.7450578 ; 0.7450578

0.7491364 = 0.7491364 ; 0.7491364

0.7531868 = 0.7531868 ; 0.7531868

0.7572088 = 0.7572088 ; 0.7572089

0.7612023 = 0.7612024 ; 0.7612024

0.7651672 = 0.7651673 ; 0.7651673

0.7691032 = 0.7691033 ; 0.7691033

0.7730104 = 0.7730105 ; 0.7730104

0.7768884 = 0.7768885 ; 0.7768885

0.7807372 = 0.7807372 ; 0.7807373

0.7845565 = 0.7845566 ; 0.7845566

0.7883464 = 0.7883464 ; 0.7883464

0.7921066 = 0.7921066 ; 0.7921066

0.7958369 = 0.7958369 ; 0.7958369

0.7995372 = 0.7995373 ; 0.7995373

0.8032075 = 0.8032075 ; 0.8032075

0.8068476 = 0.8068476 ; 0.8068476

0.8104572 = 0.8104572 ; 0.8104572

0.8140363 = 0.8140364 ; 0.8140363

0.8175848 = 0.8175848 ; 0.8175848

0.8211026 = 0.8211025 ; 0.8211026

0.8245894 = 0.8245893 ; 0.8245893

0.8280451 = 0.8280451 ; 0.8280451

0.8314696 = 0.8314696 ; 0.8314697

0.8348629 = 0.8348629 ; 0.8348629

0.8382248 = 0.8382247 ; 0.8382247

0.8415551 = 0.841555 ; 0.8415549

0.8448536 = 0.8448536 ; 0.8448536

0.8481205 = 0.8481203 ; 0.8481203

0.8513553 = 0.8513552 ; 0.8513552

0.8545578 = 0.854558 ; 0.854558

0.8577285 = 0.8577286 ; 0.8577287

0.8608669 = 0.8608669 ; 0.860867

0.8639728 = 0.8639728 ; 0.8639728

0.8670461 = 0.8670462 ; 0.8670462

0.8700869 = 0.870087 ; 0.870087

0.873095 = 0.873095 ; 0.873095

0.8760701 = 0.8760701 ; 0.8760701

0.8790122 = 0.8790122 ; 0.8790123

0.8819213 = 0.8819213 ; 0.8819213

0.8847971 = 0.8847971 ; 0.8847971

0.8876396 = 0.8876396 ; 0.8876396

0.8904487 = 0.8904487 ; 0.8904487

0.8932242 = 0.8932243 ; 0.8932243

0.8959663 = 0.8959662 ; 0.8959663

0.8986745 = 0.8986745 ; 0.8986745

0.9013488 = 0.9013488 ; 0.9013489

0.9039893 = 0.9039893 ; 0.9039893

0.9065958 = 0.9065957 ; 0.9065957

0.909168 = 0.909168 ; 0.909168

0.911706 = 0.911706 ; 0.911706

0.9142098 = 0.9142097 ; 0.9142098

0.9166791 = 0.9166791 ; 0.9166791

0.9191139 = 0.9191139 ; 0.9191139

0.9215141 = 0.921514 ; 0.921514

0.9238796 = 0.9238795 ; 0.9238795

0.9262103 = 0.9262102 ; 0.9262102

0.928506 = 0.9285061 ; 0.9285061

0.9307669 = 0.9307669 ; 0.930767

0.9329928 = 0.9329928 ; 0.9329928

0.9351835 = 0.9351835 ; 0.9351835

0.9373389 = 0.937339 ; 0.937339

0.9394591 = 0.9394592 ; 0.9394592

0.9415441 = 0.9415441 ; 0.9415441

0.9435934 = 0.9435934 ; 0.9435934

0.9456073 = 0.9456073 ; 0.9456074

0.9475855 = 0.9475856 ; 0.9475856

0.9495282 = 0.9495282 ; 0.9495282

0.951435 = 0.951435 ; 0.951435

0.953306 = 0.953306 ; 0.953306

0.9551411 = 0.9551412 ; 0.9551412

0.9569403 = 0.9569404 ; 0.9569404

0.9587035 = 0.9587035 ; 0.9587035

0.9604305 = 0.9604305 ; 0.9604306

0.9621214 = 0.9621214 ; 0.9621214

0.9637761 = 0.9637761 ; 0.9637761

0.9653945 = 0.9653944 ; 0.9653944

0.9669765 = 0.9669765 ; 0.9669765

0.9685221 = 0.9685221 ; 0.9685221

0.9700313 = 0.9700313 ; 0.9700313

0.971504 = 0.9715039 ; 0.9715039

0.97294 = 0.97294 ; 0.97294

0.9743394 = 0.9743394 ; 0.9743394

0.9757022 = 0.9757021 ; 0.9757021

0.9770282 = 0.9770281 ; 0.9770281

0.9783173 = 0.9783174 ; 0.9783174

0.9795697 = 0.9795698 ; 0.9795698

0.9807853 = 0.9807853 ; 0.9807853

0.9819639 = 0.9819639 ; 0.9819639

0.9831055 = 0.9831055 ; 0.9831055

0.9842101 = 0.9842101 ; 0.9842101

0.9852777 = 0.9852777 ; 0.9852777

0.9863081 = 0.9863081 ; 0.9863081

0.9873014 = 0.9873014 ; 0.9873014

0.9882575 = 0.9882576 ; 0.9882576

0.9891765 = 0.9891765 ; 0.9891765

0.9900582 = 0.9900582 ; 0.9900582

0.9909027 = 0.9909027 ; 0.9909027

0.9917098 = 0.9917098 ; 0.9917098

0.9924795 = 0.9924796 ; 0.9924796

0.993212 = 0.9932119 ; 0.993212

0.993907 = 0.993907 ; 0.993907

0.9945646 = 0.9945646 ; 0.9945646

0.9951847 = 0.9951847 ; 0.9951847

0.9957674 = 0.9957674 ; 0.9957674

0.9963126 = 0.9963126 ; 0.9963126

0.9968203 = 0.9968203 ; 0.9968203

0.9972904 = 0.9972904 ; 0.9972904

0.9977231 = 0.997723 ; 0.9977231

0.9981181 = 0.9981181 ; 0.9981181

0.9984756 = 0.9984756 ; 0.9984756

0.9987954 = 0.9987954 ; 0.9987954

0.9990777 = 0.9990777 ; 0.9990777

0.9993224 = 0.9993224 ; 0.9993224

0.9995294 = 0.9995294 ; 0.9995294

0.9996988 = 0.9996988 ; 0.9996988

0.9998306 = 0.9998306 ; 0.9998306

0.9999247 = 0.9999247 ; 0.9999247

0.9999812 = 0.9999812 ; 0.9999812



0 = 0 ; 0

0.006135784 = 0.006135883 ; 0.006135885

0.01227134 = 0.0122716 ; 0.01227154

0.01840647 = 0.01840678 ; 0.01840673

0.02454119 = 0.02454128 ; 0.02454123

0.0306749 = 0.03067486 ; 0.0306748

0.03680708 = 0.03680725 ; 0.03680722

0.0429382 = 0.04293833 ; 0.04293826

0.04906752 = 0.04906771 ; 0.04906768

0.05519509 = 0.05519526 ; 0.05519525

0.06132066 = 0.06132074 ; 0.06132074

0.06744379 = 0.06744399 ; 0.06744392

0.07356446 = 0.07356455 ; 0.07356457

0.07968237 = 0.07968249 ; 0.07968244

0.08579737 = 0.08579735 ; 0.08579732

0.09190872 = 0.09190897 ; 0.09190895

0.09801697 = 0.09801711 ; 0.09801714

0.1041216 = 0.1041216 ; 0.1041216

0.1102219 = 0.1102223 ; 0.1102222

0.1163183 = 0.1163187 ; 0.1163186

0.1224106 = 0.1224107 ; 0.1224107

0.1284982 = 0.1284981 ; 0.1284981

0.1345807 = 0.1345807 ; 0.1345807

0.140658 = 0.1406582 ; 0.1406582

0.1467303 = 0.1467305 ; 0.1467305

0.1527971 = 0.1527972 ; 0.1527972

0.158858 = 0.1588582 ; 0.1588582

0.1649131 = 0.1649132 ; 0.1649131

0.1709618 = 0.1709619 ; 0.1709619

0.177004 = 0.1770042 ; 0.1770042

0.1830396 = 0.1830399 ; 0.1830399

0.1890684 = 0.1890686 ; 0.1890687

0.1950903 = 0.1950904 ; 0.1950903

0.2011045 = 0.2011046 ; 0.2011046

0.2071114 = 0.2071114 ; 0.2071114

0.2131103 = 0.2131104 ; 0.2131103

0.2191012 = 0.2191012 ; 0.2191012

0.2250838 = 0.225084 ; 0.2250839

0.2310579 = 0.2310581 ; 0.2310581

0.2370234 = 0.2370237 ; 0.2370236

0.2429801 = 0.2429802 ; 0.2429802

0.2489275 = 0.2489276 ; 0.2489276

0.2548656 = 0.2548657 ; 0.2548657

0.2607938 = 0.2607942 ; 0.2607941

0.2667127 = 0.2667128 ; 0.2667128

0.2726213 = 0.2726214 ; 0.2726214

0.2785195 = 0.2785197 ; 0.2785197

0.2844074 = 0.2844076 ; 0.2844076

0.2902845 = 0.2902847 ; 0.2902847

0.2961509 = 0.2961509 ; 0.2961509

0.3020057 = 0.3020059 ; 0.3020059

0.3078495 = 0.3078496 ; 0.3078497

0.3136816 = 0.3136818 ; 0.3136818

0.3195019 = 0.3195021 ; 0.319502

0.3253103 = 0.3253103 ; 0.3253103

0.331106 = 0.3311063 ; 0.3311063

0.3368898 = 0.3368899 ; 0.3368899

0.3426606 = 0.3426608 ; 0.3426607

0.3484185 = 0.3484187 ; 0.3484187

0.3541633 = 0.3541635 ; 0.3541635

0.3598949 = 0.359895 ; 0.3598951

0.3656131 = 0.365613 ; 0.365613

0.3713171 = 0.3713172 ; 0.3713172

0.3770073 = 0.3770074 ; 0.3770074

0.3826833 = 0.3826834 ; 0.3826835

0.388345 = 0.3883451 ; 0.388345

0.3939919 = 0.393992 ; 0.3939921

0.399624 = 0.3996242 ; 0.3996242

0.4052413 = 0.4052413 ; 0.4052413

0.4108431 = 0.4108432 ; 0.4108432

0.4164295 = 0.4164296 ; 0.4164296

0.422 = 0.4220003 ; 0.4220003

0.427555 = 0.4275551 ; 0.4275551

0.4330938 = 0.4330938 ; 0.4330938

0.438616 = 0.4386162 ; 0.4386162

0.444122 = 0.4441222 ; 0.4441222

0.4496112 = 0.4496114 ; 0.4496113

0.4550835 = 0.4550836 ; 0.4550836

0.4605385 = 0.4605387 ; 0.4605387

0.4659763 = 0.4659765 ; 0.4659765

0.4713967 = 0.4713967 ; 0.4713967

0.4767992 = 0.4767992 ; 0.4767992

0.4821836 = 0.4821838 ; 0.4821838

0.4875499 = 0.4875502 ; 0.4875502

0.4928981 = 0.4928982 ; 0.4928982

0.4982276 = 0.4982277 ; 0.4982277

0.5035383 = 0.5035384 ; 0.5035384

0.5088301 = 0.5088301 ; 0.5088302

0.5141026 = 0.5141028 ; 0.5141028

0.5193559 = 0.519356 ; 0.519356

0.5245895 = 0.5245897 ; 0.5245897

0.5298036 = 0.5298036 ; 0.5298036

0.5349975 = 0.5349976 ; 0.5349976

0.5401714 = 0.5401714 ; 0.5401715

0.545325 = 0.545325 ; 0.545325

0.5504579 = 0.550458 ; 0.550458

0.5555702 = 0.5555702 ; 0.5555702

0.5606615 = 0.5606616 ; 0.5606616

0.5657318 = 0.5657318 ; 0.5657318

0.5707806 = 0.5707808 ; 0.5707808

0.5758082 = 0.5758082 ; 0.5758082

0.5808139 = 0.5808139 ; 0.580814

0.5857978 = 0.5857978 ; 0.5857979

0.5907595 = 0.5907597 ; 0.5907597

0.5956993 = 0.5956993 ; 0.5956993

0.6006165 = 0.6006165 ; 0.6006165

0.6055108 = 0.6055111 ; 0.605511

0.6103826 = 0.6103828 ; 0.6103828

0.6152316 = 0.6152316 ; 0.6152316

0.6200572 = 0.6200572 ; 0.6200572

0.6248594 = 0.6248595 ; 0.6248595

0.629638 = 0.6296383 ; 0.6296383

0.6343932 = 0.6343933 ; 0.6343933

0.6391243 = 0.6391245 ; 0.6391245

0.6438313 = 0.6438316 ; 0.6438316

0.6485143 = 0.6485144 ; 0.6485144

0.6531728 = 0.6531729 ; 0.6531729

0.6578067 = 0.6578067 ; 0.6578067

0.6624157 = 0.6624158 ; 0.6624158

0.6669998 = 0.6669999 ; 0.6669999

0.6715588 = 0.671559 ; 0.671559

0.6760926 = 0.6760927 ; 0.6760927

0.6806009 = 0.680601 ; 0.680601

0.6850834 = 0.6850837 ; 0.6850837

0.6895405 = 0.6895406 ; 0.6895406

0.6939713 = 0.6939715 ; 0.6939715

0.6983762 = 0.6983762 ; 0.6983763

0.7027546 = 0.7027547 ; 0.7027547

0.7071068 = 0.7071068 ; 0.7071068

0.7114322 = 0.7114322 ; 0.7114322

0.7157307 = 0.7157308 ; 0.7157308

0.7200024 = 0.7200025 ; 0.7200025

0.724247 = 0.7242471 ; 0.7242471

0.7284644 = 0.7284644 ; 0.7284644

0.7326542 = 0.7326543 ; 0.7326543

0.7368164 = 0.7368166 ; 0.7368166

0.7409511 = 0.7409511 ; 0.7409512

0.7450578 = 0.7450578 ; 0.7450578

0.7491363 = 0.7491364 ; 0.7491364

0.7531867 = 0.7531868 ; 0.7531868

0.7572088 = 0.7572088 ; 0.7572089

0.7612023 = 0.7612024 ; 0.7612024

0.7651672 = 0.7651673 ; 0.7651673

0.7691032 = 0.7691033 ; 0.7691033

0.7730104 = 0.7730105 ; 0.7730104

0.7768884 = 0.7768885 ; 0.7768885

0.7807373 = 0.7807372 ; 0.7807373

0.7845566 = 0.7845566 ; 0.7845566

0.7883464 = 0.7883464 ; 0.7883464

0.7921065 = 0.7921066 ; 0.7921066

0.7958369 = 0.7958369 ; 0.7958369

0.7995371 = 0.7995373 ; 0.7995373

0.8032075 = 0.8032075 ; 0.8032075

0.8068476 = 0.8068476 ; 0.8068476

0.8104572 = 0.8104572 ; 0.8104572

0.8140362 = 0.8140363 ; 0.8140363

0.8175848 = 0.8175848 ; 0.8175848

0.8211026 = 0.8211025 ; 0.8211026

0.8245893 = 0.8245893 ; 0.8245893

0.8280449 = 0.8280451 ; 0.8280451

0.8314695 = 0.8314696 ; 0.8314697

0.8348628 = 0.8348629 ; 0.8348629

0.8382246 = 0.8382247 ; 0.8382247

0.8415549 = 0.841555 ; 0.8415549

0.8448536 = 0.8448536 ; 0.8448536

0.8481203 = 0.8481203 ; 0.8481203

0.8513551 = 0.8513552 ; 0.8513552

0.8545579 = 0.854558 ; 0.854558

0.8577285 = 0.8577286 ; 0.8577287

0.860867 = 0.860867 ; 0.860867

0.8639728 = 0.8639728 ; 0.8639728

0.8670462 = 0.8670462 ; 0.8670462

0.8700869 = 0.870087 ; 0.870087

0.873095 = 0.873095 ; 0.873095

0.8760701 = 0.8760701 ; 0.8760701

0.8790122 = 0.8790122 ; 0.8790123

0.8819212 = 0.8819213 ; 0.8819213

0.8847971 = 0.8847971 ; 0.8847971

0.8876396 = 0.8876396 ; 0.8876396

0.8904487 = 0.8904487 ; 0.8904487

0.8932243 = 0.8932243 ; 0.8932243

0.8959663 = 0.8959662 ; 0.8959663

0.8986745 = 0.8986745 ; 0.8986745

0.9013488 = 0.9013488 ; 0.9013489

0.9039893 = 0.9039893 ; 0.9039893

0.9065956 = 0.9065957 ; 0.9065957

0.9091679 = 0.909168 ; 0.909168

0.9117059 = 0.911706 ; 0.911706

0.9142097 = 0.9142098 ; 0.9142098

0.9166791 = 0.9166791 ; 0.9166791

0.9191139 = 0.9191139 ; 0.9191139

0.921514 = 0.921514 ; 0.921514

0.9238795 = 0.9238795 ; 0.9238795

0.9262102 = 0.9262102 ; 0.9262102

0.928506 = 0.9285061 ; 0.9285061

0.9307669 = 0.9307669 ; 0.930767

0.9329928 = 0.9329928 ; 0.9329928

0.9351835 = 0.9351835 ; 0.9351835

0.9373391 = 0.937339 ; 0.937339

0.9394592 = 0.9394592 ; 0.9394592

0.9415441 = 0.9415441 ; 0.9415441

0.9435934 = 0.9435934 ; 0.9435934

0.9456074 = 0.9456073 ; 0.9456074

0.9475856 = 0.9475856 ; 0.9475856

0.9495282 = 0.9495282 ; 0.9495282

0.951435 = 0.951435 ; 0.951435

0.953306 = 0.953306 ; 0.953306

0.9551411 = 0.9551412 ; 0.9551412

0.9569404 = 0.9569404 ; 0.9569404

0.9587035 = 0.9587035 ; 0.9587035

0.9604305 = 0.9604305 ; 0.9604306

0.9621213 = 0.9621214 ; 0.9621214

0.9637761 = 0.9637761 ; 0.9637761

0.9653944 = 0.9653944 ; 0.9653944

0.9669764 = 0.9669765 ; 0.9669765

0.968522 = 0.9685221 ; 0.9685221

0.9700311 = 0.9700313 ; 0.9700313

0.9715039 = 0.9715039 ; 0.9715039

0.97294 = 0.97294 ; 0.97294

0.9743394 = 0.9743394 ; 0.9743394

0.9757022 = 0.9757021 ; 0.9757021

0.9770282 = 0.9770281 ; 0.9770281

0.9783174 = 0.9783174 ; 0.9783174

0.9795697 = 0.9795698 ; 0.9795698

0.9807853 = 0.9807853 ; 0.9807853

0.9819639 = 0.9819639 ; 0.9819639

0.9831055 = 0.9831055 ; 0.9831055

0.9842101 = 0.9842101 ; 0.9842101

0.9852776 = 0.9852777 ; 0.9852777

0.9863081 = 0.9863081 ; 0.9863081

0.9873014 = 0.9873014 ; 0.9873014

0.9882575 = 0.9882576 ; 0.9882576

0.9891765 = 0.9891765 ; 0.9891765

0.9900582 = 0.9900582 ; 0.9900582

0.9909026 = 0.9909027 ; 0.9909027

0.9917097 = 0.9917098 ; 0.9917098

0.9924796 = 0.9924796 ; 0.9924796

0.993212 = 0.9932119 ; 0.993212

0.993907 = 0.993907 ; 0.993907

0.9945646 = 0.9945646 ; 0.9945646

0.9951847 = 0.9951847 ; 0.9951847

0.9957675 = 0.9957674 ; 0.9957674

0.9963127 = 0.9963126 ; 0.9963126

0.9968203 = 0.9968203 ; 0.9968203

0.9972905 = 0.9972904 ; 0.9972904

0.997723 = 0.997723 ; 0.997723

0.9981181 = 0.9981181 ; 0.9981181

0.9984756 = 0.9984756 ; 0.9984756

0.9987954 = 0.9987954 ; 0.9987954

0.9990778 = 0.9990777 ; 0.9990777

0.9993225 = 0.9993224 ; 0.9993224

0.9995294 = 0.9995294 ; 0.9995294

0.9996988 = 0.9996988 ; 0.9996988

0.9998307 = 0.9998306 ; 0.9998306

0.9999248 = 0.9999247 ; 0.9999247

0.9999812 = 0.9999812 ; 0.9999812

ここからわかることは、

CPU と GPU で sin() の値が違うこと
実際に GPU によって sin() の値が変わること
CPU の Mathf.Sin() は環境によらず同じ値を返していること
NVIDIA の sin() 式はどちらとも微妙に違うこと
- 計算誤差か、 FMA か、そもそも方式が違う (テーブルルックアップとか) か
- もちろん環境依存なので、たまたま私のデバイスは違うっぽい、ということだけわかる

です。

それから、 GPU の sin() に環境間での再現性を求めてはいけない、ということもわかりますね。
例えばノイズを地形生成とかに使うと同じシードでも微妙に再現できない、みたいな可能性もあります。

おわりに

本稿では、シェーダー乱数の比較検討を行いました。
みなさん frac(sin(...)) のやつはご存知だったかもしれませんが、知らない関数も多かったのではないでしょうか。
本記事が発見のきっかけになったなら幸いです。

また、高速で頑健なシェーダー乱数 IbukiHash を提案しました。
ぜひ使ってあげてください。ライセンスは CC0 です。
どうせ乱数なんてコピペするものなのですから、これを機に frac(sin(...)) のやつから切り替えてみたりしていただければと思います！

IbukiHash

GLSL での実装が欲しい方は、以下の Shadertoy を参照してください。

https://www.shadertoy.com/view/XX3yRn

*1:Jarzynski, Mark, and Marc Olano. "Hash functions for gpu rendering." UMBC Computer Science and Electrical Engineering Department (2020). https://www.jcgt.org/published/0009/03/02/

*2:https://mrl.cs.nyu.edu/~perlin/noise/

*3:Gustavson, Stefan, and Ian McEwan. "Tiling simplex noise and flow noise in two and three dimensions." J Comput Graph Tech 11.1 (2022). https://jcgt.org/published/0011/01/02/paper.pdf

*4:Valdenegro-Toro, Matias, and Hector Pincheira. "Implementing Noise with Hash functions for Graphics Processing Units." arXiv preprint arXiv:1903.12270 (2019). https://arxiv.org/pdf/1903.12270

*5:A Hash Function for Hash Table Lookup, http://www.burtleburtle.net/bob/hash/doobs.html

*6:Blum, Lenore, Manuel Blum, and Mike Shub. "A simple unpredictable pseudo-random number generator." SIAM Journal on computing 15.2 (1986): 364-383.

*7:"google/cityhash: Automatically exported from code.google.com/p/cityhash", https://github.com/google/cityhash

*8:Schechter, Hagit, and Robert Bridson. "Evolving sub-grid turbulence for smoke animation." Proceedings of the 2008 ACM SIGGRAPH/Eurographics symposium on Computer animation. 2008. https://www.cs.ubc.ca/~rbridson/docs/schechter-sca08-turbulence.pdf

*9:Hash without Sine, https://www.shadertoy.com/view/4djSRW

*10:Howes, Lee, and David Thomas. "Efficient random number generation and application using CUDA." GPU gems 3 (2007): 805-830.

*11:Jorge Jimenez – Next Generation Post Processing in Call of Duty: Advanced Warfare, https://www.iryoku.com/next-generation-post-processing-in-call-of-duty-advanced-warfare/

*12:Integer Hash11, https://www.shadertoy.com/view/llGSzw

*13:Integer Hash33, https://www.shadertoy.com/view/XlXcW4

*14:Integer Hash21, https://www.shadertoy.com/view/4tXyWN

*15:Jones, David. "Good practice in (pseudo) random number generation for bioinformatics applications." URL http://www.cs.ucl.ac.uk/staff/d.jones/GoodPracticeRNG.pdf (2010).

*16:https://en.wikipedia.org/wiki/KISS_(algorithm)

*17:Press, William H., et al. Numerical recipes. Cambridge University Press, London, England, 1988.

*18:Rivest, Ronald. "RFC1321: The MD5 message-digest algorithm." (1992).

*19:Tzeng, Stanley, and Li-Yi Wei. "Parallel white noise generation on a GPU via cryptographic hash." Proceedings of the 2008 symposium on Interactive 3D graphics and games. 2008. ソースコードは https://github.com/1iyiwei/parallel-white-noise/blob/master/src/code/OpenGL/md5.frag

*20:O’neill, Melissa E. "PCG: A family of simple fast space-efficient statistically good algorithms for random number generation." ACM Transactions on Mathematical Software (2014).

*21:Jarzynski, Mark, and Marc Olano. "Hash functions for gpu rendering." UMBC Computer Science and Electrical Engineering Department (2020).

*22:Jarzynski, Mark, and Marc Olano. "Hash functions for gpu rendering." UMBC Computer Science and Electrical Engineering Department (2020).

*23:Press, William H. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.

*24:Hash functions., http://www.azillionmonkeys.com/qed/hash.html

*25:Wheeler, David J., and Roger M. Needham. "TEA, a tiny encryption algorithm." Fast Software Encryption: Second International Workshop Leuven, Belgium, December 14–16, 1994 Proceedings 2. Springer Berlin Heidelberg, 1995.

*26:Integer Hash Function, http://web.archive.org/web/20071223173210/http://www.concentric.net/~Ttwang/tech/inthash.htm

*27:Marsaglia, George. "Xorshift rngs." Journal of Statistical software 8 (2003): 1-6.

*28:Marsaglia, George. "Xorshift rngs." Journal of Statistical software 8 (2003): 1-6.

*29:Cyan4973/xxHash: Extremely fast non-cryptographic hash algorithm, https://github.com/Cyan4973/xxHash

*30:wangyi-fudan/wyhash: The FASTEST QUALITY hash function, random number generators (PRNG) and hash map., https://github.com/wangyi-fudan/wyhash

*31:Prospecting for Hash Functions, https://nullprogram.com/blog/2018/07/31/

*32:Prospecting for Hash Functions, https://nullprogram.com/blog/2018/07/31/

*33:Creating reliable and efficient random hash for use in GPU shaders | by Lumi | Oct, 2024 | Medium, https://medium.com/@lumi_/creating-reliable-and-efficient-random-hash-for-use-in-gpu-shaders-fe5b5c9b6b72

*34:A fast and simple 32bit floating point hash function | briansharpe, https://briansharpe.wordpress.com/2011/11/15/a-fast-and-simple-32bit-floating-point-hash-function/

*35:Salmon, John K., et al. "Parallel random numbers: as easy as 1, 2, 3." Proceedings of 2011 international conference for high performance computing, networking, storage and analysis. 2011.

*36:Rijmen, Vincent, and Joan Daemen. "Advanced encryption standard." Proceedings of federal information processing standards publications, national institute of standards and technology 19 (2001): 22.

*37:heptaplex-collapse noise, https://www.shadertoy.com/view/ms3czf

*38:Steele Jr, Guy L., and Sebastiano Vigna. "Computationally easy, spectrally good multipliers for congruential pseudorandom number generators." Software: Practice and Experience 52.2 (2022): 443-458.

2024-12-05

NVIDIA Nsight Graphics を使ってみる

はじめに

みなさんはシェーダーのプロファイリングをしたくなったとき、どうやっていますか？
FPS 測定？でもオーバヘッドが大きいし、具体的にどこが重いのかわからない……

みたいな悩みを、 NVIDIA Nsight Graphics で解決できるかもしれません。

ほかにいい手段をご存知の方はぜひご連絡ください……

手順

まずはこちらから NVIDIA Nsight Graphics をインストールします。

インストールの待ち時間に、 Unity プロジェクトの用意をします。
対象となるシェーダーを書いて、マテリアルを作成、それを適用した Plane なり Cube なりを配置します。
カメラに写ってさえいれば OK です。むしろ余計なものがあると解析の邪魔になるので最低限で大丈夫です。

シェーダーですが、以下の #pragma が必要となるので、シェーダーコードに書いておきます。

// #pragma fragment の下あたりに入れる  
  
// デバッグシンボルを含める  
#pragma enable_d3d11_debug_symbols  
  
// shader model 6.0 を使う  
#pragma use_dxc

これがないと詳細な解析ができません。

なお、 #pragma enable_d3d11_debug_symbols はパフォーマンスに若干の悪影響を与える可能性があります。 ( →参照 )
デバッグが終わったら消しておくのが無難です。
なんで use_dxc で SM6.0 が使えるのかはこちらを参照してください。

ビルドしたアプリを解析することになるので、アプリビルドの準備をします。

まずは Project Settings / Player / Other Settings / Rendering から Auto Graphics API for Windows のチェックを外し、一番上に Direct3D12 が来るようにします。

デフォルトの Direct3D11 だと解析に失敗します。
Vulkan でも動作しますが、デコンパイル結果などが異なり、見慣れない感じになります。
そちらに慣れている方は Vulkan のほうがいいかもしれませんが、私は DirectX をお勧めします。

そうしたら、 Windows 用にアプリをビルドします。
Development Build で大丈夫です。

ビルドしたら、 NVIDIA Nsight Graphics を 管理者権限で 実行します。
管理者権限でないと一部データの取得に失敗します。

起動したら、 Start Activity... を選んでください。

Application Executable: にさっきビルドしたアプリのパスを入れます。
そうしたら右下の Launch Frame Debugger を押します。

ここで Steam を起動していると Interfering Processes Detected と表示され、 Steam のアプリが干渉するとか言われるのですが、気にせず Keep を押して大丈夫です。

うまくいけば、ビルドしたアプリが起動します。

ここで Direct3D11on12 is unsupported みたいな警告が出ますが、とりあえずは無視しても大丈夫そうです。

アプリ上で F11 キーを押すとキャプチャできるので、押します。

すると、 Unity の Profiler みたいな画面が開きますので、解析をします。

まずは Events タブから、描画していそうなイベントを探します。
色がついているので、それを見ながら探してください。
私の場合は DrawIndexedInstanced というイベントがそれでした。検索フィルタもあるのでそれを活用すると見つけやすいかもです。

クリックすると、 API Inspector タブが更新されます。
PS (Pixel Shader) ボタンを押して Profile Shaders ↗ を押します。

すると Shader Profiler が開きます。
ここでは命令数やその内訳、かかった時間などを見ることができます。
カーソルを合わせると詳細な説明がありますので適宜見てください。

ここで、 Correlation 行に警告が出ていた場合はたぶんシェーダーコードの設定を書き忘れています。
#pragma をちゃんと書いたか確認してください。

詳しい説明は本家の説明書を参照してください。

ストールの種類と説明 (MIO Throttle とは何ぞや、みたいなの) はこちら。

処理にかかった時間は下段の Events / GPU Duration のほか、 Events タブの GPU ms から確認できます。
Hot Spots からは重い処理がどこかを突き止めることができます。

動的にシェーダーを編集することもできる (!) らしいですが、一度編集してしまうとソースコードと照らし合わせた解析ができなくなるっぽいので、あまり有効ではなさそうです。

おわりに

シェーダーのプロファイリングなんもわからん状態だったのでとても感動して記事にしました。
とりあえず導入してみた感じですが、プロファイリングが詳しいしいろんな機能があるようなのでいろいろ触ってみたいと思います。

参考文献

User Guide — NsightGraphics 2024.3 documentation
公式のユーザーガイドです。

はじめに

モンゴメリ乗算とは

モンゴメリリダクション

乗法逆元

拡張ユークリッドの互除法

Jeffrey Hurchalla 法 ( が 2 冪の場合)

モンゴメリ乗算

モンゴメリ加算・減算

の計算

冪剰余

発展問題

が偶数のときのモンゴメリ乗算

中国の剰余定理

Garner のアルゴリズム

モンゴメリ乗算中に 2 で割る

乗法逆元の応用

Extended Binary GCD

Hybrid Extended GCD

乗法逆元を使った倍数判定

実践: Miller-Rabin 素数判定法

余談

その他の応用

おわりに

はじめに

手順

テクスチャを描く

Unity へインポート

FontAsset 生成

テクスチャを上書き

文字情報の置換

試してみる

おわりに

はじめに

いつ 2 冪以外の配列を返すのか

.NET 9 環境

Unity 環境

余談

おわりに

はじめに

IBufferWriter<T> とは

ArrayBufferWriter<T>

PipeWriter

アロケーション問題

IBufferWriter<T> を実装してみよう

ArrayPoolBufferWriter<T> の実装

バッファ拡大時のコピーを回避する

ReadOnlySequence<T> とは

ArrayPoolSegmentedBufferWriter<T> の実装

余談

おわりに

はじめに

ランダムソート

ソートなので遅い・均等ではない

危険なランダムソート

Fisher-Yates shuffle

擬似乱数まわりの最適化

擬似乱数生成器の選定

乱数を絞りつくす

を求める

64 bit 擬似乱数 を生成する

を の範囲に均等に変換できるように調整する

を分割してインデックス 2 ～ 20 ぶんの乱数を得る

Fisher-Yates shuffle を行う

次の を求めて、必要なぶん繰り返す

結果

小手先の高速化

Next() メソッドの (手動) インライン展開

上限を削る

タプルでの交換をやめる

SIMD 化

配列アクセス時の範囲チェック削除

のキャッシュ

Fisher-Yates shuffle の誤った実装例

MergeShuffle

Feistel 構造を利用したシャッフル

おわりに

はじめに

シェーダー乱数の定義

PractRand テストについて

シェーダー乱数一覧

Jeffrey Hurchalla 法 ( $m$ が 2 冪の場合)

$R^{2}$ の計算

冪剰余 $a^{b} \bmod m$

$m$ が偶数のときのモンゴメリ乗算

`IBufferWriter<T>` とは

`ArrayBufferWriter<T>`

`PipeWriter`

`IBufferWriter<T>` を実装してみよう

`ArrayPoolBufferWriter<T>` の実装

`ReadOnlySequence<T>` とは

`ArrayPoolSegmentedBufferWriter<T>` の実装

$n$ を求める

64 bit 擬似乱数 $r$ を生成する

$r$ を $0 \le s \lt n$ の範囲に均等に変換できるように調整する

$r$ を分割してインデックス 2 ～ 20 ぶんの乱数を得る

次の $n$ を求めて、必要なぶん繰り返す

`Next()` メソッドの (手動) インライン展開

$n$ のキャッシュ

余談：GPU によって `sin()` の返り値が違う問題