2026-05-04

DUALSHOCK 4 のジャイロ・タッチパッドを Unity で使う

はじめに

Unity の InputSystem はクソ。 (挨拶)
ソースだけ欲しい方はこれを読み飛ばしてソースの項を見てください。

Unity の新しい InputSystem は、様々なデバイスを扱うことができます。たぶん。

だがしかし、デフォルト状態だと制約が結構あります。
特に今私が使っている DUALSHOCK 4 (PS4 用純正コントローラ) を PC で利用する場合を考えると、

ジャイロ (角速度センサー)：アクセス不可
加速度センサー：アクセス不可
タッチパッド：クリックのみ検出可、タッチ位置は検出不可
振動・ライトバー変更：有線 (USB) 接続時のみ可、無線 (Bluetooth) 接続時不可

と、惨憺たる有様です。本来の力を発揮できていません。

したがって、本記事ではこれらの要素を検出し、 Unity (PC) から使えるようにすることを目標とします。

検証環境は以下の通りです。

Unity 6000.3.0b4
InputSystem 1.19.0
DUALSHOCK4 CUH-ZCT2J

なお、 InputSystem の基本的な使い方やインストール方法などについてはマスターしているものとします。
それらについては他の記事を参照してください。

注意点ですが、 InputSystem は 絶望的にバグが多いことで私の中で有名 なので、最新版 (現時点で 1.19.0) を使うことを推奨します。
現在のデフォルトインストールは確か 1.14.x なので、手動更新が必要です。古いとなんかのバグを踏んでスクリプトが動かなくなったりしました。

ソース

先にソースを出したほうが分かりやすいかもしれないので先に貼ります。
適当に貼りつければ動くと思います。 MIT です。
注意点としては、 unsafe コードを含むので (ポインタ操作が不可避) 、 unsafe を許可する asmdef に置く必要があります。

Dualshock4Controller.cs

詰まった点など

無線だと、そもそもジャイロとかのデータが取れない

大前提として、データが降ってこないことにはどうしようもありません。
取れるデータの詳細は下記サイトに詳しいです：

USB 接続時: https://www.psdevwiki.com/ps4/DS4-USB
Bluetooth 接続時: https://www.psdevwiki.com/ps4/DS4-BT

問題は Bluetooth 接続時で、初期状態だと HID INPUT Reports 0x01 でデータが送られてくるのですが、これが簡易バージョンであり、普通のボタン入力しかデータを得ることができません。
(ジャイロ・加速度・タッチ入力などは詳細バージョンのデータにのみ含まれています。)

これを詳細バージョンの HID INPUT Reports 0x11 にしたいわけなのですが、

This report is sent once the GET REPORT FEATURE 0x02 is received.

ということで、 FEATURE Report 0x02 を投げる必要があります。
だがしかし、我らが InputSystem には FEATURE Report を投げる機能は存在しません (たぶん) 。したがって手詰まりです。

一応 HidD_GetFeature などを駆使すれば FEATURE Report を投げることは可能ですが、 Windows 限定になりますし取り回しも悪い (extern が必要) ですし、できれば Unity 内で完結したいですよね。

なのですが、 1 日ほど実験した結果、どうやら HID OUTPUT Report 0x11 を投げることでも詳細モードに切り替わる ことが分かりました。

HID OUTPUT Report 0x11 は、振動やライトの状態を変更させるためのレポートです。
なので、接続を検知した瞬間 Dualshock4ControllerDevice.OnAdded() に HID OUTPUT Report 0x11 を投げてやれば、晴れて無線でもジャイロなどのデータが取れるようになる、というわけです。
HID OUTPUT Report は InputDevice.ExecuteCommand() で送信することができます。

なお、このとき対象となる構造体のサイズを明示的に指定する必要がある ([StructLayout] の Size パラメータを必ず書く) ことに注意してください。 (一敗)
そうしないとパディングか何かでサイズがずれて送信に失敗します。

ちなみに、振動とライトの値を変更するレポートは有線 (0x05) と無線 (0x11) で微妙にフォーマットが違います。注意してください。

ジャイロ・加速度・タッチの座標が勝手に補正される

これは InputSystem の問題です。
タッチ座標とかを State に定義するときは以下のようにすると思うのですが、

// 抜粋  
  
[InputControl(name = "touch1/id", layout = "Integer", format = "BIT", offset = 0, bit = 0, sizeInBits = 7)]  
[InputControl(name = "touch1/flag", layout = "DiscreteButton", format = "BIT", offset = 0, bit = 7, sizeInBits = 1, parameters = "minValue=0,maxValue=0,nullValue=1")]  
[InputControl(name = "touch1", layout = "Vector2")]  
[InputControl(name = "touch1/x", layout = "Axis", format = "BIT", offset = 0, bit = 8, sizeInBits = 12, minValue = 0, maxValue = 1919)]  
[InputControl(name = "touch1/y", layout = "Axis", format = "BIT", offset = 0, bit = 20, sizeInBits = 12, minValue = 0, maxValue = 942)]  
[FieldOffset(34)]  
public readonly uint touch1;

このとき、 touch1/x などの値が 勝手に 0 ～ 1 や -1 ～ 1 に補正されてしまいます 。
ジャイロはまだともかく、タッチパッドの座標でもこれをやられてしまい、大変取り扱いづらいです。

それなら layout = "Integer" を使えばいいじゃん、と思われるかもしれませんが、それだと Vector2 や Vector3 の軸として扱うことができません。
あくまで layout = "Axis" のままやる必要があります。

じゃぁどうするのかというと、 /x や /y に parameters = "scale,scaleFactor=4095" を追加します。
4095 は 12 bit の最大値です。ビット数に合わせて変更してください。
内部的に最大値で割って補正しているようなので、逆補正をかけて生の値を取り出すイメージです。

あとついでに。
ここで touch1/flag があるのですが、これの生の値は

0 = 押している
1 = 押していない

と普通のボタンとは逆になっています。
これを普通のボタンに変換するには、安直には parameters = "invert" で行けそうな気がしますが、これを使うと -1 になってしまうのでダメです。
正解は上のように DiscreteButton を駆使することです。これは dpad でも使われていますね。

ジャイロの単位がわからん

未だにわかりません。
センサー (Bosch BMI055) のデータシートを見た限りでは deg/s っぽいのですが、正しい係数が分かりません。
分かる方は連絡してください。よろしくお願いします…

なお、加速度のほうは 1G (重力加速度; 9.80665 m/s²) で 16384 ぐらいだったので分かりやすかったです。

おわりに

Unity の InputSystem はクソ。 (結言)

でも拡張性があることだけは認めてやらんこともないです。
いや最初から対応していればこんなに苦労することはなかったのでやっぱりクソ。

ところで、この機能が完成するまで 4 日ぐらいかかったのですが、これだけ頑張っても対応できるのは DUALSHOCK 4 だけで他のデバイスは全然なんですよね……
Rewired とか買ったほうがいいのかもしれません、よ。知らんけど。

2026-01-15

浮動小数点数の《正確な》乱数生成って難しいね

はじめに

一様分布な浮動小数点数の乱数 (例えば、 $[0.12, 3.45)$ の範囲の均等な乱数) が欲しいときがしばしばあります。
皆さんもおそらく random.NextDouble() * (max - min) + min みたいなコードを書いたことがあるかと思います。
しかし、このコードが実は複数の問題をはらんでいることをご存知でしょうか？
本稿では、そういった問題とその解決策、またパフォーマンスとの兼ね合いについて紹介したいと思います。

