Форум программистов, компьютерный форум, киберфорум
palva
Войти
Регистрация
Восстановить пароль
Блоги Сообщество Поиск  

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

Запись от palva размещена 24.12.2013 в 03:11
Показов 42846 Комментарии 0

Для генерации случайных чисел в фортране имеется встроенная подпрограмма 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)
Надоела реклама? Зарегистрируйтесь и она исчезнет полностью.
Всего комментариев 0
Комментарии
 
Новые блоги и статьи
Установка нескольких штампов электронной подписи в строго определенных местах файла docx
ВладимирСамохин 19.07.2026
(В!) Работа с Электронной подписью - это неотъемлемая часть современного документооборота. Но что делать, если нужно поставить несколько штампов электронной подписи в строго определенных местах. . .
сукцессия 35. Научная статья о проделанной работе
anaschu 19.07.2026
Написал в формате латекс и пдф
Вангую, что это не пройдёт модерацию, и на неделе я запущу свой сервер.
Hrethgir 19.07.2026
Эта публикация сейчас в песочнице и ждёт приглашения. https:/ / habr. com/ ru/ sandbox/ 295048/ начало и оглавление - Как «пернатого» заставить осваивать новые горизонты опыта через масштабирование. . .
сукцессия 33. открытые вопросы от клауде
anaschu 19.07.2026
"Что накопилось за эту часть А — тринадцать правок, из которых шесть пришли из ваших вопросов и каждая оказалась реальной ошибкой, а не калибровкой: односторонний симбиоз, отсутствующий листопад,. . .
32 сукцессия
anaschu 19.07.2026
сукцессия 28‑мерное ядро стабилизировано Коллеги, фиксирую разбор инженерных правок и их изоморфную проекцию на экономику, меметику и половой отбор. Модель теперь не «подкручивает» сходимость —. . .
сукцессия 31: модель микоризы - это модель ещё нескольких явлений, социальных и экономических
anaschu 18.07.2026
Теория «Всего»: апдейт v1. 1. 2 — 28‑мерное ядро стабилизировано Коллеги, фиксирую разбор инженерных правок и их изоморфную проекцию на экономику, меметику и половой отбор. Модель теперь не. . .
сукцессия 30. Массив проверяющих друг друга моделей
anaschu 18.07.2026
Архитектура сети взаимопроверяющих моделей микоризной сукцессии (v2. 0) Развитие тензорного ОДУ-ядра и создание кросс-платформенного калибровочного полигона Уважаемые коллеги! В продолжение. . .
Грибы - это женщины, деревья - это мужчины. Анти инь янь для союза мужчины и женщины.
anaschu 18.07.2026
ГЛАВНЫЙ НАУЧНО-ФИЛОСОФСКИЙ ВЫВОД: Сексуально-Репродуктивный Капитализм против Государства Моногамии Коллеги, мы вышли на финишную прямую 20-мерного ОДУ-моделирования вековой сукцессии (ветка. . .
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2026, CyberForum.ru