Надлежащий способ генерировать случайное число с плавающей запятой, заданный двоичным генератором случайных чисел?
Допустим, у нас есть генератор двоичных случайных чисел, int r(); это вернет ноль или единицу с вероятностью 0,5.
Я посмотрел на Boost.Random, и они генерируют, скажем, 32 бита и делают что-то вроде этого (псевдокод):
x = double(rand_int32());
return min + x / (2^32) * (max - min);
У меня есть некоторые серьезные сомнения по этому поводу. Двойник имеет 53 бита мантиссы, и 32 бита не могут должным образом генерировать полностью случайную мантиссу, среди прочего, например, ошибки округления и т. Д.
Что бы быстро создать равномерно распределенный float или же double в полуоткрытом диапазоне [min, max), предполагая IEEE754? Акцент здесь лежит на правильность распределения, не скорость.
Чтобы правильно определить правильное, правильное распределение будет равно тому, которое мы получим, если мы возьмем бесконечно точный равномерно распределенный генератор случайных чисел, и для каждого числа мы округлим до ближайшего представления IEEE754, если это представление все еще будет в пределах [min, max)иначе номер не будет учитываться при распределении.
П.С .: Мне было бы интересно и правильные решения для открытых диапазонов.
Решение
Вот правильный подход без попытки эффективности.
Мы начнем с класса bignum, а затем — рациональной обертки указанных bignum.
Мы производим ассортимент «достаточно большой, чем» наш [min, max) диапазон, так что округление нашего smaller_min а также bigger_max производит значения с плавающей запятой вне этого диапазона, в нашем рациональном построении на bignum.
Теперь мы поделили диапазон на две части идеально посередине (что мы можем сделать, так как у нас есть рациональная система Bignum). Мы выбираем одну из двух частей наугад.
Если после округления верх и низ выбранного диапазона будут (A) вне [min, max) (с той же стороны, имейте в виду!) вы отклоняете и перезапускаете с самого начала.
Если (B) верх и низ вашего диапазона округляется до одинакового double (или же float если вы возвращаете число с плавающей точкой), все готово, и вы возвращаете это значение.
В противном случае (C) вы рекурсируете на этом новом, меньшем диапазоне (подразделите, выбирайте случайным образом, тестируйте).
Нет никаких гарантий, что эта процедура останавливается, потому что вы можете либо постоянно углубляться до «края» между двумя округлениями doubles, или вы можете постоянно выбирать значения за пределами [min, max) спектр. Вероятность этого (однако, никогда не останавливается) равна нулю (при условии хорошего генератора случайных чисел и [min, max) ненулевого размера).
Это также работает для (min, max)или даже выбирая число в округленном достаточно толстом наборе Кантора. Пока мера действительного диапазона реалов, округляющих до правильных значений с плавающей запятой, не равна нулю, и диапазон имеет компактную поддержку, эта процедура может быть запущена и имеет вероятность 100% завершения, но без жесткого верхнего ограничено по времени, которое может быть сделано.
Другие решения
Проблема здесь в том, что в IEEE754 двойники, которые могут быть представлены, не являются равнораспределенными. То есть, если у нас есть генератор, генерирующий действительные числа, скажем, в (0,1), и затем сопоставленный с представимыми числами IEEE754, результат не будет равномерно распределен.
Таким образом, мы должны определить «равное распределение». Тем не менее, предполагая, что каждое число IEEE754 является просто представителем для вероятности нахождения в интервале, определенном округлением IEEE754, процедура первого генерирования равнораспределенных «чисел» и округления до IEEE754 будет генерировать (по определению) an » Экви-распределение «IEEE754 номеров.
Следовательно, я полагаю, что приведенная выше формула станет произвольно близкой к такому распределению, если мы просто выберем достаточно высокую точность. Если мы ограничим задачу нахождением числа в [0,1), это означает ограничение на множество деномализованных чисел IEEE 754, которые представляют собой однозначное число до 53-разрядного целого числа. Таким образом, должно быть быстро и правильно генерировать только мантиссу с помощью 53-разрядного двоичного генератора случайных чисел.
Арифметика IEEE 754 всегда является «арифметикой с бесконечной точностью, за которой следует округление», то есть число IEEE754, представляющееб является ближайшим кb (иначе говоря, вы можете думать о a * b, рассчитанном с бесконечной точностью, затем округленном до числа IEEE754 с закрытием). Следовательно, я считаю, что min + (max-min) * x, где x — деномализованное число, — это осуществимый подход.
(Примечание: как ясно из моего комментария, я сначала не знал, что вы указываете на случай, когда min и max отличаются от 0,1. Денормализованные числа обладают свойством, что они равномерно распределены. Следовательно, вы получаете распределение по уравнению сопоставление 53 битов с мантиссой. Далее вы можете использовать арифметику с плавающей запятой, потому что она верна вплоть до точности машины. Если вы используете обратное отображение, вы восстановите равномерное распределение.
Посмотрите этот вопрос для другого аспекта этой проблемы: Масштабирование Int равномерного случайного диапазона в Double
std::uniform_real_distribution.
Есть действительно хороший разговор С.Т.Л. с конференции Going Native этого года, которая объясняет, почему вы должны использовать стандартные дистрибутивы, когда это возможно. Короче говоря, свернутый вручную код имеет тенденцию быть смехотворно низкого качества (думаю, std::rand() % 100), или имеют более тонкие недостатки однородности, такие как в (std::rand() * 1.0 / RAND_MAX) * 99, который является примером, приведенным в докладе, и является частным случаем кода, размещенного в вопросе.
РЕДАКТИРОВАТЬ: я взглянул на реализацию libstdc ++ std::uniform_real_distributionи вот что я нашел:
Реализация выдает число в диапазоне [dist_min, dist_max) с помощью простого линейного преобразования из некоторого числа, полученного в диапазоне [0, 1), Он генерирует этот номер источника, используя std::generate_canonical, реализация которого может быть найдена здесь (в конце файла). std::generate_canonical определяет количество раз (обозначается как k) диапазон распределения, выражается как целое число и обозначается здесь как r*, поместится в мантиссе целевого типа. То, что он тогда делает, по сути, генерирует одно число в [0, r) для каждого rсегмент мантиссы и, используя арифметику, заполняет каждый сегмент соответственно. Формула для полученного значения может быть выражена как
Σ(i=0, k-1, X/(r^i))
где X является стохастической переменной в [0, r), Каждое деление на диапазон эквивалентно сдвигу на количество битов, используемых для его представления (т.е. log2(r)) и так заполняет соответствующий сегмент мантиссы. Таким образом, используется вся точность целевого типа, и поскольку диапазон результата [0, 1)экспонента остается 0** (по модулю смещения), и у вас не возникает проблем с единообразием, которые возникают у вас, когда вы начинаете возиться с показателем степени.
Я бы не стал доверять простоте, что этот метод криптографически безопасен (и у меня есть подозрения о возможных ошибках при расчете размера r), но я полагаю, что это значительно более надежно с точки зрения единообразия, чем реализация Boost, которую вы опубликовали, и определенно лучше возиться с std::rand,
Возможно, стоит отметить, что код Boost на самом деле является вырожденным случаем этого алгоритма, где k = 1это означает, что это эквивалентно если входной диапазон требует как минимум 23 бита для представления его размера (стандарт IEE 754 с одинарной точностью) или как минимум 52 бита (двойная точность). Это означает, что минимальный диапазон составляет ~ 8,4 миллиона или ~ 4,5е15, соответственно. В свете этой информации, я не думаю, что если вы используете бинарный генератор, реализация Boost довольно собираюсь сократить это.
После краткого взгляда на реализация libc ++, похоже, что они используют один и тот же алгоритм, реализованный немного по-другому.
(*) r на самом деле диапазон ввода плюс один. Это позволяет использовать max значение urng в качестве действительного ввода.
(**) Строго говоря, закодированный показатель не 0, так как IEEE 754 кодирует неявное ведение 1 перед основанием значенияand. Концептуально, однако, это не имеет отношения к этому алгоритму.
AFAIK, правильный (и, вероятно, также самый быстрый) способ состоит в том, чтобы сначала создать 64-разрядное целое число без знака, где 52 дробных бита являются случайными битами, а показатель степени равен 1023, что, если тип, вставленный в двойное число (IEEE 754), будет равномерно распределено случайное значение в диапазоне [1,0, 2,0). Таким образом, последний шаг состоит в том, чтобы вычесть 1.0 из этого, что приводит к равномерно распределенному случайному двойному значению в диапазоне [0.0, 1.0).
В псевдокоде:
rndDouble = bitCastUInt64ToDouble (1023 << 52 | rndUInt64 & 0xfffffffffffff) — 1,0
Этот метод упоминается здесь:
http://xoroshiro.di.unimi.it
(См. «Создание одинаковых двойников в единичном интервале»)
РЕДАКТИРОВАТЬ: рекомендуемый метод с тех пор изменился на:
(x >> 11) * (1. / (UINT64_C (1) << 53))
См. Ссылку выше для деталей.