本稿の擬似乱数生成器は、このシグネチャで実装されていると仮定します。
(C# / .NET の Random とは違うことにご留意ください。)

public class Random  
{  
    // [0, 2^64) の範囲の一様な整数乱数を取得  
    public ulong Next();  
}

また、浮動小数点数の型は基本的に double (IEEE 754 の binary64) を前提とします。
float や Half でも同様の議論は適用可能ですので適宜読み替えてください。

本稿では、以降 double のビットパターンや内部表現を理解されている前提で話を進めていきます。
要所要所で説明は入れていますが、よくわからないという方は Wikipedia などをご参照ください。

本稿でよく取り扱う形式の早見表を以下に示します。

形式名	C# 型名	指数部bit	仮数部bit	指数部ゲタ	ラウンドトリップ桁数
binary16	`Half`	5	10	15	5
binary32	`float`	8	23	127	9
binary64	`double`	11	52	1023	17

種類	指数部	仮数部
0	0	0
非正規化数	0	>0
正規化数	1～(emax-1)	any
∞	emax	0
NaN	emax	>0

$[0.0, 1.0)$ の範囲において

一番標準となる $[0.0, 1.0)$ の範囲において一様乱数を生成する場合を考えます。

念のため、 $\lbrack a, b \rbrack$ は閉区間 (その値を含む) 、 $(a, b)$ は開区間 (その値を含まない) ことを表します。要するに 0.0 <= x && x < 1.0 です。
数学的には閉区間でも開区間でも確率は同じ (ある特定の値を取る確率は 0 なので) ですが、工学的にはある種の離散分布となるので確率が違ってきます。

hole さんのスライドに詳しいので、まだ見ていない方は一度ご覧になるとよいかと思います。

$2^{64}$ で割る手法

// 2^64  
const double POW_2_64 = 2.0 * (1ul << 63);  
  
return random.Next() / POW_2_64;  
  
// -- or --  
  
return random.Next() * (1.0 / POW_2_64);

前述のスライドで「除算法」とされていた手法です。
除算版と乗算版は本質的に同一です。 C# では今のところ最適化されませんが、 C 言語 (clang) なら掛け算のほうに統一されます。

開区間にならない (`1.0` が出る)

random.Next() が 0xffffffff_ffffffff を返した場合に 1.0 が返ります。

「誤差だよ誤差！」と思った方もいらっしゃるかもしれませんが、例えば指数分布に変換するコードを書いた場合問題になりえます。

double r = random.Next() * (1.0 / POW_2_64);  
  
// 1.0 - r とすることで (0.0, 1.0] の範囲にする意図  
// しかし実際は r == 1.0 になりうるので、  
// -log(0.0) == ∞ を返しうる  
return -Math.Log(1.0 - r);

乱数が絡む関係上、これを見逃すと「ユーザ側でごくまれに計算エラーが発生するが、開発機では全く再現できない」ということになりがちです。
未然に防ぐためにも、必ず定義域内に収まるようにすべきです。

表現可能な値すべてが出ない

例えば、この手法で出現しうる 0 より大きい最小の数は 5.421010862427522e-20 (0x3bf0000000000000) です。

一方で、本来の 0 より大きい最小の数は、正規化数の範囲では 2.225073858507201e-308 (0x0010000000000000) 、非正規化数も含めれば 4.940656458412465e-324 (0x0000000000000001) です。桁の桁が違いますね。

これまた「誤差だよ誤差！」と思われたかもしれませんが、上と同じように指数分布に変換した場合に違いが顕著となります。

// 正規化数の範囲での最大値は 708.3964185322641  
// 非正規化数も含めれば 744.4400719213812  
// しかし、この手法では 44.3614195558365 にしかならない  
return -Math.Log(r);

ちょっと遅い

詳しいパフォーマンスについては後で触れますが、この手法はちょっと遅くなります。

というのも、 long → double に数値変換する intrinsics はあります (cvtsi2sd) が、 ulong → double では AVX-512 でない限り (_mm256_cvtepu64_pd) ありません。
ということで、前処理の命令が挟まって微妙に遅くなります。

$2^{53}$ で割る手法

return (random.Next() >> 11) * (1.0 / (1ul << 53));

前述の $2^{64}$ で割る手法の欠点に対応したのが、この $2^{53}$ で割る手法です。

なぜ $2^{53}$ なのかというと、 double 型の仮数部の精度が 53 bit だからです。
実際の仮数部は 52 bit 幅ですが、暗黙的に存在する先頭の 1 (ケチ表現) の分で 53 bit 精度になります。

開区間になる

random.Next() が 0xffffffff_ffffffff を返した場合に 0.9999999999999999 (0x3fefffffffffffff) が返ります。
したがって、 1.0 以上になることはありません。

表現可能な値すべてが出ない

この問題は引き続き起こります。

この手法では、 0 より大きい最小の値は 1.110223024625157e-16 (0x3ca0000000000000) となります。

ちょっと速い

シフトがある分遅くなるのでは？と思われるかもしれませんが、 $2^{64}$ で割る手法に比べて速くなります。
具体的にアセンブリを見てみましょう。 BenchmarkDotNet の DisassemblyDiagnoser を使ってみます。

BenchmarkDotNet v0.15.8, Windows 11 (10.0.26100.7171/24H2/2024Update/HudsonValley)  
12th Gen Intel Core i7-12700F 2.10GHz, 1 CPU, 20 logical and 12 physical cores  
.NET SDK 10.0.100  
  [Host]     : .NET 10.0.0 (10.0.0, 10.0.25.52411), X64 RyuJIT x86-64-v3  
  Job-OJRDUL : .NET 10.0.0 (10.0.0, 10.0.25.52411), X64 RyuJIT x86-64-v3

まずは $2^{64}$ のほう。

; FPTest.MultiplyFullWidth()  
;         return Rng.Next() * (1.0 / (2.0 * (1ul << 63)));  
;         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  
       sub       rsp,28  
       mov       rcx,[rcx+8]  
       cmp       [rcx],ecx  
       call      qword ptr [7FFA9699FC48]; Seiran.Next()  
       vxorps    xmm0,xmm0,xmm0  
       mov       rdx,rax  
       shr       rdx,1  
       mov       ecx,eax  
       and       ecx,1  
       or        rcx,rdx  
       test      rax,rax  
       cmovns    rcx,rax  
       vcvtsi2sd xmm0,xmm0,rcx  
       jns       short M00_L00  
       vaddsd    xmm0,xmm0,xmm0  
M00_L00:  
       vmulsd    xmm0,xmm0,qword ptr [7FFA9666EC08]  
       add       rsp,28  
       ret  
; Total bytes of code 65

次に、 $2^{53}$ のほう。

; FPTest.MultiplyAppropriateWidth()  
;         return (Rng.Next() >> 11) * (1.0 / (1ul << 53));  
;         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  
       sub       rsp,28  
       mov       rcx,[rcx+8]  
       cmp       [rcx],ecx  
       call      qword ptr [7FFA9698FC48]; Seiran.Next()  
       shr       rax,0B  
       vxorps    xmm0,xmm0,xmm0  
       mov       rdx,rax  
       shr       rdx,1  
       mov       ecx,eax  
       and       ecx,1  
       or        rcx,rdx  
       test      rax,rax  
       cmovns    rcx,rax  
       vcvtsi2sd xmm0,xmm0,rcx  
       jns       short M00_L00  
       vaddsd    xmm0,xmm0,xmm0  
M00_L00:  
       vmulsd    xmm0,xmm0,qword ptr [7FFA966614F8]  
       add       rsp,28  
       ret  
; Total bytes of code 69

うーん、シフト命令が増えているので遅くなっているように見えますね。
しかし、実際のベンチマークでは、

Method	Mean	Error	StdDev	Code Size
MultiplyFullWidth	5.808 ns	0.1336 ns	0.1538 ns	115 B
MultiplyAppropriateWidth	1.928 ns	0.0613 ns	0.0776 ns	119 B

このように明らかに $2^{53}$ のほうが速くなります。

同じようなコードを C で書いて、

double mul64(uint64_t num) {  
    return (double)(num) * (1.0 / (2.0 * ((uint64_t)1 << 63)));  
}  
  
double mul53(uint64_t num) {  
    return (double)(num >> 11) * (1.0 / (1.0 * ((uint64_t)1 << 53)));  
}

Compiler Explorer にかけてみると (x86-64 clang 21.1.0; -O3)、

.LCPI0_0:  
        .long   1127219200  
        .long   1160773632  
        .long   0  
        .long   0  
.LCPI0_1:  
        .quad   0x4330000000000000  
        .quad   0x4530000000000000  
.LCPI0_2:  
        .quad   0x3bf0000000000000  
mul64:  
        movq    xmm1, rdi  
        punpckldq       xmm1, xmmword ptr [rip + .LCPI0_0]  
        subpd   xmm1, xmmword ptr [rip + .LCPI0_1]  
        movapd  xmm0, xmm1  
        unpckhpd        xmm0, xmm1  
        addsd   xmm0, xmm1  
        mulsd   xmm0, qword ptr [rip + .LCPI0_2]  
        ret  
  
.LCPI1_0:  
        .quad   0x3ca0000000000000  
mul53:  
        shr     rdi, 11  
        cvtsi2sd        xmm0, rdi  
        mulsd   xmm0, qword ptr [rip + .LCPI1_0]  
        ret

命令数が結構減っていることが分かるかと思います。

余談：閉区間にするには

閉区間、つまり $\lbrack 0.0, 1.0 \rbrack$ の範囲にしたい場合、 $2^{53}$ の代わりに $2^{53} - 1$ で割る手法があります。
この手法は Unity の Random.value (これは float 型ですが) などで使われています。
ただし、表現可能な値すべてが出ない問題は引き続き発生しますので注意が必要です。

余談： .NET での実装

現時点での .NET の Random.Shared.NextDouble() はこの手法で実装されています。

なんでわざわざ Shared を付けたのかというと、 seed を指定した Random は別の実装になっており、そちらは Next() * (1.0 / int.MaxValue) (Next() は int.MaxValue より小さい正の乱数を返す) という感じになっているためです。

仮数部ビットパターン法

ulong bin = (0x3fful << 52) | (random.Next() >> 12);  
return BitConverter.UInt64BitsToDouble(bin) - 1.0;

仮数部のビットパターンを直接埋める方式です。
BitConverter.UInt64BitsToDouble は、通常のキャスト (double)x ではなくビットパターンをそのまま移行するメソッドです。

random.Next() >> 12 が仮数部 (mantissa) です。したがって 52 bit の精度を持ちます。
0x3fful << 52 が符号部 (sign) と指数部 (exponent) です。ここでは、 $+2^{0} = 1$ になるように指定しています。
したがって、 bin は $+2^{0} \times 1.mantissa$ つまり $[1.0, 2.0)$ になります。
そこから $1.0$ を引くことで $[0.0, 1.0)$ の乱数を得ます。かしこい！

開区間になる

random.Next() が 0xffffffff_ffffffff を返した場合に 0.9999999999999998 (0x3feffffffffffffe) が返ります。
したがって、 1.0 以上になることはありません。

表現可能な値すべてが出ない

そもそも上で述べた最大値の時点で、表現可能な最大値ではないです。

0 より大きい最小の値は 2.220446049250313e-16 (0x3cb0000000000000) です。
これは、 $2^{53}$ で割る手法の 2 倍ぐらい大きい (精度が悪い) です。

また問題として、この手法によって得られる値は 仮数部の最下位ビットが常に 0 になります 。
したがって 1 bit 分精度が落ちているとも言えます。

最下位ビットが 0 になる理由を説明しましょう。
まず、仮数部をランダムな値で埋めた時点では $1.mantissa \times 2^{0}$ になりますね。
ここから $1.0$ を引くと $0.mantissa \times 2^{0}$ になります。ただこれだと小数点の左隣が $1$ でないので、正規化を行います。
そうなると、少なくとも 1 ビット以上シフトしないと小数点の左隣が $1$ になりませんので、右端に 0 がシフト数と同じぶんだけ (1 つ以上) 挿入されることになります。
したがって、右端 (最下位) ビットが常に 0 になる、というわけです。

パフォーマンス

$2^{-53}$ を掛ける手法との比較です。

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.928 ns	0.0613 ns	0.0776 ns	119 B
BitOp	1.661 ns	0.0581 ns	0.0622 ns	101 B

若干速いです。

余談: dSFMT について

double の出力に特化した擬似乱数生成器 dSFMT は、この手法を応用して $[0.0, 1.0)$ の乱数を生成しています。
内部状態を仮数部 52 bit に保持して、残り 12 bit (0x3ff) を定数としたメモリパターン (0x3ffxxxxx_xxxxxxxx) をそのまま保持することで、ダイレクトに $[1.0, 2.0)$ の乱数を吐き出す仕組みです。

ただ、そのせいで $[0.0, 1.0)$ 乱数の仮数部の最下位ビットが常に 0 になる問題が避けられません。むずかしい。

余談: 別の区間への応用

Inigo Quilez :: computer graphics, mathematics, shaders, fractals, demoscene and more

例えば、 $[-1.0, +1.0)$ の区間の値を得るために、同様の手法が適用できます。

ulong bin = (0x400ul << 52) | (random.Next() >> 12);    // [2.0, 4.0)  
return BitConverter.UInt64BitsToDouble(bin) - 3.0;      // [-1.0, 1.0)

exponent+mantissa 法

ulong rng = Rng.Next();  
int entropy = 64;  
  
int exponent = 0x3fe;  
  
do  
{  
    if ((rng & 1) != 0)  
    {  
        break;  
    }  
    else  
    {  
        exponent--;  
  
        entropy--;  
        rng >>= 1;  
        if (entropy == 0)  
        {  
            rng = Rng.Next();  
            entropy = 64;  
        }  
    }  
} while (exponent != 0);  
  
ulong mantissa = (entropy >= 52 ? rng : Rng.Next()) & ((1ul << 52) - 1);  
  
return BitConverter.UInt64BitsToDouble((ulong)exponent << 52 | mantissa);

前述した hole さんのスライドで紹介されていた方式です。

まず、指数部 (exponent) を決定します。初期値は 0x3fe つまり $2^{-1}$ (0.5) です。
そこから 1 bit 乱数を生成し、真なら確定 (ループ終了)、偽ならデクリメント (つまり 1/2) してループを続行します。
つまり、これは $p=0.5$ の幾何分布に基づく (50% の確率で成功する試行が成功するまでの試行回数に基づく) といえます。

次に仮数部 (mantissa) を 52 bit 乱数で埋めます。

あとはそのビットパターンで double 型を組み立てれば完成です。

開区間になる

最大の値は 0.9999999999999999 (0x3fefffffffffffff) ですので、開区間になります。

表現可能なすべての値が出る

ここにきてようやく、表現可能なすべての値が出せるアルゴリズムが出ました。

ただし注記しておきたいのは、 厳密に正しい確率で出るわけではありません 。
そうするには後述する Downey の補正を入れる必要があります。

遅い

コードを読んでいただければわかると思うのですが、 1 ビットずつ処理する必要があり遅いです。
また、仕方のないことですが乱数を複数消費する可能性があります。

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.928 ns	0.0613 ns	0.0776 ns	119 B
ExponentMantissa_Naive	9.996 ns	0.1565 ns	0.1307 ns	185 B

改造してみる

ulong rng = Rng.Next();  
int entropy = 64;  
  
int exponent = 0x3fe;  
  
do  
{  
    if (rng != 0)  
    {  
        int clz = BitOperations.LeadingZeroCount(rng);  
        exponent -= clz;  
        entropy -= clz;  
        break;  
    }  
    else  
    {  
        exponent -= 64;  
        rng = Rng.Next();  
    }  
} while (exponent >= 0);  
  
ulong mantissa = (entropy >= 52 ? rng : Rng.Next()) & ((1ul << 52) - 1);  
  
if (exponent < 0)  
{  
    mantissa >>= -exponent;  
    exponent = 0;  
}  
  
return BitConverter.UInt64BitsToDouble((ulong)exponent << 52 | mantissa);

おそらく 1 bit ずつ処理しているのがネックなので、そこを改造してみましょう。

bit が 0 の間ループするということは、 64 bit 乱数の先頭からの 0 の数を数えられればよさそうですね。
BitOperations.LeadingZeroCount というぴったりのメソッドがあるのでそれを使いましょう。

「先頭」は 1ul << 63 のほうです。具体的には、 BitOperations.LeadingZeroCount(1ul) == 63 になります。
毎回 TrailingZeroCount とどっちがどっちなのか分からなくなります……

仮数部の処理はほぼ同じですが、一応指数部が負になった時のことも考えています。まぁほとんどあり得ませんが……

Method	Mean	Error	StdDev	Code Size
ExponentMantissa_Naive	10.055 ns	0.1603 ns	0.1339 ns	185 B
ExponentMantissa_Clz	1.944 ns	0.0616 ns	0.0546 ns	196 B

だいぶ高速化しました。乗算法とほぼ同じぐらいです。

Abseil 実装

Google が作った C++ 効率化ライブラリである Abseil による実装 (GenerateRealFromBits) です。

ulong rng = Rng.Next();  
  
if (rng == 0)  
{  
    return 0;  
}  
  
int exp = 0x3fe;  
int clz = BitOperations.LeadingZeroCount(rng);  
rng <<= clz;  
exp -= clz;  
rng >>= 11;  
  
ulong result = (ulong)exp << 52 | (rng & ((1ul << 52) - 1));  
return BitConverter.UInt64BitsToDouble(result);

本質的には exponent+mantissa 法に近いです。

まず、 0 を引いた場合はそのまま 0 を返します。これは近似的な実装ですね ( $1/2^{64}$ の確率で発生)。

続いて、 BitOperations.LeadingZeroCount で先頭から続く 0 の数を数えます。

そして、指数部 exp から clz を引きます。要するに exponent+mantissa 法で 1 bit ずつやっていた処理をまとめてやっている感じですね。

rng をシフトしているのは、桁合わせのような感じです。
まず rng <<= clz によって、必ず rng の先頭ビット (1ul << 63 のところ) が 1 になります。
そこから rng >>= 11 することで、末尾から 53 bit 目が 1 になります。これが浮動小数点数におけるケチ表現の 1 に当たるわけですね。

あとは残りの rng を仮数部として、ビットパターンを組み立てて返します。

開区間になる

最大の値は 0.9999999999999999 (0x3fefffffffffffff) ですので、開区間になります。

表現可能なすべての値は出ない

0 より大きい最小の値は 5.421010862427522e-20 (0x3bf0000000000000) です。
とはいえ、 $2^{53}$ で割る手法に比べると大体 $2^{12}$ ぐらい分解能が高いです。

すべての値は出ないことと引き換えに、乱数消費は 64 bit 1 個で固定です。
while/for ループなどが使いにくい環境 (シェーダーとか?) で便利かもしれません。

パフォーマンス

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.927 ns	0.0317 ns	0.0265 ns	119 B
Abseil	1.668 ns	0.0384 ns	0.0340 ns	139 B

乗算法より高速です。

Downey による補正

Downey 氏が発表した手法です。 *1

ulong rng = Rng.Next();  
int entropy = 64;  
  
int exponent;  
for (exponent = 0x3fe; exponent > 0; exponent--)  
{  
    if ((rng & 1) != 0)  
    {  
        break;  
    }  
  
    entropy--;  
    rng >>= 1;  
    if (entropy == 0)  
    {  
        rng = Rng.Next();  
        entropy = 64;  
    }  
}  
  
ulong mantissa = Rng.Next() & ((1ul << 52) - 1);  
  
// ポイント！仮数部 0 のとき 50% の確率で exponent を増やす  
if (mantissa == 0 && (rng & 1) != 0)  
{  
    exponent++;  
}  
  
return BitConverter.UInt64BitsToDouble((ulong)exponent << 52 | mantissa);

前半の exponent を決めるあたりは exponent+mantissa 法と同じですね。
ポイントとなるのは mantissa == 0 だった時の処理で、 50% の確率で exponent を増やしています。
どうしてこんなことをしているのかというと、そのままだと境界部分の確率が正しくないからです。

これは Downey の補正を入れる前の、それぞれの値が得られる範囲を表した図です。
色のついたエリアの値は ▲ の値に丸められる……という図です。
一見正しそうに見えますが、本来はこのようになるべきです。

特定の値の区間は、その値から始まって次の値で終わる (切り捨てのイメージ) のではなく、特定の値を中心としてその周囲に半分ずつ ± した範囲にあります (四捨五入のイメージ)。
そして、 $2^{n}$ ちょうどの値については、前半 ( $\lt 2^{n}$ ) と後半 ( $\gt 2^{n}$ ) で範囲の大きさが異なります。一般的な区間 (オレンジ) の 3/4 のサイズですね。
この範囲の違いをどう処理するかというと、 50% の確率でそのまま (つまり後半部) 、もう 50% の確率で次の指数部の区間 (つまり $2^{n+1}$ の前半部) に送る、とします。スマートですね！

開区間にならない

指数部で $[0.5, 1.0)$ の区間 (0x3fe) を引いたうえで補正がかかる条件を満たすと 1.0 を返します。
そうなる確率は $1/2^{54}$ です。無視できるほど低いとは言えません。

そうなった場合、 2 つの手法が考えられます。

対策 1 : `clamp` する

要するに 0.9999999999999999 を返す手法です。
当然ながら確率が偏ることになるので、できればしたくありませんね。
ただ、固定時間で実行したい場合や、乱数消費の個数を固定したい場合には役立つかもしれません。

対策 2 : 再抽選する

結果を破棄して再度乱数を生成する手法です。
多少時間はかかりますが確率が偏らないので、個人的にはこちらをお勧めします。
(無限ループの可能性は擬似乱数生成器が壊れていない限り無視できるでしょう。)

表現可能なすべての値が「均等に」出現する

この手法最大の利点です。
加えて、補正をかける処理自体は別の手法と組み合わせて適用できます。

遅い

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.927 ns	0.0317 ns	0.0265 ns	119 B
Downey	9.336 ns	0.0735 ns	0.0652 ns	196 B

1 bit ずつの処理があるのでそれはそう。
このあたりは LeadingZeroCount などと組み合わせれば改善できます。この手法の本質は補正なので問題ではないです。

「さいきょうの」手法

こちらのブログで紹介されていた、「ぼくのかんがえたさいきょうの手法」です。

int exp = 0x3fe;  
ulong frac;  
  
while (true)  
{  
    ulong i = Rng.Next();  
    int l = BitOperations.TrailingZeroCount(i) + 1;  
  
    exp -= 64 - l;  
    if (exp <= 0)  
    {  
        frac = i;  
        exp = 0;  
        break;  
    }  
  
    if (l > 52)  
    {  
        frac = i >> (l - 52 - 1);  
        break;  
    }  
    else if (l > 0)  
    {  
        int s = 52 - l + 1;  
        frac = i << s;  
  
        i = Rng.Next();  
        frac |= i & ((1ul << s) - 1);  
        break;  
    }  
}  
  
return BitConverter.UInt64BitsToDouble((ulong)exp << 52 | (frac & ((1ul << 52) - 1)));

本質的には exponent+mantissa 法ですね。
TrailingZeroCount でまとめて判定しているため、安直な exponent+mantissa 法よりは速そうです。

開区間になる

最大の値は 0.9999999999999999 (0x3fefffffffffffff) ですので、開区間になります。

表現可能なすべての値が出る

表現可能なすべての値が出現します。

ただし Downey の補正がないため、厳密に正しい確率で、というわけではありません。
「さいきょう」はなかなか険しいですね。

パフォーマンス

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.927 ns	0.0317 ns	0.0265 ns	119 B
Saikyou	2.399 ns	0.0281 ns	0.0234 ns	205 B

微妙な速度です。精度が向上していることを考えると十分に速いとも言えそうです。

`random_real`

.NET の Random でも使われている xoshiro256** の発案者 Sebastiano Vigna 氏のサイトで紹介されている手法です。
執筆者さん自体は Taylor R Campbell 氏だとか。

int exponent = -64;  
ulong significand;  
int shift;  
  
// exponent+mantissa 法と同じイメージ  
while ((significand = Rng.Next()) == 0)  
{  
    exponent -= 64;  
  
    if (exponent < -1074)  
    {  
        return 0;  
    }  
}  
  
shift = BitOperations.LeadingZeroCount(significand);  
if (shift != 0)  
{  
    exponent -= shift;  
    significand <<= shift;  
    significand |= Rng.Next() >> (64 - shift);  
}  
  
// ポイント: 仮数部の最下位を 1 にする  
significand |= 1;  
  
// significand * 2^{exponent}  
return Math.ScaleB((double)significand, exponent);

手法の説明をしていきましょう。

無限のランダムなビット列からランダムに抽出し、それを $\lbrack 0.0, 1.0 \rbrack$ の実数の 2 進展開の小数部として解釈します。
要するに、乱数ビット列 01010011 が得られたら 2 進小数 0.01010011 にする、みたいな話ですね。
このあたりは exponent+mantissa 法と同じイメージです。

ビット列が最初の 1 から 53 bit めに達すると、 double の仮数部の精度の関係上これ以上ビットを追加できなくなりますね。
この時、丸めが行われます。

一般に、浮動小数点数の丸めは最近接丸め（偶数）です。

最も近い値に丸める (それはそう)
完全に二値の中間だった場合は、仮数部の最下位ビットが 0 (→偶数) になるほうを採用する

という丸め方式です。
この 2 行目が曲者で、通常の数学においては偏りを消す方向に働くのですが、今回は偏りを生む原因になります。

実例を見ていきましょう。

for (int offset = 0; offset <= 20; offset++)  
{  
    long x = 1L << 54;  
    x += offset;  
  
    Console.WriteLine($"{x:x16} {BitConverter.DoubleToUInt64Bits((double)x):x16} {(double)x:g17}");  
}

// int64_t (hex)     double (hex)    double (value)     
// ------------------------------------------------  
0040000000000000 4350000000000000 18014398509481984  
0040000000000001 4350000000000000 18014398509481984  
0040000000000002 4350000000000000 18014398509481984  
0040000000000003 4350000000000001 18014398509481988  
0040000000000004 4350000000000001 18014398509481988  
0040000000000005 4350000000000001 18014398509481988  
0040000000000006 4350000000000002 18014398509481992  
0040000000000007 4350000000000002 18014398509481992  
0040000000000008 4350000000000002 18014398509481992  
0040000000000009 4350000000000002 18014398509481992  
004000000000000a 4350000000000002 18014398509481992  
004000000000000b 4350000000000003 18014398509481996  
004000000000000c 4350000000000003 18014398509481996  
004000000000000d 4350000000000003 18014398509481996  
004000000000000e 4350000000000004 18014398509482000  
004000000000000f 4350000000000004 18014398509482000  
0040000000000010 4350000000000004 18014398509482000  
0040000000000011 4350000000000004 18014398509482000  
0040000000000012 4350000000000004 18014398509482000  
0040000000000013 4350000000000005 18014398509482004  
0040000000000014 4350000000000005 18014398509482004

$2^{54} + offset$ を double に変換してみました。
double が表現できるのは 53 bit までなので、下位 2 bit は丸めの影響を受けます。

とりあえず、 double(hex) の列について、最下位が 2 の行は 5 つあるのに対して 3 の行は 3 つしかない (つまり偏っている) ということを観察してみてください。
その理由を詳しく調べてみることにします。

ここで、 int64_t (hex) の列の 0040000000000002 の行に注目してみてください。
これを 2 進数表示すると、

              |       53 bit      |   
0b 0000 0000 0100 0000 .... 0000 0010

といった感じで、ちょうど仮数部から 0b10 がはみ出す形になっています。
53 bit の仮数部を整数部分と考えると、 0b10 は小数部分にあたります。
小数部分として解釈するとこれは 0.5 になりますね。 0.5 は「完全に二値の中間である」値ですので、偶数丸めの対象となります。
このとき一番近い二値は 4350000000000000 と 4350000000000001 なのですが、このうち偶数 (最下位ビットが 0) のほうを選択するので、最終的な結果は 4350000000000000 になります。

このとき、下位ビットと仮数部との関係は下表のようになります。

下位ビット	処理
0b00	仮数部+0
0b01	仮数部+0
0b10	偶数丸め
0b11	仮数部+1

問題となるのは、今観察したように偶数の行が増えて奇数の行が減る、つまり偏ってしまうことです。
したがって今回は偶数丸めを防ぎたいわけなのですが、 C# から (というか、 C/C++ 以外のほとんどの言語では) 丸めモードを変更するのは簡単ではありません。
じゃあどう対策するのかというと、最下位ビットを 1 にセットします。
そうすることによって「完全に二値の中間である」状態にならなくすることで、偶数丸めを防ぐことができます。

この手法に確率的な正当性があるのかというと、あります。
「完全に二値の中間となる」、つまり小数部分が 0b10000000... ちょうどになる確率は $\lim_{n \to \infty} 1/2^{n} = 0$ です。
本来絶対に起こりえない事象であることを考えると、潰してしまっても問題ない、というわけです。

さて、それでは最下位ビットを 1 にセットしたと考えて、前の表をフィルタしたものを見てみましょう：

// int64_t (hex)     double (hex)    double (value)     
// ------------------------------------------------  
0040000000000001 4350000000000000 18014398509481984  
0040000000000003 4350000000000001 18014398509481988  
0040000000000005 4350000000000001 18014398509481988  
0040000000000007 4350000000000002 18014398509481992  
0040000000000009 4350000000000002 18014398509481992  
004000000000000b 4350000000000003 18014398509481996  
004000000000000d 4350000000000003 18014398509481996  
004000000000000f 4350000000000004 18014398509482000  
0040000000000011 4350000000000004 18014398509482000  
0040000000000013 4350000000000005 18014398509482004

double(hex) の列で最下位が 2 の行は 2 つあるのに対して 3 の行も 2 つと、確率が揃ったことを観察してみてください。

浮動小数点数の丸めに詳しい方なら「要するに Sticky Bit みたいな感じ？」と思われたかもしれません。その通りです。

開区間にならない

(double)significand のところで 1.0 になる可能性があります。

もっとも、この現象についてはソースコードのコメントに記載されていて、開区間にしたければ再抽選せよ、とあります。

表現可能なすべての値が「ほぼ均等に」出現する

上述した工夫により、表現可能なすべての値が出現します。

なお、 0.0 と 1.0 の出現可能性が低い問題はあるそうです。
(ほかの値は両側から丸められるのに対し、端にあるこれらの値は片側からしか丸められないためです。)

もっとも、 0.0 は出現確率が $2^{-1074}$ ですし、 1.0 に至っては棄却対象なので、ほぼ問題にならないでしょう。

パフォーマンス

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.927 ns	0.0317 ns	0.0265 ns	119 B
RandomReal	7.146 ns	0.1659 ns	0.1775 ns	366 B

ちょっと遅いですね。
調べてみた感じ、 (double)significand が遅いらしく ((double)Rng.Next() だけでも 5.6 ns ぐらいかかる) 、小手先で多少の短縮をしたとしても大幅な改善は難しそうです。

MarcDense

ulong r = Rng.Next();  
int lzc = BitOperations.LeadingZeroCount(r);  
  
if (lzc <= 64 - 52)  
{  
    int exponent = 0x3fe - lzc;  
    ulong mantissa = r & ((1ul << 52) - 1);  
  
    return BitConverter.UInt64BitsToDouble((ulong)exponent << 52 | mantissa);  
}  
  
return (0.5 / (1ul << 63)) * r;

Higher density uniform floats

もともとは float 用だったコードを double に移植したものです。
正確さと速さのバランスを取ったような実装です。

実装の説明

LeadingZeroCount の結果が 12 以下なら仮数部ビットパターン法 (っぽいもの) 、それ以上なら $2^{-64}$ 乗算法になります。

設計者さん曰く、そもそも理想を目指したとしても「多次元に均等分布する擬似乱数生成器」がない限り絶対に発生しえないパターンが生じてしまうので、そういった極限の状況はケアしないことで現実的な実用性を重視した、とのことです。

開区間になる

最大の値は 0.9999999999999999 (0x3fefffffffffffff) ですので、開区間になります。

表現可能なすべての値は出ない

0 より大きい最小の値は 5.4210108624275222e-20 (要するに $2^{-64}$ ) です。

なお、 $[2^{-12}, 1.0)$ の区間では表現可能なすべての値が出現します。
$[0.0, 2^{-12})$ の区間では、密ではない (すべての値は出ない) ものの、等間隔に出現します。

もともとの float 用のコードでは、 $[2^{-40}, 1.0)$ という広い区間で表現可能なすべての値が出現します。
double 用に引き延ばしたので微妙になっているところはあるかと思います。

すべての値は出ないかわりに、 Abseil の実装と同様に乱数消費が 1 個固定です。

パフォーマンス

Method	Mean	Error	StdDev	Median	Code Size
MultiplyAppropriateWidth	1.873 ns	0.0109 ns	0.0091 ns	1.871 ns	119 B
MarcDense	1.649 ns	0.0530 ns	0.0470 ns	1.633 ns	168 B

それなりに速いです。

Perfect

ulong expRange = 52ul << 52;  
ulong one = BitConverter.DoubleToUInt64Bits(1.0);  
  
int tailBits = 0;  
  
ulong mantissa = Rng.Next() >> 12;  
while (mantissa == 0)  
{  
    one -= expRange;  
  
    if (one < expRange)  
    {  
        mantissa = Rng.Next() >> 29 << 17;  
        tailBits = 17;  
        break;  
    }  
  
    mantissa = Rng.Next() >> 12;  
}  
  
double num = BitConverter.UInt64BitsToDouble(one | mantissa) - BitConverter.UInt64BitsToDouble(one);  
ulong numAsInt = BitConverter.DoubleToUInt64Bits(num);  
  
  
tailBits += (int)((one >> 52) - (numAsInt >> 52));  
if (tailBits > 52)  
{  
    tailBits = 52;  
}  
  
  
ulong tail = Rng.Next() >> -tailBits;  
  
return BitConverter.UInt64BitsToDouble(numAsInt + tail);

Perfect Random Floating-Point Numbers | Speculative Branches

GitHub のほうに論文があるので、興味のある方は確認してみてください。

実装の説明

まずは $\lbrack 1.0, 2.0)$ の乱数から $1.0$ を引いて $\lbrack 0.0, 1.0)$ の乱数を得るのと同じ感じで num を求めます。
で、そのままだと「仮数部ビットパターン法」で前述したように最下位ビットが 0 になる問題が生じるので、乱数 tail を生成して穴埋めする、みたいな感じです。

