Random deviates from standard algortihm

A customer filed the following bug report on Connect:

While investigating the period of various random number generators, I found a serious bug in Microsoft .NET System.Random implementation.

Microsoft Documentation says that the Random class is based on Donald E. Knuth's subtractive random number generator algorithm. For more information, see D. E. Knuth. "The Art of Computer Programming, volume 2: Seminumerical Algorithms". Addison-Wesley, Reading, MA, second edition, 1981.

The problem was discovered when I use .NET Reflector to see what is actually implemented. Knuth was very specific about the lagged Fibonacci sequence coefficient (24 and 55). Somehow Microsoft mistyped the value to be 21 (this.inextp = 0x15 in public Random(int Seed) in source code), in place of 31 (=55-24). Due to this mistype, the random number generated no longer has the guanrantee on the period to be longer than (2^55-1). The Knuth version of lags (24, 55) has been tested by many people about the property of the generated numbers but not the one created by Microsoft.

It is very easy to identify the typo. It seems that Microsoft just copied the routine (ran3) from Numerical Recipes v.2 page. 283 verbatim (even the variable names are copied), not from Knuth book which has a complete different way of initializing. You see that Numerical Recipe version uses 31 correctly, not 21 used by Microsoft.

Looking at the sources, it seems it's still the case today.

I'm not sure what the implication of our deviation is. If we wanted to fix it, we'd likely would have to quirk the implementation on .NET Framework. For .NET Core, we can debate the merits of quirking, but it's likely we break customers with changing the seed.

Thoughts?

Related/duplicate: #18996

@Clockwork-Muse

@Clockwork-Muse

I checked earlier to see if there is an existing bug but #18996 isn't a dupe, but it's related. The customer issue is more specific to our existing algorithm.

I spent a bit of time pondering about the implications of the bug; see this gist for a write-up. TL;DR: the proof for the period in the $(k, l) = (55, 24)$ case can not be applied, but test batteries don't seem to notice, since other implementation quirks appear to matter more.

The wording in the documentation is:

The current implementation of the Random class is based on a modified version of Donald E. Knuth's subtractive random number generator algorithm. [Emphasis mine].

So on the one hand in saying it's a modified version, deviation from that in Knuth 1981 isn't a bug as such. On the other hand, the use of current suggests that no promise is made to keep to the algorithm in use, so if there's a problem in the deviation it can be changed.

It should probably be noted that the deviations come in quite different flavors, some of which are less problematic than others:

The recurrence in Next() is performed modulo $2^31 - 1$, while much of the relevant theory in Knuth is concerned with what happens when the modulus is a power of two. As far as I can tell, this actually makes the generator more robust towards conventional statistical tests, since the parity of the output is no longer determined by the parities of the relevant part of the state. (Details in the gist mentioned above.)
In Next(Int32, Int32), when the distance between the arguments is large, the generator mixes together two outputs of the generator. When this scheme is used, PractRand would not find any bias at all in 512GB of random numbers, so this deviation is actually very good (although I had to think for a moment to convince myself that it did not break the RNG entirely).
Then finally, there the one mentioned in the post, that the generator accidentally uses $(k, l) = (55, 34)$ instead of $(k, l) = (55, 24)$, meaning that one no longer knows the period of the generator. As mentioned in the gist, it seems that PractRand is unable to really tell the difference between the two, so from a practical point of view, you might still be good, but as I see it, the discrepancy will almost universally be a bad thing, as you'll want your RNG to have a solid theoretical foundation (but then again, the theory breaks already by not using a power of two as your modulus).

Now if you do opt for the breaking change, I will agree with the poster of the previous issue that it would make sense to take the opportunity to replace the algorithm entirely, rather than fixing only the value of $l$, simply because Next() does succumb to the tests for relatively small inputs.

Edit: I ended up spelling out the comment below in this blog post but let me leave the original comment here for ease of reference.

Somewhat unrelated of the above, here's a concrete example of very strong bias showing that Next(0, int.MaxValue) is completely broken:

The algorithm for Next(0, b) is roughly the following:

Generate a random integer between 0 and 2^31 - 1.
Divide by 2^31 - 1 to get a floating-point between 0 and 1.
Multiply by b and take the integer part of the result.

For now, let us assume that we are doing a good job at the first step (which might not exactly be the case due to other issues mentioned above). We are then left with the floating-point conversions, which, to me at least, have a highly opaque impact on randomness.

