New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Random deviates from standard algortihm #23298

Open
terrajobst opened this Issue Aug 16, 2017 · 12 comments

Comments

Projects
None yet
7 participants
@terrajobst
Member

terrajobst commented Aug 16, 2017

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?

@terrajobst terrajobst added the bug label Aug 16, 2017

@Clockwork-Muse

This comment has been minimized.

Show comment
Hide comment
@Clockwork-Muse

Clockwork-Muse Aug 16, 2017

Collaborator

Related/duplicate: #12746

Collaborator

Clockwork-Muse commented Aug 16, 2017

Related/duplicate: #12746

@terrajobst

This comment has been minimized.

Show comment
Hide comment
@terrajobst

terrajobst Aug 16, 2017

Member

@Clockwork-Muse

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

Member

terrajobst commented Aug 16, 2017

@Clockwork-Muse

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

@fuglede

This comment has been minimized.

Show comment
Hide comment
@fuglede

fuglede Aug 20, 2017

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.

fuglede commented Aug 20, 2017

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.

@JonHanna

This comment has been minimized.

Show comment
Hide comment
@JonHanna

JonHanna Aug 20, 2017

Collaborator

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.

Collaborator

JonHanna commented Aug 20, 2017

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.

@fuglede

This comment has been minimized.

Show comment
Hide comment
@fuglede

fuglede Aug 21, 2017

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.

fuglede commented Aug 21, 2017

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.

@fuglede

This comment has been minimized.

Show comment
Hide comment
@fuglede

fuglede Aug 27, 2017

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 commented Aug 27, 2017

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.

@joperezr joperezr added this to the Future milestone Sep 25, 2017

@colgreen

This comment has been minimized.

Show comment
Hide comment
@colgreen

colgreen May 9, 2018

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 commented May 9, 2018

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)

@fuglede

This comment has been minimized.

Show comment
Hide comment
@fuglede

fuglede May 9, 2018

@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

fuglede commented May 9, 2018

@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

@colgreen

This comment has been minimized.

Show comment
Hide comment
@colgreen

colgreen May 9, 2018

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.

colgreen commented May 9, 2018

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.

@colgreen

This comment has been minimized.

Show comment
Hide comment
@colgreen

colgreen May 10, 2018

@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 commented May 10, 2018

@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/

@fuglede

This comment has been minimized.

Show comment
Hide comment
@fuglede

fuglede May 10, 2018

@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 commented May 10, 2018

@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.

@colgreen

This comment has been minimized.

Show comment
Hide comment
@colgreen

colgreen May 10, 2018

@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)

colgreen commented May 10, 2018

@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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment