Форум программистов, компьютерный форум, киберфорум
Наши страницы
palva
Войти
Регистрация
Восстановить пароль
Оценить эту запись

Генерация случайных чисел в фортране

Запись от palva размещена 24.12.2013 в 03:11
Обновил(-а) palva 09.01.2014 в 16:36

Для генерации случайных чисел в фортране имеется встроенная подпрограмма RANDU. Эта подпрограмма предполагает, что состояние генератора хранится в двух пользовательских переменных целого типа. Их сохранение от вызова к вызову генератора, а также инициализация, а если требуется рандомизация, должны обеспечиваться пользователем. Это создает дополнительные хлопоты по сравнению с подобными программами для других языков, где состояние генератора хранится где-то во внутренней памяти. Исторически также сложилось, что для хранения требуется две переменных, а не одна. Это добавляет неудобств. Естественно, что программисты стали использовать генераторы, реализованные в разных численных библиотеках или даже писали генераторы сами. Позднее в современных версиях фортрана появились более удобные средства (RANDOM_NUMBER, RANDOM_SEED). Но мы сначала опробуем работу RANDU. Этот генератор возвращает псевдослучайное значение типа REAL*4, равномерно распределенное между нулем и единицей. (В документации сказано, что единица случиться не может, но не стоит на это надеяться.) Приведем пример:
Fortran
1
2
3
4
5
6
7
integer*2 s1,s2
real*4 r
print *,s1,s2
call randu(s1,s2,r)
print *,s1,s2
print *,r
end
Мы объявили две переменных s1 и s2 для хранения состояния и дальше обратились к RANDU, передавая эти переменные в качестве параметров. RANDU выработал нам псевдослучайное число r, которое мы напечатали. Чисто из любопытства мы напечатали значения переменных s1 и s2 до и после обращения к RANDU. Обратите внимание, что переменным s1 и s2 мы не присвоили начальные значения, и там оказался мусор. Запуская программу несколько раз и даже на разных компьютерах мы видим, что мусор там оказывается один и тот же. Поэтому и число, генерируемое программой RANDU будет одним и тем же. После работы RANDU, как мы видим, состояние генератора изменилось. Поэтому дальнейший вызов RANDU, если бы он был, дал бы нам другое случайное число. Как правило многократный вызов RANDU осуществляется с одними и теми же переменными состояния.

Теперь рассмотрим вопрос о начальном состоянии генератора. Если вы действительно хотите, чтобы случайные числа при запусках программы были одними и теми же, вы не должны надеяться на постоянство мусора. Вам нужно обязательно присвоить переменным s1 и s2 какие-нибудь начальные значения. Маленькие числа присваивать не желательно. Уж очень плохого качества получаются при этом первые случайные числа. Мне нравится начинать с таких чисел: 1234 и 4321. Если же вам нужно, чтобы от запуска к запуску программа давала разные случайные числа, тогда вам нужно обеспечить, чтобы начальное значение генератора тоже было разным. Для этого обычно используют текущее время. Можно использовать встроенную функцию TIME которая вообще-то дает символьное представление времени (8 символов), но мы можем подсунуть ей целый массив и использовать те значения массива, куда лягут минуты и секунды. Следующая программа печатает первые 20 случайных чисел, причем делает это по-разному от запуска к запуску.
Fortran
1
2
3
4
5
6
7
8
9
10
11
integer*2 s1,s2,tm(4)
real*4 r
integer i
CALL TIME(tm)
s1=tm(3)
s2=tm(4)
do i=1,20
    call randu(s1,s2,r)
    print *,r
end do
end
Теперь рассмотрим тонкий вопрос. Может ли генератор дать нам число 1.0. Документация говорит, что не может. Мы, однако, хотим проверить. Для начала напишем программу, которая будет присваивать выработанное случайное выражение целой переменной. По правилам фортрана целая часть числа будет при этом отброшена. Нам остается поймать случай, когда в целую переменную попадет значение 1, и напечатать соответствующее состояние генератора.
Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
integer*2 s1,s2,t1,t2
real*4 r
integer n
s1=1234
s2=4321
do i=1,100000000
    t1=s1
    t2=s2
    call randu(s1,s2,r)
    n=r
    if(n==1) then
        print *,t1,t2
        exit
    end if
end do
end
И что вы думаете? Я несколько раз приписывал нуль к числу итераций, пока не добрался до ста миллионов. И только тогда программа напечатала 14570 21825. Теперь мы можем воспроизвести это чудо и отдельно с ним разобраться.
Fortran
1
2
3
4
5
6
7
8
9
10
11
integer*2 s1,s2
real*4 r
integer n
s1=14570
s2=21825
call randu(s1,s2,r)
n=r
if(r.eq.1.0) then
    print *,r,n
end if
end
Как видите, это число ведет себя как единица, так что скорее всего и является единицей. Мы уж не будем анализировать его двоичный вид. Понятно, что генератор случайных чисел основан на целой арифметике, и где-то внутри целое число преобразуется в real*4. Теоретически оно получается меньше единицы, но разрядной сетки real*4 не хватает для представления всех битов из целого числа. Округление и дает нам единицу. Не подумайте, что мы нашли невинную диковинку. На самом деле это очень существенную вещь. Лично я однажды попался на этом. Программа, через полгода устойчивой работы вдруг упала. (Правда, было это на бейсике, который дал мне понятную диагностику, позволившую мне догадаться, в чем дело.) В дальнейшем мы увидим, как возможность единицы заставляет нас усложнять алгоритмы и включать дополнительные проверки.

Классический вопрос, который появляется у новичка, -- это как заполнить случайными числами массив целого типа. Иногда дополнительно указывают диапазон возможных значений. Никаких условий на характер случайности не налагается, то есть предполагается, что все числа из диапазона равновероятны. Далее обычно нужно с этими числами что-нибудь сделать. Мы вычислим размах, то есть разность максимального и минимального числа диапазона. Итак, программа:
Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
integer*2 s1,s2,tm(4)
real*4 r
integer a(100) ! Массив, который нужно заполнить
integer lb,ub ! Нижняя и верхняя граница диапазона
integer i,p,imx,imn
call time(tm)
s1=tm(3)
s2=tm(4)
lb=-100
ub=100
! Генерация
do i=1,100
    call randu(s1,s2,r)
    p = (ub-lb+1)*r
    p = lb+p
    if(p>ub) p=ub ! проверка, нужная при r=1
    a(i)=p
end do
! Считаем размах
imx=a(1)
imn=a(1)
do i=2,100
    if(a(i)>imx) imx=a(i)
    if(a(i)<imn) imn=a(i)
end do
print *,a ! печать массива
print *,imx-imn ! печать размаха
end
Если есть возможность использовать средства фортрана 90, то генерацию можно написать несколько короче. Программа RANDOM_SEED, вызванная без параметров, инициализирует генератор, используя текущее время. О переменных, хранящих состояние генератора, заботиться уже не нужно. Вычисление размаха мы опустили.
Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
integer a(100) ! Массив, который нужно заполнить
integer lb,ub ! Нижняя и верхняя граница диапазона
integer i,p
call random_seed;
lb=-100
ub=100
! Генерация
do i=1,100
    call random_number(r)
    p = (ub-lb+1)*r
    p = lb+p
    if(p>ub) p=ub ! проверка, нужная при r=1
    a(i)=p
end do
print *,a ! печать массива
end
Вторая частая задача это сгенерировать массив так, чтобы целые числа не повторялись. Мы конкретизируем: перетасовать колоду карт (числа от 1 до 36).
Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
integer*2 s1,s2,tm(4)
real*4 r
integer a(36) ! Массив, который нужно заполнить
integer lb,ub ! Нижняя и верхняя граница диапазона
integer i,p,t
call time(tm)
s1=tm(3)
s2=tm(4)
lb=1
ub=36
! Начальное заполнение (в стиле фортрана 90)
forall(i=1:36) a(i)=i
! или в стиле фортрана 77
!   do i=1,36
!       a(i)=i
!   end do
! Тасование
do i=36,2,-1
    call randu(s1,s2,r)
    p = (ub-lb+1)*r
    p = lb+p
    if(p>ub) p=ub ! проверка, нужная при r=1
    if(p.ne.i) then
        t=a(i)
        a(i)=a(p)
        a(p)=t
    end if
    ub=ub-1
end do
print *,a ! печать массива
end
Последним примером будет генерация непрерывной случайной величины, имеющей нормальное распределение. Один из алгоритмов ее порождения основан на центральной предельной теореме. Если мы сложим 12 величин, имеющих равномерное распределение, то получим величину по распределению очень близкую к нормальному. Останется подогнать математическое ожидание и стандартное отклонение. Мы сгенерируем 10000 значений такой величины и посчитаем для сравнения выборочную среднюю и выборочную дисперсию.
Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
integer*2 :: s1=1234, s2=4321
real*4 r
real*8 :: Mean=3.0, Sigma=0.5 ! Требуемые средняя и стандартное отклонение.
real*8 vm, vs ! Выборочная средняя и выборочное стандартное отклонение
real*8 s
integer i, k
integer, parameter :: n=10000 ! объем выборки
 
vm=0
vs=0
do i=1,n
    ! начало вычисления нормальной величины
    s=0.0
    do k=1,12
        call randu(s1,s2,r)
        s=s+r
    end do
    s=(s-6.0)*Sigma+Mean
    ! конец вычисления нормальной величины
    vm=vm+s   ! накопление суммы
    vs=vs+s*s ! накопление суммы квадратов
end do
vm=vm/n ! выборочная средняя
vs=sqrt((vs/n-vm*vm)*n/(n-1)) ! выборочное стандарное отклонение
print *,vm, vs
end
Выборочные параметры полученной выборки совпали с заданными с точностью до нескольких тысячных.

(Эксперименты проводились на Intel Fortran 11.1)
Размещено в Без категории
Просмотров 6043 Комментарии 0
Всего комментариев 0
Комментарии
 
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Рейтинг@Mail.ru