Генерирование псевдослучайных чисел Метод середины квадрата (стр. 2 из 3)

Построение модели

Случайная величина γ равномерно распределена на интервале [0,1], если определяются соотношениями:

Р_γ(γ)=1, γ∊[0,1]; F_γ(γ)=

Математическое ожидание

и дисперсия (среднеквадратическое отклонение) для непрерывных и дискретных случайных величин:

=

=

=

=

Теперь об упомянутом увеличении длины отрезка апериодичности. Есть вероятность, что числа, генерируемые при помощи алгоритма обнулятся, при условии х₀ < 10⁴. чтобы избежать этого нужно длину отрезка апериодичности L увеличить. Практические расчеты показали, что L (длина отрезка апериодичности) достаточно мала, и это может явиться причиной отказа от алгоритма середины квадрата. Длину отрезка апериодичности можно несколько увеличить, если при выполнении условия

в качестве следующего случайного числа взять

Теперь функция γ_nпримет следующий вид:

Φ(γ_n)=

if ,

Для реализации алгоритма нужно определить наибольшее возможное значение k(чем больше это число, тем больше число значащих разрядов, а, следовательно, и длина отрезка апериодичности).

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

Например, можно принять, что γ_n есть величина типа Single, а

типа Double. При длине слова в 32 бита (1 разряд для знака; 7 бит для характеристики СС и 24 разряда под мантиссу ffffff ) мантисса может принимать 2¹⁰ *2 ¹⁰ *2 ⁴=1024*1024*32 различных значения, что соответствует 7…8 значащим цифрам.

При длине слова в 64 бита под мантиссу отводится 24+32 разряда. При помощи 56 разрядов можно закодировать 2⁵⁶=(1024)⁵*64, что соответствует 15..16 значащим цифрам. Таким образом, при величинах такого типа максимально возможное значение k равно 3 (т. к. 2k<7; 4k<15).

Разработка алгоритма

Приступим к разработке алгоритма.

Алгоритм seredina [середина квадрата]

Шаг0 [инициализация] х:=х₀;

Шаг1 [цикл, вычисление последовательности M случайных чисел х [1..M].]

For j:=1 to M do [окончаниешаг2 и stop]

Шаг2 [генерация нового случайного числа]

Y:=X² [Y имеет 4k разрядов]

Число Х получаем удалением по k разрядов с каждого конца Y

X[j]:=X;

При k=2 и х₀=0.2134 в соответствии с первым оператором Шаг2 получим Y=

=0.04 5539 56. после применения второго оператора х₁=0.5539. если х₀<10⁴, то все числа, генерируемые при помощи данного алгоритма, будут тождественно равны нулю. Здесь нужно осуществлять описанное выше увеличение отрезка апериодичности.

Для реализации алгоритма потребуется написать процедуру, реализующую метод середины квадрата, которую назовем Rand, в которой при первом обращении задается целое число IX из [0.999999] и между вызовами оно не должно меняться. На выходе имеем числа: Rand числа типа Single из (0,1), IX числа типа LongInt из (0,999999), K числа типа LongInt из (0, К+1). Далее зададим переменные. Обнулим 7 и далее разряды числа Х, вычислим квадрат Х. Выберем 4…9 разряды числа Y. Определим новое целое случайное число IX из [0.999999] и случайную величину из [0.1]. Вычислим число К из [0.K+1].

Рисунок 1. блок-схем программы

Нет Да

Рисунок 2. блок-схема процедуры Rand

Кодирование алгоритма

Programseredina {метод середины квадрата}

Usescrt;

Var

Urand: Single;

Dm, Dd: Double;

n, i, j, k, l, ll: LongInt;

F: Text;

Procedure Rand (var Ix, k: LongInt; Urand: Single);

Var

x, z: Single;

y, z1: Double; i: LongInt;

Begin

z:=Ix; x:=z{*1.e-6}; {обнулим 7 и далее разряды числа Х.}

z1:=x+1.e-6; y:=z1*z1*1.e-12; {вычислим квадрат Х.}

z1:=y*1.e+3; y:=z1-Trunc(z1); {выберем 4…9 разряды числа Y.}

if y<1.e-6 then

Begin

z1:=y*1.e+6;

y:=z1-Trunc(z1);

End;

z1:=y*1.e+6;

y:=Trunc(z1);

Ix:=Trunc(y); y:=y*1.e-6; Urand:=y;

{определили новое целое случайное число IX из [0.999999] и случайную величину из [0.1]}

k:=Round((k+1)*y); {вычисляем число К из [0, К+1]}

End; {Rand}

Begin

Clrscr;

Assign (F, ’d:/F/Rand.dat’);

Rewrite(F);

n:=3–23; Dm:=0; Dd:=0; ll:=50000;

for l:=1 to ll do

begin

k:=100;

Rand (n, k, Urand);

Dm:=Dm+Urand/ll;

Dd:=Dd+(Urand-0.5)*(Urand-0.5)/ll;

End;

Writeln (F, ’M=’, Dm, ’Sig^2=’, Dd);

Close(F);

End.

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

Анализ сложности алгоритма

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

Пусть G(n) – алгоритм для решения некоторого класса задач, а n – размерность размерность отдельной задачи из этого класса. Определим

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

Алгоритм G(n) полиномиальный, если

(n) растет не быстрее, чем полином от n, в противном случае алгоритм экспоненциальный.

Функцию f(n) определяют как О

и говорят, что она порядка g(n) для больших n, если

Функция h(n) является О

для больших n, если

Если f(n) есть О

, то эти две функции возрастают с одинаковой скоростью при n→∞. Если f(n) есть О, то g(n) возрастает гораздо быстрее, чем f(n).

Алгоритм G(n) называется полиномиальным, если

(n)= О, в противном случае алгоритм является экспоненциальным.

Рассмотрим процедуру Rand. В ней всего 4 шага, одно сравнение. Для процедуры G(n) рабочая функция имеет вид

(n)= О(n) и она является полиномиальной.

Проверка правильности программы

Анализ структуры программы

Для этого воспользуемся управляющим графом (ориентированный граф, вершины которого – операторы, а дуги соединяют операторы, между которыми возможна передача управления).

6 7 8 9 10 11 12 14

22 21 20 19 18 17 15