JDK 17 における擬似乱数生成と LXM ファミリについて

はじめに

2021/09/14 に GA となった JDK 17 では、 擬似乱数生成まわりの強化 が行われました。

更新点を見ると、 LXM という新しい擬似乱数生成器が実装されているようです。
しかし、このアルゴリズムについて聞いたことがなかったため (実装当時) 、調べてみました。

ついでに、このアルゴリズムが導入されるきっかけとなった JEP 356: Enhanced Pseudo-Random Number Generators についてもある程度解説しようと思います。

本記事では、特に注釈がなければ下記環境を使用して検証しています。

openjdk 17.0.6 2023-01-17
OpenJDK Runtime Environment Temurin-17.0.6+10 (build 17.0.6+10)
OpenJDK 64-Bit Server VM Temurin-17.0.6+10 (build 17.0.6+10, mixed mode, sharing)

JEP 356 の要約

そもそもどういった提案で何が追加されたのか、のざっくりとしたお話です。
原文は JEP 356 を参照してください。

旧 API では、擬似乱数生成器ごとの互換性が微妙だったり (SplittableRandom は Random を継承していないので代入できないなど) 、各生成器で nextDouble() などをそれぞれが独自に (コピペで) 実装していてよろしくなかったり、といった問題がありました。

また、 SplittableRandom には "弱点" が見つかっており、新しい別のアルゴリズムを導入して問題を回避したいニーズがありました。
JEP 内ではこの弱点について具体的に触れられていませんが、おそらく こちらの文献 *1 を指しているものと思います。
問題点の説明の部分を要約します。 SplittableRandom は内部に gamma という定数を持っており、 split() で新しいインスタンスを生成した際に自身と異なる gamma を設定することで別々の乱数列を得る仕組みになっています。この gamma は乱数の品質に深くかかわっており、 "弱い" gamma を設定すると品質が下がってしまいます。 Java では弱い gamma を弾くアルゴリズムが入っているのですが、それでも処理しきれない値 ( 2 ^ {64} / k や 2 ^ {s} m + 1 など) が存在しており、それが選ばれた場合に乱数の品質が下がってしまう、という問題がありました。

問題への対応だけではなく、ポジティブな面もあります。
擬似乱数生成器によっては split() や jump() といった異なる乱数列を得る仕組みを持つものがあり、並列実行において便利に扱うことができます。前者は SplittableRandom 、後者は xoroshiro128+ アルゴリズムなどが該当します。
この機能を統一的に扱うことができる API を用意したい、というニーズがありました。

インタフェースの追加

以下のインタフェースが提供されます。箇条書きは実装関係を示しています。

これらについて、ざっくりと解説をしていきます。 (網羅しているわけではないので、詳細はそれぞれリンク先の API ドキュメントから確認してください。)

RandomGenerator

擬似乱数生成器としての基底のインタフェースです。基本的に nextLong() を実装すれば nextDouble() などはデフォルト実装で生やしてくれます。

これのおかげで、擬似乱数生成器を独自実装する際に Random を継承する必要がなくなりました。不要な内部状態などがついてこなくなります。

StreamableGenerator

異なる乱数列を出力するインスタンスを無限に生成する rngs() をサポートします。

rngs() は Stream<RandomGenerator> 、つまり擬似乱数生成器のインスタンスが得られるストリームを返します。何に使うのかというと、 (特にマルチスレッドにおける) 並列実行です。

本質的に乱数生成は内部状態の変更を伴います。したがって、スレッド間でインスタンスを共有してしまうと、内部状態が正しく更新されず出力が崩壊したり、ロックで時間と手間がかかったりと良いことがありません。
(ちなみに Random はロックをかけるほうの手法を採用しており、代償としてオーバーヘッドがかなり大きかったのです。)
そこで、各スレッドで別々のインスタンスを利用することによって、速度と品質を同時に得ることができる、というわけです。従来の ThreadLocalRandom に近い雰囲気です。

ただ、各スレッド用の個別インスタンスを作る際、単純にクローンを作ったり適当なシードを与えて初期化したりした場合、運が悪いと出力される乱数列が重複してしまい、ランダムではなくなってしまう (相関が生まれる) 可能性があります。
これを回避するアプローチは 2 つあり、それらが SplittableGenerator, JumpableGenerator にあたります。

SplittableGenerator

従来 SplittableRandom で使えた機能の一般化である split() をサポートします。
(実際に SplittableRandom もこのインタフェースを実装するようになりました。)

内部状態とは別にある種の定数を持ち、この定数をインスタンスごとに変えることで乱数列自体を別々にするしくみを持ちます。
定数が異なれば、初期の内部状態 (≒シード) が完全一致していてもそれぞれ異なる乱数列が出力されるため、乱数の重複を回避できます。

ちなみに、現在 SplittableGenerator を実装している擬似乱数生成器が split() で定数をどうやって決めているのかというと、 nextLong() ベースでランダムに決定しているようです。
つまり奇跡が起こって親子で定数が一致する可能性も 0 ではありませんが、そこは楽観的に決定しているようです。
もっとも、そうなる確率は天文学的に低く、もしそうなっても内部状態が離れていれば実用上重複しないので、心配する必要はありません。
具体例を挙げると、 L64X128MixRandom で同じ乱数列に入る確率は 2 ^ {-63} 、かつ親との距離が \pm 2 ^ {64} 以内になる確率は約 2 ^ {-190} ですので、 GUID の衝突よりもはるかに低確率です。

