В предыдущих частях работы (Часть 1, Часть 2) мы разобрали, что такое линейный конгруэнтный метод (ЛКМ), и как на его основе работает функция Rnd(), вшитая в скриптовый язык VBA, созданный Microsoft. Именно эта функция и "ответственна" за генерацию псевдослучайных чисел. Мы поняли, как ведет себя функция, если в качестве ее аргумента ввести число меньшее либо равное 0. Также мы выяснили, что эта функция работает с мнимым и реальными значениями своих аргументов, также мы поняли, как соотносятся некоторые мнимые значения аргумента функции с их реальными значениями.

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

Загадка числа 13421962

Прошлую часть работы я закончил на том, что при вводе в качестве аргумента числа (-0.1) функция Rnd() выдаст значение rnd, соответствующее начальному значению (x₀), равному 13421962: 

Мнимое значение x₀ – это то значение, которое вводится в качестве аргумента функции (в данном случае -0.1), реальное значение x₀ – это то значение, с которым работает функция для расчета остатка (x₁), и, соответственно, для расчета значения rnd. В данном случае реальное значение x₀ будет равно 13421962 (CCCD8A₁₆), а вычисленный остаток (x₁) - 5785501 (о том, как находится остаток рассказано в первой части работы). Окончательное значение rnd, которое видит пользователь, в языке VBA рассчитывается как:

Итак, стало ясно: передо мной появился новый вызов, который поставил передо мной задачу: почему -0.1, введенное в качестве аргумента функции, дает реальное x₀, равное 13421962?   

Для начала я решил попробовать «скормить» функции различные значения аргументов. Например, -0.9:

На тот момент это мало чем могло помочь… Но тогда же у меня и появилась некоторая догадка: что если попробовать ввести в качестве аргумента отрицательные значения степени 2?

Как оказалось, такое решение имело смысл, и при этом, весьма и весьма немалый:

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

Дробные значения, равные четной степени двойки, имеют шесть 16-ричных разрядов, старшие четыре из которых равны 8000:

И тут у меня появилось понимание, почему мы получаем реальное x₀, равное 8000BF₁₆, если мы берем в качестве аргумента функции число -1. Все встало на свои места – ведь -1 можно представить как -(2⁰). А это четная степень двойки! Значит, представленную выше таблицу можно изобразить и так:

Таким образом, число -1 становится не точкой отсчета, а всего лишь одним из многих чисел в ряду мнимых значений x₀. Но где находится начало этого ряда? Как уже было выяснено во второй части, число 0 стоит особняком и его вполне спокойно можно вынести за скобки данного вопроса. Где заканчивается этот ряд – мы уже знаем. Это самое большое число (по модулю), которое можно представить по правилам представления 32-разрядных чисел стандарта IEEE-754.

Начало начал

В физике существуют понятия классической и квантовой механики. Физические законы, применимые в классической механике, могут не работать в квантовой, и наоборот – законы, применимые в квантовой механике, не совсем релевантны классической механике.

Работая с числами, мы также можем сталкиваться с подобным явлением. Пока речь идет об относительно крупных величинах, ЛКМ работает вполне предсказуемо, но как только речь заходит о малых величинах, этот алгоритм может выдавать, на первый взгляд, неожиданные значения.

Попробуем спуститься вниз и продолжить «вниз» таблицу для отрицательных степеней двойки, берущихся в качестве аргумента функции Rnd(). Получим, что до 2^-126 функция будет вести себя стандартно:

Далее мы могли бы спрогнозировать, что степень -127 даст реальное x₀ равное 128 (80₁₆). Но нет – степень -127 даст значение x0 равное 400080₁₆. Степень -128 – значение 200080₁₆. В итоге получим следующее:

Итак, мы видим, что минимальное число,отличное от нуля, которому VBA Excel присваивает реальное x₀ – это -2^-149. Это и есть тот самый квант, та самая альфа, начало начал, которое мы и искали! Значение 2^-150 VBA уже считает достаточно близким к нулю, и соответственно, дает для него реальное x₀ такое же, какое программа дает и для нуля:

 Может появиться вопрос: как можно вводить в качестве аргумента столь малые значения? Очень просто: при помощи символа степени (^). То есть, чтобы ввести в качестве аргумента функции число 2^-149, нам надо набрать следующее:

round = Rnd(-(2 ^ -149))