開区間になる

最大の値は 0.99999999999999989 (0x3fefffffffffffff) ですので、開区間になります。

表現可能なすべての値が出る

穴埋めのおかげで表現可能なすべての値が出ます。
ただし、 Downey の補正がないため、厳密に均等な確率ではありません。

ところで、ブログと論文をよく読むと、丸めを Round To Nearest (RTN; 最近接丸め) にするモードがあります。
これを使うと Downey の補正と同じ計算ができそうです。

ブログのほうの実装は誤っているようです。論文の Algorithm 3 の実装を参考にしました。

// tail の定義の行を以下に置換する  
ulong tail = ((Rng.Next() >> (-tailBits - 1)) + 1) >> 1;

これで名実ともに Perfect になりますね。
なお、この変更を入れると例によって開区間ではなくなる ( $1.0$ が出る) ので注意が必要です。

パフォーマンス

Method	Mean	Error	StdDev	Code Size
MultiplyAppropriateWidth	1.673 ns	0.0440 ns	0.0390 ns	119 B
Perfect	2.596 ns	0.0426 ns	0.0356 ns	228 B
Perfect_Rtn	2.579 ns	0.0583 ns	0.0487 ns	240 B

乗算法に比べると少し遅いですが、完璧な分布が得られることを考えると速いほうとも言えます。
また、切り捨て (Perfect) と最近接丸め (Perfect_Rtn) では速度差がほとんどありませんでした。

まとめ

私見も交えた各手法のまとめはこんな感じです。星の数は多いほど良いです。

手法	開区間	表現可能精度	パフォーマンス
$2^{64}$ で割る	✖️	🌟🌟🌟	🌟🌟🌟
$2^{-64}$ を掛ける	✖️	🌟🌟🌟	🌟🌟🌟
$2^{-53}$ を掛ける	✅	🌟🌟	🌟🌟🌟🌟🌟
仮数部ビットパターン法	✅	🌟	🌟🌟🌟🌟🌟
ExponentMantissa_Naive	✅	🌟🌟🌟🌟	🌟
ExponentMantissa_Clz	✅	🌟🌟🌟🌟	🌟🌟🌟🌟🌟
Abseil	✅	🌟🌟🌟	🌟🌟🌟🌟🌟
Downey	✖️	🌟🌟🌟🌟🌟	🌟
「さいきょう」	✅	🌟🌟🌟🌟	🌟🌟🌟🌟
RandomReal	✖️	🌟🌟🌟🌟🌟	🌟🌟
MarcDense	✅	🌟🌟🌟	🌟🌟🌟🌟🌟
Perfect	✅	🌟🌟🌟🌟	🌟🌟🌟🌟
Perfect_Rtn	✖️	🌟🌟🌟🌟🌟	🌟🌟🌟🌟

詳細なパフォーマンスのまとめはこんな感じです。

Method	Mean	Error	StdDev	Code Size
DivideFullWidth	4.994 ns	0.0473 ns	0.0442 ns	115 B
MultiplyFullWidth	4.915 ns	0.0680 ns	0.0636 ns	115 B
MultiplyAppropriateWidth	1.673 ns	0.0440 ns	0.0390 ns	119 B
BitOp	1.462 ns	0.0433 ns	0.0384 ns	101 B
ExponentMantissa_Naive	8.729 ns	0.1143 ns	0.1013 ns	185 B
ExponentMantissa_Clz	1.687 ns	0.0398 ns	0.0332 ns	196 B
Abseil	1.469 ns	0.0243 ns	0.0203 ns	139 B
Downey	8.262 ns	0.0672 ns	0.0562 ns	196 B
Saikyou	2.662 ns	0.0734 ns	0.0651 ns	205 B
RandomReal	6.294 ns	0.1448 ns	0.1284 ns	366 B
MarcDense	1.471 ns	0.0315 ns	0.0279 ns	168 B
Perfect	2.596 ns	0.0426 ns	0.0356 ns	228 B
Perfect_Rtn	2.579 ns	0.0583 ns	0.0487 ns	240 B

速い順に並べ替えるとこんな感じです。

速度的には仮数部ビットパターン法や Abseil ・ MarcDense、精度的には Downey や RandomReal などがよさそうです。
トータルの星の数 (バランス) 的には ExponentMantissa_Clz や Perfect_Rtn がよさげです。

$[min, max)$ の範囲において

さて、ここまでは基本となる $[0.0, 1.0)$ の範囲でしたが、それでも十分ややこしいことは理解いただけたかと思います。
以降は任意の範囲 $[min, max)$ を扱います。もっとひどいことになります。

なお、引数の例外処理は省略するものとします。
つまり、常に $min \lt max$ でかつ $min, max$ ともに有限の数を表す (±∞ や NaN ではない) ものとします。

普通のやつ

// NextDouble() returns [0.0, 1.0)  
return min + Rng.NextDouble() * (max - min);

普通の人はこう書くのではないかと思います。
お察しかと思いますが、これは複数の問題をはらんでいます。これで解決するなら本記事は存在しません。

開区間にならない

丸めの問題で max と同じ値が出る可能性があります。

