Одномерный генератор случайных действительных чисел

Генерация непрерывной равномерной случайной величины как оказалась задача не из простых, может потому что сама постановка задачи слегка абсурдна ведь логически получить случайность это значит не найти решения. Однако оставив простейшее "как будет на английском - я не знаю?", в статьях посвященных алгоритмам вы обнаружите работы по созданию только последовательностей случайных чисел.

Класс генераторов "True Random" использует физические явления и внешние связи, так например для генерации десятичного случайного числа вы можете найти рекомендацию использования "атмосферного датчика". Естественно как любителю программирования такое положение дел мне показалось несправедливым, и довольно долгое время "задача созревала". Вариант решения, как и следовало ожидать из постановки задачи, появился случайно, как дополнение к задаче компрессии для упаковки твердых сфер. Задача не нашла аналитического решения как и нет пока доказательств его отсутствия, соответственно источник по внешним признакам вполне подходящий.

Идея не нова, использовалась например в UNIX системах, но с причиной по которой я не могу приведенный алгоритм использовать как функцию столкнулся только при тщательном изучении его работы. Математически обеспечить конечную флуктуацию бесконечного числа параметров невозможно, потому если генератор действительно непрерывный, то в отличие от арифметического количество его значений бесконечно. На практике я не сталкивался со сбоями в его работе, но в сумме это не более чем месяцы работы, хоть "второй закон термодинамики" тоже на моей стороне, но строгой логической надежности он не обеспечивает. Потому на использование алгоритма в качестве функции для систем с высокой надежностью я не претендую, но допускаю что вместе с дополнительными доработками формальную надежность можно существенно повысить.

Источник непрерывной равномерно распределенной случайной величины это взаимодействие элемента с границей в модели твердых сфер. Внутренние взаимодействие элементов имеет отличное от равномерного распределение, это было проверено на тестах Diehard и потому было исключено из алгоритма, но оставлено в упрощенном виде для дополнительного "перемешивания".

Равномерное распределение на поверхности сферы это специальное распределение, и для перехода к плоскому варианту требуется дополнительное преобразование. Это преобразование обнаружил экспериментально, помощи математиков для оформления гипотезы в теорему я не получил, ответы варьировались от "это же очевидно" до "я не занимаюсь этой областью". Звучит эта гипотеза следующим образом:

Равномерное случайное распределение на поверхности D - мерной сферы без двух произвольно выбранных координат равнозначно равномерному случайному распределению в объеме D-2 мерного шара.

Двумерная проекция случайных точек с поверхности четырех мерной сферы.

Биения возле равномерного распределения связаны с качеством перехода от одномерного генератора математического пакета к многомерному. Применялся достаточно известный способ проецирования на поверхность сферы объемного распределения в кубе со случайной перестановкой координат.

RandomSphere[Rn_: 2, Pn_: 1, Rb_: 1] :=  Module[{i, j, m, p, r, s, S, X, Xi, Xj, Pm},  X = Array[0 &, {Pn, Rn}]; Pm = Rn; s = 1/Sqrt[2];  For[p = 1, p <= Pn, p++, i = RandomInteger[{1, Rn}]; S = 0;    For[r = 1, r <= Rn, r++,    X[[p, r]] =      If[r != i, RandomReal[{-1, 1}], RandomChoice[{-1, 1}]];     S += X[[p, r]]^2];   X[[p]] *= Rb/Sqrt[S];   For[m = 1, m <= Pm, m++,    i = RandomInteger[{1, Rn}]; j = i;     While[i == j, j = RandomInteger[{1, Rn}]];    Xi = X[[p, i]];    Xj = X[[p, j]];    X[[p, i]] = s (Xj - Xi);    X[[p, j]] = s (Xj + Xi)]]; Return[X]]

Использование потока из алгоритма этой публикации я посчитал недостаточно доказательным, хоть гипотеза была обнаружена именно на нем. Меня очень удивило что запыленный на поверхности четырех мерный шарик на просвет оказался плоским, в прочем пяти мерный вместо затемнения у краев как у трехмерного показывает концентрацию в центре. Это еще раз убедило меня в том что представить многомерные пространства "на пальцах" не получится.

Практически гипотезу можно использовать для получения равномерного распределения на отрезке, ведь он является одномерной проекцией распределения на поверхности трех мерной сферы. Использование такого простого перехода позволяет так же избавится от "проклятия границ" ведь дополнительные преобразования над величинами отсутствуют а потому брать больший объем сферы не нужно.

Численное моделирование динамики элементов допускает оптимизации:
Первая оптимизация связана с возможностью исключить расчет на промежутках отсутствия взаимодействия элементов. И даже более получать точное значение времени взаимодействия, что фактически означало бы для случая шага по времени его значение равным пределу машинной точности. Так как для алгоритма существенны только взаимодействия элементов в границей, расчет времени приведен только для этого случая:

