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

はじめに

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

ランダムソート

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

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

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

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

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

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

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

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

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

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

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

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

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


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

危険なランダムソート


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

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

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

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

Fisher-Yates shuffle

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

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

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

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

擬似乱数まわりの最適化

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

擬似乱数生成器の選定

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

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


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

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

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

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

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

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

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

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

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


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

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

乱数を絞りつくす

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

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

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

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

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

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

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

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

n を求める

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

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

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


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

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

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

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

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

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

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


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

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

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

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

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


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

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

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

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

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

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

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

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

Fisher-Yates shuffle を行う

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

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

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

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

結果

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

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

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

小手先の高速化

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

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

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

例えば、

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

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

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

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

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

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

上限を削る

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

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

タプルでの交換をやめる

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

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

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

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

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

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

SIMD 化

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

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

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

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

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

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

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

n のキャッシュ

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

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

Fisher-Yates shuffle の誤った実装例


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

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

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

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

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

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

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

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

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

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

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

MergeShuffle

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

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

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

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

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

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

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

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

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

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

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

具体的な流れとしては、

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

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

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

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

おわりに

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

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

Fast shuffle by batching

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

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