JumpableGenerator

固定長ジャンプを高速に行う jump() をサポートします。

jump() は、 nextLong() を大量に ( 2 ^ {64} や 2 ^ {128} といったオーダーで) 呼び出したのと同じ状態遷移を高速に得られる機能です。
これによって、そのジャンプ距離と同じ回数 nextLong() などを呼び出さない限りは乱数列が重複しない保証を得る (そんな回数は現実的に呼び出せないので事実上独立した乱数列を得る) ことができます。

通常の乱数列を 1 本の長いテープに例えると、 SplittableGenerator は「異なるテープを何本も用意する」、 JumpableGenerator は「1 本のテープを等間隔に分割して使う」といったイメージです。
こう書くと JumpableGenerator のほうは 1 本あたりの長さ (≒周期; 使える量) が減ってしまって良くないように思えるかもしれませんが、切ったとしてもそれぞれ 2 ^ {64} 個分の長さはあって現実的に使いきれる量ではないので、心配は無用です。
また SplittableGenerator に比べると品質のばらつきが比較的少ない利点もあります。前述した「 SplittableRandom で gamma 値によっては品質が下がってしまう問題」のようなことが起こりにくいです。

LeapableGenerator / ArbitrarilyJumpableGenerator

JumpableGenerator の拡張として、 LeapableGenerator と ArbitrarilyJumpableGenerator も存在します。

LeapableGenerator は、 jump() よりも長距離の固定長ジャンプを行う leap() もサポートします。
ArbitrarilyJumpableGenerator は、さらに任意長のジャンプを行える jump(distance) もサポートします。

これらの見分けが難しいかもしれませんが、基本的にはジャンプ距離の違いと、距離を指定できるかどうかです。

  • JumpableGenerator : 中距離の固定長ジャンプ (例: 2 ^ {64})
  • LeapableGenerator : ＋長距離の固定長ジャンプ (例: 2 ^ {128})
  • ArbitrarilyJumpableGenerator : ＋任意長ジャンプ

LeapableGenerator の存在意義は、擬似乱数生成器の分割自体を並列化することにあります。
jumps() で n 個のインスタンスを生成するとき、 jump() を n 回連続で呼ぶ必要があり、 O(n) で時間がかかることになります。この時、最初に leap() でいくつかの子インスタンスを作り、そこから jump() で孫インスタンスを切り出すようにすれば処理時間が (n / leap() で作ったインスタンスの数) となり、最善を尽くせば O(\sqrt{n}) に短縮することができます。先のテープのたとえで言えば、最初にテープを何本かに切り分けて、それらを重ねてまとめて切っていくイメージです。

ちなみに、固定長のジャンプ距離については厳密な指定はなく、周期を p として jump() が 2 ^ {64} または \sqrt{p} ぐらい、 leap() が 2 ^ {128} ぐらい、といった目安が参考として記載されています。実用上は特に気にする必要はありませんが、もし擬似乱数生成器を実装する側になった場合は参考にしてください。

なお、 2023/03/17 時点では ArbitrarilyJumpableGenerator を実装している擬似乱数生成器はありません。

JumpableGenerator シリーズの注意点

注意点として、 JumpableGenerator において rngs() や jumps() で作成したインスタンスからさらに rngs() などで孫を作ってはいけないことが挙げられます。
jumps() は、内部的には

  1. 自身のコピーを作る
  2. 自身を jump() させる
  3. コピーを返す

といったことを繰り返しています。
ここで、子でも rngs() してしまうと最初の親と子の jump() 呼び出し回数 (距離) が一致し、擬似乱数生成器の内部状態がかなり近い (最悪一致する) 状態となり、結果として親子でほぼ同一の乱数列を永遠に出力し続けるようになってしまいます。
jumps() の戻り値が Stream<JumpableGenerator> ではなく Stream<RandomGenerator> になっているのはそのためで、孫の作成を自然と防止するようにできています。

孫が欲しい場合 (例えば、単一スレッドで複数の乱数列が必要なものをさらにマルチスレッド処理する場合など) は LeapableGenerator の leaps() を使いましょう。大本 → leaps() → 子供 → jumps() → 孫 、といった感じで生成できます。もちろん、 leaps() を繰り返して子孫を作ると同じ問題が発生します。 leaps() と jumps() を使い分けるのがポイントです。

なお、 SplittableGenerator の splits() では子孫で衝突する問題は事実上発生しません。一つの乱数列を共有している JumpableGenerator とは異なり、乱数列自体がインスタンスごとに異なるためです。

実装例

インタフェースが分かったところで、自分でアルゴリズムを実装する際の例を紹介します。
今回は Seiran128 アルゴリズムを実装してみました。
このアルゴリズムは LeapableGenerator の要件を満たしているので、それでざっくり実装してみます。

package Java17Test.src;  

import java.nio.ByteBuffer;  
import java.security.SecureRandom;  
import java.util.concurrent.atomic.AtomicLong;  
import java.util.random.RandomGenerator.LeapableGenerator;  

public class Seiran128 implements LeapableGenerator {  

    private long state0, state1;  

    private static final AtomicLong SEEDER = new AtomicLong(ByteBuffer.wrap(SecureRandom.getSeed(8)).getLong());  

    public Seiran128() {  
        this(SEEDER.getAndAdd(0x9e3779b97f4a7c15L));  
    }  

