Swift
乱数
swtws

Swift4.2のRandom APIについて

この記事はSwift Tweets 2018 Summerで発表したLTのもうちょっと踏み込んだになります。

まだ変更の多い部分ですので、最新情報はswiftリポジトリをご参照ください。

Swift4.2にて導入されるRandom APIについて

なぜRandom APIが必要だったのか

iOS/macOS開発の場合、これまではarc4random_uniformを使っていたと思います。
ところで、クロスプラットフォーム開発にて乱択アルゴリズムを使いたい〜等となると、Linux環境にはarc4random_uniformが存在しないため自分で用意する必要があります。

Linux環境向けの実装として、下請けにはrandomrand 等が使えますが、これらはUpper Boundをとらない関数になっています。
arc4random_uniformに合わせるために引数upperBoundで剰余を取るというのが思いつきますが、その方法だとModulo Biasがかかってしまいます。ほとんどの場合Modulo Biasは微小なので気付かないし影響もないですが、最悪のケースでは要素間の選択確率が2倍開きます。ので一般的な実装においては使うべきではないです。

強引な最悪のケース

クロスプラットフォームでの確認も面倒だし、標準で用意しておいてほしい!

といった内容がプロポーザルのMotivationに書いてありました。
https://github.com/apple/swift-evolution/blob/master/proposals/0202-random-unification.md

RandomNumberGenerator

乱数生成器のインターフェースとして用意されているのがRandomNumberGeneratorプロトコルです。

public protocol RandomNumberGenerator {
  // This determines the functionality for producing a random number.
  // Required to implement by all RNGs.
  mutating func next() -> UInt64
}

この他にnext<T : FixedWidthInteger & UnsignedInteger>() -> Tnext<T : FixedWidthInteger & UnsignedInteger>(upperBound: T) -> Tのデフォルト実装が用意されており、UInt64を生成する部分を実装するだけでそれらも使えます。特に後者のupperBound付きの関数はModulo Biasを考慮した実装になっており、安全に使えます。

Random

そして公式に用意されているRNGの実装がRandomです(が、名前は変更予定のようです)。

Randomはプラットフォームごとに適切な下請けの関数を使うようになっており、これを使って書けばクロスプラットフォームでの実装が簡単になります。
Randomのコメントによると初期化いらずでスレッドセーフというところまでは共通ですが、cryptographic qualityが下請けによって変わると書いてあるのでそこだけ注意が必要です(後述)。

Randomはstructになっていますが、これ自体は内部状態を持たないため、Random.defaultは毎回新しいインスタンスを返すようになっています。defaultが用意されているのは次に書く型ごとの実装のためのようです。

各型の乱数生成

staticメソッドとしてFIxedWidthInteger.random(in: Range<Self>, using: inout RandomNumberGenerator)BinaryFloatingPoint.random(in: Range<Self>, using: inout RandomNumberGenerator)等が追加されました。
数値以外ではCollectionについてもCollection.randomElement(using: inout RandomNumberGenerator)が追加されています。

usingを省略した場合Random.defaultが使われるようになっており、基本的に背後のRandomNumberGeneratorを意識せずに使用できるAPIになっています。
このデザインだとRandom()で作ってusingに渡す処理を一行でかけないので、そのためにRandom.defaultを用意しているとのことです(参考)。

まとめ

乱数周りはarc4random_uniformをメインに自分で色々と用意していくスタイルだったので、Random APIの実装によって大いに時短できそうです。

Swift4.2からはRandom APIをガンガン使ってクロスプラットフォームで再利用しやすいプログラムを書いていきましょう。

補:RandomはCryptographically secureなのか?

Randomのコメントより、

/// While the Random.default generator is automatically seeded and
/// thread-safe on every platform, the cryptographic quality of the stream of
/// random data produced by the generator may vary. For more detail, see the
/// documentation for the APIs used by each platform.
///
/// - Apple platforms use arc4random_buf(3).
/// - Linux platforms use getrandom(2) when available; otherwise, they read
/// from /dev/urandom.

またSE-0202 Random Unification#Random Number Generatorより、

The aspiration is that this RNG should be cryptographically secure, provide reasonable performance, and should be thread safe. If a vendor is unable to provide these goals, they should document it clearly.

後者は項名がRandom Number Generatorになっていてややこしいですが、ここで対象としているのはRandomNumberGeneratorプロトコルでなくRandom(default RNG)とその下請けの実装のようです。

現在のRandomの3つの下請けはどれもcryptographically secureなようですが、cryptographically secureな実装を用意できないプラットフォームのことを考慮してRandom自体のsecurityは保証しないことになっているようです。
安全性が要求される箇所にはcryptographically secureを保証できる手段を用いるのが良さそうです(といってunsecureな実装が足されて、かつそのプラットフォームを使うはめになることはそうそう無さそうですが)。

ところでthread safetyについてはコメントでは断言されててプロポーザルでは曖昧になってるんですがどっちが正しいんでしょう?

おまけ:浮動小数点数の乱数生成について

浮動小数点数の生成についても調べていたので、ついでにまとめます。

浮動小数点数について

浮動小数点数は$(符号)*2^{(指数部)} * 1.(仮数部)$という形で表現されます。$1.(仮数部)$は2進数で1.0..<10.0 (=2)という範囲を取ります。

浮動小数点数1..<2の範囲には$2^{(仮数部ビット数)}$だけの点があり、この区間で指数部は$0$で一定なので仮数部のインクリメントでの増加量が点間の間隔となり、Doubleの場合$2^{-52}$となります。これはSwiftではBinaryFloatingPoint.ulpOfOneとして定義されています(1.0.nextUp - 1 == .ulpOfOne)。

