I (think I've) implemented the method described in the paper I mentioned in the OP:
func pseudorandomFloatInClosedUnitRange<R>(using prng: inout R) -> Float
where R: PseudorandomNumberGenerator
{
// Start by choosing an exponent (with correct probability):
var rnd64 = prng.next()
var i: UInt8 = 0
while (rnd64 >> i) & 1 != 0 && i < 126 {
i += 1
if i & 63 == 0 { rnd64 = prng.next() }
}
var exp = UInt8(126) &- i
// Note: Can be at most 126, not 127 (as in exp of 1). This is because
// we have to increment it with 50% probability if the random
// significand happens to be zero (see below).
// Choose a random 23-bit mantissa (and waste some rnd bits for now):
rnd64 = prng.next()
let significand = UInt32(truncatingIfNeeded: rnd64) & 0x7FFFFF
// If it's zero, then with 50% probability we should move to
// the next exponent range:
if significand == 0 && (rnd64 >> 23) & 1 != 0 { exp &+= 1 }
return Float(bitPattern: (UInt32(exp) << 23) | significand)
}
But this is only for the closed range 0 ... 1.
Assuming this works as intended, I wonder how it can be extended to any range
(including -Float.greatestFiniteMagnitude ... Float.greatestFiniteMagnitude).
Can it somehow be used (I guess it can't be just scaled and translated) for other ranges or do we have to use a completely different algorithm?
Will there be a problem with -0 and 0 when using ranges that span from - to + etc ...
And I have no idea how "good" this method is compared to others (as I haven't really found that many, let alone a comparison between them).
Here's the function again, in a complete program that can be compiled (swiftc -O test.swift) and run.
import func QuartzCore.CABase.CACurrentMediaTime
//-----------------------------------------------------------------------------
// Some stuff that it uses:
//-----------------------------------------------------------------------------
protocol PseudorandomNumberGenerator {
associatedtype State
var state: State { get }
/// Creates a new instance of this generator with a state that is
/// completely determined by the given seed. Any `UInt64` is allowed and
/// maps to a valid state. Each seed may or may not map to a unique state
/// depending on the generator.
init(seed: UInt64)
/// Creates a new instance of this generator with the given state.
/// The initializer fails if the given state is invalid for this generator.
init?(state: State)
/// Advances the state of the generator and returns the next random bit
/// pattern.
mutating func next() -> UInt64
}
extension PseudorandomNumberGenerator {
/// Creates a new instance of this generator by calling the required seeded
/// initializer `init(seed: UInt64)` with a truly random seed.
init() {
self.init(seed: UInt64.random(in: .min ... .max))
}
/// Creates a new instance of this generator by calling the required seeded
/// initializer with the given `seed` or a truly random seed if it's nil.
init(seed: UInt64?) {
if let seed = seed {
self.init(seed: seed)
} else {
self.init()
}
}
}
/// The wyrand generator, translated from:
/// https://github.com/wangyi-fudan/wyhash
struct WyRand : PseudorandomNumberGenerator {
typealias State = UInt64
private (set) var state: State
init(seed: UInt64) {
//print("WyRand(seed: \(seed))")
self.state = seed
}
init(state: State) {
// Every UInt64 value is a valid WeRand state, so no need for this
// initializer to be failable, it will satisfy the requirement anyway.
self.state = state
}
mutating func next() -> UInt64 {
state &+= 0xa0761d6478bd642f
let mul = state.multipliedFullWidth(by: state ^ 0xe7037ed1a0b428db)
return mul.high ^ mul.low
}
}
//-----------------------------------------------------------------------------
// The function
//-----------------------------------------------------------------------------
func pseudorandomFloatInClosedUnitRange<R>(using prng: inout R) -> Float
where R: PseudorandomNumberGenerator
{
// Choose an exponent (with correct probability):
var rnd64 = prng.next()
var i: UInt8 = 0
while (rnd64 >> i) & 1 != 0 && i < 126 {
i += 1
if i & 63 == 0 { rnd64 = prng.next() }
}
var exp = UInt8(126) &- i
// Note can be at most 126, not 127 (as in exp of 1). This is because
// we have to increment it with 50% probability if the random
// significand happens to be zero (see below).
// Choose a random 23-bit mantissa (and waste some rnd bits for now):
rnd64 = prng.next()
let significand = UInt32(truncatingIfNeeded: rnd64) & 0x7FFFFF
// If it's zero, then with 50% probability we should move to
// the next exponent range:
if significand == 0 && (rnd64 >> 23) & 1 != 0 { exp &+= 1 }
return Float(bitPattern: (UInt32(exp) << 23) | significand)
}
//-----------------------------------------------------------------------------
// Demo
//-----------------------------------------------------------------------------
func test() {
var prng = WyRand()
let numRuns = 10
let floatsPerRun = 10_000_000
var meanSum = Double(0)
for _ in 0 ..< numRuns {
let t0 = CACurrentMediaTime()
var sum = Double(0)
for _ in 0 ..< floatsPerRun {
let v = pseudorandomFloatInClosedUnitRange(using: &prng)
sum += Double(v)
}
let t1 = CACurrentMediaTime()
print(t1 - t0, "s | mean: \(sum / Double(floatsPerRun))")
meanSum += sum
}
print("mean mean", meanSum / Double(floatsPerRun * numRuns))
}
test()
I just copied that from the paper without thinking ...