Let us consider the special case where b = int.MaxValue. If i is the outcome of the first step, the last two steps boil down to x = (int)((i * (1.0/int.MaxValue)) * int.MaxValue). It is not hard to show that x == i - 1 whenever $i = 2^{n+22} + 2^{n+2}k + 2^n$ for some $n \in { 0, \dots, 8 }$, $k \in { 0, \dots, 2^{21} - 2^{20} - 1 }$, and that x == i otherwise. That is, the program below outputs no errors:

#include <stdio.h>

int main()
{
    int maxInt = 2147483647;
    double maxIntRec = 1.0 / maxInt;
    int n = 1;  // Index of most significant bit
    int modulus;
    int period;
    for (int i = 1; i < maxInt; m++)
    {
        // When i == 2^n, update criterion for i == x vs. i == x - 1
        if (i == 1 << n) {
            printf("i = 2^%d\n", n);
            modulus = 1 << n - 22;
            period = (1 << n - 20) - 1;
            n++;
        }
        int x = (int)(i * maxIntRec * maxInt);
        // Test i % 2^(n-20) == 2^(n-22)
        if ((i & period) == modulus) x++;
        if (x != i) printf("Mistake at i: %d\n", i);
    }
    return 0;
}

Under the assumption that i is taken uniformly at random, this means that the probability that x is odd is about 0.5034, which in turn means that Next(0, int.MaxValue) suffers from strong bias in its least significant bit. To see what this means in practice, let us do a quick coin flip simulation;

public static void Main()
{
    var rng = new Random(0);
    var sampleSize = 1000000;
    var numbers = Enumerable.Repeat(0, sampleSize).Select(_ => rng.Next(0, int.MaxValue));
    var heads = numbers.Count(x => x % 2 == 1);
    var tails = sampleSize - heads;
    Console.WriteLine($"Heads: {heads}, tails: {tails}");
    // Output: Heads: 503291, tails: 496709
}

Now on the other hand, the likelihood of seeing such an outcome, or something more extreme, under the hypothesis of fairness of the coin is about 4.67e-11:

In [1]: import scipy.stats as st

In [2]: 2*st.binom(1000000, 0.5).cdf(496709)
Out[2]: 4.6722165389139607e-11

Now you might argue that using Next(0, int.MaxValue) is a bit silly when we already have Next(), but that does not keep anybody from doing it; moreover, it is not clear to me that the floating-point conversions will not also be problematic for other values of b.

@fuglede

I've been working in the area of PRNGs recently and looking at different methods and implementations, including System.Random, so will summarise my thoughts here while they are fresh.

NextDouble()

The IEEE754 double precision format has 53 bits of precision (52 explicit, and one implicit) and as such is capable of representing 2^53 distinct values in the interval [0,1). Random.NextDouble() will generate only 2^31 of those possible values, i.e. only one in every 2^22 possible values are generated.

There are two aspects of this class that produce this behaviour:

(1) The underlying PRNG method generates 32 random bits per iteration. It could be called twice to obtain more bits, but the current implementation doesn't do this.

(2) One of the 32 random bits is always discarded.

InternalSample() contains the PRNG logic, and all other methods call this method directly or indirectly as the source of random bits. However, InternalSample() returns a non-negative Int32, i.e. the sign bit is always zero. Thus all other code in this class is based on working with a 31 bit random source.

Next(int minValue, int maxValue)

This method must handle a special case where the range of possible values is larger than Int32.MaxValue (2^31-1). E.g. in the extreme case:

minValue = Int32.MinValue
maxValue = Int32.MaxValue

And therefore:

range = Int32.MaxValue - Int32.MinValue
= 2^32-1

I.e. in order to be able to generate all possible 2^32-1 values, 32 bits of randomness are required. Since InternalSample() generates only 31 bits it is called twice, and one extra bit is obtained from the second call and merged with the other 31 bits. This all happens within GetSampleForLargeRange() using a rather convoluted and inelegant approach. However, the whole approach would be unnecessary if a PRNG source that generates more bits was used.

NextBytes(byte[])

Invokes the PRNG to obtain 31 bits, and uses only 8 of them (the comment on the method says 7 bits, but this is incorrect). This is fine, but could be a lot faster if we e.g. generated 32 bits and filled 4 bytes per PRNG iteration.

Construction

The PRNG in use appears to have a very computationally expensive initialisation routine in which it constructs a 56 byte table. Thus if you want to rapidly initialise lots of instances (e.g. with fixed seeds to give multiple deterministic random sequences), System.Random will be a very slow way of achieving this.

Overall Performance

Below are some numbers from a BenchmarkDotNet based comparison between System.Random and xoroshiro256**. The specific implementation used is:

Xoshiro256StarStarRandom.cs

This is a recently published PRNG method that passes all known statistical tests, generates 64 bits per iteration, and is very fast to compute and initialise. I believe the above implementation addresses all of the issues I have mentioned in this post, and also the concerns of @fuglede regarding bias.

Method	System.Random	Xoroshiro256StarStar
Next() [10M]	86.63 ms	25.23 ms
Next(int) [10M]	111.61 ms	41.31 ms
Next(int,int) [10M]	116.98 ms	47.21 ms
NextDouble() [10M]	94.97 ms	23.29 ms
NextBytes() [100M]	816.58 ms	15.43 ms

The timings are for 10 million calls (e.g. 20ms/10,000,000 = 2ns per call) apart from NextBytes() which is for 100 calls of filling a byte array of length 1,000,000 (xoshiro256** is about 50x faster mostly because it fills 8 bytes per iteration instead of 1, and also uses unsafe pointers to assign the 8 bytes in one operation)

@colgreen

@colgreen: You're not saying anything suggesting the contrary but you may want to note that (as mentioned also in the gist referenced above) the fact that two runs are used to generate the full 32 bits actually appears to make the existing algorithm for Next(Int32.MinValue, Int32.MaxValue) much stronger, at least in the eyes of PractRand.

I didn't go through your implementation, but you're right, the concrete example of bias that System.Random suffers from does not appear to be an issue here. (Concretely, in System.Random it is caused by the division by int.MaxValue which you do not need.)

Regarding xoshiro256** in particular, the author of PCG does not seem overly impressed:

But in fact we don't have to fully invert the output function to see problems. I tried just multiplying by 0x8e38e38e38e38e39 and found that xoshiro256** failed statistical tests. Then I tried just multiplying by 0x8e39 and that failed too. Finally, I tied multiplying by 0x39 (a.k.a. 57) and that failed as well—in only 9.5 seconds.

http://www.pcg-random.org/posts/a-quick-look-at-xoshiro256.html

the fact that two runs are used to generate the full 32 bits actually appears to make the existing algorithm for Next(Int32.MinValue, Int32.MaxValue) much stronger, at least in the eyes of Practrand.

Yeh, superficially it makes sense to me that sampling an extra bit from one more iteration (which also moves all further samples on by one iter) would likely improve the quality of the generated randomness. If however the underlying PRNG was generating 32bits of good quality pseudo-randomness then all this would be unnecessary and we could obtain a significant performance boost.

caused by the division by int.MaxValue which you do not need.

Yeh, it just seems like problems that shouldn't exist have been addressed by additional layers of complexity with their own problems. So at some level the source I posted is an attempt to go right back to basics and avoid the unnecessary complexity.

Regarding xoshiro256** in particular, the author of PCG does not seem overly impressed

Interesting (and concerning!). Thanks for the link. Concerns about the underlying PRNG itself aside, the rest of the logic in the class I linked, i.e. for converting 64 random bits to doubles, integer ranges etc. should all be good and independent of the PRNG concerns. But yeh, will look into switching to PCG.

Thanks.

@fuglede

@fuglede

Regarding xoshiro256** in particular, the author of PCG does not seem overly impressed

Here's a response from Sebastiano Vigna:

https://www.reddit.com/r/programming/comments/8gx2d3/new_line_of_fast_prngs_released_by_the_author_of/dyqe1u4/

@colgreen

@colgreen:

Here's a response from Sebastiano Vigna:

Thanks for the reference. I will have to keep an eye on that one; looks like it might take a bit of popcorn though.

Concerns about the underlying PRNG itself aside, the rest of the logic in the class I linked, i.e. for converting 64 random bits to doubles, integer ranges etc. should all be good and independent of the PRNG concerns.

One thing you may want to pay attention to -- just in case you aren't already -- is the expected behavior of NextDoubleInner; and here we are getting slightly off-topic for the GitHub issue but let me leave the comment regardless as it will be something for the .NET teams to also take into account if they decide to fix System.Random. By limiting yourself to 2^53 outputs you get the nice property that every point in the image appears with equal probability (under the assumption of uniform NextInnerULong) but the method will only output a small fraction of the possible doubles in [0, 1). This also means that I don't quite follow this part of the first comment:

The IEEE754 double precision format has 53 bits of precision (52 explicit, and one implicit) and as such is capable of representing 2^53 distinct values in the interval [0,1).

They're annoying to count, but the number of distinct values, let's call it $N$, is much greater; somewhere below 2^62; see e.g. this StackOverflow answer.

