SE-0202: Random Unification

Would it be possible to also supply a generator that can be seeded that is fast and for the same seed produces a known sequence. I am not proposing this as the default but as part of the library for people that want to do Monte Carlo type simulations etc.

It is safe to read the same descriptor from different threads, however we should probably use a mutex to prevent reading at the same time. We can move this discussion to the implementation pr for specifics.

I came into this with the idea that stateful generators are generally classes and non-stateful ones are structures (Random).

@Alejandro, is this the semantics you're going for? Could you clarify the proposal?

I’ll let Alejandro speak for his own intent, but the implementation in the proposal does have this semantics (possibly up to rounding error of a few ulps, perhaps magnified by the issue Jens flagged).

This line of discussion started with this remark:

“Uniform in representable values” would make floating point randomness unusable for almost every numerical application I can think of: generating other distributions, Monte Carlo methods, game environment seeding, generating audio noise, etc. I’ve never heard of a random number library working that way. (Does such a library exist?) I don’t think the proposal’s choice of “numerically uniform” really requires further clarification or justification; this is pretty standard.

1 Like

A couple of folks have referred to PRNG selection as if it were much too complex a problem to be tackled by mere mortals. This is largely FUD. Any one of several common PRNGs would make a fine default and could easily stand as the sole PRNG included in the stdlib.

PRNGs do not in fact differ all that much in their performance. The modern ones deliver vastly better randomness quality than the old C library standards but are still pretty darn quick. The Wikipedia page for Xorshift, which is in the same family as (and runs at about the same speed as) the Mersenne Twister, claims that it can generate a random 32-bit value in 10 clock cycles on x86 architectures. That is fast enough that performance simply isn't a relevant concern when building a general-purpose library that wants to default to generating high-quality output. Anyone that needs a random number in less than 10 clock cycles is welcome to port their own favorite PRNG hack.

Just to put some actual numbers on this, the hoary old C library rand, a linear congruential generator, generates a billion random ints in 6.75 seconds on my desktop. Xorshift compiled from C code does it in 9.07 seconds and a C implementation of the Mersenne Twister does it in 10.36. So were talking about roughly a 35% performance penalty for modern generators that leave essentially nothing to be desired in terms of period and randomness quality. (Interestingly, the random library routine, which is a "nonlinear additive feedback" generator that is supposedly significantly better than the traditional linear congruential generators, is even faster than rand. I'm not sure why GameplayKit opted to include an LCG instead of this.)

Similarly, memory consumption is also a complete nonissue. I couldn't even identify a modern algorithm that uses more than a handful of bytes in its state. You could instantiate hundreds of thousands of RNG instances and no phone is going to bat an eyelash. (But in practical terms, the number of separate generators that you'd want to actually use is most likely less than 10.)

Quality and algorithm selection are exquisitely difficult problems for cryptographically secure RNGs because the stakes are standards are both very high. But there's no reason to let this free-floating anxiety slop over into PRNG selection. You don't have to select an absolutely perfect PRNG for all applications. You just have to do better than the crap that's included in the C library and follow generally agreed-upon industry standards.

The world has already selected a standard PRNG, namely, the Mersenne Twister. It's the default for SPSS, SAS, Microsoft Excel, GAUSS, GLib, GNU Multiple Precision Arithmetic Library, GNU Octave, GNU Scientific Library, gretl, IDL, Julia, CMU Common Lisp, Embeddable Common Lisp, Steel Bank Common Lisp, Maple, MATLAB, Free Pascal, PHP, Python, R, Ruby, SageMath, Scilab, and Stata. Just in case they didn't jump out at you: SPSS. SAS. OCTAVE. R. MATLAB. STATA.

Personally, I think there's some value in investigating a few of the newer generators such as Xorshift, but if the prime concern is that people are going to point and make faces at the Swift stdlib because it uses an unusual PRNG, by all means just take the Mersenne Twister.

8 Likes

If I understand it correctly, that code in fact chooses uniformly among some number (252 for Double, 223 for Float) of evenly-spaced values that span the range. Thus in general there will be many representable values in the range which are impossible to produce. For example, the range 0..<1 contains approximately 262 Double and 230 Float values.

The suggested alternative, namely to operate as though a real number were chosen uniformly from the range and then rounded to a representable value, is more subtle. It would make every representable value in the range possible to produce, with probability proportional to the gap between it and the next.

I have made an attempt at implementing such a function, and it works almost-correctly. There are some thorny rounding issues that I haven’t quite sorted out, though I’m sure @scanon knows how to do it right.

A sketch of one possible implementation:
extension BinaryFloatingPoint {
    static func random<T: RandomNumberGenerator>(in range: Range<Self>, using generator: T) -> Self {
        // Handle the easy cases
        precondition(!range.isEmpty)
        guard range.lowerBound.isFinite && range.upperBound.isFinite else {
            return .nan
        }
        if (range.lowerBound.nextUp == range.upperBound) {
            return range.lowerBound
        }
        
        // Subdivide the interval into 2^64 sections, choose one at random,
        // and repeat recursively until the subinterval has only one value.
        func randomHelper(_ subrange: Range<Self>) -> Self {
            let n: UInt64 = generator.next()
            
            // FIXME - This will underflow if Self.exponentBitCount is small
            let scale: Self = 0x1p-64
            
            // To avoid overflow here, we multiply before subtracting.
            // FIXME - Need to handle underflow
            let stepSize = (subrange.upperBound * scale) - (subrange.lowerBound * scale)
            
            // FIXME - Need this to be rounded down, not rounded to nearest
            let a = subrange.lowerBound + Self(n) * stepSize
            
            let b = a + stepSize
            if (a == b) || (a.nextUp == b) {
                return a
            }
            return randomHelper(a ..< b)
        }
        
        // Do the subdividing thing, and if any weird rounding issues make the
        // value fall outside the target range, then try again until it works.
        var result = randomHelper(range)
        while !range.contains(result) {
            result = randomHelper(range)
        }
        return result
    }
}

Note that the vast majority of the time the randomHelper() function returns a value on its first pass, with no recursion necessary. For Double it will need a second pass approximately 1 in 1024 times, and for Float only 1 in 240.

Below is a possible seeded generator that the library could provide. Whilst implementing this I noticed that next in RandomNumberGenerator should be mutating.

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

/// Swift version of Xorshift+ from: https://en.wikipedia.org/wiki/Xorshift
struct Xorshift128Plus: RandomNumberGenerator {
    private var xS: UInt64
    private var yS: UInt64
    
    /// Two seeds, `x` and `y`, are required for the random number generator (default values are provided for both).
    init(xSeed: UInt64 = 0, ySeed:  UInt64 = UInt64.max) {
        xS = xSeed == 0 && ySeed == 0 ? UInt64.max : xSeed // Seed cannot be all zeros.
        yS = ySeed
    }
    
    mutating func next() -> UInt64 {
        var x = xS
        let y = yS
        xS = y
        x ^= x << 23 // a
        yS = x ^ y ^ (x >> 17) ^ (y >> 26) // b, c
        return yS &+ y
    }
}

var rng = Xorshift128Plus()
print(rng.next()) // 18446743798831644671
print(rng.next()) // 18446743798823260096
print(rng.next()) // 16140918656667222015
print(rng.next()) // 13835128698895859711
print(rng.next()) // 16140954115928756175

The generator is fast, passes BigCrush even when reversed, and takes up little memory.

2 Likes

It's probably just the bad example seed, but those are some very suspect looking pseudorandom numbers.

You may be misunderstanding. What I said, and I think others as well, is that PRNG selection is actually that. There is a selection to be made. It‘s not insanely difficult, probably not difficult at all once you understand the basics, but there are tradeoffs to be made.

But the main argument for me is that whatever would get chosen could not be changed ever. This could leave us in a situation where next year a way better algorithm (better magic numbers more likely, but bear with me) might be found, but Swift would forever be left with a subpar old generator, which eventually nobody should use anymore. I don‘t want to end up in a place where newcommers ask for seeded PRNG and the answer is „yeah, there is this old thing in the stdlib, but don‘t ever use that, it‘s horribly outdated“.

1 Like

If we can't even change an implementation detail like the choice which RNG should be used by default, I predict Swift will become really ugly quite soon:
After all, the ideal generator would be completely unpredictable, so if someone actually relies on deterministism, he should use a custom generator in the first place.

Nonetheless, it might be convenient to be able to specify a custom generator to be used as default (but besides testing, I see little reason to do so).
Imho the same is true for the cryptography context: When someone without knowledge of the problems with randomness designs encryption, the default PRNG won't make a difference for the quality of the results.

The ideal PRNG, in the context of the post you're replying to, is not completely unpredictable. It's deliberately entirely predictable given the seed. This would be in addition to the (unseedable) secure RNG that is the default outlined in the proposal.