例えば、 min == 3.082039625533209e-09, max == 0.00024414275416177517, rng == 0.99999999999999989 のときに 0.00024414275416177517 が返ります。

NextDouble() が 1.0 を返すバグを含んでいた場合、 max を超える可能性すらあります。

実例として、 min == 5.84856512116677e-08, max == 0.00010836065264006332 のときに 0.00010836065264006333 が返ります。

この問題は実際に C++ で発生しています。
std::generate_canonical ( $\lbrack 0.0, 1.0)$ の範囲の乱数を生成する関数) が 1.0 を返すバグと std::uniform_real_distribution ( $\lbrack min, max )$ の範囲の乱数を生成する関数) がこの式で実装されていることの合わせ技によって、 $max$ が含まれるどころか $max$ より大きい値が返る可能性まであります。

望ましくない結果を返す可能性がある

ここでの「望ましくない結果」というのは、「正規化数を入力したにもかかわらず、 ±∞ や NaN などの値が返る状態」を指します。

例えば、 min == -1e+308, max == 1e+308 の場合は max - min の計算でオーバーフローするため ∞ を返します。

確率が均等にならない

出やすい値や出にくい値、最悪の場合全く出ない値が生じる可能性があります。

簡単な集計プログラムを書いてみましょう。
全集計を簡単にするために Half 型でやります。
乱数の代わりに $[0.0, 1.0)$ の範囲にあるすべての値をこの式に入れて、その結果の分布を集計するプログラムです。
まずは、 $[0.0, \pi)$ の範囲、つまり result = r * Half.Pi でやってみましょう。

Half min = Half.Zero;  
Half max = Half.Pi;  
  
Dictionary<Half, int> bucket = new();  
for (Half i = min; i < max; i = Half.BitIncrement(i))  
{  
    bucket[i] = 0;  
}  
  
  
for (Half i = Half.Zero; i < Half.One; i = Half.BitIncrement(i))  
{  
    Half r = min + (max - min) * i;  
  
    if (bucket.TryGetValue(r, out var c))  
    {  
        bucket[r] = c + 1;  
    }  
    else  
    {  
        bucket[r] = 1;  
    }  
}  
  
foreach (var pair in bucket.OrderBy(pair => pair.Key))  
{  
    Console.WriteLine($"{pair.Key:g}\t{pair.Value}");  
}

結果をグラフにするとこんな感じです。

理想的にはすべての数が 1 回ずつ均等に出現するはずですが、ところどころ 0 回や 2 回出現している数値があります。

次に、 $[\pi, 2 \pi)$ の範囲でやってみましょう。 result = Half.Pi + r * Half.Pi です。

縦軸は対数です。
もともと一様分布ではない (小さいほど多い) ので左側が天元突破しているのは置いておいて、特に右半分に注目してみてください。
やはり 1 ではない (理想よりも多い) 確率で出現している値がたくさんあります。

原因として考えられるのは、まずは鳩の巣原理です。
Half 型において $[0.0, 1.0)$ の範囲には 15360 個の数値がありますが、 $[\pi, 2 \pi)$ の範囲には 1024 個の数値しかありません。
これをマッピングするとなると、どうしても同じ値になってしまう値が発生します。

加えて、浮動小数点数の丸めによる問題もあるかと思います。
$[0.0, \pi)$ の範囲には 16968 個の数値があるため、 $[0.0, 1.0)$ からきちんとマッピングされていれば 0 ～ 1 個で済むはずが、 2 個になっている箇所が結構あります。
丸めによって同じ値にマッピングされてしまった値があるのだと思います。

以上は精度の低い Half 型での問題ということで強調されている面もあるかと思いますが、本質的に float や double でも同様の問題が発生します。

パフォーマンス

今後のことを考えて、複数の min, max の値についてパフォーマンスを測ることにします。
今回は、以下の五種類について測定を行います。

$-\pi, \pi$ : 異符号で指数が同じ
$-0, 2^{52}+1$ : 異符号で、終点の仮数部が 1 (エッジケース)
$0, \pi$ : 同符号で 0 から
$1, \pi$ : 同符号で指数が 1 異なる
$e, \pi$ : 同符号で指数が同じ

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B

どれも同じぐらいですね。それはそう。

標準ライブラリにおける実装について

こちらの文献 *2 によると、 C++, Fortran, Java, Julia, Matlab, Octave, Python, R, Rust, Scilab, Swift の標準ライブラリはこの式で実装されています。

なお、我らが C# にはそもそも NextDouble(min, max) が存在しません。かなしい。
ただ一応、こちらの文献にはこの式を使うとよい、との記述があります。

Lerp 式

double r = Rng.NextDouble();  
return min * (1.0 - r) + max * r;

Lerp (線形補間) のような式です。数学的には「普通のやつ」と同一ですが、実際にはいくつか差異があります。

開区間にならない

例えば、 min == 0.32975114157467966, max == 0.34467819389987275, rng == 0.99999999999999989 のときに 0.34467819389987275 が返ります。

望ましくない結果を返さない

「普通のやつ」とは違い、オーバーフローによって望ましくない結果を返すことはありません。
（もちろん直接 +∞ とかを食わせた場合は別ですが）

確率が均等にならない

「普通のやつ」と同様に $[0.0, 1.0)$ の乱数を引き延ばしている以上、出やすい値や出にくい値、全く出ない値が生じる可能性があります。

単調増加にならないことがある

驚くべきことに、この式では r が増えたときに必ずしも結果が増加するとは限りません。

具体例を挙げると、 min == 0.31236300804549622, max == 0.66538509490050768, rA == 0.7628541362561102, rB == 0.76285413625611032 のとき lerpA == 0.58166736719260559 > lerpB == 0.58166736719260548 となります。

ただ、乱数生成においてこの性質が問題になるかどうかは微妙なところです。線形補間として利用する場合には注意しておきたいポイントかもですね。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
Lerp	-3.141592653589793	3.141592653589793	1.710 ns	0.0192 ns	0.0179 ns	1.705 ns	163 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
Lerp	0	3.141592653589793	1.930 ns	0.0216 ns	0.0192 ns	1.927 ns	163 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
Lerp	-0	4503599627370497	1.756 ns	0.0361 ns	0.0337 ns	1.748 ns	163 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
Lerp	1	3.141592653589793	1.757 ns	0.0459 ns	0.0430 ns	1.746 ns	163 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
Lerp	2.718281828459045	3.141592653589793	1.716 ns	0.0143 ns	0.0127 ns	1.717 ns	163 B

ほとんど変わりませんが、若干速い傾向にあるぐらいです。

FMA Lerp

double r = rng.NextDouble();  
return Math.FusedMultiplyAdd(r, max - min, min);  
  
// ---- OR ----  
  
return Math.FusedMultiplyAdd(r, max, (1.0 - r) * min);  
  
// ---- OR ----  
  
return Math.FusedMultiplyAdd(r, max, Math.FusedMultiplyAdd(-r, min, min));

「普通のやつ」や Lerp 式に FMA を適用したバリエーションです。

計算結果の誤差が多少小さくなる可能性はありますが、本質的には同じなので、同様の問題が引き続き発生します。

これらの細かい違いについては、こちらの記事で調査しています。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
NormalFma	-3.141592653589793	3.141592653589793	1.943 ns	0.0263 ns	0.0246 ns	1.936 ns	152 B
LerpFma1	-3.141592653589793	3.141592653589793	1.732 ns	0.0227 ns	0.0212 ns	1.736 ns	168 B
LerpFma2	-3.141592653589793	3.141592653589793	1.944 ns	0.0211 ns	0.0197 ns	1.943 ns	161 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
NormalFma	0	3.141592653589793	1.727 ns	0.0271 ns	0.0254 ns	1.725 ns	152 B
LerpFma1	0	3.141592653589793	1.926 ns	0.0219 ns	0.0194 ns	1.927 ns	168 B
LerpFma2	0	3.141592653589793	1.766 ns	0.0613 ns	0.0656 ns	1.758 ns	161 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
NormalFma	-0	4503599627370497	1.918 ns	0.0078 ns	0.0073 ns	1.918 ns	152 B
LerpFma1	-0	4503599627370497	1.736 ns	0.0239 ns	0.0224 ns	1.741 ns	168 B
LerpFma2	-0	4503599627370497	1.937 ns	0.0112 ns	0.0094 ns	1.933 ns	161 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
NormalFma	1	3.141592653589793	1.928 ns	0.0138 ns	0.0116 ns	1.926 ns	152 B
LerpFma1	1	3.141592653589793	1.697 ns	0.0076 ns	0.0071 ns	1.698 ns	168 B
LerpFma2	1	3.141592653589793	1.934 ns	0.0195 ns	0.0182 ns	1.932 ns	161 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
NormalFma	2.718281828459045	3.141592653589793	1.916 ns	0.0059 ns	0.0049 ns	1.916 ns	152 B
LerpFma1	2.718281828459045	3.141592653589793	1.710 ns	0.0086 ns	0.0076 ns	1.711 ns	168 B
LerpFma2	2.718281828459045	3.141592653589793	1.933 ns	0.0103 ns	0.0096 ns	1.929 ns	161 B

多少のブレはあるものの、大きな違いはなさそうです。

HalfLerp

double r = Rng.NextDouble();  
return 2 * (min / 2 + (max / 2 - min / 2) * r);

「普通のやつ」に似ているのですが、いったん半分にして計算してから 2 倍しています。
何が違うのかというと、半分になったおかげでオーバーフローしません。あと 2 冪倍は非正規化数でない限り精度を失わないので同じ精度で計算できます。

開区間にならない

例えば、 min == 0.57526014729317487, max == 0.89380543555156478, rng == 0.99999999999999989 のときに 0.89380543555156478 を返します。

望ましくない結果は返さない

前述したように、オーバーフローしないため ±∞ を返したりはしません。

確率が均等にならない

「普通のやつ」と同様に $[0.0, 1.0)$ の乱数を引き延ばしている以上、出やすい値や出にくい値、全く出ない値が生じる可能性があります。

なお、「非正規化数でない限り」と以前書いたように、非正規化数だった場合は精度を失って確率の偏りが顕著になる場合があります。
極端な例を挙げれば、 halfLerp(-double.Epsilon, double.Epsilon) は常に 0 を返します。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
HalfLerp	-3.141592653589793	3.141592653589793	1.927 ns	0.0206 ns	0.0193 ns	1.923 ns	171 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
HalfLerp	0	3.141592653589793	1.729 ns	0.0258 ns	0.0215 ns	1.725 ns	171 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
HalfLerp	-0	4503599627370497	2.004 ns	0.0669 ns	0.0771 ns	1.972 ns	171 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
HalfLerp	1	3.141592653589793	1.906 ns	0.0112 ns	0.0100 ns	1.907 ns	171 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
HalfLerp	2.718281828459045	3.141592653589793	1.923 ns	0.0113 ns	0.0100 ns	1.921 ns	171 B

これもほとんど変わらないですね。

Matt 式

int SampleToPowerOfTwoExponent(int exponent)  
{  
    while (exponent > -0x3fe)  
    {  
        int lz = BitOperations.LeadingZeroCount(Rng.Next());  
        if (lz == 64)  
        {  
            exponent -= 64;  
        }  
        else  
        {  
            return Math.Max(-0x3ff, exponent - 1 - lz);  
        }  
    }  
    return -0x3ff;  
}  
  
int SampleExponent(int emin, int emax)  
{  
    int c = 0;  
    while (true)  
    {  
        int lz = BitOperations.LeadingZeroCount(Rng.Next());  
        if (lz == 64)  
        {  
            c += 64;  
        }  
        else  
        {  
            return emax - 1 - ((c + lz) % (emax - emin));  
        }  
    }  
}  
  
// assumes a >= 0, b >= 0  
double SampleRange(double a, double b)  
{  
    ulong aa = BitConverter.DoubleToUInt64Bits(a);  
    ulong bb = BitConverter.DoubleToUInt64Bits(b);  
  
    int ea = (int)(aa >> 52) - 0x3ff;  
    int eb = (int)(bb >> 52) - 0x3ff;  
    ulong ma = aa & ((1ul << 52) - 1);  
    ulong mb = bb & ((1ul << 52) - 1);  
  
    if (mb != 0)  
    {  
        eb++;  
    }  
  
    while (true)  
    {  
        int e;  
        if (ea == -0x3ff)  
        {  
            e = SampleToPowerOfTwoExponent(eb);  
        }  
        else  
        {  
            e = SampleExponent(ea, eb);  
        }  
  
        double v = BitConverter.UInt64BitsToDouble((ulong)(e + 0x3ff) << 52 | (Rng.Next() & ((1ul << 52) - 1)));  
        if (a <= v && v < b)  
        {  
            return v;  
        }  
    }  
}  
  
if (min >= 0.0)  
{  
    return SampleRange(min, max);  
}  
if (max < 0.0)  
{  
    return -SampleRange(-max, -min);        // TODO: it will be (min, max]  
}  
  
double absmax = Math.Max(Math.Abs(min), Math.Abs(max));  
double result;  
do  
{  
    result = SampleRange(0, absmax);  
    if ((Rng.Next() & 1) != 0)  
    {  
        result = -result;  
    }  
} while (result < min || max <= result);  
  
return result;

以下のブログ記事のコードをベースに、 float 用だったのを double に改造した実装です。

実装の説明

メインとなるのは SampleRange() です。

まず、 b (絶対値が大きいほう) の仮数部が 0 でなければ、指数部を 1 増やします。
逆に言えば、仮数部が 0 のときは 1 小さくなるということです。 $[a, b)$ の $b$ のほうに当たるわけですね。

次に、 a (絶対値が小さいほう) の指数部が 0 (非正規化数領域) かどうかで処理を分岐します。

SampleToPowerOfTwoExponent と SampleExponent では、似たようなことをやっています。
LeadingZeroCount による効率化を除けば、やっていることは「確率 50% で exponent を減算、もう 50% で処理を終了」といった感じで、要するに幾何分布に従う乱数で埋めているわけなのですが、末端まで行った時の処理が異なります。
SampleToPowerOfTwoExponent では Math.Max でクランプしていますが、 SampleExponent では % (emax - emin) でループさせています。
こうしている理由は、非正規化数エリアの確率の扱いに由来します。

非正規化数 (指数部 0x000) がカバーするエリアの重み (出現確率) は、その一つ上 (指数部 0x001) のエリアの重みと等しいです。
対して、それ以外のエリアでは、一つ上のエリアの重みは一つ下のエリアの重みの 2 倍になります。

SampleToPowerOfTwoExponent のほうは非正規化数を含むため、重みのトータルが $1 + \sum_{i=0}^{e-1} 2^{i} = 2^{e}$ になります。
一方、 SampleExponent のほうは重みのトータルが $\sum_{i=0}^{e-1} 2^{i} = 2^{e} - 1$ になります。あふれた $1$ に当たった場合は再抽選とみなしてもう一度最初からやる……のですが、これはつまり % でループさせてしまっても問題ないということです。

幾何分布には「無記憶性」という性質があります。
コイントスをしていて $n$ 回失敗している状態でさらに $k$ 回失敗する確率と、単に最初から $k$ 回失敗する確率は同じ (つまり、過去のことは未来の確率に関係しない) ということです。

なので、わざわざ再抽選しなくとも、過去のことはなかったことにして (emax - emin を引いて) よい、それを繰り返しても構わない (% (emax - emin) としてよい) ということになります。

以上から指数部が求められたら、仮数部を 52 bit 一様乱数で埋めます。
あとは範囲チェックして OK ならそれを返し、 NG なら再生成するという流れです。

以上で説明したのは a, b が同符号だった場合でしたが、異符号だった場合はちょっと処理が増えます。
$[0.0, \mathbb{max}(|a|, |b|))$ として上記の生成処理をやった後、ランダムに符号をつけます。それが範囲内なら返して範囲外なら再生成、という流れです。