If you want to include all of the missing doubles as possible outputs, assigning probabilities becomes trickier, as you probably don't want every double with probability $1/N$ as that would cause the expected value of NextDoubleInner to become much smaller than 1/2 since they cluster close to 0. Rather, a given positive double close to 0 would have to appear with smaller probability than one in the interval [1/2, 1). See this implementation of integer-to-double conversion for more information.

@fuglede

@fuglede

looks like it might take a bit of popcorn though.

Indeed.

By limiting yourself to 2^53 outputs you get the nice property that every point in the image appears with equal probability (under the assumption of uniform NextInnerULong) but the method will only output a small fraction of the possible doubles in [0, 1)

Yes I see. And thanks for the links - I will need some time to digest the info.

So at the very least I need to correct my comment. It is very convenient to calc N * 1/(2^53), and it is an improvement on N * 1/(2^31) (as per System.Random).

Perhaps we should move our discussion to email, and if we get anywhere then post a summary here when we're done? (if so then go so my github profile -> homepage link -> email addr at bottom of page)

Move files to shared CoreLib partition (#10830)

System.Random instance created through its no-argument constructor could theoretically use a different PRNG algorithm than the one it uses through its int constructor, while preserving backward compatibility

That might be the way forward. The parameter less ctor is the commonly used one per https://apisof.net/catalog/System.Random..ctor(Int32)

How about:

cstr() - legacy PRNG
cstr(int seed) - legacy PRNG
cstr(RandomAlgorithm prng) - A specified PRNG
cstr(RandomAlgorithm prng, int seed) - A specified PRNG, with a given seed

And:

enum RandomAlgorithm 
{
    Legacy,  // or MersenneTwister, or LegacyMersenneTwister
    Foo,     // E.g. xoshiro256**
    Bar,
    // etc.
}

Where on day one there may be only a single new PRNG in the enum, but we have the freedome in the future to add new PRNGs as necessary, e.g. as and when problems are reported with the FooPrng.

That might be the way forward.

I don't think so. The arguments given in this thread (ML scenarios, etc.) revolve around being able to reproduce the random sequence. The Random class itself might not be serializable in .NET Core, but the way you'd make this work is to serialize just the seed, then to have all applications instantiate a random instance from the persisted seed. This very strongly implies that the desired behavior from this thread is both to use a different algorithm and to have a way to reproduce the output sequence, which means we can't key the improved behavior off of the parameterless ctor.

I further don't think we should have an enum value that controls the algorithm. The two immediate concerns that come to mind are: (a) how does the caller determine which algorithm to use?; and (b) what's the bar for determining which algorithms to include in-box? For most applications we say that if you need a specialized algorithm or data structure that we don't provide in-box, pull in an existing third-party implementation. Or, for frameworks that require some functionality, the frameworks can provide the implementation. In the case of ML, this means that the ML package can provide Random-derived types for all the different random number generators that are commonly used by ML scenarios.

If the intent is that the PRNG algorithm used by System.Random will never change, even in future versions of .NET (and perhaps even in the case when a System.Random instance is initialized through its no-argument constructor), then the System.Random class should stop stating that the algorithm is subject to changes and should document the current algorithm in full, flaws and warts and all, just as Java does with java.util.Random.

For most applications we say that if you need a specialized algorithm or data structure that we don't provide in-box, pull in an existing third-party implementation.

In light of this, one wonders if we should stop providing all the random "helpers" (as in, NextInt, etc), and just provide random bytes from the OS, expecting third party libraries to handle this sort of thing.

This is one of the times I wish we'd made packages for almost everything in the standard library, since that would allow us to "fix" this easier.

@GrabYourPitchforks

@GrabYourPitchforks

We can't have it both ways. Either Random utilizes a stable algorithm under the covers (as has always been the case in .NET Framework, and has historically been the reason we've not made changes in .NET Core, and has also been the motivating factor for issues like #27429 and #18996 (comment) and #6203 (comment)). Or we're free to improve the algorithm (which is hinted at on MSDN, but in a very loose fashion). But we can't exist in both worlds concurrently.

So... which world do we want to be in? :)

In a world of in-place updates changing the algorithm isn't doable, but that's not the world we live in with .NET Core anymore, so I'd argue it's false dichotomy :-)

We could decide to introduce a new constructor/static factory method on the same type or we could decide to make a breaking change between major releases (potentially with an app compat switch). These are all better outcomes than a new type.

We could decide to introduce a new constructor/static factory method on the same type or we could decide to make a breaking change between major releases (potentially with an app compat switch). These are all better outcomes than a new type.