В итоге -2^-149 даст значение rnd, соответствующее реальному значению x₀, равному 129. Я думаю, что многие уже догадались, но все-таки для формальности нужно сказать, что -2^-149 – это самое малое отрицательное число (по модулю), которое можно получить, используя 32-битное представление чисел по стандарту IEEE-754 (80000001₁₆):

Более подробно формат представления IEEE-754 был разобран во второй части. Но, напомню, что розовый блок – это блок знака, в котором содержится число 1, что означает, что мы имеем дело с отрицательным числом. Лиловый блок – блок экспоненты. В данном случае значение в лиловом блоке равно 0, но мы должны вычесть из него 127, и в итоге получим значение -127. Соответственно, самый младший разряд зеленого блока мантиссы соответствует степени 2^-149, самый старший – 2^-126. 

Небольшой оффтоп: интересно, что число (80000000₁₆) в формате IEEE-754 будет соответствовать т. н. отрицательному нулю. Таким образом, в указанном стандарте ноль может быть представлен в двух видах: в отрицательном и положительном:

Не совсем понятно, почему разработчики решили назначить для начальной точки именно число 129, а не 0, или, хотя бы единицу. Но что есть, то есть – минимальное число, воспринимаемое VBA, соответствует реальному значению x₀ 129. Таков выбор разработчиков и с этим приходится жить.

Квантовые вычисления

Итак, согласно приведенной выше таблице, мы видим, что число -2^-147 соответствует реальному x₀, равному 132 (84₁₆). А число -2^-148 – реальному числу 130 (82₁₆). А есть ли такое число, которое соответствует реальному 131 (81₁₆)?

Для наглядности изобразим все эти числа на числовой оси:

Будем считать, что число 2^-149 является одним шагом, который для большего пафоса можно назвать квантом. Понятно, что расстояние от 0 до -2^-149 равняется одному кванту. Далее:

2^-148 = 2∙2^-149 = 2^-149 + 2^-149

Следовательно, расстояние от -2^-149 до -2^-148 также равняется одному кванту. Далее: 

2^-147 = 2∙2^-148 = 4∙2^-149  = 2^-149 + 2^-149 + 2^-149 + 2^-149 = 2^-148 + 2^-149 + 2^-149

Отсюда делаем вывод, что расстояние от -2^-148 до -2^-147 будет равняться двум квантам.  

Отсюда напрашивается довольной простой вывод – реальному числу 131 будет соответствовать число, находящееся посередине между -2^-147 и -2^-148. Исходя из всех вышеприведенных данных, понятно, что между двумя этими числами находится число -2^-148 – 2^-149. Вот его-то мы и введем в качестве аргумента функции Rnd():

Rnd(-(2 ^ -148 + 2 ^ -149))

Наша гипотеза полностью подтвердилась! В итоге мы можем представить числовой интервал от -2^-147 до 0 в следующем виде:

Какой вывод из следует всего этого? Мы знаем что числу -2^-126 соответствует реальное x₀, равное 8388736 (800080₁₆). А числу -2^-149 – реальное число 129 (81₁₆). Соответственно, можно предположить, что интервал от -2^-126 до -2^-149 разделен на 8388607 кванта (8388736 - 129).

Проверим это утверждение. Число -2^-126 + 2^-149, согласно этой логике, должно соответствовать реальному числу, на единицу меньшему, чем реальное число, соответствующее аргументу 2^-126 - это число 8388735 (80007F₁₆). Итак, вводим следующий аргумент:

Rnd(-(2 ^ -126 - 2 ^ -149))

И получаем:

Все сходится! А что относительно числа -2^-126 - 2^-149? Проверяем:

Следовательно, можно сделать вывод, что квант, равный 2^-149, будет присутствовать и на интервале от 2^-126 до 2^-125. 

Следим за руками

А теперь я всех призываю внимательно следить за манипуляциями с числами! Мы знаем, что числу -2^-125 соответствует реальное число 129 (81₁₆). Соответственно, мы могли бы спрогнозировать, что числу -2^-125 + 2^-149 будет соответствовать число 128 (80₁₆). Проверяем и получаем следующую картину:

Упс! Число 128 куда-то испарилось… Может быть, мы пропустили какую-то «ступень», которой соответствует то самое неуловимое число 128? Возможно, чтобы найти эту ступень, мы должны прибавить к -2^-125 не целый квант, а только его половину (2^-150)? Но, как мы уже знаем, число 2^-150 в VBA соответствует нулю, поэтому если мы введем в качестве аргумента число -2^-125 + 2^-150, то мы получим картину, соответствующую аргументу 2^-125 без каких-либо дополнительных слагаемых:

Далее все идет без изменений и число -2^-125 + 2^-148 дает реальное x₀, равное 126 (7E₁₆).

Число 128 пока что является для нас весьма трудноуловимым – своеобразным цифровым нейтрино, скрывающимся в потоке других числовых частиц, которые генерирует разобранный нами уже, фактически, на атомы алгоритм ЛКМ.

Впрочем, в отсутствии этого числа есть и свои причины. Как мы уже стало понятно, на интервале от -2^-149 до -2^-126 имеется 8388607 кванта. А сколько квантов имеется по другую сторону от -2^-126 – на отрезке от 2^-126 до -2^-125 + 2^-149?

Во второй части исследования мы уже выяснили, что алгоритм даст одинаковые значения для чисел, сравнимых по модулю с числом 16277216 (1000000₁₆). Итак:

7F₁₆ ≡ 100007F₁₆ (mod 1000000₁₆)  

То есть, для алгоритма нет никакой разницы между числами 7F₁₆ и 100007F₁₆, и если взять два этих числа в качестве реального x₀, на выходе функции мы получим одно и то же число (0,407612).

100007F₁₆  = 16777343₁₀, 800080₁₆ = 8388736₁₀.

Таким образом, на исследуемом нами интервале мы будем иметь 8388607 квантов (16777343 - 8388736). Ровно столько же квантов мы имеем и на противоположном интервале. Итак, мы получили симметричную картину – по обе стороны от числа -2^-126 находится ровно по 8388607 квантов. А кванты с номером 8388608 являются своеобразными переходами, если так можно можно выразиться, на новые энергетические уровни.

Выходим на новый уровень

Итак, мы разобрались с тем, что происходит справа от числа 2^-125. Теперь нужно понять, что происходит слева от этого числа. 

Для этого протестируем некоторые значения мнимых аргументов функции Rnd() и получим следующую картину:

В итоге делаем вывод, что в интервале от -2^-125 до -2^-124 мы будем иметь дело не с найденным квантом, а с шагом, равным двум квантам (2^-148).

Для удобства условимся, что квантом мы будем называть самый малый возможный шаг - в VBA это 2^-149. А шагом мы будем называть любой шаг, размер которого будет больше кванта.

Соответственно, на участке от -2^-124 до 2^-123 шаг снова увеличится в два раза и уже будет составлять уже 2^-147 (4 кванта).

В отличие от предыдущего интервала картина уже не будет симметричной – в правой части мы получим 8388608 шагов (8388737 – 129). А в левой – 8388607 шагов (16777344 - 8388737). Что интересно, мы смогли поймать «нейтрино» предыдущего интервала – числу 128 (80₁₆) соответствует точка, являющая правой границей для шага перехода к следующему числовому интервалу. Мнимое значение x0 в этом месте равно -2^-123 + 2^-147. В то же время в левой части интервала исчезло число 129 (81₁₆).

Та же самая закономерность будет наблюдаться и далее для любой нечетной степени двойки, большей или равной -125. Итак, если n – нечетное число, n ≥ -125, q – размер шага, то мы получим следующее:

Вычислить, какое реальное значение x₀ соответствует каждой степени -2 не составит особого труда. Таблицы соответствия для некоторых значений степеней также были приведены во второй части данной работы.

Возвращаясь к -0.1

Теперь у нас имеется вся теоретическая база для того, чтобы с успехом решить заявленную в самом начале этого текста загадку числа 13421962. Итак, -0,1 лежит между числами-1/8 (-(2^-3)) и -1/16 (–(2^-4)). Соответственно, n для этого числа будет равен -5, а показатель степени размера кванта на данном промежутке составит n-22 = -5 -22 = -27.

Мы уже знаем, что аргументу -1/16 соответствует реальное значение 8388797 (8000BD₁₆). Итак, для того, чтобы найти реальное число, соответствующее аргументу -0.1, нам нужно узнать, сколько нужно сделать шагов размером 2^-27 от числа -1/16 до числа -1/10. Затем это число мы прибавляем к значению 8388797 и, таким образом, получаем реальное значение, соответствующее мнимому значению аргумента -0.1.