3 Likes

@Alejandro, here's another alternative which hides the default generator.

public protocol RandomNumberGenerator : AnyObject {
  func next() -> UInt64
}

extension Bool {
  public static func random(
    using generator: RandomNumberGenerator? = nil
  ) -> Bool {
    var number: UInt64 = 0
    if let generator = generator {
      number = generator.next()
    } else {
      withUnsafeMutableBytes(of: &number) { bytes in
        _stdlib_random(bytes)
      }
    }
    return number % 2 == 0
  }
}
  • The size of the RandomNumberGenerator? parameter is two words.

  • There's no need for the other next() and next(upperBound:) methods.

  • Access to the default generator would only be via Bool, Double, Int, etc.

All of the real subtlety is in what you do at the endpoints. Under the interpretation where you uniformly choose a value with fixed granularity, this is straightforward (since it's exactly the same as the integer case).

Under the "sample a uniform real number on the interval, and round to the nearest representable value" interpretation, it's considerably more subtle. The first issue: there are a bunch of real numbers in [0,1) that round to 1.0, but most people will not expect (0..<1.0).random() to return 1.0. Truncating instead of rounding addresses this, but then you give double-weight to zero in any interval that spans it, so you actually probably want to round towards minus infinity.

Rounding this way also has the nice feature that (1..<2).random() produces all values in the binade with equal probability.

3 Likes

Yes, that is what the 3rd FIXME comment in the code under the disclosure triangle in my previous post was getting at. I want it to round down, I just don’t know how to actually *make it* round down. (I am aware of the existence of IEEE 754 floating-point rounding modes, but I don’t know how to change them in Swift.)

The 2nd FIXME deals with situations like (-1.0)..<(0.5.ulp.nextDown) where, as written, positive numbers will never occur because the upper bound is too small (relative to lower bound) to affect stepSize. Also, really small bounds that underflow to 0 when multiplied by 2−64, like Double.leastNormalMagnitude.

The 1st is just the implicit assumption that 2−64 is representable, which it might not be for some super-small floating point type.

Why? By this logic, no software system could ever provide a PRNG, because they're all operating under these same constraints. And yet they do...

It seems to me the solution to this issue is entirely trivial. Part of the benefit of putting RNGs under a common API is that it makes them easily substitutable. Just as Random.default makes no specific claims about exactly what RNG you'll be connecting with, neither should PseudoRandom.default or PseudoRandom(seed:).

If the default PRNG were to change, clients could easily instantiate, e.g., MersenneTwisterPseudoRandom(seed:) to recover the prior behavior. Or if they're particularly forward-looking and anticipate a need for long-term reproducibility, they can just bypass the default and instantiate that specific PRNG in their original code.

2 Likes

(hope this isn't too off-topic but here I go)

The rationale behind min() and max() being functions is that they take linear time to compute (have to scan the entire collection), while first and last only take constant time. While random() does take constant time too, it returns different values even when the collection remains constant. In conclusion, there are two different requirements at work here:

  • Computed properties have to return in constant time.
  • Computed properties can only change if the receiver changes.
5 Likes

They do, and most use, like you said, Mersenne Twister. Now have a look at this CS Stackexchange question and answers about MT. Mersenne Twister is not a great PRNG by todays standards. It performs worse in BigCrush than Xorshift variants (which have already been mentioned here), while being slower, completely predictable after only a few outputs and using a huge state array of 2.5KiB.

Yet most languages ship with this as their PRNG. Why? Because they can't change to a better one without breaking compatibility. MT being the PRNG in all those languages you mentioned is a perfect example of the pitfall I want Swift to avoid.

This would defeat the purpose of using a seeded PRNG where you absolutely need the promise that output for one seed will stay the same.

Changing the PRNG would thus be a change that is source-breaking (at run-time no less), so that's unlikely to happen.

1 Like

I cannot see why the library can't provide one seeded random number generator (rng) at first and then at some point in the future provide others. I don't think you can retire an old rng but you can add new ones, adding new ones is not a breaking change and therefore no problem. I would suggest using the name of the algorithm rather than 'default', because you generally want to know which generator you are using.

3 Likes

This sounds like a perfect role for one of the “non-standard libraries” that people keep talking about. A well-curated statistics module would be a welcome addition to the ecosystem. There is no need to burden the standard library with multiple random number generators, at least not right away.

4 Likes