Вторая оптимизация связана с тем что взаимодействие двух элементов никак не влияет на остальные. Так как взаимодействие элемент элемент как источник случайной величины было исключено, то необходимо знать только времена до взаимодействия с границей, что пропорционально количеству элементов N. Это означает существенную экономию в вычислениях так как меж элементное взаимодействие в простейшем случае означает сложность пропорциональную квадрату количества элементов O(N(N-N1)/2) или O(N^2). В более сложной реализации с использованием списков для хранения времени возможны и лучшие значения для меж элементного взаимодействия, не выводящие за пределы N^(1+1/D) даже на плотностях близким к плотной упаковке, но эта статья не посвящена алгоритму компрессии.

Проверку качества распределения производил на тестах Diehard, но для этой статьи я выбрал модификацию Parking Test, потому как доказательность "Numerous experiments prove" я лично не смог идентифицировать, так же оригинальный тест не позволял прямо проверять многомерное распределение и предназначен для средних величин заполнения, что не позволяет его использовать для оценки качества непрерывности.

Описание теста: броски монеток на поднос и расчет вероятности того, что монетки наложатся.

Математика теста:

Вероятность попадания сферического элемента единичного объема, размером

$r = {\left( {\frac{{Gamma\left[ {\frac{M}{2} + 1} \right]}}{{{\pi ^{\frac{M}{2}}}}}} \right)^{\frac{1}{M}}}$

где M - мерность пространства, Gamma- Гамма функция. В область линейного размера R-r, где

$R = {\left( {V\frac{{Gamma\left[ {\frac{M}{2} + 1} \right]}}{{{\pi ^{\frac{M}{2}}}}}} \right)^{\frac{1}{M}}}$

а V ее объем, равна единице. Используем равномерное распределение по объему, тогда

$P\left( x \right) = \begin{cases} {{{\left( {\frac{x}{{R - r}}} \right)}^M}} &\text{${0 \le x \le R - r}$}\\ 1&\text{${x > R - r}$} \end{cases}$

вероятность единичного конфликта X, попадания с пересечением двух элементов с размером r, т.е. попадания в область с внешним размером 2r равна:

$X = \frac{{{{\left( {2r} \right)}^M}}}{{{{\left( {R - r} \right)}^M}}} = \frac{{{2^M}}}{{{{\left( {\frac{R}{r} - 1} \right)}^M}}}$

так как (T(T-1))/2-количество расстояний для Т участников, то, вероятность отсутствия конфликта Y соответственно равна:

$Y = {\left( {1 - X} \right)^{\frac{{T\left( {T - 1} \right)}}{2}}} = {\left( {1 - \frac{{{2^M}}}{{{{\left( {\frac{R}{r} - 1} \right)}^M}}}} \right)^{\frac{{T\left( {T - 1} \right)}}{2}}} = {{\rm{e}}^{ - \frac{{{2^M}}}{{{{\left( {\frac{R}{r} - 1} \right)}^M}}}\frac{{T\left( {T - 1} \right)}}{2}}}$

последнее равенство использует предельный переход, а результат в целом не учитывает влияние границ, но все же формула достаточно точна и это можно увидеть на сравнениях с генератором математического пакета. Использование арифметических генераторов для размерностей выше шести в литературе не рекомендуется, потому я не стал разбираться кто более прав в данном случае, и ограничился тем что на размерностях меньше четырех я имею достаточно точный инструмент проверки распределения.

Сравнение аналитики и получаемых значений в математическом пакете. По осям графика вероятность наложения в процентах и коэффициент заполнения в процентах (100*T/V). Параметры расчета M=3, V=8000

Результаты сравнения с C#

Результаты сравнения показывают что на высоких размерностях поведение алгоритма немного но лучше соответствуют аналитике чем библиотечный C# вариант арифметического генератора Кнудта. Для одномерного случая, с той точностью что считал, методика Diehard отличия не проявила.

Многомерный поток - заполнение объема единичной сферы, показал "память" от высокой размерности и для использования оказался не пригоден. Визуально отличить распределение от равномерного невозможно, как и некоторые тесты его тоже распознают таковым, но именно ParkingTest и его модификация результаты выбраковывает. Так же причиной несоответствия может являться обобщение гипотезы на высокие размерности. Внимательно я этот вопрос не исследовал, а ограничился только проверкой действия всего алгоритма стандартными тестами.

Производительность, замеры производительности производил оценочно, к сожалению алгоритм не смог догнать простейшие реализации библиотечных вариантов но все же показал сравнимое по порядку с ними время исполнения. Значительно перегнал вариант из математического пакета, но это сравнение вряд ли заслуживает внимания.

Выводы:

Генератор показал свойства равномерного непрерывного распределения для одномерного случая. В генераторе нет недостатков его арифметических аналогов, так например отсутствует период повтора, теоретически это обеспечивается влиянием иррациональности при отражении элементов, плотность иррациональных чисел выше. Практически же это реализуется благодаря случайной величине на хвостах значений времени до столкновения, получающейся при операциях сложения и вычитания, потому можно сказать что данный генератор получает свойства класса TrueRandom благодаря дополнительной внешней связи с памятью. Генератор не имеет "времени разогрева" хотя нормализация распределения элемент-элемент требует для сотни элементов более чем миллиона значений, события создаваемые на поверхности сферы показывают необходимые свойства сразу, не смотря на то что стартовые значения задаются с помощью арифметического ГСПЧ.

Реализация алгоритма на C#