2026-07-31

System.Random を使うな～ Slay the Spire 2 の乱数偏りバグから得られる教訓

はじめに

ある日、この記事が目に留まりました。

『Slay the Spire 2』にて数学者が「乱数の偏りバグ」をガチ分析し修正に至る。絶対に入手できないカードもあったバグ、“不運”は気のせいじゃなかった - AUTOMATON

曰く、 Slay the Spire 2 という有名なローグライクカードゲームにおいて、

本来均等にランダムな確率で得られるはずのカードが、偏った確率で得られる
- よりによって悪い効果のカードが多く出現する状態だった
結果として入手確率が 0% (入手不能) になったカードすら存在する

というものです。
しかもこれらは「実験・体感による偏り」ではなく、「理論的に求められた結果」だというのです。

気になったので、少し調べてみました。
なお、 (私が未プレイということもあり) ゲーム自体の内容についてはほとんど触れません。擬似乱数生成器の性質のほうに焦点を当てていますのでご了承ください。

Slay the Spire 2 乱数偏りバグの概要

ここでは Slay the Spire 2 の乱数偏りバグの概要を説明します。
詳細は元の記事を参照してください。

Correlated randomness in Slay the Spire 2 - Andy Tockman

Slay the Spire 2 において、乱数が関わるイベントは以下のような感じで初期化されているそうです：

int seed = /* 乱数 */;  
  
// 各イベント用の乱数生成器を作成  
var eventRng1 = new Random(seed + hash("event1"));  
var eventRng2 = new Random(seed + hash("event2"));  
var eventRng3 = new Random(seed + hash("event3"));  
// ...

ここで、 hash() はハッシュ関数 (入力に依存して何らかの int 値が決定的に得られる関数、 GetHashCode() みたいなもの) です。
また、擬似乱数生成器 Random は、標準の System.Random を使っているそうです。

それで、これらの eventRng* から得られる乱数を使ってランダムなイベントが発生するわけなのですが、その内容が偏っているというのです。
先ほどの日本語記事から引用すると、

地下水路ルートで「ネオーの骨」を取得した場合、プレイヤーは54.25%の確率で「負債」を、40.32%の確率で「腐敗」を得る。
また、繁茂の地ルートの場合は73.74%の確率で「苦悩」、18.85%の確率で「羞恥」を得る。
提供される呪いカードのプールは10種類ある

10 種類の呪いカードから均等に選ぶ場合、ある特定のカードが選ばれる確率は当然 10% になるべきですが、なんと 54.25% (!) の確率で特定のカードが選ばれていたそうです。
しかもよりによって強力なデバフを持つカードが付与されていたらしく、プレイヤーが本来の想定よりも不利になっていたのだとか。

また、同様の確率の偏りによって抽選確率が 0% となり、実質入手不可能となっていたカードもあったそうです。

こちらの提供プールにも確率の偏りがあり、「リバウンド」を取得できる確率は0%となっている。

なお、オリジナルの記事が公開されてから 4 日後に修正パッチが配信され、別の擬似乱数生成器 (xoshiro256**) に置き換えられる形でこれらの問題は修正されました。対応が早くて素晴らしい！

`System.Random` の概要

System.Random は皆さんご存知、 C# (.NET) 標準搭載の擬似乱数生成器です。
System.Random の内部実装については以下に詳しいです。（手前味噌）

かなり長いので要約すると、

アルゴリズムは古く、得られる乱数の品質が低い
- シードなしのコンストラクタおよび Shared については .NET 6 で改善されたが、 シードありのほうは依然古く悪い実装
絶対値が同じシードは同じ乱数列を生成する
初期状態のパターン数が少ない
遷移関数の実装ミスにより周期が計算できない (理論が破綻している)
アルゴリズムが安直で Next(max) ・ Next(min, max) の結果が偏る
NextDouble() の精度が低い
NextBytes() の結果が偏る・遅い
シリアライズできないため、乱数列の再現が難しい
拡張も難しい
乱数性テストで速やかに失敗する
得られた乱数から内部状態が復元可能で、過去の逆算・未来の予測が容易である
全体的に遅い
消費メモリも比較的多い

といった感じです。ひどいありさま。
標準搭載の実装ならある程度まともであるはずだ、という先入観を大きく裏切る状態となっています。

`System.Random` シードの相関

今回問題になったのは、「 System.Random に異なるシードを与えたにもかかわらず、結果に相関が生まれる」という現象らしいです。

論より証拠、まずは以下のような簡単なプログラムを実行してみましょう。

// 適当な定数  
int hashA = 1357924680;  
int hashB = 1470369258;  
  
for (int i = 0; i < 1024; i++)  
{  
    // マスターとなるシード値 (適当に決める)  
    int seed = Random.Shared.Next();  
  
    // 子となる乱数生成器 A / B (マスターシード値に定数を足す)  
    var rngA = new Random(seed + hashA);  
    var rngB = new Random(seed + hashB);  
  
    // A, B の最初の乱数を出力  
    Console.WriteLine($"{rngA.NextDouble()}\t{rngB.NextDouble()}");  
}

シード値が異なる (Random のコンストラクタに異なる値を渡している) わけですから、出力はそれぞれランダムに、つまり相関がないことが期待されますね。

実行すると、以下のようなテキストが出力されます。

0.04242573727035231  0.032940200079670266  
0.16686262617207254 0.1573770889813905  
0.17914017950144603 0.18862571669212808  
0.10575448773137969 0.09626895054069765  
0.09201346435212227 0.08252792716144022  
0.24937726894830226 0.2398917317576202  
0.6580073254453053  0.6674928626359873  
0.39119978826083235 0.40068532545151436  
0.43487830107793135 0.4253927638872493  
0.4942762146211584  0.5037617518118405  
……

なんとなく比例しているような傾向がみられますね。 A と B が似たような値になっています。
これをグラフに描画してみましょう。横軸が A 、縦軸が B です。

完璧に比例関係にあるわけではないですが、何本かの線の上に乗るように点群があることが分かります。

適当な定数 (hashA, hashB) を変えても同じようになるでしょうか。定数にランダムな値を使ってグラフを描画すると以下のようになります。

やはり、何本かの線の上に乗るように点群がありますね。
加えて、 $y=x$ を軸として左右対称になっていそうです。

乱数生成にちょっと詳しい方であれば、「乱数の初期化直後の値は品質が低いから、何個か乱数を消費して《暖めて》おく」みたいな噂を聞いたことがあるかもしれません。
これを試してみましょう。

int hashA = Random.Shared.Next();  
int hashB = Random.Shared.Next();  
  
for (int i = 0; i < 1024; i++)  
{  
    int seed = Random.Shared.Next();  
    var rngA = new Random(seed + hashA);  
    var rngB = new Random(seed + hashB);  
  
    // 乱数を 100 個消費する  
    for (int warmup = 0; warmup < 100; warmup++)  
    {  
        rngA.NextDouble();  
        rngB.NextDouble();  
    }  
  
    // 101 個めの乱数を表示  
    Console.WriteLine($"{rngA.NextDouble()}\t{rngB.NextDouble()}");  
}

驚くべきことに、傾向は変わりません。
つまり、初期化時に生まれた相関は 101 回目でも引き続き生きている、ということになります。
試してみた感じ、 warmup が 1000 でも同じ雰囲気でした。こわいね。

なお、《暖めて》おくとよい、という噂については、一昔前の内部状態の大きい擬似乱数生成器 (メルセンヌツイスタとか) はともかく、現代に設計された擬似乱数生成器 (xoshiro256** など) では当てはまりません。内部状態すべてをまんべんなく乱数で埋めることができれば、最初の 1 回めから正しく動作することが期待されます。

一応ですが、まともな乱数だとどうなるかを確認してみましょう。
Seiran128 の簡易実装を使ってみましょう。

// Seiran128 擬似乱数生成器  
public class Seiran  
{  
    private ulong state0, state1;  
  
    // ほんとうはよくない (入力が 64 bit に対して内部状態が 128 bit あるため；詳細は後述)  
    // シグネチャを合わせるため致し方なし  
    public Seiran(ulong seed)  
    {  
        state0 = seed = seed * 6364136223846793005 + 1442695040888963407;  
        state1 = seed = seed * 6364136223846793005 + 1442695040888963407;  
  
        Next();  
    }  
  
    public ulong Next()  
    {  
        ulong s0 = state0, s1 = state1;  
  
        ulong result = BitOperations.RotateLeft((s0 + s1), 29) + s0;  
  
        state0 = s0 ^ BitOperations.RotateLeft(s1, 29);  
        state1 = s0 ^ s1 << 9;  
  
        return result;  
    }  
  
    public double NextDouble()  
    {  
        return (Next() >> 11) * (1.0 / (1ul << 53));  
    }  
}

これで上述のプログラムの Random を置き換えてみると、以下のグラフが得られます。

相関はなく、十分にランダムになっているように見えますね。

どうしてそうなるのか

どうしてそうなるのかを考えてみましょう。