開閉区間の実装が不完全

オリジナルの実装では $[a, b)$ の区間を想定していますが、 $a, b$ がともに正の場合のみ正しいです。
ともに負の場合は $(a, b \rbrack$ 、正負に分かれている場合は $\lbrack a, b \rbrack$ になります。

望ましくない結果を返さない

∞ や NaN といった望ましくない結果を返しません。

確率が「ほぼ」均等になる

範囲内の値すべてがほぼ正しい確率で出現します。

「ほぼ」とつけたのは Downey の補正がないからです。
これについても改造するのは簡単でしょう。

あと、符号を跨いだ場合に ±0 それぞれのぶんで実質 0 が 2 倍出やすくなる問題もあります。
これについても Downey の補正をかければ確率がそれぞれ 1/2 になるので解決するはずです。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
Matt	-3.141592653589793	3.141592653589793	12.517 ns	0.2780 ns	0.4644 ns	12.495 ns	882 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
Matt	0	3.141592653589793	6.399 ns	0.0416 ns	0.0389 ns	6.377 ns	861 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
Matt	-0	4503599627370497	18.473 ns	0.1117 ns	0.0990 ns	18.482 ns	786 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
Matt	1	3.141592653589793	9.061 ns	0.0412 ns	0.0386 ns	9.061 ns	786 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
Matt	2.718281828459045	3.141592653589793	50.068 ns	0.1239 ns	0.1159 ns	50.067 ns	770 B

結構重いです。最悪の場合は「普通のやつ」の 25 倍程度の時間がかかっています。
また、入力に依存してかなり速度が異なることも分かります。最速は 6 ns ですが最遅は 50 ns です。

GammaCO

double gamma = Math.Max(Math.BitIncrement(min) - min, max - Math.BitDecrement(max));  
  
// ceiling(b / g - a / g)  
static ulong CeilInt(double a, double b, double g)  
{  
    double s = b / g - a / g;  
    double epsilon;  
    if (Math.Abs(a) <= Math.Abs(b))  
    {  
        epsilon = -a / g - (s - b / g);  
    }  
    else  
    {  
        epsilon = b / g - (s + a / g);  
    }  
  
    ulong si = (ulong)Math.Ceiling(s);  
    return s != si ? si : (si + (epsilon > 0 ? 1ul : 0ul));  
}  
  
ulong hi = CeilInt(min, max, gamma);  
ulong k = 1 + Rng.NextULong(hi - 1);    // [1, hi]  
ulong khi = k >> 2;  
ulong klo = k & 0x3;  
  
if (Math.Abs(min) <= Math.Abs(max))  
{  
    // k == hi ? min : max - k * gamma;  
    return k == hi ? min : 4 * (max / 4 - khi * gamma) - klo * gamma;  
}  
else  
{  
    // min + (k - 1) * gamma;  
    return 4 * (min / 4 + khi * gamma) + (klo - 1) * gamma;  
}

文献 "Drawing random floating-point numbers from an interval" *3 に載っていた手法です。

実装の説明

まず gamma が何かというと、大雑把に言えば ulp を取得する操作です。
ulp は連続する 2 つの浮動小数点数の間隔を表します。簡単に言えば、仮数部の最下位ビットが 0/1 に変わった時の数値の差分です。例えば、 $\mathrm{ulp}(1.0) = 2^{-52}$ になります。

次に、 CeilInt は、コメントにもあるように Ceiling(b / g - a / g) を正確に計算するためのメソッドです。
g (gamma) は ulp なので 2 冪です。ということは除算は基本的には正確に実行することができます。
問題は引き算のほうなのですが、うまいこと処理することで (Dekker's exact summation algorithm だそうです) 誤差なく計算できるようにしています。

次に、 k に $\lbrack 1, hi \rbrack$ の一様分布整数乱数を入れます。
NextULong(max) の実装には、もちろん Rng.Next() % max などではなく Lemire 式を使うとよいでしょう。

// [0, max) の一様分布整数乱数を取得  
public ulong NextULong(ulong max)  
{  
    ulong hi = Math.BigMul(Next(), max, out var lo);  
  
    if (lo < max)  
    {  
        ulong mod = (0ul - max) % max;  
        while (lo < mod)  
        {  
            hi = Math.BigMul(Next(), max, out lo);  
        }  
    }  
  
    return hi;  
}

Rng.Next() % max ではダメな理由は、確率が偏るためです。
max が 2 冪でない限り、必ず 0 が出る確率＞ max-1 が出る確率になります。

それ以降の計算は、雑に言えば a + (k - 1) * g みたいな感じです。コメントアウトされているコードがあるかと思いますが、これをオーバーフローしないように注意深く実装するとこういう感じになります。

なので、一言で言えば $min + \mathrm{rand}((max - min)/ulp) * ulp$ という感じです。伝わりましたでしょうか……？

この手法の利点としては、ビット操作を必要としないこと、ループを必要としないこと (ただし、 NextULong() 内でループは必要です) でしょうか。

開区間になる

なお、論文のオリジナル実装には $\lbrack a, b \rbrack, (a, b), \lbrack a, b), (a, b \rbrack$ すべての実装が載っています。

望ましくない結果を返さない

論文によれば、 a + (b - a) * x のようにオーバーフロー起因の望ましくない結果は返さないとされています。

論文オリジナルの実装では、一部の計算でオーバーフローを起こして正しくない結果を返す場合があったそうです。
本項冒頭のコードは、その問題が修正された正誤表のほうのコードを使用しています。

確率は「ほぼ」均等になるが、すべての値は出ない

「ほぼ」とつけたのは、空間的に一様分布にならない可能性があるためです。
論文曰く、 max - min が gamma の倍数でない場合にひとつだけ距離間隔が縮まる場合があるとのことです。

また、すべての値は出ません。
gamma は min と max でより大きいほうの ulp に依存するので、小さい ulp の区間では飛び飛びにしか値が出現しません。

上の図で説明しましょう。
浮動小数点数は通常上側の目盛りのように左 (小さいほう) に行けば行くほど間隔が密になるわけなのですが、この手法だと下側の目盛りのように一番大きな間隔に依存していて左に行っても等間隔になります。これを論文では "Spatial equidistribution" と表現しています。
なので、一番大きな目盛りの区間 (4 ～ 8) ではすべての値が出現しますが、それより小さな区間 (1 ～ 4) では飛び飛びに出現しない値 (例えば、 2.5 や 1.25) が出てきます。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
GammaCO	-3.141592653589793	3.141592653589793	9.146 ns	0.0238 ns	0.0211 ns	9.151 ns	1,127 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
GammaCO	0	3.141592653589793	9.272 ns	0.1411 ns	0.1319 ns	9.243 ns	1,127 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
GammaCO	-0	4503599627370497	9.091 ns	0.2059 ns	0.2022 ns	8.980 ns	1,127 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
GammaCO	1	3.141592653589793	9.119 ns	0.0355 ns	0.0332 ns	9.125 ns	1,127 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
GammaCO	2.718281828459045	3.141592653589793	9.125 ns	0.0359 ns	0.0319 ns	9.124 ns	1,127 B

「普通のやつ」の 4.5 倍程度の時間がかかっていますが、 Matt 法のように入力値依存で速度が変わることはなさそうです。

余談: PHP での実装について

PHP の Random\Randomizer::getFloat は、このアルゴリズムで実装されているそうです。

Cauldron

int sign;  
  
ulong minExponent, minMantissa, maxExponent, maxMantissa;  
  
// bit 表現の取得  
{  
    double a, b;  
  
    sign = (min < 0.0 ? 1 : 0) + (max < 0.0 ? 1 : 0);  
    if (sign == 0)  
    {  
        a = min;  
        b = max;  
    }  
    else if (sign == 1)  
    {  
        a = 0;  
        b = (max > -min) ? max : -min;  
    }  
    else  
    {  
        a = max;  
        b = min;  
    }  
  
  
    ulong aBits = BitConverter.DoubleToUInt64Bits(a);  
    ulong bBits = BitConverter.DoubleToUInt64Bits(b);  
  
    minExponent = (aBits >> 52) & 0x7ff;  
    minMantissa = aBits & ((1ul << 52) - 1);  
    maxExponent = (bBits >> 52) & 0x7ff;  
    maxMantissa = bBits & ((1ul << 52) - 1);  
}  
  
// 指数部が等しいとき  
if (minExponent == maxExponent)  
{  
    ulong result = (minExponent << 52) | (Rng.NextULong(maxMantissa - minMantissa + 1) + minMantissa);  
  
    if (sign == 1)  
    {  
        result |= Rng.Next() << 63;  
    }  
    else if (sign == 2)  
    {  
        result |= 1ul << 63;  
    }  
  
    return BitConverter.UInt64BitsToDouble(result);  
}  
  
// 指数部が 1 だけ異なるとき  
if (minExponent + 1 == maxExponent && minExponent > 0)  
{  
    ulong invMinMantissa = ((1ul << 52) - 1) - minMantissa;  
    ulong range = invMinMantissa + maxMantissa + 1;  
  
    ulong exponent, mantissa;  
    ulong x = 0;  
    int i = 0;  
  
    while (true)  
    {  
        if (i <= 3)  
        {  
            x = Rng.Next();  
            i = 64;  
        }  
  
        if ((x & 1) != 0)  
        {  
            i--;  
            x >>= 1;  
  
            exponent = maxExponent;  
            mantissa = Rng.NextULong(range);  
            if (mantissa <= maxMantissa)  
            {  
                break;  
            }  
        }  
        else if ((x & 2) != 0)  
        {  
            i -= 2;  
            x >>= 2;  
  
            exponent = minExponent;  
            mantissa = Rng.NextULong(range);  
            if (mantissa <= invMinMantissa)  
            {  
                mantissa = ((1ul << 52) - 1) - mantissa;  
                break;  
            }  
        }  
        else  
        {  
            i -= 2;  
            x >>= 2;  
        }  
    }  
  
    ulong result = exponent << 52 | mantissa;  
  
    if (sign == 1)  
    {  
        result |= x << 63;  
    }  
    else if (sign == 2)  
    {  
        result |= 1ul << 63;  
    }  
  
    return BitConverter.UInt64BitsToDouble(result);  
}  
  
// それ以外  
while (true)  
{  
    ulong exponent = maxExponent;  
  
    ulong x;  
    while ((x = Rng.Next()) == 0)  
    {  
        exponent -= 64;  
    }  
  
    exponent -= (ulong)BitOperations.TrailingZeroCount(x);  
  
    if (exponent < minExponent || exponent > maxExponent)  
    {  
        exponent = 0;  
    }  
  
    x = Rng.Next();  
    ulong result = exponent << 52 | x >> 12;  
  
    if (sign == 1)  
    {  
        result |= x << 63;  
    }  
    else if (sign == 2)  
    {  
        result |= 1ul << 63;  
    }  
  
    double doubleResult = BitConverter.UInt64BitsToDouble(result);  
    if (min <= doubleResult && doubleResult <= max)  
    {  
        return doubleResult;  
    }  
}

https://github.com/camel-cdr/cauldron/blob/main/cauldron/random.h

実装の説明

ながい！けど長さには理由があります。

最初のブロックでは初期化を行っています。
sign は 0 なら min, max がどちらも正の場合、 1 なら正負を跨ぐ場合、 2 ならどちらも負の場合となります。

次のブロックでは、 min, max の指数部が同じだった場合の特殊処理を行っています。
この場合は仮数部を $\lbrack min, max \rbrack$ の一様分布整数乱数で埋めればよいですね。

その次のブロックでは、 min, max の指数部が 1 異なる場合 (かつ、非正規化数でない場合) の特殊処理を行っています。
1/2 の確率で maxExponent 側の指数部の処理を、 1/4 の確率で minExponent 側の指数部の処理を、残りの 1/4 は再抽選としています。
このブロックがなぜ存在するのかというと、最後のブロックの処理にすべてやらせると、仮数部がランダム生成な都合上 $\lbrack 1 - 2^{-52}, 1 + 2^{-52} \rbrack$ みたいなケースで採択率が $1/2^{52}$ のようにひどいことになるためです。

それ以外の場合は、最後のブロックに来ます。
最後のブロックでは、まず exponent を求めます。
exponent は $p=0.5$ の幾何分布に基づく (50% の確率で成功する試行が成功するまでの試行回数に基づく) ため、連続する 0 のビット数を BitOperations.TrailingZeroCount を使って効率よく数えます。
そうしたら符号部と仮数部をランダムに生成して、範囲内ならそれを返し、範囲外なら再抽選するようになっています。

閉区間になる

このアルゴリズムでは閉区間 $\lbrack min, max \rbrack$ になります。

望ましくない結果を返さない

全編通してビット操作で組み立てているので、オーバーフローなどの問題は発生しません。

確率は「ほぼ」均等・すべての値を返しうる

ほぼ、とつけたのは例によって Downey の補正がないためです。
これも補正の実装はすぐできるでしょう。

非正規化数まわりの確率が正しくないかも、といったコメントがオリジナルのソースコードにありましたが、私が調べた限りでは問題なさそうでした。
別件の問題はありましたが……これについては後述します。

半無限ループの可能性

例えば min == 2.2250738585072004e-308 (0x000ffffffffffffe), max == 2.2250738585072019e-308 (0x0010000000000001) を与えた場合に、半無限ループに陥ります。
というのも、片方が非正規化数なので最後のブロックに行くのですが、仮数部の関係でどちらの指数部でも採択率が $2/2^{52}$ となり、ほぼ採択されず半無限にリトライし続ける状態になります。

「min, max の指数部が 1 異なる場合」の処理をうまく修正する必要がありそうですね。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
Cauldron	-3.141592653589793	3.141592653589793	7.011 ns	0.0212 ns	0.0188 ns	7.016 ns	855 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
Cauldron	0	3.141592653589793	6.831 ns	0.1054 ns	0.0935 ns	6.793 ns	839 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
Cauldron	-0	4503599627370497	15.236 ns	0.1784 ns	0.1669 ns	15.225 ns	839 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
Cauldron	1	3.141592653589793	24.717 ns	0.0667 ns	0.0624 ns	24.742 ns	1,036 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
Cauldron	2.718281828459045	3.141592653589793	2.298 ns	0.0111 ns	0.0093 ns	2.293 ns	1,007 B

入力値にもよりますが、速いときは速いです。ただ、入力値によっては「普通のやつ」の 12 倍になるなど、改善の余地はありそうです。

パフォーマンスを向上できそうな点としては、以下が考えられます。

あふれた exponent を棄却するのではなく Matt 式のように % でループさせる
非正規化数をうまく扱えるように頑張る

オレオレ手法

さて、ここで自分で実装するならどうするかを検討してみましょう。

まず、低精度なやつならいくらでも作れるので、高精度な (つまり、範囲内のすべての表現可能な数が正しい確率で出現する) 手法を考えるものとします。
そのうえで、できる限り速いと嬉しいですね。
浮動小数点演算は何が起こるか分からないので、ビットパターンを直接構築する手法で考えてみましょう。

その前に - 確率のおさらい

一般に、指数部は幾何分布乱数で、仮数部は一様分布乱数で埋めればよいです。
ただ、一部例外があるので、おさらいしておきましょう。

$2^{n}$ (仮数部が 0) のとき

Downey の補正の項で前述したように、一般に $2^{n}$ (仮数部が 0) の場合の確率は通常 (仮数部が非 0 の場合) の $3/4$ になります。
なお、「通常」と書いたように例外もあって、指数部が 0 (非正規化数) の場合は $1/2$ (-0 を含める場合) または $1$ (-0 を除外する場合) 、指数部が 1 (非正規化数のひとつ上) の場合は $1$ となります。

-0 の取り扱い

そのまま実装すると「0」の確率が 2 倍になってしまうので、 -0 を棄却 (再抽選) するか、 +0 ・ -0 の確率をそれぞれ 1/2 にするか、を行う必要があります。

今回は後者の実装を採用します。

開閉区間について

閉区間の場合、その端点 (min か max そのもの) に当たる確率は、端点の仮数部が 0 でない場合通常の $1/2$ になります。
なぜかというと、無限精度の乱数を浮動小数点数に丸めることを考えたとき、端点が占める面積はそれ以外の場所の半分だけになるからです。

この図は、仮数部が 2 bit の浮動小数点数で $\lbrack 1.5, 3.0 \rbrack$ の範囲の乱数を生成した場合の範囲の図です。
上段が実数 $\mathbb{R}$ (無限精度) 、下段が浮動小数点数 FP を表します。
このとき、下段の 1.5 と 3 は 1/2 だけが範囲内になっているのがおわかりいただけますでしょうか。

なお、仮数部が 0 の場合は 1/2 になるとは限らず、もうちょっとややこしくなります。上図を参考に考えてみてください。

また、開区間の場合は、当然端点に当たる確率は $0$ になります。

ところで、開閉区間の実装において、「閉区間 $\lbrack min, max \rbrack$ だけ実装すれば $\lbrack \mathrm{next}(min), \mathrm{prev}(max) \rbrack$ で開区間 $(min, max)$ を表せるから、実装を使いまわせるのでは？」という考えが浮かぶかもしれません。

$\mathrm{prev}(x)$ は $x$ より前 (負方向) に隣接する浮動小数点数 (つまり Math.BitDecrement(x)) 、 $\mathrm{next}(x)$ は $x$ の後 (正方向) に隣接する浮動小数点数 (つまり Math.BitIncrement(x)) を表します。

しかし、これは誤りです。
理由を説明しましょう。

上図は正しい $(1.5, 3.0)$ です。
端点の隣 (1.75 と 2.5) の領域は 100% 塗られています。

一方、上図は誤った実装 $\lbrack \mathrm{next}(1.5), \mathrm{prev}(3.0) \rbrack$ です。
端点の隣 (1.75 と 2.5) の領域が 50% だけ塗られている、つまり半分の確率でしか出現しなくなっています。

対処としては、もちろん専用の実装をするのが最善ではありますが、次善の策として $\lbrack min, max \rbrack$ で生成してから $min, max$ のどちらかに等しければ再生成、という安直な処理もあります。

場合分け

直接コードを示す前に、どうやって設計するか考えてみましょう。
まずは場合分けです。

符号部が同じ場合
- 指数部が同じ場合
- 指数部が 1 異なる場合
  - min が非正規化数の場合
  - min が正規化数の場合
- min が非正規化数の場合
- それ以外 (指数部が 2 以上異なり、 min が正規化数の場合)
符号部が異なる場合
- 指数部が同じ場合
  - min と max が非正規化数の場合
  - それ以外 (min と max が正規化数の場合)
- それ以外 (指数部が異なる場合)

符号部が同じ - 指数部が同じ

この場合は、仮数部を一様分布整数乱数を使って minMantissa + nextULong(maxMantissa - minMantissa) で埋めればよいです。
前述したように端点の確率は 1/2 にすべきなので、これが minMantissa に等しいときは 50% の確率で maxMantissa に振り替える処理が必要です。

符号部が同じ - それ以外

こっちを先に説明したほうが分かりやすいと思うので先に書きます。

仮数部は 52 bit 一様乱数で埋めます。

次に、指数部は幾何分布乱数をもとに決定します。
Downey の補正のため、仮数部が 0 の場合は 50% の確率で指数部を +1 します。
溢れた指数部は exponent % (maxExponent - minExponent + 1) でループさせます。理由は Matt 式で述べているのと同じです。
最後に、生成した乱数が min か max に等しい場合、 50% の確率で再抽選を行います。

この時点で得られるのは $\lbrack 2^{emin}, 2^{emax+1} \rbrack$ の乱数ですので、 $\lbrack min, max \rbrack$ の範囲からはみ出す場合があります。
なので範囲チェックをして OK ならそれを返し、 NG なら再生成します。

符号部が同じ - `min` が非正規化数の場合

上記の処理とほぼ同じですが、溢れた指数部を Math.Max(exponent, 0) でクランプする点が異なります。
この理由も Matt 式で述べたとおりですね。

符号部が同じ - 指数部が 1 異なる - `min` が正規化数の場合

なんで「指数部が 1 異なる場合」を特別扱いする必要があるのかというと、そのままだと $\lbrack \mathrm{prev}(2.0), \mathrm{next}(2.0) \rbrack$ みたいな入力を与えた場合に採択率が $1/2^{52}$ になってしまい、処理に異常に時間がかかるからです。

説明を簡単にするため、各場所に名前を付けましょう：

α: 指数部が小さいほうで、仮数部が 0
β: 指数部が小さいほうで、仮数部が非 0
γ: 指数部が大きいほうで、仮数部が 0
δ: 指数部が大きいほうで、仮数部が非 0

この場合、各エリアの確率の倍率 (上図で β が出る確率を 1 としたときの倍率) は以下のようになります。

場所	α	β	γ	δ
始点	1/2	1/2	✖	✖
中間	✖	1	3/2	2
終点	✖	✖	1/2	1

基本的には、一様分布整数乱数を使って仮数部を生成します。
r = minMantissa + nextULong((1ul << 52) - minMantissa + ((maxMantissa + 1) * 2)) みたいな感じですね。
それで、 r が 1ul << 52 より小さければ minExponent 側で、以上なら maxExponent 側にするみたいな感じです。
あとは上記の確率に則って再抽選を行えばよいです。

符号部が同じ - 指数部が 1 異なる - `min` が非正規化数の場合

概ね min が正規化数の場合と同様ですが、確率の倍率が異なります。

場所	α	β	γ	δ
始点	1/2	1/2	✖	✖
中間	✖	1	1	1
終点	✖	✖	1/2	1/2

こっちは簡単ですね。始点か終点なら 1/2 の確率で再抽選すればよいです。

符号部が異なる - 指数部が同じ - `min` と `max` が非正規化数

この場合も、「指数部が 1 異なる場合」と同様に特別扱いする必要があります。
例えば、 $\lbrack \mathrm{prev}(-0.0), \mathrm{next}(0.0) \rbrack$ を与えた場合を考えてみましょう。仮数部をランダムに発生させる手法だと採択率が $2/2^{52}$ になり、非常に時間がかかってしまいます。

この場合の手法としても「指数部が 1 異なる場合」と概ね同様です。
r = nextULong(minMantissa + MaxMantissa + 2) で一様分布整数乱数を得て、 minMantissa 以下ならそれを、より大きければ r - minMantissa - 1 を仮数部として構成し、あとはそれぞれマイナス・プラスの符号をつければよいです。

なお、 r が 0 になった時か、 minMantissa または maxMantissa に等しくなった時は 1/2 の確率で再抽選しましょう。

符号部が異なる - 指数部が同じ - `min` と `max` が正規化数

指数部が同じで正規化数な場合は、符号部が + になる確率と - になる確率はほぼ均等だと思われるので、 1 bit の乱数を使ってランダムに決めます。
仮数部も 52 bit 乱数でランダムに決めます。
指数部は例によって幾何分布乱数で決めます。 Math.Max でクランプするのを忘れずに。

あとは前述の表を見ながら再抽選するわけですね。
過程を省くとこんな感じになります。

if (x == min) {  
    // 1/2 の確率で再抽選  
} else if (x == max) {  
    if (mantissa == 0 && exponent >= 2) {  
        // 3/4 の確率で再抽選  
    } else {  
        // 1/2 の確率で再抽選  
    }  
} else if (mantissa == 0 && exponent >= 2) {  
    // 1/4 の確率で再抽選  
}

あとは範囲チェックして範囲内なら採択、範囲外なら再生成します。

符号部が異なる - 指数部が異なる

指数部が異なる場合、符号部が + になる確率と - になる確率は大幅に偏ります。そのため、単純にランダムに決定すると棄却率が 1/2 に近づいてしまいます。
そのため、以下のようなアルゴリズムである程度偏りを再現します。

// 指数部の差  
int exponentDiff = Math.Min(Math.Abs((int)(minExponent >> 52) - (int)(maxExponent >> 52)), 63);  
  
// [0, 2^exponentDiff] の乱数を生成して、それが 0 なら  
if (Rng.NextULong((1ul << exponentDiff) + 1) == 0)  
{  
    // 符号と指数シフトを設定  
    sign = minExponent == Math.Min(minExponent, maxExponent) ? (1ul << 63) : 0ul;  
    exponentShift = exponentDiff;  
}  
else  
{  
    sign = minExponent == Math.Max(minExponent, maxExponent) ? (1ul << 63) : 0ul;  
    exponentShift = 0;  
}

exponentShift は、この後に続く指数部の処理で最初から指数部を exponentShift だけ減らしておくための変数です。
結果を $2^{-shift}$ している感じです。

後の処理は指数部が同じ場合と同じです。

なんで「ある程度」偏らせる（＝完全に同じ確率に偏らせるわけではない）ので良いのかというと、棄却採択法だからです。
本来の確率分布より多少大きくなっても結局棄却されるので問題ないのです。

コード

以上をコードに落とし込むとこんな感じになります。

ulong minBits = BitConverter.DoubleToUInt64Bits(min);  
ulong maxBits = BitConverter.DoubleToUInt64Bits(max);  
  
ulong minSign = minBits & (1ul << 63);  
ulong maxSign = maxBits & (1ul << 63);  
ulong minExponent = minBits & (0x7fful << 52);  
ulong maxExponent = maxBits & (0x7fful << 52);  
ulong minMantissa = minBits & ((1ul << 52) - 1);  
ulong maxMantissa = maxBits & ((1ul << 52) - 1);  
  
ulong r = 0;  
int entropy = 0;  
  
  
if (minSign == maxSign)  
{  
    if (minExponent == maxExponent)  
    {  
        ulong mantissa = minMantissa + Rng.NextULong(maxMantissa - minMantissa);  
  
        if (mantissa == minMantissa)  
        {  
            if ((Rng.Next() & 1) != 0)  
            {  
                mantissa = maxMantissa;  
            }  
        }  
  
        return BitConverter.UInt64BitsToDouble(minSign | minExponent | mantissa);  
    }  
    else if (minExponent + (1ul << 52) == maxExponent)  
    {  
        if (minExponent == 0)  
        {  
            while (true)  
            {  
                ulong mantissa;  
                ulong exponent;  
                ulong mantissaRange = ((1ul << 52) - minMantissa) + ((maxMantissa + 1));  
                ulong rr = minMantissa + Rng.NextULong(mantissaRange);  
                if (rr < (1ul << 52))  
                {  
                    mantissa = rr;  
                    exponent = minExponent;  
                }  
                else  
                {  
                    mantissa = rr - (1ul << 52);  
                    exponent = maxExponent;  
                }  
  
                ulong x = minSign | exponent | mantissa;  
                if (x == minBits || x == maxBits)  
                {  
                    if (entropy < 1)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    bool flag = (r & 1) != 0;  
                    r >>= 1;  
                    entropy--;  
  
                    if (flag)  
                    {  
                        continue;  
                    }  
                }  
  
                return BitConverter.UInt64BitsToDouble(x);  
            }  
        }  
        else  
        {  
            while (true)  
            {  
                ulong mantissa;  
                ulong exponent;  
                ulong mantissaRange = ((1ul << 52) - minMantissa) + ((maxMantissa + 1) << 1);  
                ulong rr = minMantissa + Rng.NextULong(mantissaRange);  
                if (rr < (1ul << 52))  
                {  
                    mantissa = rr;  
                    exponent = minExponent;  
  
                    if (mantissa == minMantissa)  
                    {  
                        if (entropy < 1)  
                        {  
                            r = Rng.Next();  
                            entropy = 64;  
                        }  
  
                        bool flag = (r & 1) != 0;  
                        r >>= 1;  
                        entropy--;  
  
                        if (flag)  
                        {  
                            continue;  
                        }  
                    }  
                }  
                else  
                {  
                    mantissa = (rr - (1ul << 52)) >> 1;  
                    exponent = maxExponent;  
  
                    if (mantissa == 0)  
                    {  
                        if (mantissa == maxMantissa)  
                        {  
                            if (entropy < 2)  
                            {  
                                r = Rng.Next();  
                                entropy = 64;  
                            }  
  
                            bool flag = (r & 3) != 0;  
                            r >>= 2;  
                            entropy -= 2;  
  
                            if (flag)  
                            {  
                                continue;  
                            }  
                        }  
                        else  
                        {  
                            if (entropy < 2)  
                            {  
                                r = Rng.Next();  
                                entropy = 64;  
                            }  
  
                            bool flag = (r & 3) == 0;  
                            r >>= 2;  
                            entropy -= 2;  
  
                            if (flag)  
                            {  
                                continue;  
                            }  
                        }  
                    }  
                    else if (mantissa == maxMantissa)  
                    {  
                        if (entropy < 1)  
                        {  
                            r = Rng.Next();  
                            entropy = 64;  
                        }  
  
                        bool flag = (r & 1) != 0;  
                        r >>= 1;  
                        entropy--;  
  
                        if (flag)  
                        {  
                            continue;  
                        }  
                    }  
                }  
  
                ulong x = minSign | exponent | mantissa;  
                return BitConverter.UInt64BitsToDouble(x);  
            }  
        }  
    }  
    else if (minExponent == 0)  
    {  
        while (true)  
        {  
            ulong mantissa;  
            {  
                if (entropy < 52)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                mantissa = r & ((1ul << 52) - 1);  
                r >>= 52;  
                entropy -= 52;  
            }  
  
            int exponentShift = 0;  
            do  
            {  
                if (entropy <= 0)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                int ctz = Math.Min(BitOperations.TrailingZeroCount(r), entropy);  
                r >>= ctz + 1;  
                entropy -= ctz + 1;  
                exponentShift += ctz;  
  
            } while (entropy == -1);  
  
            if (mantissa == 0)  
            {  
                if (entropy < 1)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                bool flag = (r & 1) != 0;  
                r >>= 1;  
                entropy--;  
  
                if (flag)  
                {  
                    exponentShift--;  
  
                    if (exponentShift < 0)  
                    {  
                        continue;  
                    }  
                }  
            }  
  
            ulong exponent = (ulong)Math.Max(((long)maxExponent >> 52) - exponentShift, 0) << 52;  
  
            ulong x = minSign | exponent | mantissa;  
            if (x == minBits || x == maxBits)  
            {  
                if (entropy < 1)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                bool flag = (r & 1) != 0;  
                r >>= 1;  
                entropy--;  
  
                if (flag)  
                {  
                    continue;  
                }  
            }  
  
            double result = BitConverter.UInt64BitsToDouble(x);  
            if (min <= result && result <= max)  
            {  
                return result;  
            }  
        }  
    }  
    else  
    {  
        while (true)  
        {  
            ulong mantissa;  
            {  
                if (entropy < 52)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                mantissa = r & ((1ul << 52) - 1);  
                r >>= 52;  
                entropy -= 52;  
            }  
  
            int exponentShift = 0;  
            do  
            {  
                if (entropy <= 0)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                int ctz = Math.Min(BitOperations.TrailingZeroCount(r), entropy);  
                r >>= ctz + 1;  
                entropy -= ctz + 1;  
                exponentShift += ctz;  
  
            } while (entropy == -1);  
  
            if (mantissa == 0)  
            {  
                if (entropy < 1)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                bool flag = (r & 1) != 0;  
                r >>= 1;  
                entropy--;  
  
                if (flag)  
                {  
                    exponentShift--;  
  
                    if (exponentShift < 0)  
                    {  
                        continue;  
                    }  
                }  
            }  
  
            ulong exponent;  
            {  
                ulong exponentRange = ((maxExponent - minExponent) >> 52) + 1;  
                ulong exponentSub = (ulong)exponentShift < exponentRange ? (ulong)exponentShift : ((ulong)exponentShift % exponentRange);  
                exponent = maxExponent - (exponentSub << 52);  
            }  
  
            ulong x = minSign | exponent | mantissa;  
            if (x == minBits || x == maxBits)  
            {  
                if (entropy < 1)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                bool flag = (r & 1) != 0;  
                r >>= 1;  
                entropy--;  
  
                if (flag)  
                {  
                    continue;  
                }  
            }  
  
            double result = BitConverter.UInt64BitsToDouble(x);  
            if (min <= result && result <= max)  
            {  
                return result;  
            }  
        }  
    }  
}  
else  
{  
    if (minExponent == maxExponent)  
    {  
        if (minExponent == 0)  
        {  
            while (true)  
            {  
                ulong rr = Rng.NextULong(minMantissa + maxMantissa + 2);  
                if (rr <= minMantissa)  
                {  
                    if ((rr == 0 || rr == minMantissa) && (Rng.Next() & 1) != 0)  
                    {  
                        continue;  
                    }  
  
                    return BitConverter.UInt64BitsToDouble(minSign | rr);  
                }  
                else  
                {  
                    rr -= minMantissa + 1;  
  
                    if ((rr == 0 || rr == maxMantissa) && (Rng.Next() & 1) != 0)  
                    {  
                        continue;  
                    }  
  
                    return BitConverter.UInt64BitsToDouble(rr);  
                }  
            }  
        }  
        else  
        {  
            while (true)  
            {  
                ulong sign;  
                {  
                    if (entropy <= 0)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    sign = r << 63;  
                    r >>= 1;  
                    entropy--;  
                }  
  
                ulong mantissa;  
                {  
                    if (entropy < 52)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    mantissa = r & ((1ul << 52) - 1);  
                    r >>= 52;  
                    entropy -= 52;  
                }  
  
                int exponentShift = 0;  
                do  
                {  
                    if (entropy <= 0)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    int clz = Math.Min(BitOperations.TrailingZeroCount(r), entropy);  
                    r >>= clz + 1;  
                    entropy -= clz + 1;  
                    exponentShift += clz;  
  
                } while (entropy == -1);  
  
                ulong exponent = (ulong)Math.Max(((long)maxExponent >> 52) - exponentShift, 0) << 52;  
  
                ulong x = sign | exponent | mantissa;  
                if (x == minBits)  
                {  
                    if (entropy < 1)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    bool flag = (r & 1) != 0;  
                    r >>= 1;  
                    entropy--;  
  
                    if (flag)  
                    {  
                        continue;  
                    }  
                }  
                else if (x == maxBits)  
                {  
                    if (mantissa == 0 && exponent >= 2)  
                    {  
                        if (entropy < 2)  
                        {  
                            r = Rng.Next();  
                            entropy = 64;  
                        }  
  
                        bool flag = (r & 3) != 0;  
                        r >>= 2;  
                        entropy -= 2;  
  
                        if (flag)  
                        {  
                            continue;  
                        }  
                    }  
                    else  
                    {  
                        if (entropy < 1)  
                        {  
                            r = Rng.Next();  
                            entropy = 64;  
                        }  
  
                        bool flag = (r & 1) != 0;  
                        r >>= 1;  
                        entropy--;  
  
                        if (flag)  
                        {  
                            continue;  
                        }  
                    }  
                }  
                else if (mantissa == 0)  
                {  
                    if (exponent >= 2)  
                    {  
                        if (entropy < 2)  
                        {  
                            r = Rng.Next();  
                            entropy = 64;  
                        }  
  
                        bool flag = (r & 3) == 0;  
                        r >>= 2;  
                        entropy -= 2;  
  
                        if (flag)  
                        {  
                            continue;  
                        }  
                    }  
                    else if (exponent == 0)  
                    {  
                        if (entropy < 1)  
                        {  
                            r = Rng.Next();  
                            entropy = 64;  
                        }  
  
                        bool flag = (r & 1) != 0;  
                        r >>= 1;  
                        entropy--;  
  
                        if (flag)  
                        {  
                            continue;  
                        }  
                    }  
                }  
  
                double result = BitConverter.UInt64BitsToDouble(x);  
                if (min <= result && result <= max)  
                {  
                    return result;  
                }  
            }  
        }  
    }  
    else  
    {  
        while (true)  
        {  
            ulong sign;  
            int exponentShift;  
            {  
                ulong lesserExponent = Math.Min(minExponent, maxExponent);  
                ulong greaterExponent = Math.Max(minExponent, maxExponent);  
  
                int exponentDiff = Math.Min((int)((greaterExponent - lesserExponent) >> 52), 63);  
  
                if (Rng.NextULong((1ul << exponentDiff) + 1) == 0)  
                {  
                    sign = lesserExponent == minExponent ? (1ul << 63) : 0ul;  
                    exponentShift = exponentDiff;  
                }  
                else  
                {  
                    sign = greaterExponent == minExponent ? (1ul << 63) : 0ul;  
                    exponentShift = 0;  
                }  
            }  
  
            ulong mantissa;  
            {  
                if (entropy < 52)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                mantissa = r & ((1ul << 52) - 1);  
                r >>= 52;  
                entropy -= 52;  
            }  
  
            do  
            {  
                if (entropy <= 0)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                int ctz = Math.Min(BitOperations.TrailingZeroCount(r), entropy);  
                r >>= ctz + 1;  
                entropy -= ctz + 1;  
                exponentShift += ctz;  
  
            } while (entropy == -1);  
  
            ulong exponent;  
            {  
                ulong greaterExponent = Math.Max(minExponent, maxExponent);  
  
                exponent = (ulong)Math.Max(((long)greaterExponent >> 52) - exponentShift, 0) << 52;  
            }  
  
            ulong x = sign | exponent | mantissa;  
            if (x == minBits)  
            {  
                if (entropy < 1)  
                {  
                    r = Rng.Next();  
                    entropy = 64;  
                }  
  
                bool flag = (r & 1) != 0;  
                r >>= 1;  
                entropy--;  
  
                if (flag)  
                {  
                    continue;  
                }  
            }  
            else if (x == maxBits)  
            {  
                if (mantissa == 0 && exponent >= 2)  
                {  
                    if (entropy < 2)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    bool flag = (r & 3) != 0;  
                    r >>= 2;  
                    entropy -= 2;  
  
                    if (flag)  
                    {  
                        continue;  
                    }  
                }  
                else  
                {  
                    if (entropy < 1)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    bool flag = (r & 1) != 0;  
                    r >>= 1;  
                    entropy--;  
  
                    if (flag)  
                    {  
                        continue;  
                    }  
                }  
            }  
            else if (mantissa == 0)  
            {  
                if (exponent >= 2)  
                {  
                    if (entropy < 2)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    bool flag = (r & 3) == 0;  
                    r >>= 2;  
                    entropy -= 2;  
  
                    if (flag)  
                    {  
                        continue;  
                    }  
                }  
                else if (exponent == 0)  
                {  
                    if (entropy < 1)  
                    {  
                        r = Rng.Next();  
                        entropy = 64;  
                    }  
  
                    bool flag = (r & 1) != 0;  
                    r >>= 1;  
                    entropy--;  
  
                    if (flag)  
                    {  
                        continue;  
                    }  
                }  
            }  
  
            double result = BitConverter.UInt64BitsToDouble(x);  
            if (min <= result && result <= max)  
            {  
                return result;  
            }  
        }  
    }  
}

大長編です。ただ、場合分けが結構多いので実行時に実際に通るパスはそれほど長くはないです。

閉区間になる

上に述べたように、この手法は閉区間 $\lbrack min, max \rbrack$ になります。
開区間 (例えば $\lbrack min, max)$ ) にしたい場合は、一番安直な方法としてはこの手法で生成した後 max と同じなら再抽選するコードを入れることです。
本当にパフォーマンスを追求するなら、きちんとコードを全部修正する必要がありますが、まぁ面倒くさいですね……

望ましくない結果を返さない

ビットパターンから構築しているので、オーバーフロー起因の ∞ や NaN を返すことはありません。

パフォーマンス

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
MyMethod	-3.141592653589793	3.141592653589793	7.077 ns	0.1655 ns	0.1840 ns	7.086 ns	3,144 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
MyMethod	0	3.141592653589793	6.461 ns	0.1199 ns	0.1001 ns	6.428 ns	3,137 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
MyMethod	-0	4503599627370497	38.659 ns	0.7987 ns	1.2898 ns	38.805 ns	3,453 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
MyMethod	1	3.141592653589793	7.273 ns	0.0195 ns	0.0183 ns	7.275 ns	3,340 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
MyMethod	2.718281828459045	3.141592653589793	2.318 ns	0.0295 ns	0.0261 ns	2.310 ns	3,328 B

最速の場合 (同符号・同指数の場合) は「普通のやつ」と遜色ないレベルです。
しかし、最悪の場合 (異符号・仮数部 1 の場合) は 17 倍程度遅くなってしまっています。

改善の余地

このアルゴリズムではケアできていませんが、「指数部の差が 2 以上あって、仮数部が 1」みたいなパターンのとき、指数部が maxExponent になる確率が $1/2$ で仮数部が maxMantissa を超えない確率が $2/2^{52} \approx 0$ とかなので、再抽選率が約 $1/2$ ぐらいになって時間がかかってしまいます。
うまいことこういうパターンを処理できればより高速になりそうですね。棄却採択法のミソはどれだけ元の分布に近くて簡単な関数を作れるか、です。

まとめ

パフォーマンスを全部まとめた表を以下に示します。

Method	min	max	Mean	Error	StdDev	Median	Code Size
Normal	-3.141592653589793	3.141592653589793	1.918 ns	0.0169 ns	0.0158 ns	1.920 ns	155 B
Lerp	-3.141592653589793	3.141592653589793	1.710 ns	0.0192 ns	0.0179 ns	1.705 ns	163 B
NormalFma	-3.141592653589793	3.141592653589793	1.943 ns	0.0263 ns	0.0246 ns	1.936 ns	152 B
LerpFma1	-3.141592653589793	3.141592653589793	1.732 ns	0.0227 ns	0.0212 ns	1.736 ns	168 B
LerpFma2	-3.141592653589793	3.141592653589793	1.944 ns	0.0211 ns	0.0197 ns	1.943 ns	161 B
HalfLerp	-3.141592653589793	3.141592653589793	1.927 ns	0.0206 ns	0.0193 ns	1.923 ns	171 B
Matt	-3.141592653589793	3.141592653589793	12.517 ns	0.2780 ns	0.4644 ns	12.495 ns	882 B
GammaCO	-3.141592653589793	3.141592653589793	9.146 ns	0.0238 ns	0.0211 ns	9.151 ns	1,127 B
Cauldron	-3.141592653589793	3.141592653589793	7.011 ns	0.0212 ns	0.0188 ns	7.016 ns	855 B
MyMethod	-3.141592653589793	3.141592653589793	7.077 ns	0.1655 ns	0.1840 ns	7.086 ns	3,144 B
Normal	0	3.141592653589793	1.743 ns	0.0217 ns	0.0203 ns	1.740 ns	155 B
Lerp	0	3.141592653589793	1.930 ns	0.0216 ns	0.0192 ns	1.927 ns	163 B
NormalFma	0	3.141592653589793	1.727 ns	0.0271 ns	0.0254 ns	1.725 ns	152 B
LerpFma1	0	3.141592653589793	1.926 ns	0.0219 ns	0.0194 ns	1.927 ns	168 B
LerpFma2	0	3.141592653589793	1.766 ns	0.0613 ns	0.0656 ns	1.758 ns	161 B
HalfLerp	0	3.141592653589793	1.729 ns	0.0258 ns	0.0215 ns	1.725 ns	171 B
Matt	0	3.141592653589793	6.399 ns	0.0416 ns	0.0389 ns	6.377 ns	861 B
GammaCO	0	3.141592653589793	9.272 ns	0.1411 ns	0.1319 ns	9.243 ns	1,127 B
Cauldron	0	3.141592653589793	6.831 ns	0.1054 ns	0.0935 ns	6.793 ns	839 B
MyMethod	0	3.141592653589793	6.461 ns	0.1199 ns	0.1001 ns	6.428 ns	3,137 B
Normal	-0	4503599627370497	1.935 ns	0.0624 ns	0.0553 ns	1.898 ns	155 B
Lerp	-0	4503599627370497	1.756 ns	0.0361 ns	0.0337 ns	1.748 ns	163 B
NormalFma	-0	4503599627370497	1.918 ns	0.0078 ns	0.0073 ns	1.918 ns	152 B
LerpFma1	-0	4503599627370497	1.736 ns	0.0239 ns	0.0224 ns	1.741 ns	168 B
LerpFma2	-0	4503599627370497	1.937 ns	0.0112 ns	0.0094 ns	1.933 ns	161 B
HalfLerp	-0	4503599627370497	2.004 ns	0.0669 ns	0.0771 ns	1.972 ns	171 B
Matt	-0	4503599627370497	18.473 ns	0.1117 ns	0.0990 ns	18.482 ns	786 B
GammaCO	-0	4503599627370497	9.091 ns	0.2059 ns	0.2022 ns	8.980 ns	1,127 B
Cauldron	-0	4503599627370497	15.236 ns	0.1784 ns	0.1669 ns	15.225 ns	839 B
MyMethod	-0	4503599627370497	38.659 ns	0.7987 ns	1.2898 ns	38.805 ns	3,453 B
Normal	1	3.141592653589793	1.902 ns	0.0107 ns	0.0095 ns	1.903 ns	155 B
Lerp	1	3.141592653589793	1.757 ns	0.0459 ns	0.0430 ns	1.746 ns	163 B
NormalFma	1	3.141592653589793	1.928 ns	0.0138 ns	0.0116 ns	1.926 ns	152 B
LerpFma1	1	3.141592653589793	1.697 ns	0.0076 ns	0.0071 ns	1.698 ns	168 B
LerpFma2	1	3.141592653589793	1.934 ns	0.0195 ns	0.0182 ns	1.932 ns	161 B
HalfLerp	1	3.141592653589793	1.906 ns	0.0112 ns	0.0100 ns	1.907 ns	171 B
Matt	1	3.141592653589793	9.061 ns	0.0412 ns	0.0386 ns	9.061 ns	786 B
GammaCO	1	3.141592653589793	9.119 ns	0.0355 ns	0.0332 ns	9.125 ns	1,127 B
Cauldron	1	3.141592653589793	24.717 ns	0.0667 ns	0.0624 ns	24.742 ns	1,036 B
MyMethod	1	3.141592653589793	7.273 ns	0.0195 ns	0.0183 ns	7.275 ns	3,340 B
Normal	2.718281828459045	3.141592653589793	1.925 ns	0.0209 ns	0.0196 ns	1.928 ns	155 B
Lerp	2.718281828459045	3.141592653589793	1.716 ns	0.0143 ns	0.0127 ns	1.717 ns	163 B
NormalFma	2.718281828459045	3.141592653589793	1.916 ns	0.0059 ns	0.0049 ns	1.916 ns	152 B
LerpFma1	2.718281828459045	3.141592653589793	1.710 ns	0.0086 ns	0.0076 ns	1.711 ns	168 B
LerpFma2	2.718281828459045	3.141592653589793	1.933 ns	0.0103 ns	0.0096 ns	1.929 ns	161 B
HalfLerp	2.718281828459045	3.141592653589793	1.923 ns	0.0113 ns	0.0100 ns	1.921 ns	171 B
Matt	2.718281828459045	3.141592653589793	50.068 ns	0.1239 ns	0.1159 ns	50.067 ns	770 B
GammaCO	2.718281828459045	3.141592653589793	9.125 ns	0.0359 ns	0.0319 ns	9.124 ns	1,127 B
Cauldron	2.718281828459045	3.141592653589793	2.298 ns	0.0111 ns	0.0093 ns	2.293 ns	1,007 B
MyMethod	2.718281828459045	3.141592653589793	2.318 ns	0.0295 ns	0.0261 ns	2.310 ns	3,328 B

グラフにするとこんな感じです。

Normal とか Lerp 系列はほとんど差がない一方で、正確な計算をしようとすると結構時間がかかることが分かります。
また、正確な計算をしようとすると GammaCO 以外は引数によってかなり速度にばらつきがあることが分かります。

例によって私見を交えて 🌟 にまとめるとこんな感じです。

手法	開区間	オーバーフローしないか	表現可能精度	パフォーマンス
普通のやつ	✖️	✖️	🌟	🌟🌟🌟🌟🌟
Lerp 式	✖️	✅	🌟	🌟🌟🌟🌟🌟
FMA Lerp (1)	✖️	✖️	🌟	🌟🌟🌟🌟🌟
FMA Lerp (2)	✖️	✅	🌟	🌟🌟🌟🌟🌟
FMA Lerp (3)	✖️	✅	🌟	🌟🌟🌟🌟🌟
Half Lerp	✖️	✅	🌟	🌟🌟🌟🌟🌟
Matt 式	✖️	✅	🌟🌟🌟🌟	🌟
GammaCO	✅	✅	🌟🌟🌟	🌟🌟
Cauldron	✖️	✅	🌟🌟🌟🌟	🌟🌟🌟
オレオレ手法	✖️	✅	🌟🌟🌟🌟🌟	🌟🌟

パフォーマンス第一の場合は Lerp 式にしておくのがよさそうです。ただ、そのままだと意図せず開区間にならない問題があるので、範囲をはみ出したら再抽選する処理も入れておくとよいと思います。
精度がほしい場合はオレオレ手法も検討してみてください。

余談

擬似乱数生成器の選定について

本稿では擬似乱数生成器部分は理想的な実装であるものとして詳しくは触れませんでしたが、ここで挙げた「正しい」手法を実行するためには、かなりの「均等分布次元」が必要になります。

「均等分布次元」というのは、ここでは任意のビット列が出力される保証のあるビット長を表します。
例えば、メルセンヌツイスタ (32 bit 版) は 623 次元均等分布であることが有名ですが、これは 32 bit ワード単位での話なので、ビット長的には 623 x 32 = 19936 bit までなら任意のビット列を生成できます。
(おそらく、 bit 単位なら 19937 bit まで行けるとは思うのですが、自信がないです)

逆に言えば、均等分布次元を超えたビット列を出力した場合、絶対に出現しえないビット列が存在することになります。
実例を挙げると、線形合同法 (64 bit) は 64 bit ワード単位で 1 次元に均等分布するので、任意の 64 bit ワード、例えば 0ul は出力できますが、それ以上の長さを持つビット列は出現しない可能性があります。例えば、 [0ul, 0ul] (128 bit) は理論上絶対に出現しないと断言できます。

$[0.0, 1.0)$ の範囲でさえ、最小値が $2^{-1074}$ である、つまり 0.000～(1069個の0)～001 を生成できなければいけない以上、 1074 bit もの均等分布次元が要求されます。
そうなると xoshiro256++ (192 bit 均等分布) や xorshift1024* (1024 bit) などでは対応できず、メルセンヌツイスタ (mt19937_64) などの巨大な擬似乱数生成器を利用するか、事実上均等分布次元が ∞ である CSPRNG (暗号論的擬似乱数生成器) を使う必要が出てくるでしょう。
しかし、その場合当然生成自体が遅くなるという問題が出てきます。難しいですね。

逆に考えて、例えば xoshiro256++ を使うなら 192 bit ぶんの均等分布しか得られない前提に立って出力関数を設計する (つまり、 $2^{-192}$ までしか出ない前提で高速化パスを作る。非正規化数にならない前提で組むなど)、というのもあります。 Abseil の実装はこのあたり考えられていそうですね。
要はバランス、です。

並列化の可能性

float 型のある範囲のそれぞれ独立した一様乱数が 2 つほしいシチュエーションが時々あります。
具体的には、単位円内・単位球面上に一様分布する座標の取得だったり、ボックス＝ミュラー法による正規分布乱数への変換などがあります。

並列化と言えば SIMD ですが、ややこしいのでいったん置いておいて、もっと簡単な手法を考えます。

例えば、 $[0.0, 1.0)$ の乱数が 2 個ほしいとき、 $2^{-24}$ を掛ける手法なら 1 つあたり 24 bit の乱数があればよいのですから、 Next() が一度に 64 bit の乱数を生成できる以上、一度に 2 個分作れるわけですね。

ulong r = Rng.Next();  
float f1 = (r & ((1ul << 24) - 1)) * (1f / (1 << 24));  
float f2 = (r >> 40) * (1f / (1 << 24));  
return (f1, f2);

次に、さっき置いておいた SIMD です。
C# なら System.Numerics.Vector256 などからアクセスできますね。

// 4 つの [min, max) double 値を得る  
var vec = Vector256.Create(Rng.Next(), Rng.Next(), Rng.Next(), Rng.Next());  
var zeroOne = Vector256.ConvertToDouble(vec >> 11) * (1.0 / (1ul << 53));  
var result = Vector256.Create(min) + (max - min) * zeroOne;

ただ問題があるとすれば、高精度な手法は必然的に「再抽選が必要になる可能性」をはらんでおり、そうなると並列化が至難というところです。
なので、 Abseil の手法のように乱数固定消費で現実的な精度を出せる手法を使うことになるでしょう。

再現性について

シミュレーションやゲームにおいて、リプレイのために再現性が必要になる場合があります。
整数乱数の場合は完全に再現可能と言えますが、浮動小数点数乱数の場合はどうでしょうか？

IEEE 754 に準拠している場合 (現代における大半のプログラミング言語は該当します) 、以下の演算はどの環境でもビット単位で同じ結果を返す保証があります：

加減乗除 (+, -, *, /)
剰余 (Math.IEEERemainder) ※ % とは異なるので注意
平方根 (Math.Sqrt)
FMA (Math.FusedMultiplyAdd)
整数との相互変換
CopySign ・単項マイナス演算子・絶対値など、数値演算を含まない関数
比較 (==, < など)

さて、逆に言えばこれ以外の演算 (特に Math.Sin や Math.Log といった数学関数) を使用した場合は完全に同じ値が返る保証はないということになります。
なぜかというと、無限に正確に計算しようとすると計算コストが馬鹿にならないためです。多少の誤差と引き換えに高速に計算できるようにしてあります。
(大抵の場合、誤差は 1 ulp 以下です。したがって、仮数部の最下位ビットに ±1 を足した程度のわずかな誤差ではあります。)

C# の場合、例えば Math.Sin にはこのように書いてあります：

This method calls into the underlying C runtime, and the exact result or valid input range may differ between different operating systems or architectures.
(訳) このメソッドは基盤となる C ランタイムを呼び出すため、正確な結果や有効な入力範囲はオペレーティングシステムやアーキテクチャによって異なる場合があります。

で、その C ランタイムのほうはどうなのかというと、こちらに記述があります：

In most cases, the result produced is within +/-1 ULP (unit of least precision) of the correctly rounded result, though there may be cases where there's greater inaccuracy.
(訳) ほとんどの場合、生成される結果は、正しく丸められた結果の +/-1 ULP (最小精度単位) 以内ですが、不正確さが大きくなる場合もあります。

話を乱数に戻しましょう。
ビットパターンを利用して (整数演算のみで) 組み立てた乱数については、再現性が担保されています。
ただ、例えばこれを指数分布に加工するために -Math.Log(r) などとしてしまった場合、必ずしも同じ値が得られるとは限りません。
Java には異なる環境でも同じ値を返す保証がある StrictMath がありますが、 C# には今のところありません。なので、もし必要なら数学ライブラリを自力で実装する必要があります。険しい！

これは個人の考えにはなりますが、楽観的過ぎても悲観的過ぎても難しいので、「要はバランス」ということにしておくとよいかと思います。
例えば、乱数のシード値だけではなくて定期的に実際の値を送ったり保存したりして同期させる、など。

余談ですが、 Unity の BurstCompile 属性では、数学関数の精度を指定できます。
ここで低精度 Low を指定すると、なんと誤差が最大 350.0 ulp (!) になります。

またまた余談なのですが、 C# の仕様に恐ろしげな文章を見つけました：

浮動小数点演算は、演算の結果の型よりも高い精度で実行される場合があります。

一部のハードウェアアーキテクチャでは、double 型よりも範囲と精度が広い「extended」または「long double」浮動小数点型をサポートしており、すべての浮動小数点演算をこのより高い精度の型を使用して暗黙的に実行します。

この文章が正しいとすると、 x87 FPUの呪いのような事象が発生しうる (異なる環境で浮動小数点数演算の再現性が担保されない) ということになります。
幸いにして x64 だと発生しないっぽい？ので、現実的な範囲で問題になることは少ない……かも……
ちなみに、 .NET Core 2.0 以降なら常に RyuJIT が使われて SSE(2) 命令を使うので呪われないらしい、です。

Unity? IL2CPP? Burst? なんもわかりません。これだけで記事ができそう。

$\pi$ テスト - 精度の必要性を検証する

さて、ここまで高精度な一様分布浮動小数点数乱数の実装について説明してきましたが、この精度は本当に必要なのでしょうか？
別に「普通のやつ」でいいじゃん速いし、と言われれば、実例がないと言い返すのは難しそうです。
なので、実際に差が出るのかどうか、確かめてみましょう。

ここでは、円周率 $\pi$ の値を古典的なモンテカルロ法で求めてみることにしましょう。
手法をざっくり説明すると、

一辺 $2l$ の長さの正方形の中にランダムに点を打つ
それが半径 $l$ の円の中に入っていればカウントする
以上を繰り返すと、 (円の中の点数 / すべての点数) は $\pi / 4$ に近づく

というものです。
コードで説明するとこんな感じです。

// 誤差が分かりやすいように、乱数生成は Half でやるものとする  
  
// とりあえず l = ネイピア数とする  
// (なんでもよいが、無理数だと分かりやすいので)   
Half l = Half.E;  
  
long trial = 1ul << 32;  
long count = 0;  
  
for (long i = 0; i < trial; i++)  
{  
    // お好みの手法で [-l, l] の範囲の乱数を生成  
    Half x = Rng.NextHalf(-l, l);  
    Half y = Rng.NextHalf(-l, l);  
  
    // 円の中に入っていればカウント  
    if ((double)x * (double)x + (double)y * (double)y <= (double)l * (double)l)  
    {  
        count++;  
    }  
}  
  
// うまくいけば PI (に近い値) が出るはず  
// なお Math.PI == 3.1415926535897931  
Console.WriteLine($"pi = {4.0 * count / trial:g17}");

まず、オレオレ手法 MyMethod では pi = 3.140543058514595 を得ました。小数点以下 2 桁まで合っています。
一方、普通のやつ Normal では pi = 3.1403892487287521 を得ました。小数点以下 2 桁まで合っています。
……全然変わりませんね……先行研究にはこれで差が出るとあったのですが……
先行研究での試行回数はたった 8000 回らしいので、乱数の上振れ下振れの可能性がまぁまぁあります。闇。

ここで示されたこととしては、精度を上げても結局そんなに変わりはないということです。かなしい。
もし差が出そうなセンシティブな実験をされている方は、ぜひオレオレ手法のほうを検討してみてください。

ただこれは私見ですが、モンテカルロ法で必要になるのは精度よりスピードなのでは？というのがあり、そうなると精度って実はいらないんじゃね？という悲しい結末を迎えそうです。はい。

おわりに

「 完璧な 一様分布の浮動小数点乱数が欲しい」というただそれだけのためにどれだけ頑張る必要があるか、伝わっていれば幸いです。
一様分布ですらこれなので、ほかの分布はもっと大変でしょうね……今は考えたくもないです。

「オレオレ手法」は使っても使わなくてもいいです。
ただ Rng.NextDouble() * (max - min) + min するときは範囲外にはみ出す可能性があることを忘れないように。それだけは心に置いておいてください。

*1:Downey, A. B. (2007). Generating Pseudo-random Floating-Point Values. cit. on, 90. https://allendowney.com/research/rand/downey07randfloat.pdf

*2:Goualard, F. (2022). Drawing random floating-point numbers from an interval. ACM Transactions on Modeling and Computer Simulation (TOMACS), 32(3), 1-24. https://hal.science/hal-03282794v2/file/rand-in-range.pdf

*3:Goualard, F. (2022). Drawing random floating-point numbers from an interval. ACM Transactions on Modeling and Computer Simulation (TOMACS), 32(3), 1-24. https://hal.science/hal-03282794/file/rand-in-range.pdf ; https://frederic.goualard.net/publications/corrigendum.pdf

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).

はじめに

ソース

詰まった点など

無線だと、そもそもジャイロとかのデータが取れない

ジャイロ・加速度・タッチの座標が勝手に補正される

ジャイロの単位がわからん

おわりに

はじめに

の範囲において

で割る手法

開区間にならない (1.0 が出る)

表現可能な値すべてが出ない

ちょっと遅い

で割る手法

開区間になる

表現可能な値すべてが出ない

ちょっと速い

余談：閉区間にするには

余談： .NET での実装

仮数部ビットパターン法

開区間になる

表現可能な値すべてが出ない

パフォーマンス

余談: dSFMT について

余談: 別の区間への応用

exponent+mantissa 法

開区間になる

表現可能なすべての値が出る

遅い

改造してみる

Abseil 実装

開区間になる

表現可能なすべての値は出ない

パフォーマンス

Downey による補正

開区間にならない

対策 1 : clamp する

対策 2 : 再抽選する

表現可能なすべての値が「均等に」出現する

遅い

「さいきょうの」手法

開区間になる

表現可能なすべての値が出る

パフォーマンス

random_real

開区間にならない

表現可能なすべての値が「ほぼ均等に」出現する

パフォーマンス

MarcDense

実装の説明

開区間になる

表現可能なすべての値は出ない

パフォーマンス

Perfect

実装の説明

開区間になる

表現可能なすべての値が出る

パフォーマンス

まとめ

の範囲において

普通のやつ

開区間にならない

望ましくない結果を返す可能性がある

確率が均等にならない

パフォーマンス

標準ライブラリにおける実装について

Lerp 式

開区間にならない

望ましくない結果を返さない

確率が均等にならない

単調増加にならないことがある

パフォーマンス

FMA Lerp

パフォーマンス

HalfLerp

開区間にならない

望ましくない結果は返さない

確率が均等にならない

パフォーマンス

Matt 式

$[0.0, 1.0)$ の範囲において

$2^{64}$ で割る手法

開区間にならない (`1.0` が出る)

$2^{53}$ で割る手法

対策 1 : `clamp` する

`random_real`

$[min, max)$ の範囲において

$2^{n}$ (仮数部が 0) のとき

符号部が同じ - `min` が非正規化数の場合

符号部が同じ - 指数部が 1 異なる - `min` が正規化数の場合

符号部が同じ - 指数部が 1 異なる - `min` が非正規化数の場合

符号部が異なる - 指数部が同じ - `min` と `max` が非正規化数

符号部が異なる - 指数部が同じ - `min` と `max` が正規化数

$\pi$ テスト - 精度の必要性を検証する

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

$R^{2}$ の計算

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

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