На самом деле мы уже знаем количество шагов (N) – подсказка у нас уже есть. Мы знаем, что аргументу -0.1 соответствует реальное x₀, равное 13421962. Соответственно, N = 13421962 - 8388797 = 5033165.

Проверяем нашу гипотезу. Можно вычислить N, как говорится в лоб. Расстояние  от  -0.1 до -1/16  (S) равно 3/80.

Теперь мы можем вычислить N:

Понятно, что число 5033164,8 округляется до 5033165, что дает нам N, совпадающее с ранее полученными данными. В принципе, гипотеза является доказанной, но доказанной, если так можно выразиться, постфактум.

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

Итак, в данном случае точность у нас равняется 2^-27. Это означает, что нам необходимо представить число -0.1 в двоичном виде с точностью до 27 знака после запятой. Сделать это можно в любом калькуляторе, который переводит десятичные числа в двоичные числа. Переводим десятичную дробь -0.1 в двоичную:

-0.1₁₀ = -0.0001 1001 1001 1001 1001 1001 1001 1001 …₂ 

Мы видим, что в 27 разряде после запятой в двоичной дроби имеется 0. Но в 28 разряде имеется единица. Соответственно, по правилам двоичного округления при округлении до 27 разрядов последний разряд числа у нас будет равняться 1.

Следовательно:      

-0.1₁₀ ≈ -0.0001 1001 1001 1001 1001 1001 101₂ 

В итоге получим, что:

-0.1₁₀ ≈ -(2^-4) – 2^-5 – 2^-8 – 2^-9 – 2^-12 – 2^-13 – 2^-16 – 2^-17 – 2^-20 – 2^-21 – 2^-24 – 2^-25 – 2^-27

Теперь нам следует убрать из этой суммы первое слагаемое и выяснить, сколько в получившемся значении содержится шагов размеров 2^-27. Надеюсь понятно, почему нам не следует трогать число -(2^-4) – оно является своеобразной вехой, от которого идет отсчет. Мы знаем, что числу -(2^-4) соответствует реальное реальное x₀, равное 8388797 и для того, чтобы получить реальное x₀ для -0.1 нам нужно будет прибавить к числу 8388797 искомое число N.

Итак, 2^-5  = 2²² ∙ 2^-27, 2^-8  = 2¹⁹ ∙ 2^-27 и т. д. В итоге имеем:

-0.1₁₀ ≈ -(2^-4) – 2^-27 ∙ (2²² + 2¹⁹ + 2¹⁸ + 2¹⁵ + 2¹⁴ + 2¹¹ + 2¹⁰ + 2⁷ + 2⁶ + 2³ + 2² + 1) ≈ -(2^-4) – 2^-27 ∙ ∙ (4194304 + 524288 + 262144 + 32768 +16384+ 2048 + 1024 + 128 + 64 + 8 + 4 + 1) ≈ ≈ -(2^-4) – 2^-27 ∙ 5033165

Вау!!! Хотя, такой результат и был ожидаем, но все же, на мой взгляд, он достоин восхищения. Мы вместе прошли все этапы очень долгого пути, чтобы, наконец, с полным правом сказать: мы знаем о функции Rnd() почти все…    

Небольшая формальность: введем в качестве аргумента функции Rnd() следующие числа:

-0.1
a = -(2 ^ -4 + 2 ^ -5 + 2 ^ -8 + 2 ^ -9 + 2 ^ -12 + 2 ^ -13 + 2 ^ -16 + 2 ^ -17 + 2 ^ -20 + 2 ^ -21 + 2 ^ -24 + 2 ^ -25 + 2 ^ -27)

Получаем:

Итак, тайна числа 13421962 полностью раскрыта! Для любителей поломать себе мозги предлагаю по этому же алгоритму «разложить» любую другую десятичную дробь и сравнить полученные результаты с теми, которые выдаст алгоритм. Далеко ходить не надо, можно взять уже упомянутое число -0.9:

Еще не конец

Внимательные читатели должны были заметить, что выше я сказал о том, что теперь мы знаем о функции Rnd() почти всё. Но, наверное, если бы все было так просто, как кажется на первый взгляд, вряд ли мы могли говорить о том, что имеем дело с продуктом корпорации из списка Big Five. А закладывать "сюрпризы" в свои детища они весьма и весьма любят - вот и здесь, как оказалось, не обошлось без такого сюрприза. Итак, поехали!