初期化シーケンスの問題 - 線形性

System.Random の初期化処理 (new Random(Seed)) は GitHub で確認することができます。
必要な部分だけ抜粋してコメント付きで解説します。

// 内部状態  
private int[]? _seedArray;  
private int _inext;  
private int _inextp;  
  
private void Initialize(int seed)  
{  
    Debug.Assert(_seedArray is null);  
  
    int[] seedArray = new int[56];  
  
    // シード値が負の場合、絶対値を取る  
    // (→ 絶対値が同じシードは同じ乱数列を生成する)  
    int subtraction = (seed == int.MinValue) ? int.MaxValue : Math.Abs(seed);  
    int mj = 161803398 - subtraction;  
    seedArray[55] = mj;  
    int mk = 1;  
  
    // ラグ可変の Fibonacci 法らしきもの  
    int ii = 0;  
    for (int i = 1; i < 55; i++)  
    {  
        // [1..56] の 55 要素を使う ([0] は使わない)   
        if ((ii += 21) >= 55)  
        {  
            ii -= 55;  
        }  
  
        seedArray[ii] = mk;  
  
        // mk = (mj - mk) mod (2^31 - 1)  
        mk = mj - mk;  
  
        if (mk < 0)  
        {  
            mk += int.MaxValue;  
        }  
  
        mj = seedArray[ii];  
    }  
  
    // Next とはラグの異なる Fibonacci 法で 4 周する  
    for (int k = 1; k < 5; k++)  
    {  
        for (int i = 1; i < 56; i++)  
        {  
            // n = (i + 30) mod 55  
            int n = i + 30;  
            if (n >= 55)  
            {  
                n -= 55;  
            }  
  
            // seedArray[i] = (seedArray[i] - seedArray[((i + 30) mod 55) + 1]) mod (2^31 - 1)   
            seedArray[i] -= seedArray[1 + n];  
            if (seedArray[i] < 0)  
            {  
                seedArray[i] += int.MaxValue;  
            }  
        }  
    }  
  
    // この時点で seedArray[1..56] は 0 <= x <= 2^31 - 2 の範囲  
    _seedArray = seedArray;  
  
    // ラグの初期化 (※初期化バグの原因、別件)  
    _inext = 0;  
    _inextp = 21;  
}

非常に雑に言うと、初期化処理は

seedArray[i] = (a - b) mod (2^31 - 1)

を繰り返している感じです。

英語のソース記事曰く、「出力は abs(Seed) に対して線形である」ということになります。
線形である、というのは、 x * S + y (x と y は整数定数、 S はシード) で表現できるということです。

まず最初の mj = 161803398 - subtraction は線形です。 -1 * subtraction + 161803398 と言い換えることができます。
次に、 seedArray[i] = (a - b) mod (2^31 - 1) について、これも線形です。
a と b が線形のとき、その引き算は (x1 * S + y1) - (x2 * S + y2) = (x1 - x2) * S + (y1 - y2) となり、線形になるためです。
数学的帰納法みたいですね。

初期化バグの原因、と書いたのは、ここでは触れない別のバグの原因となっているコードです。
詳しくはこちらを参照してください。
雑に説明すると、これのせいで理論が破綻して、周期計算ができなくなっていて、品質が保証できない状態になっています。

遷移関数の問題 - やっぱり線形性

続いて、 Next() にあたる処理 (InternalSample()) はこのようになります。

internal int InternalSample()  
{  
    Debug.Assert(_seedArray is not null);  
  
    // locINext = (_inext + 1) % 55 みたいな感じ (※ 1 基点)  
    int locINext = _inext;  
    if (++locINext >= 56)  
    {  
        locINext = 1;  
    }  
  
    // locINextp = (_inextp + 1) % 55 みたいな感じ (※ 1 基点)  
    int locINextp = _inextp;  
    if (++locINextp >= 56)  
    {  
        locINextp = 1;  
    }  
  
    int[] seedArray = _seedArray;  
    // 結果の乱数  
    // (seedArray[locINext] - seedArray[locINextp]) mod (2^31 - 1)  
    int retVal = seedArray[locINext] - seedArray[locINextp];  
  
    if (retVal == int.MaxValue)  
    {  
        retVal--;  
    }  
    if (retVal < 0)  
    {  
        retVal += int.MaxValue;  
    }  
  
    seedArray[locINext] = retVal;  
    _inext = locINext;  
    _inextp = locINextp;  
  
    return retVal;  
}

すごく簡略化すると、

_inext = (_inext + 1) mod 55;  
_inextp = (_inextp + 1) mod 55;  
  
return seedArray[_inext] = (seedArray[_inext] - seedArray[_inextp]) mod (2^31 - 1);

みたいな感じです。要するに引き算なので、線形です。

というわけで、初期化関数および遷移関数は線形であることが分かりました。
では、線形だと何が悪いのでしょうか？

線形な演算は x * S + y で得られると書きました。
ここで、 S を d 変えた (S = S + d した) 場合を考えてみると、 x * (S + d) + y になりますね。これは、もともとの値から x * d だけ大きい値になります。 (mod による影響は考えないものとします。)
それで何が問題かというと、何度遷移した (Next() を呼んだ) としても常にこの関係性が維持されることです。
つまり、片方の乱数生成器から a が得られたとき、もう片方の乱数生成器からは a + x * d が必ず得られます。
これをグラフに描くと直線 $y = x + z$ ( $z$ は定数) みたいになるので、それであのグラフが得られた、というわけですね。

相関すると何がいけないのか

相関が起きてしまうとどうなるのかというと、「特定のイベントが起こった時、別のイベントが起こるかどうかが推測できてしまう」といった実害が生じます。

例えば、以下のようなコードがあったとき、

var rngA = new Random(seed + 123456789);  
var rngB = new Random(seed + 987654321);  
  
// イベント A の判定 (10%)  
if (rngA.NextDouble() < 0.1)  
{  
    Console.WriteLine("Event A");  
}  
  
// イベント B の判定 (10%)  
if (rngB.NextDouble() < 0.1)  
{  
    Console.WriteLine("Event B");  
}

本来ランダムな事象は相関が全くないはずなので、 Event A が出力されるかどうかと Event B が出力されるかどうかは全くランダムで互いに関係しないはずです。
しかし今回のように相関が起こると、「Event A が出たときは Event B が出やすい」とか、あるいはその逆のように、偏りが生じてしまうことになります。
また、制作者が意図しない形でプレイヤーにヒントを与えてしまうことにもなります。

別件 1 - シードの数が少ない

これは別件ですが、シード値の空間 (とりうるパターンの数) は高々 31 bit ぶん ( $2^{31}$ ) しかありません。
これによって、一般のご家庭の PC でさえ数秒で全探索 (すべてのシードについてシミュレーションして、都合のいいシードを探すこと) が可能になっています。
一般にこれは擬似乱数生成器として (あるいはゲームとして) 望ましい性質ではないでしょう。

現代において、シード値は最低最悪でも 64bit 、理想を言えば擬似乱数生成器の内部状態のサイズと同じだけ必要です。 Seiran128 なら 128bit 、 xoshiro256** なら 256bit ですね。
シード値としてよく使われる日付や時刻は 32 bit にも満たない場合が多い (値自身が 64 bit だったとしても、そのうち有効な情報量はだいぶ小さいです - 例えば年のビットは 1 年に 1 回しか変わりませんよね？) ので不適切です。
ちゃんと CSPRNG (RandomNumberGenerator など) を使って、擬似乱数生成器の内部状態すべてをまんべんなくランダムに初期化する必要があります。

別件 2 - 再開時の挙動について

これもまた別件ですが、 Slay the Spire 2 では、乱数のシードとその消費個数を記憶していて、シードで初期化したあと消費個数ぶんだけ乱数を捨てる (消費する) することでプレイの再開を行っていたそうです。消費個数が少ないうちはなんとかなりますが、多くなると破綻するのでゲーム性によりそうな方式ですね。
これは System.Random がシリアライズできない (状態の保存ができない) ため生じている問題です。

得られる教訓

System.Random 、特にシード付きコンストラクタは使ってはいけません (MUST NOT) 。
乱数の品質が低く、上述のような意図しない問題を生じる可能性が高いです。

.NET 6 以降については、シードなしのコンストラクタ new Random() や、 Random.Shared を使う分には上記の問題は発生しないので大丈夫です。
ただし、それ以前 - 具体的には現行の Unity 環境とか .NET Framework とか - ではシードなしのほうも古い実装になるので危険です。
ただし (2) 、仮に Random.Shared を使っていたとしてもその実装が正しい・高速である保証はありません。「上記の問題は発生しない」と書いたように、別件の問題はいくつか残っています。
あと別件ですが API がしんどいです。 Next() が int.MaxValue を返さないとか、 NextUInt64() 相当のメソッドがないとか、 IRandom インタフェースがないとか……