I don't think this puts us in a better world. This discounts that libraries might be using the Random type and might be expecting it to have a given algorithm. An app-wide switch isn't workable for this scenario since it removes the ability of the individual callers to have control over their logic.

Edit: the reason this matters is that Random has always been treated as a serialization equivalent, even if the type isn't marked [Serializable]. If a client on Full Framework and a server on Core agree on the seed, the guarantee we have always made is that their outputs are the same. We even have a unit test for this.

FWIW here's where I see things landing:

If the caller currently provides a seed, they require consistency with previous versions. There's quite literally no reason for them to have provided a seed other than to have reproducible results. Us changing the algorithm out from under them would break these scenarios.
If the caller doesn't provide a seed, no back-compat is needed. By definition they didn't require reproducible output, nor did they rely on any specific algorithm.
If the caller needs high-quality random data, their scenario should mandate a tolerance that they're willing to accept. This implies that the scenario itself will dictate what algorithm is to be used, or it's up to the individual caller to ensure that they're using an algorithm that fulfills those requirements. In neither case is a general "we give you some random data but we don't tell you what the algorithm is, and BTW the algorithm is subject to change" API sufficient for this purpose.

My argument hinges on that it's a bad idea for us to change (1) and silently break consumers; and that for (3) callers should be calling specific APIs which are contracted to fulfill their specific requirements, but it's not up to the Framework to provide such algorithms in-box.

@GrabYourPitchforks

I broadly agree with @GrabYourPitchforks, but bear in mind that the PRNG is only part of the problem. There are other issues related to how the various methods (Next, NextDouble...) convert the random bits into the required return types.

It would still be necessary to document well what is going on.

Another quirk of System.Random (that I once encountered causing a problem in the wild), is that the parameterless constructor uses the system clock as the basis for the seed, with the effect that if you construct multiple Random(s) in quick succession, they can all have the same seed! In the Redzen.Random classes I used another PRNG as the source of seed values, which itself is seeded by RNGCryptoServiceProvider - which I believe is seeded by system entropy. I also uses a low locking strategy to ensure this would scale reasonably well if anyone ever did want to construct instances at a high rate (which is questionable, but real world code gets messy sometimes).

See:

https://github.com/colgreen/Redzen/blob/master/Redzen/Random/DefaultRandomSeedSource.cs

Another quirk of System.Random (that I once encountered causing a problem in the wild), is that the parameterless constructor uses the system clock as the basis for the seed, with the effect that if you construct multiple Random(s) in quick succession, they can all have the same seed!

This is true in .NET Framework, but not .NET Core. The latter uses a private Random based on OS entropy to initialize a per-thread Random that generates seeds. I actually ran into this a few days ago, and .NET Framework was having exactly this problem, but .NET Core did not.

@GrabYourPitchforks

I broadly agree with @GrabYourPitchforks, but bear in mind that the PRNG is only part of the problem. There are other issues related to how the various methods (Next, NextDouble...) convert the random bits into the required return types.

It would still be necessary to document well what is going on.

In this sense, see again how this is done, for example, in the documentation for java.util.Random, especially in the nextInt and nextDouble methods.

Navigation Menu

Repository navigation

Random deviates from standard algortihm #23198

NextDouble()

Next(int minValue, int maxValue)

NextBytes(byte[])

Construction

Overall Performance

13 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Random deviates from standard algortihm #23198

Description

Activity

Clockwork-Muse commented on Aug 16, 2017

terrajobst commented on Aug 16, 2017

fuglede commented on Aug 20, 2017

JonHanna commented on Aug 20, 2017

fuglede commented on Aug 21, 2017

fuglede commented on Aug 27, 2017

colgreen commented on May 9, 2018

NextDouble()

Next(int minValue, int maxValue)

NextBytes(byte[])

Construction

Overall Performance

fuglede commented on May 9, 2018

colgreen commented on May 9, 2018

colgreen commented on May 10, 2018

fuglede commented on May 10, 2018

colgreen commented on May 10, 2018

13 remaining items

danmoseley commented on Aug 1, 2020

colgreen commented on Aug 1, 2020

GrabYourPitchforks commented on Aug 2, 2020

peteroupc commented on Aug 2, 2020

Clockwork-Muse commented on Aug 3, 2020

terrajobst commented on Aug 3, 2020

GrabYourPitchforks commented on Aug 3, 2020

GrabYourPitchforks commented on Aug 3, 2020

colgreen commented on Aug 3, 2020

danmoseley commented on Aug 3, 2020

peteroupc commented on Aug 3, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions