Ah, of course! 
I just copied that from the paper without thinking ...
Here's the improved function:
func pseudorandomFloatInClosedUnitRange<R>(using prng: inout R) -> Float
where R: PseudorandomNumberGenerator
{
// Start by choosing an exponent (with correct probability):
var exp = 126 &- prng.next().trailingZeroBitCount
if exp == 62 { exp &-= prng.next().trailingZeroBitCount }
if exp < 0 { exp = 0 }
// Note that `exp` 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):
let 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)
}
The times per run reported by the test program (new version below) changed from about 0.11 seconds to 0.04 s (on my old MBP, one run consists of generating and adding 10 million Float values).
I added two simpler functions for comparison (that do not do what we want)
- One that uses mimics what's currently in
Float.random(in:using:)but simplified to work only for0 ..< 1(while wasting 40 random bits per generated float) - One that just returns
0or1(while wasting 63 random bits).
If we use the zero-or-one function as a baseline, the improvement you mentioned made it (0.11 - 0.011) / (0.04 - 0.011) = 3.4 times faster.
Here's the complete updated program.
import func QuartzCore.CABase.CACurrentMediaTime
//-----------------------------------------------------------------------------
// MARK: - Some scaffolding and the PRNG we'll use
//-----------------------------------------------------------------------------
/// A type that can produce a reproducible stream of uniformly distributed
/// pseudorandom bits.
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 64 uniformly
/// distributed pseudorandom bits.
mutating func next() -> UInt64
}
extension PseudorandomNumberGenerator {
/// Creates a new instance of this generator by calling the required seeded
/// initializer with the given `seed` or a truly random seed if given `nil`.
init(seed: UInt64? = nil) {
let seed = seed ?? UInt64.random(in: .min ... .max)
self.init(seed: seed)
}
}
/// 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) {
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
}
}
//-----------------------------------------------------------------------------
// MARK: - The function we're trying to write
//-----------------------------------------------------------------------------
/// This function chooses a pseudorandom value from a continuous uniform
/// distribution within the closed range `0 ... 1` and then convert that value
/// to the nearest representable Float value.
///
/// It uses the method described in:
/// http://allendowney.com/research/rand/downey07randfloat.pdf
///
/// # TODO:
///
/// * Check that it work as intended. I have not analyzed it more than
/// by using the simple `test()` function at the end of this program.
///
/// * Make it work for any other range, eg
/// `0 ..< 1` (half open instead of closed) and
/// `-Float.greatestFiniteMagnitude ... Float.greatestFiniteMagnitude`.
///
/// - Parameters:
/// - prng: The pseudorandom number generator to use.
///
func pseudorandomFloatInClosedUnitRange<R>(using prng: inout R) -> Float
where R: PseudorandomNumberGenerator
{
// Start by choosing an exponent (with correct probability):
var exp = 126 &- prng.next().trailingZeroBitCount
if exp == 62 { exp &-= prng.next().trailingZeroBitCount }
if exp < 0 { exp = 0 }
// Note that `exp` 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):
let 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)
}
//-----------------------------------------------------------------------------
// MARK: - Two functions that are similar but not what we want, for comparison
//-----------------------------------------------------------------------------
/// This function returns a value in the half open range `0 ..< 1` using the
/// method that is currently (2020-01-30) in the stdlib, ie the one which will
/// only choose between 2**24 values evenly spaced across the range, 0.5.ulp
/// apart, meaning many representable values will never be returned.
///
/// This is probably sufficient for most purposes (might even be preferrable
/// since it's more efficient). But there are cases where it's nice to be able
/// to produce *any* representable `Float` value within any finite range.
///
/// - Parameters:
/// - prng: The pseudorandom number generator to use.
///
func pseudorandomFloatInHalfOpenUnitInterval<R>(using prng: inout R) -> Float
where R: PseudorandomNumberGenerator
{
return Float(prng.next() &>> 40) * (Float.ulpOfOne/2)
}
/// This function simply returns 0 or 1 with equal probability.
///
/// - Parameters:
/// - prng: The pseudorandom number generator to use.
///
func pseudorandomFloatZeroOrOne<R>(using prng: inout R) -> Float
where R: PseudorandomNumberGenerator
{
return Float(prng.next() & 1)
}
//-----------------------------------------------------------------------------
// MARK: - Demo
//-----------------------------------------------------------------------------
func test() {
var prng = WyRand()
let numRuns = 10
let floatsPerRun = 10_000_000
var meanSum = Double(0)
var minTimePerRun = Double.infinity
for _ in 0 ..< numRuns {
let t0 = CACurrentMediaTime()
var sum = Double(0)
for _ in 0 ..< floatsPerRun {
// Uncomment the function to test here:
// The function that returns *any* representable value in [0...1]:
let v = pseudorandomFloatInClosedUnitRange(using: &prng)
// 0.036 seconds per run on my machine.
// The function that returns *some* representable value in [0..<1]:
// let v = pseudorandomFloatInHalfOpenUnitInterval(using: &prng)
// 0.013 seconds per run on my machine.
// The function that returns 0 or 1, used as a baseline:
// let v = pseudorandomFloatZeroOrOne(using: &prng)
// 0.011 seconds per run on my machine.
sum += Double(v)
}
let t1 = CACurrentMediaTime()
let timeForThisRun = t1 - t0
minTimePerRun = min(minTimePerRun, timeForThisRun)
print(t1 - t0, "s | Mean float value: \(sum / Double(floatsPerRun))")
meanSum += sum
}
print("Total mean float value:", meanSum / Double(floatsPerRun * numRuns))
print("Minimum time per run:", minTimePerRun)
}
test()