23.04.2021 04:03
Desireless 27 2200 Источник

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

Класс генераторов "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]]

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

Иллюстрация двумерной проекции "запыленной поверхности на просвет" для прозрачных 7,4,3 мерных шаров.
Иллюстрация двумерной проекции "запыленной поверхности на просвет" для прозрачных 7,4,3 мерных шаров.

Что еще раз убедило меня в том, что представить многомерные пространства "на пальцах" не получится.

Численное моделирование динамики элементов допускает оптимизации.

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

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

Вторая оптимизация связана с тем что взаимодействие двух элементов никак не влияет на остальные. Так как взаимодействие элемент элемент как источник случайной величины было исключено, то необходимо знать только времена до взаимодействия с границей, что пропорционально количеству элементов 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
Сравнение аналитики и получаемых значений в математическом пакете. По осям графика вероятность наложения в процентах и коэффициент заполнения в процентах (100*T/V). Параметры расчета M=3, V=8000
Результаты сравнения с C#
Результаты сравнения с C#

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

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

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

В алгоритме меж элементное взаимодействие предельно упрощено, можно даже сказать исключено, оставлено только его подобие в виде выбора для обмена импульсом соседей по номеру у взаимодействующего с границей элемента:

//file CRMainSect.cs, function  void Next(), line 66

//register event
EmtsRegEvt();
//Select two elements
int Ei = Rs.E < En - 1 ? Rs.E + 1 : 0;
int Ej = 0 < Rs.E ? Rs.E - 1 : En - 1;
//Interact selected elements
EmtsColSel(Ei, Ej);

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

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

Выводы

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

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

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