Понятно, что размер шага с каждым новым интервалом увеличивается в два раза, пока в один прекрасный момент не станет равным двум. Это произойдет на интервале от -2²⁴ (-16777216) до -2²⁵ (-33554432).

Итак, по логике, если числу -16777216 соответствует реальное x₀, равное 8388811 (8000CB₁₆), то мнимому аргументу -16777218 должно соответствовать реальное число 8388812 (8000CC₁₆), поскольку шаг на этом интервале у нас стал равным 2. Проверяем:  

  Да, все сходится. А что насчет числа -16777217? Выше уже говорилось, что любое число на интервале, меньшее шага этого интервала, VBA считает нулем. Соответственно, и число 1 на указанном интервале должно уже считаться нулем. Итак, -16777216 – 1 должно давать реальное x₀, соответствующее числу -16777216, что мы и видим:

   Казалось бы – вот и всё! Какие еще могут быть сюрпризы? По этой же логике число -16777219 должно давать реальное x₀ 8388812:  

И здесь нас ждет очередной «упс». Почему-то число -16777219 выдает реальное x₀ 8388813, словно у нас неожиданно уменьшился шаг… Что же, надо идти дальше, в итоге получаем следующую картину:

Оказывается, после первых двух чисел интервала, которые будут соответствовать одному и тому же реальному x₀, мы будем иметь дело с чередованием серий чисел длиной 1 и 3 соответственно. При этом каждое число серии будет соответствовать одному и тому же реальному x₀.

А как обстоит дело в конце интервала?

Как видно, последнее целое число интервала соответствует тому же реальному x₀, которому соответствует и первое число следующего интервала (сделаем допущение, что число -2²⁵ в интервал не входит). Далее все стандартно – чередуются серии из одного и трех целых чисел. Таким образом, мы можем записать следующую формулу для интервала от -2²⁴ до -2²⁵:

-2²⁵ = -2²⁴ – 2 – 1 – 3 – 1 – 3 – 1 …– 3 – 1 – 3 – 1 – 1

  Не вдаваясь в подробности, запишем формулы для нескольких следующих промежутков (получены опытным путем):

-2²⁶ = -2²⁵ – 3 – 3 – 5 – 3 – 5 – 3 … – 5 – 3 – 5 – 3 – 2   

-2²⁷ = -2²⁶ – 5 – 7 – 9 – 7 – 9 – 7… – 9 – 7 – 9 – 7 – 4

-2²⁸ = -2²⁷ – 9 – 15 – 17 – 15 – 17 – 15 … – 17 – 15 – 17 – 15 – 8

Понятно, что в каждой из формул каждое слагаемое (кроме первого) указывает на серию из целых чисел, дающих одинаковые x₀.

Как можно получить формулу для следующего интервала (от -2²⁸ до -2²⁹)? Думаю, внимательные читатели уже могли заметить закономерность.

Понятно, что первым числом формулы будет -2²⁸. Второе число мы получим, взяв четвертое число из формулы для предыдущего интервала – оно же будет равно третьему с краю числу из этой же формулы. Вот такая рекурсия. В данном случае – это число -17. Теперь нам нужно получить числа для повторяющихся серий. В данном случае шаг q у нас будет равен 32. Соответственно, сумма чисел в сериях, как видно из предыдущих формул, должна будет равняться -2q. Обозначим эти числа как a и b:

Из предыдущих формул мы видим, что:

Таким образом, в данном случае a будет равно -31, b – -33. Последнее число суммы будет равно удвоенному последнему числу формулы для предыдущего интервала.    

Итак, теперь у нас есть формула для еще одного интервала:

-2²⁹ = -2²⁸ – 17 – 31 – 33 – 31 – 33 – 31 … – 33 – 31 – 33 – 31 – 16

По такой же схеме мы сможем вычислить формулы для любого интервала, как теоретического, так и практически доступного для работы функции Rnd() языка VBA.

На этом месте я позволю себе закончить третью часть данного материала. В четвертой части мы попытаемся понять, почему числа на интервалах с шагом большим, чем единица, ведут себя таким странным образом. Не пропустите новый выпуск! До новых встреч на Хабре!