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

はじめに

シェーダー (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 216

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

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 216

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

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 216

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

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 221

しかし、テストには失敗しています。
こういう目視では確認できないような品質がチェックできるのが 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 216

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

加えて、定数を 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 216

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

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 241

ここにきてようやく 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 240

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 216

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

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 216

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

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 216

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

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 217

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

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 242

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

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 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 216

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

ただし、 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 218

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

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 216

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

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 > 238

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

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 241

テストも問題なく、 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 238

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

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 227

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

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 242

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

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 225

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

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 242

それなりに速く、かつ 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 216

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

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 228

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

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 219

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

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 221

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

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 216

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

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 235

健闘しましたが、テスト結果は 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 216

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

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 227

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

license xxHash Library
Copyright (c) 2012-2021 Yann Collet
All rights reserved.

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 242

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 242

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

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 239

なぜか 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 216

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


ちなみに、 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 216

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

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 242

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

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 > 235

あまりにも重い (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 219

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

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 241

設計思想

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

PC


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.