title: Bug in .Net's `Random` class?
body: <p>I was looking at a question that was talking about a bad implementation of the Fisher-Yates shuffling algorithm and I was perplexed that there was a bias when implemented incorrectly.</p>
  <p>The two algorithms are these:</p>
  <pre><code>private Random _random = new Random();
  
  public int[] FisherYates(int[] source)
  {
      int[] output = source.ToArray();
      for (var i = 0; i &lt; output.Length; i++)
      {
          var j = _random.Next(i, output.Length);
          (output[i], output[j]) = (output[j], output[i]);
      }
      return output;
  }
  
  public int[] FisherYatesBad(int[] source)
  {
      int[] output = source.ToArray();
      for (var i = 0; i &lt; output.Length; i++)
      {
          var j = _random.Next(0, output.Length);
          (output[i], output[j]) = (output[j], output[i]);
      }
      return output;
  }
  </code></pre>
  <p>A really subtle different, but enough to cause a massive bias.</p>
  <p>Good implementation:</p>
  <p><a href="https://i.stack.imgur.com/zW5VN.png" rel="noreferrer"><img src="https://i.stack.imgur.com/zW5VN.png" alt="Good Fisher-Yates" /></a></p>
  <p>Bad implementation:</p>
  <p><a href="https://i.stack.imgur.com/ogde9.png" rel="noreferrer"><img src="https://i.stack.imgur.com/ogde9.png" alt="Bad Fisher-Yates" /></a></p>
  <p>To be clear about these plots, I start with the numbers 0 to 99, create 10_000_000 shuffles using whichever algorithm and then I average the values in each of the shuffles to get a single set of figures. If the shuffle is trying random then all 100 figures would belong to the same normal distribution.</p>
  <p>Now, that's all fine, but I thought I'd check to see if these methods produce valid results:</p>
  <pre><code>public int[] OrderByRandomNext(int[] source) =&gt; source.OrderBy(x =&gt; _random.Next()).ToArray();
  
  public int[] OrderByRandomNextDouble(int[] source) =&gt; source.OrderBy(x =&gt; _random.NextDouble()).ToArray();
  </code></pre>
  <p>Both nice a neat, but are they a fair shuffle?</p>
  <p><code>OrderByRandomNext</code>:</p>
  <p><a href="https://i.stack.imgur.com/9XyOX.png" rel="noreferrer"><img src="https://i.stack.imgur.com/9XyOX.png" alt="OrderByRandomNext" /></a></p>
  <p><code>OrderByRandomNextDouble</code>:</p>
  <p><a href="https://i.stack.imgur.com/a6ETv.png" rel="noreferrer"><img src="https://i.stack.imgur.com/a6ETv.png" alt="OrderByRandomNextDouble" /></a></p>
  <p>Notice that the <code>1</code> and the <code>100</code> figures are significantly lower in each?</p>
  <p>Well, I thought it might be an artefact of how <code>OrderBy</code> works. So I tested it with another random number generator - one brought to use from Eric Lippert in his improving Random series.</p>
  <pre><code>public int[] OrderByBetterRandomNextDouble(int[] source) =&gt; source.OrderBy(x =&gt; BetterRandom.NextDouble()).ToArray();
  
  public static class BetterRandom
  {
      private static readonly ThreadLocal&lt;RandomNumberGenerator&gt; crng =
          new ThreadLocal&lt;RandomNumberGenerator&gt;(RandomNumberGenerator.Create);
  
      private static readonly ThreadLocal&lt;byte[]&gt; bytes =
          new ThreadLocal&lt;byte[]&gt;(() =&gt; new byte[sizeof(int)]);
  
      public static int NextInt()
      {
          crng.Value.GetBytes(bytes.Value);
          return BitConverter.ToInt32(bytes.Value, 0) &amp; int.MaxValue;
      }
  
      public static double NextDouble()
      {
          while (true)
          {
              long x = NextInt() &amp; 0x001FFFFF;
              x &lt;&lt;= 31;
              x |= (long)NextInt();
              double n = x;
              const double d = 1L &lt;&lt; 52;
              double q = n / d;
              if (q != 1.0)
                  return q;
          }
      }
  }
  </code></pre>
  <p>Well, here's the chart:</p>
  <p><a href="https://i.stack.imgur.com/EAdi0.png" rel="noreferrer"><img src="https://i.stack.imgur.com/EAdi0.png" alt="BetterRandom" /></a></p>
  <p>There's no bias!</p>
  <p>Here's my code to generate the data (run in LINQPad):</p>
  <pre><code>void Main()
  {
      var n = 100;
      var s = 1000000;
  
      var numbers = Enumerable.Range(0, n).ToArray();
  
      var algorithms = new Func&lt;int[], int[]&gt;[]
      {
          FisherYates,
          OrderByRandomNext,
          OrderByRandomNextDouble,
          OrderByBetterRandomNextDouble,
      };
  
      var averages =
          algorithms
              .Select(algorithm =&gt;
                  Enumerable
                      .Range(0, numbers.Length)
                      .Select(x =&gt;
                          Enumerable
                              .Range(0, s)
                              .Select(y =&gt; algorithm(numbers))
                              .Aggregate(0.0, (a, v) =&gt; a + (double)v[x] / s))
                      .ToArray())
              .Select(x =&gt; new
              {
                  averages = x,
                  distribution = Accord.Statistics.Distributions.Univariate.NormalDistribution.Estimate(x.Skip(1).SkipLast(1).ToArray()),
                  first = x.First(),
                  last = x.Last(),
              })
              .Select(x =&gt; new
              {
                  x.averages,
                  x.distribution,
                  x.first,
                  x.last,
                  first_prob =x.distribution.DistributionFunction(x.first),
                  last_prob = x.distribution.DistributionFunction(x.last),
              })
              .ToArray();
  
      var d = 
  
      averages.Dump();
  }
  
  private Random _random = new Random();
  
      public int[] FisherYates(int[] source)
      {
          int[] output = source.ToArray();
          for (var i = 0; i &lt; output.Length; i++)
          {
              var j = _random.Next(i, output.Length);
              (output[i], output[j]) = (output[j], output[i]);
          }
          return output;
      }
  
  public int[] OrderByRandomNext(int[] source) =&gt; source.OrderBy(x =&gt; _random.Next()).ToArray();
  
  public int[] OrderByRandomNextDouble(int[] source) =&gt; source.OrderBy(x =&gt; _random.NextDouble()).ToArray();
  
      public int[] OrderByBetterRandomNextDouble(int[] source) =&gt; source.OrderBy(x =&gt; BetterRandom.NextDouble()).ToArray();
  
      public static class BetterRandom
      {
          private static readonly ThreadLocal&lt;RandomNumberGenerator&gt; crng =
              new ThreadLocal&lt;RandomNumberGenerator&gt;(RandomNumberGenerator.Create);
  
          private static readonly ThreadLocal&lt;byte[]&gt; bytes =
              new ThreadLocal&lt;byte[]&gt;(() =&gt; new byte[sizeof(int)]);
  
          public static int NextInt()
          {
              crng.Value.GetBytes(bytes.Value);
              return BitConverter.ToInt32(bytes.Value, 0) &amp; int.MaxValue;
          }
  
          public static double NextDouble()
          {
              while (true)
              {
                  long x = NextInt() &amp; 0x001FFFFF;
                  x &lt;&lt;= 31;
                  x |= (long)NextInt();
                  double n = x;
                  const double d = 1L &lt;&lt; 52;
                  double q = n / d;
                  if (q != 1.0)
                      return q;
              }
          }
      }
  </code></pre>
  <p>Here's the data that I generated:</p>
  <pre>
  distribution                                             | first              | last               | first_prob             | last_prob            
  -------------------------------------------------------- | ------------------ | ------------------ | ---------------------- | ---------------------
  N(x; μ = 49.50267467345823, σ² = 0.0008896228453062147)  | 49.505465999987585 | 49.49833699998965  | 0.5372807100387846     | 0.44218570467529394  
  N(x; μ = 49.50503062243786, σ² = 0.0009954477334487531)  | 49.36330799998817  | 49.37124399998651  | 3.529550818615057E-06  | 1.115772521409486E-05
  N(x; μ = 49.505720877539765, σ² = 0.0008257970106087029) | 49.37231699998847  | 49.386660999990106 | 1.7228855271333998E-06 | 1.712972513601141E-05
  N(x; μ = 49.49994663264188, σ² = 0.0007518765247716318)  | 49.50191999998847  | 49.474235999989205 | 0.5286859991636343     | 0.17421285127499514  
  </pre>
  <p>Here's my question. What's up with <code>System.Random</code> and the bias it introduces?</p>