[stdlib] Floating-point random-number improvements by NevinBR · Pull Request #33560 · swiftlang/swift

Overview

This patch resolves multiple issues with generating random floating-point numbers.

The existing random methods on BinaryFloatingPoint will crash for some valid ranges (SR-8798), cannot produce all values in some ranges (SR-12765), and do not follow the proposed and documented semantics. This patch solves these problems:

  • SR-8798: BinaryFloatingPoint.random(in:) crashes on some valid ranges
  • SR-12765: BinaryFloatingPoint.random(in:) cannot produce all values in range

Summary of changes

i) Finite ranges where the distance between the bounds exceeds greatestFiniteMagnitude previously caused a trap at run-time. Now (with this patch) they are handled correctly to produce a uniform value in the range.

ii) Generating a random floating-point value in -1..<1 left the low-bit always 0. If the magnitude was less than 1/2, then the lowest 2 bits would be 0. If less than 1/4, 3 bits, and so forth. Other ranges which spanned multiple binades had similar issues. Now all values in the input range are produced with correct probability.

iii) The proposal (SE-0202: Random Unification) which added random number generation to Swift was quite sparse in its mention of floating-point semantics. However, during the review thread, discussion about the intended semantics arose with comments like this:

I think it makes sense to have the floating-point implementation behave as though a real number were chosen uniformly from the appropriate interval and then rounded.

Agreed. I would consider any implementation that does not have this behavior (within an ulp or two) to be incorrect.

To which the author of the proposal responded:

Yes, these are the semantics that I’m trying achieve here (those being uniform in real numbers).

Concordantly, the documentation comments for the floating-point random methods state:

  /// The `random(in:using:)` static method chooses a random value from a
  /// continuous uniform distribution in `range`, and then converts that value
  /// to the nearest representable value in this type.

However, prior to this patch, the implementation did not match that behavior, and many representable values could not be produced in certain ranges.

Mathematical details

In order to achieve the desired semantics, it is necessary to define precisely what “converts that value to the nearest representable value” should mean. This patch takes the following axiomatic approach:

Range

Single-value ranges: random(in: x ..< x.nextUp) always produces x.

Adjacent intervals: If x < y < z, then random(in: x ..< z) is equivalent to generating random(in: x ..< y) with probability (y-x)/(z-x), and otherwise random(in: y ..< z) with the remaining probability (z-y)/(z-x).

In order to satisfy these two principles, random(in: x ..< y) must behave as if a real number r were generated uniformly in [x, y), then rounded down to the nearest representable value. Note that the rounding must be downward, as any other choice would violate one of the two principles.

ClosedRange

Subintervals: If x <= y < z, then repeatedly generating random(in: x ..< z) until the result lands in x ... y is equivalent to generating random(in: x ... y).

This rule ensures consistency of results produced by the Range and ClosedRange versions of random(in:). As a result, it also guarantees that partitioning a closed interval into disjoint closed subintervals is consistent as well.

In order to satisfy this principle, random(in: x ... y) must be equivalent to random(in: x ..< y.nextUp) if the latter is finite.

In the edge-case that y == .greatestFiniteMagnitude, we utilize the adjacent intervals principle on [x, y) and [y, y + y.ulp). Although the latter endpoint is not representable as a finite floating-point value, the conceptual idea still holds, and the probability of producing y is proportional to y.ulp just as it is for all other values in the same binade.

This patch implements those semantics.

Similarity with random integer methods

It is interesting to note that the strategy of generating a uniform real number in a half-open interval then rounding down, is equivalent to how the random methods work for integer types. That is, T.random(in: x ..< y) behaves as if choosing a real number r uniformly in [x, y) then rounding down to the next representable value, regardless of whether T is an integer or (with this patch) a floating-point type.

Similarly for closed ranges, T.random(in: x ... y) behaves as if extending to a half-open interval bounded above by the next representable value larger than y (or in case of overflow, then where that value would be as a real number), generating a random value in the new range, and rounding down.