乱数の再現性が必要な場合は、擬似乱数生成器を自分で実装するのが良いでしょう。
シリアライズできない問題も自前実装なら解決できるでしょう。

擬似乱数生成器のおすすめは Seiran128 ……と言いたいところですが、採用実績を考えると xoshiro256++ あたりがおすすめです。
xoshiro256** (シードなしのコンストラクタで使われているアルゴリズム) でもいいですが、 ++ のほうが新しくて品質がより良く、内部状態復元攻撃にも強いとされています。

xoshiro256** の出力関数は可逆のため、出力から内部状態が簡単に逆算できてしまいます。
xoshiro256++ の出力関数は不可逆 (2 つの状態から 1 つの出力を生成しています) なので逆算が簡単ではなく、いわゆる「乱数調整」への対策ができます。
(簡単ではないだけで不可能ではない、という点に注意は必要ですが……完全に乱数調整を防止したい場合は CSPRNG (C# であれば RandomNumberGenerator クラス) を使ったほうが良いでしょう。)

ただ、「確率的に正しい」実装をするのはなかなか困難ではあります。注意しましょう。
例えば、以下のコードでは確率的に偏りが生じるのはご存知ですか？

int NextInt(int max)  
{  
    return Next() % max;  
}

正しい実装の一例はこんな感じです。

int NextInt(int max)  
{  
    ulong r = Next();  
    ulong maxUlong = (ulong)max;  
  
    if (r < maxUlong)  
    {  
        ulong mod = (0 - maxUlong) % maxUlong;  
        while (r < mod)  
        {  
            r = Next();  
        }  
    }  
  
    return (int)((r * maxUlong) >> 32);  
}

詳細は "Lemire bound method" とかでググってみてください。

本当は自前実装にもいろいろ罠があって難しいのですが、ここでは省略します。

じゃあなんかいい感じのライブラリはあるのか、というと、難しいです。
Math.NET Numerics はアルゴリズムがあまりに古すぎ、 NRandom は速さの為か正確性を犠牲にしている実装になっているので……
なにか良い実装のライブラリがあれば教えてください。文句をつけたんだから自分で作れ？それはそう。今作ってます。ですがいつできあがるか分からないので……

おわりに

標準ライブラリにあるからといって品質が高いとは限らないの、 C 言語の rand() を思い出すようですね。
あと「令和版のカルドセプトサーガのダイスのバグみたいだな」なんてことも思いました。

とにかく new Random(seed) は使わないこと (もっと言えば System.Random 自体を使わないこと) 、もし使いたくなったら自前実装すること。これだけ覚えておいてください。

2026-07-28

まもけんR 村長の試練 RTA メモ

はじめに

#まもけんR pic.twitter.com/7iOXL1uxJ8
— Andante (@andanteyk) 2026年7月27日

例によってチャートを忘れそうなので書いておく。
この記事に書いてあるものはすべて体感であり、実際の設定とは異なる可能性があることに注意。

攻略スプレッドシートへのリンクはこちら。

基本戦略

作戦

とにかく即降り。
近くに民家があったとかでない限り、階段を見つけ次第降りる。

40F あたりまでにコンボパーツが揃うことを祈る。

強化の杖
一時しのぎの杖 (できれば 2 本)
氷柱の杖
カカオ菓子
それなりに強い剣 (2 印)
それなりに強い盾 (炎印)

操作系

直線ダッシュと i ダッシュ、階段ダッシュをそれぞれ使えるようにする。
基本は i ダッシュが速いが、制御が難しい場合があるのでそういう時は直線ダッシュを使う。
階段ダッシュは特に便利、明かりの本や天照でマップを読んで階段直行、が手動より速くできる。

まめちしき

焚火の部屋・店内には基本的に罠がない。「罠増しの本」を読まない限りは。

食事

「ダイエットバーガー」がおすすめ。
(12000T のあいだ SP 減少 2 倍・最大 HP +30 ・攻撃力 +15)
最大 HP も上がるので安定感がだいぶ上がる。加えて、「生きているエビ」より効果ターンが長いので深層でも切れない。
また、「生きているエビ」よりもだいぶ安い、という微妙な利点がある。本作は銀行のキャパが低いので……
ただし、パンを一つも無駄にしないで進行する勢いでないと餓死する。
（空腹状態だとダッシュできなくなるので死ななかったとしてもタイムが死ぬ）

他の択としては定番の「生きているエビ」か。
(10000T のあいだ攻撃力 +20)
こちらは拾い食いでよくなるので対処アイテムを多く持てる利点がある。

武器印・称号

武器本体は印枠が 2 個以上あれば基本なんでもいい (基本値が高いに越したことはないが) 。重要なのは印と称号。

理想を言えば「2 金回4 回4」、最低限「2 回4」とかか。

「2 回攻撃」 (x0.75 x 2回) は言わずもがな最優先。火力が 1.5 倍になる。これを入れるためにほかの赤印を入れないよう注意が必要。
「かなてこ」 (HP 満タンの敵に x1.2) も一撃で倒せるラインが格段に上がるので良い。ただ「2 回攻撃」とのシナジーは微妙。
「会心」 (確率で x1.5?) は運任せにはなるが一撃で倒せる敵が増えるので良い。
「回復 4」 (攻撃時 HP+5) は回復リソースが不足しがちな火山地帯で輝く。「2 回攻撃」とシナジーがある。もし「2 回回」できれば殴るたびに 20 回復でき、なんなら被ダメージより上回る場合も。

「錆よけ」は貴重な印枠を使うほどではないが、こんぼう +3 が手に入った場合は修正値目当てで入れるか……？ぐらい。印枠 6 とかなら。
「種族特効」は基本的には不要な (あえて入れるほどではない) イメージ。記憶が正しければ倍率が x1.25 とそこまで高くないため。基本値が強い「ドラゴンキラーィ」や「砕きのハンマー」のついでで入るぐらいか。印枠が余っていれば。
状態異常の果実が識別できていれば「混乱攻撃」「睡眠攻撃」あたりもうれしい。

まもけんシリーズは基本的に修正値 (加算) よりも特効倍率 (乗算) が強いので、修正値目当てで要らない印を入れるのは本来あまりよくない。

「疾風剣ハルピュイア」「バールのようなもの」「回復の剣」「騎士の槍」など、つよそうな武器は 26F 以降？にしか落ちていない。ただしそれ以前に店売り・変化のポーチで出る可能性はある。

称号は、

最高
- 孤高の
- 魔術師の
効果あり
- やせ我慢の
- 英雄の
- 火事場の
- 小食の
- 杖宿りの
- 殲滅の
- 杖絞りの
- 冒険者の
- 熱を帯びた
効果なし (活かせる機会が少ない)
- カリスマな
- 活食の
- 恵みの
- 成長の
- 凍てついた
- 霧を斬る
- 傭兵の
- 淀みの
- 結晶砥石の
- 守銭奴の
- 旅人の
デメリットあり・癖がある
- 一点集中の (盾なしで x2.0)
- 金溶の (100G 消費で +10)
- 大金溶の (1000G 消費で +50)
- 剛投の (投擲がダメージに)
- 博打の (x0.5 or x2.0 、期待値的には上だが…)
- 波動の (HP 20% 消費で x0.75 の波動)

「魔術師の」 (杖ごとに火力増加) が安定して強い。保存があまり拾えず杖を大量に持ち歩くことが多いので。
「金溶の」 (100G 消費してダメージ +10) も強いが、盾と比べて消費が激しいのがネック。腕輪を大量に売れば 30000G ぐらいにはなるのでなんとかなるか？
「孤高の」 (ひとりのとき x1.2) はほぼ常に効果が発動するのでありがたい。実質かなてこ。

「熱を帯びた」 (暑いエリアで火力増加・炎無効) は火山エリアが結構長いので役立つほか、盾の炎印を省略できるのでありがたい。

「一点集中の」はワンチャンチャートある……か……？そもそもそんなに見たことない。
もし引いた場合は火力は十分なので槍・灯印も視野に入る。
終盤のドラゴンに焼かれて死ぬ未来が見えるので、そのへんは長期休暇でカバーできれば。

ボーナス効果は「錆よけ」「状態異常攻撃」あたりがあるとうれしい。

民家がそんなに存在しないうえ釜がないこともしばしばあるので、合成のチャンスは限られている。できれば無駄にしたくないところだが、実際そううまくはいかない……

武器スキルは以下のスキルだと嬉しい。

「天照」は文句なしに強い（タイム削減に寄与する）が、チャージが遅く夜間は使えない問題がある。
マナ充填の魔法書があれば実質明かりの本がわりになる。

「一文字斬り」はモンハウからの脱出やダメージゾーンの突破、水路に囲まれたアイテムの回収など応用が幅広い。
射程範囲が全部壁や敵で埋まっていると不発になるので注意。

「絶剣」は一応ダメージが高い (攻撃力 2 倍かつダメージ倍率 x2 倍? 要検証) ので使えなくもない。ただ外したら終わり。
ワンチャン大好きなあなたに。予防状態なら気兼ねなく撃てる。

「装魂解放」は「水晶の剣」さえ引ければ強力なスキル。 11x11 範囲内の敵に消費した剣の攻撃力の 5 倍のダメージ (水晶の剣なら攻撃力 30 なので 150 ダメージ) 。盾の場合はそれが回復に変わる。
それ以外だと次点の「ドラゴンキラーィ」でさえ 60 ダメージ程度なので、アイテムの引きに依存する。
なお、消費する武器選択をキャンセルするとスキルゲージが初期化されるので注意。

「爆炎剣」は敵を吹き飛ばせるが自身にもダメージがあるのがネック。爆発耐性と組み合わせたい。
アルファフレイムリリィなど対策しようがない敵を実質無料で消去できるのはえらい。

「氷結剣」は凍結が便利。水晶がわりにしたり足止めしたり凍結破壊で大ダメージを狙ったり……
なお、斬撃部分はミスする可能性があるが、凍結効果は必中。

「魔迅剣」は射程無限で通常攻撃の x1.5 のダメージ。印効果や称号効果は乗らない。
チャージも早く使いやすい。

「無明」はまもけん2 と違って耐性で無効化できるようになった。したがってスキュラ系やワーバット系などには効かないが、自分も目つぶし耐性のボーナス効果を得ていれば一方的に有利になるようになった。
モンハウで同士討ちを誘うことができる。ただし戦闘演出にちょっと時間がかかるのがネック。

盾印・称号

理想的には「炎爆マヒ」とかか。最低でも「炎」は必要。
盾は印枠が少ないので厳選が大事。変な印を入れている場合ではない。本作では盾スキルが印枠を食うようになったので、かなりカツカツ。

「炎」印がアイテム焼失を防ぐために絶対に必要。ないとドラゴン系に焼かれて死ぬ。ほかにはレッドファングの「炎の息」や変化キツネ七尾以降の「浮遊する狐火の矢」も炎属性。
「爆」は大型地雷や大爆発の魔法書による事故を防げるほか、「爆発のポーチ」「爆炎剣」で自爆して敵を吹っ飛ばす戦術が使いやすくなる。ただ、こっちは最悪なくてもいい。

「風の加護」はあると強力だが SP がゴリゴリ持っていかれるので餓死の危険性が見える。ダイエットバーガー採用時はなかなか厳しい。
「重」は空腹が加速するほか、基本的に一人旅なので投げ技で仲魔と分断される問題が起きにくくあまり活用できないかも。
「ダメージ反射」は結構役立つこともあるが、それで倒すと経験値が得られなくなるので一長一短。

「錆よけ」は入れられる枠がないので、神器のボーナス効果についていることを祈る。
「鈍足無効」は前作より使い手が少ないので重要度が下がったが、ボーナス効果であればありがたい。
「混乱無効」のほうが嬉しいかも？
「毒無効」も結構役に立つが、印として入れるほどではないのでボーナス効果でついていればありがたいぐらい。

意外なところでは「麻痺無効」が有効。強化いちしのチャートをやるなら必須。
ただ異種合成か神器のボーナス効果にしかないので入手が難しい。麻痺の果実を識別できたら持っておきたい。

「強護の盾」「風の加護の盾」「魔法反射の盾」「回復の盾」は 26F 以降？にしか落ちていない。
「強護」は風の加護に比べると消費コストがなく、エルフ系の矢やダメージ床が軽減できるので便利っちゃ便利。入れ得ではある。

称号は、

最高
- 孤高の
- 金溶の (100G 消費で -10)
- 魔術師の
効果あり
- ダメ押しの
- 火事場の
- 活力の
- 錯乱の
- 小食の
- 熱を帯びた
- 本宿りの
- 枕靴下の
- 万全の
- 冷静な
- 杖宿りの
- 爆森人の
- 英雄の
- 殲滅の
- 冒険者の
- 七転八起の
- 斥候の
効果なし
- カリスマな
- 快感の
- 空腹硬化の
- 守銭奴の
- 聖なる
- 淀みの
- 旅人の
- 二人三脚な
- お祝いの
- 危機管理
- 快眠枕の
- 逆境の
- 一点集中の
- 傭兵の
- 活食の
デメリットあり・癖がある
- 臆病者の (被弾時に攻撃ボタン長押しでテレポート、 SP -10)
- 剛投の (投擲がダメージに)
- 正面装甲の (正面からの被ダメ軽減、それ以外からの被ダメ増加)
- 博打の (x0.5 or x2.0 、期待値が通常より悪くなる)
- 爆弾魔
- 氷結導師の
- 雷神導師の
- 風塵導師の
- 回復導師の
- 感染の (被ダメ時自分自身の状態異常をコピー；倍速や予防も移してしまう)
- 不動の (移動していなければ被ダメ軽減、移動していれば被ダメ増加)
- 大金溶の (1000G 消費でダメージ -50)

「金溶の」 (100G 消費でダメージ -10) がかなり強力。剣と違って被弾することはそう多くないので消費も少なめ。
「孤高の」「魔術師の」も剣と同様条件を満たしやすく強力。

「冷静な」 (混乱魅了無効) は安心できる。催眠の魔眼も無効。
「枕靴下の」はリソース純増。開幕モンハウで足元の罠を起動しても泣かない。
「熱を帯びた」なら炎印を省略可能。

「○○導師の」をメイン盾とは別に持っておいて未識別本の漢識別用に使うのもアリ。特に回復が不足するので「回復導師の」だとうれしい。
ただ、メイン盾をそれにするのはさすがにやめたほうがよさそう。

「臆病者の」 (被弾時テレポ) は理論上強い、強いのだが、実際は誤爆が多いので微妙。
あと SP を消費するのでさらに満腹度管理が大変になる。

ボーナス効果は「錆よけ」「状態異常無効」あたりだと嬉しい。

盾スキルは基本的にあまりこだわりがないが、しいて言えば「加速の風・霧」「回復の風・霧」あたりだと結構使える。

「癒しの風」「癒しの霧」は状態異常回復。ただ、やばい状態異常は得てして行動不能なため、使いどころが限られる。
呪いや火傷の回復に使えるかも。

「おみやげ転送」は「カカオ菓子」や「青い森にんにく」などを生成できるチャンス。まぁおまけ程度に。

「回復の風」「回復の霧」は HP+50 。火山地帯あたりで地味ながら役立つ。

「加速の風」「加速の霧」は倍速効果。モンハウに突っ込んだ時に時間稼ぎ (スキル演出中に次の手を考える) として。

「狐様の祝福」は HP 回復や満腹度回復に地味に役立つ。

「全員転移」は瞬間移動の実がわりになる。ただテレポート後隣接して殴られる可能性も結構あるので祈る。

「石壁の守り」は自分で発動しないといけないので面倒だが、一発でも食らうとアウトな終盤で地味に役立つ。

腕輪

識別できないことが多いのであまり頼っていない。資金源。
序盤 (～30F) は優先的に拾って換金し 20000～30000G 程度持つようにしているが、金溶でもない限りはそこまで稼ぐ必要はないかも。
呪いが怖くなければとりあえず装備して「千里眼の腕輪」か「盗賊の腕輪」を狙うのも手か。デメリットのある腕輪は「易疲の腕輪」「爆弾魔の腕輪」ぐらいなのでワンチャンに賭けてもいいかもしれない。私はやってないけど……
あとアイテム枠が足りないので、例えば「ちからの腕輪」は「麻痺の杖 [5]」よりも効果があるか？は考えたほうがいいかも。

「ちからの腕輪」は唯一 2000G なので識別しやすく、誰にでも効果があるので良い。ただ、与ダメージ増加は +2.25 程度なので役立つかどうかは微妙。
印を重ねられれば +4.5 ダメージなので強い……？

「千里眼の腕輪」は敵が見えるほう。当然ながら強い。ただ識別リソースが足りない…… 漢識別でもわかりやすいので試してみる？
おそらく 26F 以降にしか落ちていないので、序盤に拾った腕輪は適当に名前を付けておくと絞り込みやすい……かもしれない。
「盗賊の腕輪」はアイテムが見えるほう。モンハウが見えるのは強い。識別は上に同じ。

「錆よけの腕輪」は宿敵ヌメックの印消し対策ができる。あとスラグの -3 も痛いのでうれしい。ヌメックゾーンが終われば基本的には不要。

「瞬間移動の腕輪」は出たらかなり便利そう。ただし使うたびに SP-5 されるので無限に使えるわけではない。
「爆発の腕輪」も同様。

価格帯としては (太字は確定)

2000: ちから
3000: 状態異常よけ・魔物寄せ
5000: 綺麗・持久・易疲・爆弾魔・会心・回復・値切り・呪いよけ・アイテム保護・盗賊・爆発・瞬間移動・通過
10000: 命中・しあわせ・遠投・へた投げ・魅了よけ・錆よけ・千里眼・手放さず
15000: 娘守り

矢

あまり使わないが、あればそれなりに便利。
状態異常矢は杖の代わりに使える。まとめることができるので便利。
常に識別状態で祝呪もなく、キツネが化けていることもない。信頼できる。（命中不安を除けば）

道中のハイエルフが何本かドロップするので、鈍足の矢などを拾っておけると良いかも。

「爆発の矢」をモンハウに撃ち込むと全員目を覚ます。
本作では爆発が起きると部屋内の仮眠状態の敵が起きる仕様があるため。

食料

「ダイエットバーガー」を使うチャートの場合、絶対に一個も無駄にしない戦略が求められる。
「腐ったパン」だろうと男らしく食していくしかない。焚火に投げて焼くと美味しい。

「巨大なパン」は SP+110 。トルネコとかより回復量が低いので注意。

なんか深層ではパンの床落ち率が低い？気がするので、拾い食いに頼らずできるだけ持っていく。
祝福されていれば空間を圧縮できる。

果実

「復活の実」があればワンチャンできるので加速になる。効果発動よりも呪われたり燃やされたりすることのほうが多いかも。
1F から落ちていることもある。

「瞬間移動の果実」はモンハウ脱出用に。
「麻痺の果実」は盾に麻痺耐性をつけるため。果実拡散のポーチを持っていればモンハウ対策にもなる。

民家の青釜で合成可能なものはほぼ状態異常系の実。未識別なら武器に合成するのが無難。

「物忘れの果実」 (500G) を絶対に食べないようにしたいところ。
基本的には 26F あたりから出始めるが、低層の店売りで存在することがあるので注意。 500G だったら投げて識別するべき。
万が一食べてしまった場合、それ自体を食べた記憶も忘れるのでログから名づけをしておかないともう一度食べることになる。

価格帯は、

2: しなびた
10: 腐った
100: 小さい回復
200: 回復・瞬間移動・ちから毒消し・ちから
300: 倍速
400: 大きい回復
500: 状態異常系・物忘れ
1000: 特大の回復・命・封印・予防
1500: 復活・癒し・しあわせ・ふこう・無敵

200G の実は瞬間移動の可能性があるため店内で食べないこと。
あと地味なところで、憤怒の果実を民家で食べてしまってそのまま泥棒、という事故もしばしばある。もらってから食べよう。

杖

要らない杖は罠チェックにも使える。 (魔法弾の軌道上の罠は発見状態になるため。炎の魔法書の火炎弾なども同様。)
また、ほとんどの杖は壁に斜めに当てると 2 回まで屈折するのでトリックプレイも可能。あまりやったことはないけど……

以下の杖はアイテムに振った時にも効果がある。ただし店の商品には無効化されることに注意。

癒しの杖 (呪い・封印回復)
一時しのぎの杖 (階段付近にテレポート)
トンネルの杖 (パンに振ると穴あきパンに変化)
吹き飛ばしの杖 (12 マス吹き飛ばし)
場所替えの杖 (場所替え)
瞬間移動の杖 (フロアのどこかにテレポート)
飛びつきの杖 (飛びつき)
引き寄せの杖 (引き寄せ)
物知りの杖 (識別)
強化の杖 (確率で武器・盾・杖・本の修正値 +1 。杖は [30] まで、本は [3] まで)
弱化の杖 (確率で武器・盾・杖・本の修正値 -1 。 [0] まで)

「眠りの杖」「麻痺の杖」など、敵を確実に止められる杖は強い。
その意味で「混乱の杖」は安定感に欠ける (本作の混乱はスキルも撃つので信頼するのは危険) 。ただマーメイド系を混乱させると回復の旋律が自分にも当たるようになる。

明らかに安い (100G 以下の) 杖は「道しるべの杖」。敵に当てた場合何のメッセージも表示されない。
階段の方向を知れるほか、経路上の罠を発見できる。

「癒しの杖」は主に呪いを解くために使用する。アイテムが呪われている確率は結構高めなので 1 本は持っておきたい。

「一時しのぎの杖」は当然ながら非常に強力、 1F スキップできる。「氷柱の杖」や麻痺耐性と組み合わせると即座に飛ばせる。
また、敵に対して使用するとテレポート先が一瞬マップに映るため、階段の位置を知ることができる。
おそらく 26F 以降にしか落ちていない。

「トンネルの杖」は行き止まりの通路が多くなる終盤に欲しくなる場合が増えるが、そこまで温存しておくかというと微妙な位置。
終盤で拾えたらうれしい、ぐらい。
あと敵に振ると 20 ダメージとなる。トルネコとかだと 10 ダメージなので一応注意。

「回復の杖」は絶対にキープしておくべき。特に火山地帯で回復リソースがいくらあっても足りないため。
状態異常系の杖は合成してしまうと使える回数が 1 回減る (投げてぶつける分) が、回復の杖は投げることはないので合成してもよい。合成のポーチが余っているなら。

移動系では「場所替えの杖」「飛びつきの杖」ももちろん強い。矢を撃って場所替えすることで遠距離移動が可能。あとは水路のアイテム救出など。
ただ「引き寄せの杖」はほとんど使いどころがない。水路に囲まれたアイテムの地形が存在しない？ため。かなしいね。
隣接して振った時に自分の 5 ダメージのログが早ければ飛びつき、相手の 5 ダメージが早ければ引き寄せ。

「物知りの杖」はアイテムを識別できる。店の商品には無効化されるため、一旦購入する必要がある。

「強化の杖」は強化状態ではなくアイテム強化に使う。「一時しのぎの杖」に使うのが最適だが、所持品枠を圧迫するようなら「回復の魔法書」やメイン剣などに使うのもよい。
コンボは後述。 1 本は [0] でもいいので持っておきたい。

「そんたくの杖」「いじわるの杖」は、「攻撃のみが強化される」「雷撃の麻痺が即座に切れる」「Lv1 の敵に振った時何もメッセージが表示されなかった」あたりで識別が可能。

「氷柱の杖」は通路を塞ぐ・凍結破壊によるダメージ増・魔法弾反射など、活用が幅広い。「一時しのぎの杖」と組み合わせると即座に階段に飛べる。
水路に落ちたアイテム救出にも使える。

26F 以降？、「体力交換の杖」が落ちている可能性があり、未識別の杖を適当に振ったら自分の HP が激減した……という事故もまれにある。
弱い敵にも殴り殺される可能性がなくはないので、回復手段があるときに振ること。
なお、体力交換の杖は中立 NPC には無効化される。

価格帯としては (太字は確定、実際の価格は単価 x (回数+1)) 、

10: 道しるべ
100: 癒し・隠れ身・石化・一時しのぎ・トンネル
120: 回復・毒・雷電
140: 混乱・封印
150: 移動系・体力交換・物知り
160: 強化・弱化・倍速・鈍足・麻痺・憤怒
180: 眠り
200: そんたく・いじわる・衰弱
250: 空振り・誘い・柱系

ベース価格方式ではないので価格計算が難しい。振って識別したほうが早いかも……
安い杖は狙い目。明らかに安いのは道しるべ確定だし、 500G 以下の場合一時しのぎ [2 ～ 4] の可能性がある。

本

「識別の魔法書」はできれば識別が困難で貴重なポーチに使えると嬉しい。
価格が安い (単価 100G) のでそれ自体の識別は簡単。

「氷の魔法書」は「氷柱の杖」同様、凍結させて魔法弾を反射する使い道がある。民家の釜や住人に向けて撃たないように注意。
水路を凍結させることもできる。

「油揚げ化の魔法書」は貴重な食料兼回復アイテムである「極上油揚げ」を生成できる。ただしメイン剣盾などに読まないよう注意。

「呪いの本」は 1000G の価格帯で唯一「どれを？」系なので識別可能。

「解呪の魔導書」は全アイテムの呪いを解くほか、呪い状態 (バステトと変化キツネ九尾が使う、確率で与ダメージ x0.5 ・被ダメージ x2.0) の解除にも使える。非常に危険な状態異常なので覚えておくとよい。

「回復の魔法書」は言わずもがな強力な回復リソース。封印状態でさえなければいつでも使える。
床落ちのは [2 ～ 4] だが強化する場合は [3] までしか増えないことに注意。

「大爆発の魔法書」は炎水晶や風水晶の近くで読むと即死するので一応注意 (5 敗) 。
爆発耐性さえあれば「爆発のポーチ」と同様の使い方ができる。

「明かりの本」はもちろん強力。階段ダッシュを使えば高速に移動できる。

「武器強化の本」「防具強化の本」は呪いを解く効果もある。特に序盤の武具や神器を適当に装備する際にあると安心なので、即座に読まずに持っておく手もある。

「メッキの本」は地雷。意図せず錆よけ印をつけてしまうとほかの印の枠がなくなってしまう。
しかも呪いを解く効果はない。厳しい。あと錆よけ印は複数つけることもでき、見事無駄になる。

「ポーチ強化の魔法書」を読むべきかは微妙。呪われた際に割って取り出せなくなるため。
強化されたポーチを敵にぶつけると足元に落ちて罠を起動するため、トリックプレイが可能。まぁやらないとは思うが……

「流星群落としの禁書」があればモンハウ (とくに 99F のやつ) を一発で打開できるが、終盤以降にしか落ちていない？ので、難しい。
グレムリンに読まれても泣かないこと。 (2 敗)

「嫉妬の禁書」は 26F 以降？から出現し始める。
これがあるので 26F 以降の本は迂闊に読めない……が、未識別の本を読まずにクリアするのは困難だと思われるので祈るしかない。
盾が育っていれば一発は耐えられることが多いので、階段が見えていない状態で読んでも諦めない。

「聖域の禁書」はトルネコなどと異なり、敵も乗れたり貫通する特技が多かったりであまり信用ならない。
貼りつく本はこれしかないので識別は簡単。一応「引き寄せの杖」など移動系の杖で動かすことができる。

「長期休暇の禁書」は個人的にはドラゴン系に読んでいる。ダークドラゴンを封じられるのは大きい。
他の候補としては変化キツネ七尾系、なめくじ系などか？
同じ部屋の敵が怒り状態になるので、一対一ではない場合は注意。
一応だが、投げても効果はないのでトルネコ・シレン経験者は注意。

「在宅勤務の禁書」は部屋内の敵とアイテムが消滅する。得られるものはなくなるが速やかに処理できるので便利。
タイムアタック的には 99F で読めると最高。

「白紙の本」は未識別状態で登場するうえ、出現確率が低いので 2 枚目が拾えることはほとんどない。
変化候補としては「長期休暇の禁書」「在宅勤務の禁書」「女神の気まぐれの禁書」あたりか。
なお、祝福状態の白紙の本に書き込むと変化後の本が祝福状態で得られる。別の本を書き込めたりはしないので注意。
また、「書き込む」を選択してキャンセルすると祝呪が分かるので、入手し次第やっておくとよい。

「女神の気まぐれの禁書」は「困った時の巻物」と言えば伝わるだろうか。とりあえず HP を全回復 (料理による増加分は含まない？) できる。
敵に囲まれていれば周囲 1 マスに永続の麻痺を与えられるほか (耐性を貫通するためケツァコアトルも麻痺させられる) 、 SP が 10 以下？なら SP 回復もできる。
何もないときに読むとお金が 2000-5000G もらえるが、普通に売却したほうが高い。
10000G の本はこれだけなのですぐ識別できる。

価格帯は (回数制の本は単価 x 回数)

100: 識別
500: 罠見え・召喚
600: 属性系
750: 油揚げ化
1000: 罠増し・罠消し・呪い・解呪・回復・大爆発・状態異常系
1200: 明かり・迷い
1500: 武器強化・防具強化・メッキ
2000: 予防・千里眼・鑑定・小流星・目つぶし・日照・集団転移・招集・アイテム寄せ
2500: マナ充填
3000: 吸出し・幻揚げ
5000: ポーチ拡張・ポーチ強化・流星群落とし・マナ乱れ・大雨・嫉妬・聖域・長期休暇・在宅勤務・白紙
10000: 女神の気まぐれ

ポーチ

「変化のポーチ」は変化先のアイテムテーブルが店売りと同じ？らしく、低層でもハルピュイアが出る場合がある。
これが未識別の場合は大事なアイテムを初手で入れないこと。
「転送のポーチ」はおそらく出現しないため、単価 250G のポーチは変化確定。

「保存のポーチ」は言わずもがな重要。
8F までに拾える「入れる」ポーチは保存のポーチ確定。
完走までに拾える保存は 1-3 個程度なので、大切に使うこと。
なお、ポーチ容量の上限は [6] 。ポーチ拡張が拾えた場合は無駄にしないように。

「回復のポーチ」はもちろん重要。呪いが解けない場合でも最悪叩き割って足元にマナを落とせば回復可能。
これも保存と同じく、 8F までに拾える「使う」ポーチは回復のポーチ確定。

「果実拡散のポーチ」は売却アイテムを詰めておいたり、パンを一時的に詰めておくのにも使える。
空きをつくっておいて「麻痺の果実」を入れればモンハウを制圧可能。

「爆発のポーチ」は使うと HP 半減の爆発を起こす。火薬壺と違って何度も使えて便利。
敵に囲まれたときに使って吹っ飛ばす、アイテム損失覚悟でモンハウに投げ込むなど……
投げるとマナが落ちた場所を中心に爆発が発生する。
爆発耐性があれば気兼ねなく使えるようになる。使用直後はワーラビットやドラゴンなどの遠隔攻撃に注意。

「識別のポーチ」は腕輪を入れるほか、「一時しのぎの杖」「氷柱の杖」など回数を把握しておきたい杖などを入れる。
単価 500G で容量 [2] だった場合はこれ確定。

「呪いのポーチ」「御狐様のポーチ」は単価 1000G 。どちらにせよアイテムロストする。
「御狐様のポーチ」は食料生産には使える。

「強化のポーチ」は概ね強力。 10F 以上持つこともあれば、 1F で即座に割れることもある。祈ろう。
入れる候補としては「一時しのぎの杖」かサブ武器・盾。回復の魔法書は入れても上限が [3] なのでもったいない。
割れた際に水路に落ちても泣かないこと。
ちなみに、ポーチ強化していてもフロア移動時の破壊は免れない。
未識別でフロア移動時に突然割れた場合は強化か弱化なのでログから名づけておくこと。

「弱化のポーチ」は弱化するわ突然割れるわでいいことがない。即座に破棄すべき。
杖を入れていた場合はむなしく空を切るかどうかで推測できたりもする。

「祝福のポーチ」に入れるのは理論上は回数制ではない本が良い。「眠りの魔法書」や「集団転移の魔法書」、禁書類、あと「復活の実」など。
食糧難が予想される場合はパンを入れる手もある。
呪われたアイテムを入れると祝福で上書きされるので、呪いを解くためにも使えなくもない。ちょっともったいないかも。

「合成のポーチ」は異種合成できるので、異種合成 (麻痺耐性とか) 優先で入れる。通常合成は民家の釜で。
ただし果実は後で入れること。先に入れると合成されない場合がある。
単価 2000G で [5] だった場合は基本これ。

価格帯は (実際の値段は単価 x 容量)

250: 変化
500: 保存・果実拡散・識別・使う系 (マナ込みなら 600)
1000: 呪い・誰か・御狐様
2000: 強化・弱化・祝福・合成

その他のアイテム

ギフトアイテムは投擲すると状態異常を与える。

「濁酒・村一番」は混乱 (10T) 。序盤のアルファ狩りに。
確率で殴られるほかスキルを撃ってくる可能性もあるので一応注意。

「天然またたび」は魅了 (10T) 。「誘いの杖」より効果時間が長い。
アルファに効かないのであまり活用できるところは少ないか。
一応モンハウなどで囮を作って逃げるのに使えなくもない。

「極上油揚げ」は HP+100・SP+10・空振り (10T) 。空振り効果が長くなり、アルファにぶつける戦略も取りやすくなった。
敵が周囲にいるときに食べないように注意。

「幻の油揚げ」は HP+200・SP+20・空振り (10T) 。
終盤の変化キツネ七尾系に握ってもらって食糧難を回避する限界プレイもある。

「青い森にんにく」は、ちから +1 。うれしい。一応封印状態だと食べられない。

「すべすべリンス」は空振り (10T) 。アルファ狩りに便利。

「ネイルセット」は目つぶし (10T) 。通路にいるアルファに使えるが、殴られる可能性もなくはないので注意。
5000G なので序盤の金策にも使える。

「カカオ菓子」は HP+50・SP+20 。
後述するコンボのキーアイテム。 1F から落ちていることがある。
食糧難や命の危機に食べることもできる。

「魔法のカギ」は 15000G で売れる。一応 90F のカギ部屋で本来の使い方をすることもできるが、それは通常プレイでの話。
あと変化キツネ七尾系を仲魔にして敵を握ってもらい、食糧難を回避するのにも使える。

「変化キツネ」は冒険者の不倶戴天の敵。
祝呪されず、合成できず、濡れず、凍らず、燃えない。水路に投げると消滅する。
敵が呪い攻撃・アイテム封印攻撃をしてきたときに何も起こらなかった場合は所持品にキツネがいる可能性がある。
フロアに敵が多い場合、姿を表そうとしたときに「まぼろしだった」と表示されて消滅する。モンハウを疑うこと。
矢やポーチ・ギフトには化けない。あと「復活の実」に化けていることはないはず。
一応「見渡す」で注視すると呼吸するように動いているため分かる。ただ RTA では悠長に観察している場合ではないので大抵の場合騙される。

「お金」はダッシュで乗ってから拾うとアイテムとして入手でき、投げ当てると金額の 1/10 のダメージ。
一応序盤の対アルファ攻撃用に取っておく手はなくはないが、本作の村長の試練では結構対処アイテムが落ちていることが多いのでお世話になることは少ないか。
むしろノーティ系に投げられて即死する事故のほうが多いかもしれない。
ちなみに最低 10G ・最高 1000G 。

強化杖使いコンボ

「強化の杖」をアイテムに振ると、確率で修正値が +1 される。 (武器・盾・杖・本が対象。杖は [30] まで、本は [3] まで。)
また、杖使い系が「杖拾い」・「杖使い」で振った杖は回数が減らない。
以上を組み合わせてアイテムを無限に強化するコンボ。

必要なアイテムは「強化の杖」・「カカオ菓子」・強化したいアイテム。
対象としては当然「一時しのぎの杖」。 [30] まで強化できれば 30F スキップできる。 2 本あれば 60F なので実質杖使いが出る階層以降すべてスキップ可能 (実際に 2 本引いたことはあるので不可能ではない……)
また「氷柱の杖」も重要。一時しのぎの反射のほか、深層の敵を一時的に足止めするのにも使える。
あとは「物知りの杖」や手持ちの武器盾か。 +30 の剣盾はやはり安定感が変わるが、時間との相談。

杖使いが出現するのは 40-44F (青杖使い) ・ 48-53F ・ 63-69F (紫杖使い) ・ 80-83F (赤杖使い) 。
「倍速魔法」を使える紫杖使いが戦力としても便利。火山地帯に入る前にコンボが実行できると比較的安全。
青杖使いでもよいが仲魔にした瞬間「怒り魔法」をこちらに撃ってくる可能性があるので危険。周囲に危険がないことを確認してから。即座に好感度を上げてスキルを封印する必要があるため、いらないアイテムが落ちている場所で勧誘するとよい。
赤杖使いがいるフロアでは安全性というものが皆無なのであまりお勧めできない。ただ物資が揃うのはこの辺かもしれないのでつらい。

なお、杖使いは呪われていたり [0] だったり変化キツネだったりする杖でも問題なく振れる。

余談だが、最近作者のねことかげ先生が「強化の杖でアイテムも強化できる」とつぶやくまで気づかなかった。私が最初に振った時は確率に負けて何も起こらなかったらしく、完全に意識の外にあった……

【まもけんR Tips2】
『まもけん』シリーズでは、杖などから射出される魔法弾はアイテムに当たり、特別な効果を与える事があります。
『まもけんR』で「強化の杖」は装備品の強化値を稀に上昇させる場合がありますので、
誰かにプレゼントする前に使い切ってしまいましょう！#まもけんR pic.twitter.com/gPgA9LVmnq
— ねことかげ@ゲーム開発中 (@nekotokage9653) 2026年7月18日

階層別戦略

1-4F

建物エリア。
特筆すべき敵も存在しない。アイテムの識別を進めよう。

このあたり (8F まで？) に落ちているポーチは「保存のポーチ」 (入れる) か「回復のポーチ」 (使う) 確定なので、拾えたら即名前をつけておく。

まもけん2 の「村長の試練2」と違って序盤のアルファがそこまで強くないので、お金は財布に入れて大丈夫。

「復活の実」が落ちていることがある。大抵胃袋に入ることになるが、祝福されていればラッキー、そうでなくても識別できてラッキー、ぐらい。

料理効果で十分戦えるので、武器や防具は現時点では装備しない選択もあり。呪いケアのため。

最序盤に落ちている「どれを？」の本は「識別の魔法書」「武器強化の本」「防具強化の本」「メッキの本」のどれかのため、適当な武具に読むと無駄がない。

5-8F

鉱山エリア。

5F からアルファが出現し始める。 8F は確定出現？
ここで拾った武具が案外最後まで活躍することもあるので倒しておきたい。
ただ、基本的に初期位置から離れたところに湧く（開幕詰み防止のため？）のでたどり着くまで時間がかかりがち。階段とアルファの距離が遠ければスキップもありか。

しいて言えば攻撃力が高めなワーウルフとアイテム弾きでポーチを破壊してくる可能性があるゴブリン系に注意、ぐらいか。
ワーウルフは少なくとも攻撃力 31 以下だと一発で倒せない。

このあたりから店が出現する。値段識別と腕輪の売却をしておく。
ちからの腕輪 (2000G) が識別しやすくて無難に効果がある。
店売りアイテムはテーブルが異なるらしく、 11F あたりで床落ちがまだないはずの「疾風剣ハルピュイア」「物忘れの果実」 (500G) が売っていることも。注意が必要。
店売りの 500G の果実は名前を付けておき、投げて識別する。店員に投げても即座に怒られることはない。本作は投げてもログから命名できるのでありがたい。

9-12F

赤い洞窟エリア。

保存と回復以外のポーチが出始める。変化のポーチの変化先もおそらく店売りテーブルなのでいいアイテム (または危険なアイテム) が出るかもしれない。

また、民家も出現し始める。今後合成ができるとは限らないので、できる限り使えると嬉しい。
念のため説明すると、住民の好感度が ❤3 以上であれば本棚やフルーツバスケットからアイテムがもらえる。アロマポットで合成も可能 (茶色は通常合成、青が異種合成対応) 。
使い終わった釜に火炎弾 (炎の魔法書・炎剣イフリートなど) を当てると使用回数が復活する。逆に氷の魔法書を当てると鎮火してしまうので注意。
あとで捨てる前提で適当に武器を束ねて修正値を集めておくのも手。
「大爆発の魔法書」「風の魔法書」などのケアのため、住民から離れて読むこと。あと「集団転移の魔法書」ケアのためちゃんともらってから読むこと。

13-16F

崖のある洞窟エリア。

しいて言えばラミアの「叩きつけ」に注意。火力が高いうえ実質 2 回攻撃なので、油断していると落とされることもしばしば。
16F から出現するハニービーの「脱力針」にも注意。永続的に火力を減らされるので痛い。
ちから 1 = 武器修正値 1.5 ≒ 0.75 ダメージと考えよう。

ノーティは地味に床落ちアイテムを投げて無駄に消費するので、余った場所替えの杖などがあればそれで対処していく。
混乱状態だと適当な場所に投げるので地形によってはロストする場合がある。
お金の上に乗っているやつは金額によっては即死しかねないので (額面の 1/10 のダメージ) 一応注意。

17-20F

遺跡エリア。

ピンクスライムの耐久力と魅了、フェンリルの火力と耐久力が若干脅威か。

21-25F

水路のある遺跡エリア。

このあたりから？モンスターハウスが出現し始める。ダッシュの勢いで突っ込んでしまうことも多い。
まぁ対処アイテムがなくても物理で殴って解決できることが多いので、むしろアイテム回収ができてありがたいかもしれない。

26F-30F

水の洞窟エリア。ここは外周階段フロアなので基本的に外周を回るように探索する。
このあたりから「疾風剣ハルピュイア」「バールのようなもの」といった強力なアイテムや、「嫉妬の禁書」といった危険なアイテムが出現し始める。

デルフィネが普通に強いほか、毒で HP を削られがち。ダッシュ時に衝突しないように気を付ける。

31-34F

暗い沼エリア。黒い霧がかかっていて見通しが悪くてつらい。外周階段フロア。
34F は高確率でアルファが出現する？

ここは民家とお店が多い印象がある（体感）。できれば探したいところだが、外周階段なので階段を見つけやすく、スキップしてしまうこともしばしば。
前のゾーンでハルピュイアを拾ったりしている場合は探すのも手かもしれない。

35-38F

キノコの森エリア。ここも外周階段フロア。

アマゾネスが強い。見通しの悪い通路から高火力の槍でちくちくしてくる。不慮の死亡事故が起こりがちなので注意。
あとメデューサの叩きつけでポーチを破壊されることもあるので、一撃で倒せる火力がないときは対処アイテムを切ったほうがいい。

39-42F

マナ乱れの森エリア。外周階段フロア。

マナが乱れているので杖識別がやりにくい。うっかり振って無駄にする事故が多発するので注意、ただ注意してても間違えるときは間違えるので一手早い回復を心がける。

スラグに注意。妙にタフなうえ、修正値を -3 されるとプレイヤーのやる気ごと -3 される。
ここだけのためにメッキするのはもったいない気がするので、錆よけの腕輪か神器のボーナス効果にあるとうれしい。

40F から青杖使いが出現し始める。上述のコンボの起点になる、が、ここはマナ乱れなので難しいかも。
というかここまでに資材が揃ったことがない。
ドラゴンとかスラグとかを怒り状態にしてくるのでしんどい。

ドラゴンがタフなうえ火炎ブレスで燃やしてくる。火炎耐性がない場合は注意。

43F-48F

要塞フロア。

フレイムリリィが攻撃力が高くて脅威。一撃で倒せるとうれしい。
もしアルファ化している場合は杖無効・投擲無効なので、眠りの魔法書などを切るか爆発で吹っ飛ばすしかない。

ブラウニーは床落ちの本を読んで事故を起こしてくる。モンハウに陣取ってアイテムをぽんぽん消費してしまうのでかなしい。

グリフォンは一撃で倒せないことが多いうえ、モンハウに連れ去られるなどの被害報告が相次ぐ。

49-52F

空中の要塞フロア。壁は掘れない。
火山地帯手前の良心のようなフロア。

紫杖使いが 48F-53F で登場する。ここでコンボしておきたいところ。
出現する敵の関係上ここでやるのが一番安定する。行き止まりの部屋でやると敵が侵入しづらくてよい。

変化キツネ出現フロアなので、アイテムを信用しないこと。投げて識別するのを怠らない。

コカトリスはさらにタフで、大抵の場合 2-3 発当てないと倒せない。例によってモンハウに誘拐されることがあるので注意。

53-59F

火山地帯。暑いので HP の自然回復が下がる。

ヘカーテに要注意。上位淫魔なので魅了解除まで 5 回行動が必要。魅了耐性があるなら仲魔にし、魅了祭りをして評判を上げる手もなくはないが時間がかかる。

ヴァンパイアは催眠の魔眼を撃ってくるので絶対に行動させてはいけない。幸い、ここまで来れるなら大抵の場合一撃で倒せる。

同じく変化キツネ出現フロア。アイテム (特に主力の剣盾・回復のポーチ) を封印されないように。

ベノムスライムがやばい。投石だけでも毒になるうえ耐久が高く倒しづらく HP を削られがち。
暑さもあり HP がゴリゴリ持っていかれるので、ここまででどれだけ回復リソースが拾えたかが生存のカギとなる。

ネモフィラにはできるだけちょっかいを出さないこと。睡眠どころか永眠させられる可能性がある。

60-72F

火山地帯深部。引き続き暑い。外周階段フロア。
引きが悪い冒険者は大体ここで死ぬ。

ここの序盤でもベノムスライム (56-63F) が厳しい。

ヌメックに本当に注意。 2 回攻撃や火炎耐性を消されて死んだ冒険者は多い。
かといって通路から突然出てくることもあるので本当につらい。錆よけの腕輪があれば安心できる。
しかもタフなので少なくとも 2 回殴らないと倒せない。混乱でごまかそうとするとおもむろにこちらを向いて印消ししてくるので注意 (n 敗) 。
空振りならスキルも外すので完封できる。

63-69F には紫杖使いが出現し、倍速魔法を撃ってくる。
倍速ヌメックとか倍速サキュバスライムとか考えたくない敵を量産してくるので可及的速やかに降りたい。
理論上はここで仲魔にしてコンボできなくはないが、敵が危険すぎて悠長に強化している場合ではない。
あとダークエルフがアイテムを吹っ飛ばす可能性がある。

後半はサキュバスライムが出てくる。 (66-72F)
ただでさえ投石が痛いのに魅了してくるので無理。復活があると安心。
火山地帯で通路が入り組んでいることも多く、敵が溜まっていてそこから延々投げてくる場合がある。

ハイエルフ (66-70F) は距離を開けて状態異常矢を撃ってくるため地味に厄介。
ただ、倒せば装填していた矢を落とすので利用できなくはない。対処アイテムを切って (あるいは燃やされて) 持ち物枠が空いていることもあるので、ありがたく使わせてもらおう。

ブルードラゴン (70-74F) はタフ。火炎耐性がない場合無限にアイテムが燃やされて死ぬことになる。
耐性さえあればまだマシな相手。

グレムリンに「流星群落としの禁書」を読まれて即死した走者がいるらしい。 (2 敗)
「本大好き小悪魔」で「罠増しの本」を読まれると厄介。せめてログをさかのぼって名前を付けておこう。

73-76F

要塞深部フロア。ようやく暑くなくなるがモンスターは激アツ (悪いほうの意味で)。

クラーケンが地味に物理が強いうえ、スプラッシュで本を濡らしてくる。ちゃんとポーチに入れよう。

このあたりから変化キツネ七尾 (75-79F) が出てくる。ポーチの外のアイテム (とくに復活の実) を握ってくるほか、「浮遊する狐火の矢」はなにげに火属性なのでアイテムを燃やしてくる。
限界冒険者はわざと不要アイテムを握ってもらい食料を確保するとかなんとか。
逆に、アイテムに変化するキツネは以降出現しなくなるので、床落ちアイテムの信頼度は上がる。拾っている暇はあるのか？

77-82F

青い要塞フロア。水路がある。

なんといってもバステト (78-84F) がヤバい。アイテムを次々呪ううえ、人間まで呪ってくる。
絶対に行動させてはならないが、通路で出会い頭に呪われることも多い。
人間が呪われると倒しきれなかったり突然二倍撃を受けたりして即死につながることもある。しかも効果時間が結構長い。

ゴブリンロード (80-85F) も危険。絶対に水路を背にして戦ってはいけない。メイン剣を失う羽目になる。

加えて、赤杖使い (80-83F) がヤバさを引き上げてくる。以上の魔物が透明になって襲いかかってくる。
ここでコンボしている暇はないと思われるので、カカオ菓子をもし温存している場合は HP 回復アイテムとして使ってしまってもいいかもしれない。

83-89F

黄色い要塞フロア。

カーバンクル (83-87F) が危険と言えば危険。例によって視界外から槍で突いてくる。
水晶代わりに使えるので一時しのぎや瞬間移動を振って逃げにも使える。

ニンフ (86-89F) は基本的に一発で倒せるのでそこまで脅威ではない。アルファにならないことを祈る。
一応周囲 1 マスにも杖無効効果があるのでそこだけ忘れないように。

ケツァコアトルを殴ると麻痺させられる。モンハウで「目つぶしの魔法書」を読むと敵が勝手に殴ってこちらが麻痺してそのまま永眠、という事故もある (1 敗) 。
また、例によってモンハウに拉致されることがある。早めの回復を心がけよう。

90-99F

帰ってきた建物フロア。
99F は固定大部屋モンハウ。

序盤のエキドナ (88-95F) が厳しい。特に通路で接敵すると「叩きつけ」でアイテムをばらまいたあげくロストすることがある。
保存のポーチが消えるとシャレにならないので気をつけよう。

ハイリザード (92-96F) ・ダークバット (90-95F) ・セイレーン (89-94F) の状態異常も危険。
本人の戦闘力はそこまでではないが、ほかの強力な魔物に絡まれたり囲まれたりすると危ない。

最後のエリアにはダークドラゴン (97-99F) が出る。説明不要のフロア射程 40 ダメージにアイテム燃焼つき。
なんなら 3 体ぐらいいることもある。寝ていることを祈ろう。
ここまで来て死にたくはないが火炎耐性がないと回復リソースごと燃やされて死ぬ。ヌメックの印消しがここで効いてくる。

ジャーゴイルは物理に強いほか、サンダーブレードで麻痺させてくる。麻痺している間にほかの敵が寄ってきてタコ殴りにされる。

あとは地味ながらケルベロス・変化キツネ九尾が呪ってくる場合があるので注意。

99F のモンハウ対策のため、部屋に効果がある本 (眠りの魔法書や混乱の魔法書など；目つぶしの魔法書は効かない敵がいる) or 果実拡散のポーチ＋麻痺の果実、あと移動系の杖を持っておく必要がある。
別に全滅させる必要はないのでさっさと降りてしまおう。

おわりに

個人的な難易度としてはまもけん2 の村長の試練＜まもけんR の村長の試練＜＜まもけん2 の村長の試練2 という感じ。
ちょうど良くて遊びやすいので、皆さんも試してみてください。

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 送りされるので、実行が絶望的に遅くなります。

おわりに

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

はじめに

Slay the Spire 2 乱数偏りバグの概要

System.Random の概要

System.Random シードの相関

どうしてそうなるのか

初期化シーケンスの問題 - 線形性

遷移関数の問題 - やっぱり線形性

相関すると何がいけないのか

別件 1 - シードの数が少ない

別件 2 - 再開時の挙動について

得られる教訓

おわりに

はじめに

基本戦略

作戦

操作系

まめちしき

食事

武器印・称号

盾印・称号

腕輪

矢

食料

果実

杖

本

ポーチ

その他のアイテム

強化杖使いコンボ

階層別戦略

1-4F

5-8F

9-12F

13-16F

17-20F

21-25F

26F-30F

31-34F

35-38F

39-42F

43F-48F

49-52F

53-59F

60-72F

73-76F

77-82F

83-89F

90-99F

おわりに

はじめに

ソース

詰まった点など

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

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

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

おわりに

はじめに

の範囲において

で割る手法

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

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

ちょっと遅い

で割る手法

開区間になる

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

ちょっと速い

余談：閉区間にするには

余談： .NET での実装

仮数部ビットパターン法

開区間になる

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

パフォーマンス

余談: dSFMT について

余談: 別の区間への応用

exponent+mantissa 法

開区間になる

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

遅い

改造してみる

Abseil 実装

`System.Random` の概要

`System.Random` シードの相関

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

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

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

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

対策 1 : `clamp` する

`random_real`

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