    public Seiran128(long seed) {  
        state0 = seed = seed * 6364136223846793005L + 1442695040888963407L;  
        state1 = seed = seed * 6364136223846793005L + 1442695040888963407L;  
    }  

    public Seiran128(long s0, long s1) {  
        if (s0 == 0 && s1 == 0) {  
            long seed = SEEDER.getAndAdd(0x9e3779b97f4a7c15L);  
            s0 = seed = seed * 6364136223846793005L + 1442695040888963407L;  
            s1 = seed = seed * 6364136223846793005L + 1442695040888963407L;  
        }  

        state0 = s0;  
        state1 = s1;  
    }  

    @Override  
    public void jump() {  
        internalJump(new long[] { 0xF4DF34E424CA5C56L, 0x2FE2DE5C2E12F601L });  
    }  

    @Override  
    public double jumpDistance() {  
        return 0x1p64;  
    }  

    @Override  
    public long nextLong() {  
        long s0 = state0, s1 = state1;  
        long result = Long.rotateLeft((s0 + s1) * 9, 29) + s0;  

        state0 = s0 ^ Long.rotateLeft(s1, 29);  
        state1 = s0 ^ s1 << 9;  

        return result;  
    }  

    @Override  
    public LeapableGenerator copy() {  
        return new Seiran128(state0, state1);  
    }  

    @Override  
    public void leap() {  
        internalJump(new long[] { 0x185F4DF8B7634607L, 0x95A98C7025F908B2L });  
    }  

    @Override  
    public double leapDistance() {  
        return 0x1p96;  
    }  

    private void internalJump(long[] jumpTable) {  
        long t0 = 0, t1 = 0;  

        for (int i = 0; i < 2; i++) {  
            for (int b = 0; b < 64; b++) {  
                if ((jumpTable[i] & (1L << b)) != 0) {  
                    t0 ^= state0;  
                    t1 ^= state1;  
                }  
                nextLong();  
            }  
        }  

        state0 = t0;  
        state1 = t1;  
    }  
}

nextDouble() などは RandomGenerator のデフォルト実装としてついてくるので基本的には nextLong() だけ、 LeapableGenerator を実装する場合は対応する jump() ・ leap() などを実装すればよいです。

ただ、 RandomGenerator.of("Seiran128") でインスタンス化できるか試してみましたが、どうもうまくいきませんでした。
内部実装を見た限りでは internal なアノテーションを設定する必要がありそうで、今のところ外部からの追加は難しそうです。
(もしうまくやる方法をご存じの方がいらっしゃればぜひ…)

LXM アルゴリズムについて

さて、いよいよ本題に入ります。
JDK 17 では以下の擬似乱数生成器 (アルゴリズム) が追加されています。 ソースコードはこちらを参照してください。

  • Xoroshiro128PlusPlus
  • Xoshiro256PlusPlus
  • L32X64MixRandom
  • L64X128MixRandom
  • L64X128StarStarRandom
  • L64X256MixRandom
  • L64X1024MixRandom
  • L128X128MixRandom
  • L128X256MixRandom
  • L128X1024MixRandom

RandomGenerator.of("Xoshiro256PlusPlus") といった形でアルゴリズム名を指定することでインスタンスを生成することができます。

xoroshiro128++ と xoshiro256++ については 作者のサイトでの説明 や論文 *2 が存在しますが、 LXM ファミリは初出当時は jdk リポジトリや API ドキュメント以外では記述のないアルゴリズムでした。 2023 年現在は論文 *3 が発表されており、詳細を知ることができます。

以降は、 java.util.random パッケージの "The LXM Family of Random Number Generator Algorithms" の項 をベースに紹介していきます。

LXM は、

  • L : LCG → Linear Congruential Generator; 線形合同法
  • X : XBG → Xor-Based Generator; xoshiro/xoroshiro アルゴリズム
  • M : Mixer → 上の 2 つを混ぜること (いわゆる出力関数)

という略語です。
名前の通り、内部に線形合同法と xo(ro)shiro の 2 つのアルゴリズムを独立して持ち、出力の際にそれらを "mix" することによって高品質な乱数出力を実現する設計となっています。

「あの悪名高い線形合同法？怪しい…」と思われた方もいらっしゃるかもしれませんが、線形合同法も使い方を工夫したりこのように他の擬似乱数生成器と組み合わせたりすることで高品質な乱数を生成することができます。使い方を工夫した例でいえば PCG が有名です。

L64X128MixRandom という暗号めいた名前は、それぞれどんなアルゴリズムを組み合わせているかを示しています。

name LCG (bits) XBG (algorithm) Mixer (algorithm)
L32X64MixRandom 32 xoroshiro64 (32-bit) mixLea32(s+x0)
L64X128MixRandom 64 xoroshiro128 mixLea64(s+x0)
L64X128StarStarRandom 64 xoroshiro128 starstar(s+x0)
L64X256MixRandom 64 xoshiro256 mixLea64(s+x0)
L64X1024MixRandom 64 xoroshiro1024 mixLea64(s+x0)
L128X128MixRandom 128 xoroshiro128 mixLea64(sh+x0)
L128X256MixRandom 128 xoshiro256 mixLea64(sh+x0)
L128X1024MixRandom 128 xoroshiro1024 mixLea64(sh+x0)

2021/09/12 当時は Package Summary の表 "LXM Multipliers" には概ね mixLea32 が使われていると書かれていますが、 これは誤りでした。 32 bit ベースの L32X64MixRandom だけが本当に mixLea32 で、それ以外の mixLea32 の行は実際は mixLea64 を使用しています。 2023/03/17 現在のドキュメントでは修正されています。