さて、1/2..<1の範囲でも同じことを考えてみます。指数部は1..<2範囲の指数部を-1したものになり、一方仮数部のビット数は変わらない=点数は変わらないため、点の間隔は$2^{-53}$、Swiftでは.ulpOfOne/2になります。
同様に1/4..<1/2, 1/8..<1/4範囲で間隔はどんどん半分になっていきます。逆に2..<4, 4..<8の範囲では倍々になっていくので、浮動小数点数はその絶対値が小さいほど密に、大きいほど疎に存在することになります。

[0, 1)区間からの一様サンプリング

浮動小数点数の生成についてはフォーラムでも参照されていた以下の内容が参考になります。

http://xoshiro.di.unimi.it/ の Generating uniform doubles in the unit intervalの項

曰く、Doubleを生成する場合は一様乱数r: UInt64を用いてDouble(r >> 11) * (.ulpOfOne / 2)すればいいということになっています。

(.ulpOfOne / 2)1/2..<1区間の点の間隔で、さらにこれより下の区間においては点間隔が半分ずつになっていくだけなので、0..<1全範囲でこの刻みを使うことができます。
Double(r >> 11)は53ビット乱数、すなわち0..<2**53をとりますが、前述の通り(.ulpOfOne / 2)は$2^{-53}$なので掛け合わせることで0..<1の乱数を生成できます。

これより小さい刻みを使う場合は、例えば.ulpOfOne/4なら1/2..<1では丸めが生じてしまい、最大値として1.0を生成し得るようになってしまいます。

// Floatの場合
print(Float(UInt32.max >> 8) * (.ulpOfOne / 2) < 1)   // 2^-23刻み true
print(Float(UInt32.max >> 7) * (.ulpOfOne / 4) < 1)   // 2^-24刻み false

// Doubleの場合
print(Double(UInt64.max >> 11) * (.ulpOfOne / 2) < 1) // 2^-53刻み true
print(Double(UInt64.max >> 10) * (.ulpOfOne / 4) < 1) // 2^-54刻み false

逆に大きい刻みを使う場合は、例えば.ulpOfOneなら均一にサンプリングできますが、間隔が大きくなる分生成し得る値のバリエーションは半分になってしまいます。
というわけで0..<1区間からサンプリングする場合、.ulpOfOne/2刻みを使うのが最も良いということになります。

SwiftBinaryFloatingPoint.random(in: range)は以上説明した方法で[0, 1)を生成しているので、今後は気をつけなくても自然に良い実装が使えます。

https://github.com/apple/swift/blob/a6952decab6f918a9df3c6fa342153a9f9204f8e/stdlib/public/core/FloatingPoint.swift.gyb#L2391-L2458

[low, high)区間からの一様サンプリング

$[0, 1)$区間が取れるのでそれをrとして(high - low) * r + lowとすれば$[low, high)$になりそうです。
ところで1..<2から取る場合を考えると、0..<1の最大値が1.nextDownで、1..<2の最大値は2.nextDownです。
1..<2区間なので0..<1を1倍して1を足すのですが、浮動小数点数の仕様的に1.nextDonw * 2 == 2.nextDonwであり、1.nextDonw + 12.nextDown2の間に来てしまいます。
実際に計算してみると結果は丸められて2となり、範囲から外れてしまいます。

let r = Float(UInt32.max >> 8) * (.ulpOfOne / 2)
print(r == 1)    // false
print(r+1 == 2)  // true

このような計算誤差で範囲から外れるケースの対処を調べたのですが、計算の結果がupperBound以上になったら棄却するという方法が見られました。
Swiftにもこのバグが入っていたのですが同じ方法で修正されました(https://github.com/apple/swift/pull/17794)。

標準正規分布からのサンプリング

一様分布からのサンプリングができれば以下のような方法で標準正規分布からのサンプリングができます。

この3つについてSwiftで速度比較用に書いたのが以下です(Zigguratの実装に自信がないですが……)。
https://github.com/t-ae/ZigguratSampler/blob/5c1d6d391c34c471f8f23b9e729f689f84ff0cad/Tests/ZigguratSamplerTests/ZigguratSamplerTests.swift#L35-L86

Zigguratが高速との話でしたが特に最適化を入れていない状態で一番遅いです。ほか2つはこのへんに書いてあるとおりBox-Mullerの方が高速でした。

前述した.ulpOfOne/2刻みの実装を一様分布の生成に用いると、Box-Mullerよりpolar methodで生成するほうが広い範囲の値を生成できるようです。

import Foundation

func boxmuller(u1: Double, u2: Double) -> Double {
    return sqrt(-2 * log(u1)) * sin(2 * .pi * u2)
}
func polar(u1: Double, u2: Double) -> Double {
    let r = u1*u1 + u2*u2
    assert(0 < r && r < 1)
    return sqrt(-2 * log(r) / r) * u1
}

// 最大値
print(boxmuller(u1: .ulpOfOne/2, u2: 0.25))   // 8.571674348652905
print(polar(u1: (0.5.nextUp)*2-1, u2: 0))     // 12.00727336061225

// 0に近い値
print(boxmuller(u1: 1.0.nextDown, u2: Double.ulpOfOne/2))   // 1.0394658142353987e-23
print(polar(u1: (0.5.nextUp)*2-1, u2: 1.0.nextDown*2-1))  // 6.617444900424223e-24

といってこちらも正規分布が生成しうるもっと大きな値を生成できないし、どの手法を使うのが最適なのか分からないので、このあたりも公式に用意してほしいところですね。