Реализация генератора случайных чисел Бокса-Мюллера на C #

Из этот вопрос: Генератор случайных чисел, который притягивает числа к любому заданному числу в диапазоне? Я провел небольшое исследование, так как раньше сталкивался с таким генератором случайных чисел. Все, что я помню, было имя «Мюллер», так что я думаю, что нашел его здесь:

• Преобразование Бокса-Мюллера

Я могу найти множество его реализаций на других языках, но не могу правильно реализовать его на C #.

Эта страница, например, The Box-Muller Метод генерации гауссовских случайных чисел говорит, что код должен выглядеть следующим образом (это не C #):

#include 
#include 
#include 
#include 

double gaussian(void)
{
   static double v, fac;
   static int phase = 0;
   double S, Z, U1, U2, u;

   if (phase)
      Z = v * fac;
   else
   {
      do
      {
         U1 = (double)rand() / RAND_MAX;
         U2 = (double)rand() / RAND_MAX;

         u = 2. * U1 - 1.;
         v = 2. * U2 - 1.;
         S = u * u + v * v;
      } while (S >= 1);

      fac = sqrt (-2. * log(S) / S);
      Z = u * fac;
   }

   phase = 1 - phase;

   return Z;
}

А вот и моя реализация вышеизложенного на C #. Обратите внимание, что преобразование дает 2 числа, отсюда и трюк с «фазой» выше. Я просто отбрасываю второе значение и возвращаю первое.

public static double NextGaussianDouble(this Random r)
{
    double u, v, S;

    do
    {
        u = 2.0 * r.NextDouble() - 1.0;
        v = 2.0 * r.NextDouble() - 1.0;
        S = u * u + v * v;
    }
    while (S >= 1.0);

    double fac = Math.Sqrt(-2.0 * Math.Log(S) / S);
    return u * fac;
}

Мой вопрос связан со следующим конкретным сценарием, когда мой код не возвращает значение в диапазоне 0-1, и я не могу понять, как исходный код может либо.

• u = 0,5, v = 0,1
• S становится 0,5 * 0,5 + 0,1 * 0,1 = 0,26
• fac становится ~ 3. 22
• возвращаемое значение, таким образом, ~ 0,5 * 3,22 или ~ 1,6

Это не в пределах 0 .. 1 .

Что я делаю не так / не понимаю?

Если я изменю свой код так, что вместо умножения fac на u , я умножу на S , я получу значение в диапазоне от 0 до 1, но имеет неправильное распределение (кажется, имеет максимальное распределение около 0,7-0,8, а затем сужается в обоих направлениях).