ちなみに、 LXM の xoroshiro128 のパラメータ (遷移関数のシフト/ローテート定数) は (24, 16, 37) であり、 Xoroshiro128PlusPlus の (49, 21, 28) とは異なります。
LXM のパラメータは公式実装における xoroshiro128** 、 Xoroshiro128PlusPlus のは xoroshiro128++ に基づいているようです。 (ちなみに xoroshiro128++ のほうが新型です。)
解析を行う場合に混同するとおかしくなるので注意が必要です。

L64X128MixRandom の必要最小限の実装は以下のようになります。なお、異常系や付加機能 (split() など) は省いています。

package Java17Test.src;  

import java.util.random.RandomGenerator;  

public class L64X128MixRandomClone implements RandomGenerator {  

    private final long lcgAdditive;  
    private long lcgState, xbgState0, xbgState1;  

    public L64X128MixRandomClone(  
        long lcgAdditive, long lcgState, long xbgState0, long xbgState1) {  
        this.lcgAdditive = lcgAdditive | 1;  
        this.lcgState = lcgState;  
        this.xbgState0 = xbgState0;  
        this.xbgState1 = xbgState1;  
    }  

    @Override  
    public long nextLong() {  

        // Mix  
        long r = lcgState + xbgState0;  
        r = (r ^ r >>> 32) * 0xdaba0b6eb09322e3L;  
        r = (r ^ r >>> 32) * 0xdaba0b6eb09322e3L;  
        r = (r ^ r >>> 32);  

        // LCG  
        lcgState = lcgState * 0xd1342543de82ef95L + lcgAdditive;  

        // XBG (xoroshiro)  
        long s0 = xbgState0, s1 = xbgState1;  
        s1 ^= s0;  
        s0 = Long.rotateLeft(s0, 24) ^ s1 ^ s1 << 16;  
        s1 = Long.rotateLeft(s1, 37);  
        xbgState0 = s0;  
        xbgState1 = s1;  

        return r;  
    }  
}

LCG パートでは、加算する定数 (lcgAdditive) をインスタンスごとに可変にすることで、 PCG と同様に異なる乱数列の生成 (split()) を可能にしています。
この部分での周期は 2 ^ {l} ( l は lcgState のビット数) です。

ちなみに、 L128 のバリアントでは 128 bit 乗算を行っている…　わけではなく、工夫によって 64 bit 乗算と 64 bit x 64 bit → 128 bit の乗算のみで計算できるようになっています。
まず、乗算する定数 0x1d605bbb58c8abbfd は 65 bit 、つまり上位 64 bit が 0x1 になっています。これによって、乗算結果の下位 64 ビット (lo) は lo * 0xd605bbb58c8abbfdL 、上位 64 bit (hi) は Math.unsignedMultiplyHigh(lo, 0xd605bbb58c8abbfdL) + hi * 0xd605bbb58c8abbfdL + lo として計算することができます。 64 bit ごとに区切った筆算をイメージすると分かりやすいです。

XBG パートは、既存の xo(ro)shiro 系アルゴリズムから出力関数 (+, **, ++) を除いた状態になっています。この部分での周期は 2 ^ {x} - 1 ( x は xbgState の合計ビット数) です。
-1 されているのは、 xbgState がすべて 0 の場合は何回遷移しても 0 のままになるので除外しているためです。実際のコンストラクタではすべて 0 で初期化されないように補正処理が行われています。

そして、全体としての周期は 2 ^ {l} \cdot (2 ^ {x} - 1) = 2 ^ {l + x} - 2 ^ {l} となります。具体的には、 L64X128MixRandom では 2 ^ {192} - 2 ^ {64} です。各パートの周期が互いに素なため、単に掛け合わせたぶんだけ周期になります。

Mix パートが出力関数です。これを状態更新前にやっているのは、内部状態の更新と並列に計算できるようにして高速化を狙っているためです。
この mixLea64 実装は SplittableRandom の出力関数に似ていますが、定数が異なっています。
シフトを 32 bit 単位に揃えたり乗算の定数を同じにすることでより高速化を図っているようです。
定数設定は ソースコードのコメント を見た限りでは Doug Lea 氏 (SplitMix の論文の共著者) によるものらしいです。詳しくは LXM の論文を参照して下さい。

ところで、現状の JDK 17 では LXM ファミリは SplittableGenerator だけ実装してありますが、 LeapableGenerator も実装することもできます。
例えば、 L64X128MixRandom における jump() と leap() はこのように実装できます。

    @Override  
    public void jump() {  
        internalJump(new long[] { 0x2bd7a6a6e99c2ddcL, 0x0992ccaf6a6fca05L });  
    }  

    @Override  
    public double jumpDistance() {  
        return 0x1p64;  
    }  

    @Override  
    public LeapableGenerator copy() {  
        return new L64X128MixRandomClone(lcgAdditive, lcgState, xbgState0, xbgState1);  
    }  

    @Override  
    public void leap() {  
        internalJump(new long[] { 0x360fd5f2cf8d5d99L, 0x9c6e6877736c46e3L });  
    }  

    @Override  
    public double leapDistance() {  
        return 0x1p96;  
    }  

    private void internalJump(long[] jumpTable) {  
        long x0 = 0, x1 = 0, s = lcgState;  

        for (int i = 0; i < 2; i++) {  
            for (long b = 1; b != 0; b <<= 1) {  
                if ((jumpTable[i] & b) != 0) {  
                    x0 ^= xbgState0;  
                    x1 ^= xbgState1;  
                }  
                nextLong();  
            }  
        }  

        lcgState = s;  
        xbgState0 = x0;  
        xbgState1 = x1;  
    }

基本的には xo(ro)shiro の ジャンプ関数実装 を使うだけです。
LCG パートに変更がないのは、ジャンプ距離が 2 ^ {64} や 2 ^ {96} と LCG の周期 ( 2 ^ {64}) の倍数になっており、 n 周して元の位置に戻ってきているためです。
理論上は任意長ジャンプの ArbitrarilyJumpableGenerator を実装することも可能ですが、そこそこ手間がかかるので割愛します。

これは完全に余談ですが、 L128X1024MixRandom のコンストラクタがあまりにもいかついので笑ってしまいました。

LXM ファミリの品質について

JEP 356 によれば、 LXM ファミリは TestU01 や PractRand といった著名な乱数性テストに合格するとのことです。
ただ、具体的にどう合格しているのか (TestU01 なら p-value の許容範囲はどこまでか・何個のシード値で試したか、 PractRand ならどのモードで何 TB まで確認したか) は書かれていませんでした。

そこで論文を見ると、TestU01 は BigCrush (最も多くのテストを行うモード) において、現行の L64X128MixRandom 相当の実装では複数の初期化方法について各 84 回ものテストを行い問題ないことを確認した、とあるようです。
また、 PractRand では 4 TB までテストにパスすることを確認した、ともありました。

ただ、確認したとされている PractRand のバージョンは 0.90 とのことで、現行の最新版 (0.94 及び pre-0.95) より前のものになっています。 v0.94 でいくつかのテストが追加されており、それで FAIL した擬似乱数生成器もそこそこあったので、懸念材料となります。

最も、単体でもほとんどの (線形性テストを除く) テストを突破しうる xo(ro)shiro に LCG を組み合わせたうえで mix までかけているので、大丈夫ではあるだろう、とは推測できます。
参考までに、 SplittableRandom は雑に言えば state += gamma; という足し算だけのシーケンスを mix を通して出力するだけのものです。それでも PractRand にパスしている ので、それより内部が複雑な LXM でも通るだろう、というのが理由です。

実際に確認してみましょう。 LXM ファミリの中で最も内部状態が小さい (つまり、問題が顕在化しやすいと思われる) L32X64MixRandom について、 PractRand v0.94 でテストを行いました。引数は -multithreaded -te 1 -tf 2 -tlmaxonly と、最も多くのテストをかけられるモードを指定しています。

RNG_test using PractRand version 0.94
RNG = L32X64MixRandomJava64, seed = 0x19aa8a23
test set = expanded, folding = extra

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 512 megabytes (2^29 bytes), time= 2.9 seconds
  no anomalies in 1220 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 1 gigabyte (2^30 bytes), time= 7.7 seconds
  no anomalies in 1293 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 2 gigabytes (2^31 bytes), time= 15.4 seconds
  no anomalies in 1368 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 4 gigabytes (2^32 bytes), time= 29.2 seconds
  no anomalies in 1448 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 8 gigabytes (2^33 bytes), time= 60.3 seconds
  no anomalies in 1542 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 16 gigabytes (2^34 bytes), time= 117 seconds
  no anomalies in 1636 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 32 gigabytes (2^35 bytes), time= 213 seconds
  no anomalies in 1716 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 64 gigabytes (2^36 bytes), time= 455 seconds
  no anomalies in 1807 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 128 gigabytes (2^37 bytes), time= 914 seconds
  no anomalies in 1902 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 256 gigabytes (2^38 bytes), time= 1667 seconds
  no anomalies in 1983 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 512 gigabytes (2^39 bytes), time= 3541 seconds
  no anomalies in 2058 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 1 terabyte (2^40 bytes), time= 7137 seconds
  no anomalies in 2132 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 2 terabytes (2^41 bytes), time= 14111 seconds
  no anomalies in 2194 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 4 terabytes (2^42 bytes), time= 29294 seconds
  no anomalies in 2255 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 8 terabytes (2^43 bytes), time= 59262 seconds
  no anomalies in 2313 test result(s)

rng=L32X64MixRandomJava64, seed=0x19aa8a23
length= 16 terabytes (2^44 bytes), time= 119867 seconds
  no anomalies in 2366 test result(s)

時間の都合上 16 TB までの範囲ではありますが、異常は検出されませんでした。
(本来は 32 TB まで検定します。)
私の実験は完全ではありませんが、この範囲では少なくとも大丈夫だと言えそうです。

他の L64X128MixRandom などについては L32X64MixRandom よりも内部状態が大きくなっているので、それ以上の品質を確保できていると推測できます。

要約すれば、現時点では問題は検出されない、高品質な擬似乱数生成器である、ということです。

速度について

LCG と xo(ro)shiro を組み合わせて mix している関係上、それぞれ単体の xo(ro)shiro や SplittableRandom より遅くなっているのではないか、と推測できます。

実際に測ってみましょう。 1 GB (nextLong() 2 ^ {27} 回分) の乱数を生成するのにかかった時間 (実測値, ms) と、そこから計算した速度 (GB/s) を速い順に示します。環境は以下の通りです。

12th Gen Intel(R) Core(TM) i7-12700F

openjdk 17.0.6 2023-01-17
OpenJDK Runtime Environment Temurin-17.0.6+10 (build 17.0.6+10)
OpenJDK 64-Bit Server VM Temurin-17.0.6+10 (build 17.0.6+10, mixed mode, sharing)
Algorithm ms GB/s
SplittableRandom 226.00 4.42
Xoroshiro128PlusPlus 240.75 4.15
Xoshiro256PlusPlus 260.68 3.84
L64X128StarStarRandom 291.99 3.42
L64X128MixRandom 313.23 3.19
L64X256MixRandom 340.49 2.94
L64X1024MixRandom 362.83 2.76
L32X64MixRandom 386.71 2.59
L128X256MixRandom 575.84 1.74
L128X128MixRandom 578.84 1.73
L128X1024MixRandom 922.06 1.08
Random 2189.57 0.46
SecureRandom 75161.27 0.01

実測値のため環境や試行によってぶれる可能性がありますが、この結果を見ると、

  • ざっくりした傾向としては SplittableRandom > Xo(ro)shiro > LXM (L64 > L32 > L128) >> Random >>(越えられない壁)>> SecureRandom といった感じです。
    • 品質を犠牲にしてでも速度第一にする場合は SplittableRandom 、品質も大事な場合は Xoroshiro128PlusPlus がよさそうです。
  • LXM シリーズは、 L が小さい > X が小さい順に速い傾向がみられます。
    • 長い乗算が必要になるとさすがに遅くなるようです。
  • L32X64MixRandom は 32 bit ベースのため nextLong() 1 回あたり 2 回の遷移が必要になって結果として遅くなりそう、と予想していましたが、思っていたよりも速いようです。
  • 新しいアルゴリズムは、遅いものでも Random の 2.3 倍、速いものでは 9.0 倍ほど速くなっています。
    • Random が遅いのは、マルチスレッド対応のため都度ロックをかけているせいと推測しています。
  • (当然ながら) SecureRandom は非常に低速なため、必要な用途にだけ使用すべきでしょう。

といった感じです。

余談ですが、 C# に移植したバージョンで BenchmarkDotNet を使って nextLong() を測定してみると、以下のようになりました。

BenchmarkDotNet=v0.13.2, OS=Windows 11 (10.0.22000.1574/21H2)
12th Gen Intel Core i7-12700F, 1 CPU, 20 logical and 12 physical cores
.NET SDK=7.0.200
  [Host]     : .NET 7.0.3 (7.0.323.6910), X64 RyuJIT AVX2
  DefaultJob : .NET 7.0.3 (7.0.323.6910), X64 RyuJIT AVX2
Method Mean Error StdDev Ratio RatioSD
Xoroshiro128PlusPlus 0.5713 ns 0.0114 ns 0.0107 ns 0.22 0.01
Shioi 0.5972 ns 0.0127 ns 0.0113 ns 0.23 0.00
L64X128StarStar 0.6012 ns 0.0145 ns 0.0136 ns 0.23 0.01
SplitMix 0.6013 ns 0.0168 ns 0.0157 ns 0.23 0.01
Seiran 0.6019 ns 0.0206 ns 0.0183 ns 0.23 0.01
Xoshiro256PlusPlus 0.7908 ns 0.0123 ns 0.0103 ns 0.30 0.00
L64X256Mix 0.8066 ns 0.0126 ns 0.0118 ns 0.31 0.00
L64X128Mix 0.8069 ns 0.0160 ns 0.0150 ns 0.31 0.01
L64X1024Mix 1.3266 ns 0.0239 ns 0.0212 ns 0.51 0.01
L128X128Mix 1.9699 ns 0.0253 ns 0.0225 ns 0.76 0.01
L128X256Mix 2.0528 ns 0.0215 ns 0.0202 ns 0.79 0.01
SystemRandomUnseeded 2.5940 ns 0.0259 ns 0.0242 ns 1.00 0.00
L32X64Mix 2.7956 ns 0.0280 ns 0.0262 ns 1.08 0.02
L128X1024Mix 4.5566 ns 0.0337 ns 0.0299 ns 1.76 0.03
SystemRandomSeeded 50.8033 ns 0.2700 ns 0.2525 ns 19.59 0.19

補足すると、 SplitMix は SplittableRandom 相当、 SystemRandom(Un)Seeded は .NET の System.Random にシードを与えた/与えなかったもの (それぞれ内部アルゴリズムが変化します) 、 Shioi と Seiran はいつもの私の自作アルゴリズムです。

微妙に結果は異なりますが、全体の印象としては同じ雰囲気です。

内部状態復元攻撃

これだけで記事 1 本文の分量になりそうだったので、後日別の記事に…できれば…いいな…。
L32 ・ L64 シリーズについては攻撃の目処が立っています。

その他の改善点

全体の話に戻します。

乱数を加工して使いやすくするメソッドもいくつか追加・改善されています。

nextDouble(double bound) / nextDouble(double origin, double bound)

nextDouble(double bound) と nextDouble(double origin, double bound) が追加され、一定範囲への変換が楽になりました。 nextFloat() も同様です。
今まではストリームを返すほう doubles() でだけ使えていたのが nextDouble() でも使えるようになった形になります。
「そんなの nextDouble() * (bound - origin) + origin すればいいだけなのでは？」と思われた方もいると思いますが、この方法には落とし穴があります。このメソッドの戻り値は origin <= value < bound を意図していますが、前述の式では浮動小数点数の丸め誤差によってごくまれに value == bound となってしまう場合があります。今回追加されたメソッドでは、それをよしなに制御して範囲内になるように保証してくれます。

nextExponential()

指数分布 に従う乱数を生成する nextExponential() が追加されています。

nextGaussian(double mean, double stddev)

また、 正規分布 に従う乱数を生成する nextGaussian() に nextGaussian(double mean, double stddev) が追加され、平均と標準偏差を自分で計算する必要がなくなりました。
さらに、内部実装が変わってパフォーマンスも改善しています。以前は Polar Method で実装されていましたが、 改良 Ziggurat アルゴリズム に変更されてより高速化しています。
2 ^ {24} 回の呼び出しを実測してみたところ、 OpenJDK 64-bit 16.0.2 では 918 ms, 17 では 932 ms となっており、手元の環境ではなぜかあまり変化がなかったようです。あれ……？
(ちなみに、 nextExponential() も同じアルゴリズムで実装されています。)

RandomGeneratorFactory

RandomGeneratorFactory というファクトリも追加されています。これを使うと、アルゴリズム名からインスタンスを生成したり、性質に応じて (例えば内部状態が 256 bit 以下のものを) 抽出したりすることができます。
また、 RandomGeneratorFactory.getDefault().create() でアルゴリズムを明示せずにインスタンスを生成することもできます。現行実装では L32X64MixRandom が選ばれるようです。

なぜ L32X64MixRandom がデフォルトになったのかは正直よくわかりません……
L32X64MixRandom は 32 bit 環境に特化したものであり周期が短めで、一度に 64 bit 必要な乱数生成 (nextDouble() や nextLong() など) や大規模な並列実行には向いていません。
L64X128MixRandom が無難でよかったのでは、と思ってしまいます。

余談・リリース以降の更新点について

L32X64MixRandom#nextLong() の不思議な実装については、前回の記事をご覧ください。

コミットログを見ると、一時期 L32X64MixRandom と L128X256MixRandom の初期化周りで問題があったらしいですが、現在は解消されています。

また、 nextGaussian() と nextExponential() の分布が正しくなかった 問題もあったようですが、こちらも解消されています。

アルゴリズムの選び方

LXM アルゴリズムが大量に追加されていてどれをどう使えばいいか分からなくなりそうだったので、 API ドキュメントの内容をざっくり私見も含めつつまとめると、優先度順に

  1. セキュリティ関連・シード値など 予測不可能性が必須 → SecureRandom
  2. 32 bit 環境 → L32X64MixRandom
  3. 高速性 が重要 → Xoroshiro128PlusPlus
  4. 小～中規模の並列処理 がしたい → L64X256MixRandom
  5. 超大規模の並列処理 がしたい → L128X256MixRandom
  6. 多次元均等分布性 が重要 → L64X1024MixRandom
  7. とくに要件がない、 とりあえず動けばいい → L64X128MixRandom
    (または RandomGenerator.getDefault())

といった感じです。

「多次元均等分布性が重要」というのは、身近な例ではトランプなどのシャッフルです。
n 要素のシャッフルは n! 通りの結果を持ちますが、小さな擬似乱数生成器だとこれだけのパターンを生成することができない場合があり、そうなると「絶対に出現しない盤面」が出てしまいます。
例えば、トランプ (ジョーカー 2 枚) なら 54! \approx 2 ^ {237} 通り、麻雀牌なら 136! \approx 2 ^ {773} 通り *4 になるので、最低でも (2冪の指数) ビット分の均等分布性を持つ擬似乱数生成器を選ぶ必要があります。
今回のアルゴリズムでいえば、トランプなら最低でも L64X256MixRandom を、麻雀牌なら L64X1024MixRandom を、ということになります。

ただこれは理想的な処理を行った場合の話で、例えば安直に nextLong(54) などとして乱数を生成してしまうとそれだけで 64 bit 分の情報量を消費してしまうので、トランプの例では 64 \times 54 = 3456 ビット必要になってしまいます。それだと L64X1024MixRandom でも足りないのでビット数を節約するよう工夫しなければなりません。 (byte 単位に分割して自分で nextLong(max) 相当のコードを書く、など。)
さらに言えば、現時点で LXM ファミリを引数なしの RandomGenerator.of(algorithm).create() から生成した場合、内部的には long のシード値を用いて初期化されるため、初期状態は 2 ^ {64} パターンに限定されてしまいます。もし全パターンを均等に得たい場合は create(byte[] seed) に十分な情報量を持つソース (SecureRandom の出力など) を渡して初期化する必要があります。

均等分布性について

ちょうど均等分布性の話が出たので、 Package Summary の "The LXM Family of Random Number Generator Algorithms" の均等分布性に関連する部分の解説をしたいと思います。

exactly equidistributed とは

"exactly equidistributed" (勝手訳: 厳密に均等分布する) とは、 nextLong() で得られる値について、すべての値が厳密に等確率で出現させられる、ということです。

L32X64MixRandom 以外の LXM ファミリについては、この性質を満たします。
L32X64MixRandom は、 32 bit 精度 (nextInt()) では厳密に均等分布しますが、 64 bit 精度 (nextLong()) では厳密に均等分布しません。
詳しく説明すると、 L32X64MixRandom#nextLong() は ((long)nextInt() << 32) ^ (long)nextInt() として実装されているため、 nextLong() で厳密に均等分布させるためには nextInt() 2 回分 (= 2 次元に) 厳密に均等分布させる必要がありますが、 L32X64MixRandom#nextInt() は 1 次元までしか厳密に均等分布しません。そのため、 nextLong() では厳密に均等分布しないことになります。

また、 Xoroshiro128PlusPlus と Xoshiro256PlusPlus については記載がありませんが、厳密に均等分布しないと思われます。内部状態がすべて 0 にはならないので、その時に出力されるはずだった 0 の出現個数が他の値より 1 少なくなるためです。
具体例としては、 Xoroshiro128PlusPlus では 0 が出る確率が \frac{2 ^ {64} - 1}{2 ^ {128} - 1} 、それ以外の値はそれぞれ \frac{2 ^ {64}}{2 ^ {128} - 1} となります。

ちなみに、 LCG (線形合同法) は基本的に厳密に 1 次元均等分布します。
具体例で説明します。 64 bit の線形合同法 (L64) は周期が 2 ^ {64} であるので、その周期中に long がとりうる全ての値を各 1 回だけ、すなわち厳密に等確率 2 ^ {-64} で出力します。したがって厳密に均等分布するといえます。

k-equidistributed とは

"k-equidistributed" (勝手訳: k 次元均等分布する) とは、要素数 k の long[] を連続する nextLong() で埋めたとき、すべての値の組み合わせが「ほぼ」等確率で出現させられる、ということです。
こちらは "exactly equidistributed" よりは緩めで、「厳密に」等確率でなくてもよいようです。上に挙げた Xoroshiro128PlusPlus のように 2 ^ {-128} オーダーでの違いは許容しよう、といった雰囲気です。

この性質は、従来のアルゴリズムでは Mersenne Twister の「623 次元均等分布」が有名でしょう。 (ただしこれは 32 bit 精度 (int[]) での値です。)

具体例として、 2 次元均等分布する L64X128MixRandom で考えてみましょう。
まず、内部の xoroshiro128 は {0, 0} 以外の任意の 2 次元ベクトルを生成することができます。 (つまり、 xoroshiro128 は単体で 2 次元均等分布します。)
ここに LCG パートを合わせるので、 XBG パートで {0, 0} だった場合の (つまり発生しえない) {mix(s), mix(s * m + a)} (s, m, a は LCG の内部状態と乗算定数・加算定数を示します) のベクトルが生成されるパターンが 1 つ少なくなって全周期中にそれぞれ 2 ^ {64} - 1 個、それ以外のベクトルはそれぞれ 2 ^ {64} 個ずつ生成されることになります。分母 (全周期) は 2 ^ {192} - 2 ^ {64} なので、 2 ^ {-128} オーダーでの違いとなります。

なお、 RandomGeneratorFactory.of(algorithm).equidistribution() で得られる値は、こちらの k 次元均等分布のほうです。

以下に 64 bit 精度での均等分布次元をまとめた表を示します。

アルゴリズム名 均等分布次元 厳密に均等分布するか
L32X64MixRandom 1 (2 @ 32bit)
L64X128MixRandom 2 ☑
L64X128StarStarRandom 2 ☑
L64X256MixRandom 4 ☑
L64X1024MixRandom 16 ☑
L128X128MixRandom 2 ? ☑
L128X256MixRandom 4 ? ☑
L128X1024MixRandom 16 ? ☑
Xoroshiro128PlusPlus 1
Xoshiro256PlusPlus 3

? をつけた箇所は API ドキュメント及び equidistribution() ではなぜか 1 となっていますが、実際はそれぞれ 2, 4, 16 が正しいのではないかと思います。
L128X128MixRandom ・ L128256MixRandom ・ L1281024MixRandom は、内部の XBG パート (xoroshiro128 ・ xoshiro256 ・ xoroshiro1024) がそれぞれ 64 bit 精度で 2 / 4 / 16 次元均等分布します。
これと LCG パートを組み合わせて最終的な出力が得られるのですが、この時 LCG パートの出力がどのようなものであれ XBG パートが均等分布する = 任意の値 (k 次元ベクトル) をとることができるのですから、組み合わせた出力も (XBG パートでよしなにすることで) 任意の値をとることができる = 均等分布する、と考えられます。

おわりに

.NET といい Java といい、最近標準ライブラリの擬似乱数生成器の改善がアツいです。(2021年当時の感想)

全体の所感としては、高速性をとった .NET と拡張性をとった Java という印象を受けました。
また両者とも xo(ro)shiro 系アルゴリズムを採用しており、浸透してきたようで勝手に嬉しく思っています。

ただ、他に採用実績がない初出のアルゴリズムを標準ライブラリにぶち込んでいくのはなかなか挑戦的だな、という気持ちがあります。
もちろん内部構造はどれもよく研究されたものですので無謀というわけではないですが、公開されてから欠点が見つかることはよくあるので…

*1:Anonymous Author(s). 2017. Better Splittable Pseudorandom Number Generators (and Almost As Fast). PACM Progr. Lang. 1, 1, Article 1 (January 2017), 41 pages.

*2:David Blackman and Sebastiano Vigna. Scrambled linear pseudorandom number generators, 2019. To appear in ACM Trans. Math. Softw..

*3:Guy L. Steele Jr. and Sebastiano Vigna. 2021. LXM: better splittable pseudorandom number generators (and almost as fast). Proc. ACM Program. Lang. 5, OOPSLA, Article 148 (October 2021), 31 pages. https://doi.org/10.1145/3485525

*4:同種の牌も区別しているものとします。区別しなければ約 2 ^ {617} 通りになります。