Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум программистов > Fortran > [General] Random numbers in F77


Автор: valvliv 2.3.2006, 17:52
Тут нужно нагенерить случайных чисел, а код на 77-м фортране, поэтому функцию из 90-го не вплетешь. Я посмотрела в сети - кто только как не предлагает, и каждый код - по странице (причем не всегда компеляемый).

Есть ли у кого опыт в этом деле, чтобы просто и понятно??

Л.

Автор: Cr@$h 2.3.2006, 19:35
Лично я доверяю процедуре из книги Форсайта "Машинные методы математических вычислений". Вот код процедуры на F66 (77?).
Код

      REAL FUNCTION URAND(IY)
      INTEGER IY
C
C       URAND-ЭTO ДATЧИK PABHOMEPHO PACПPEДEЛEHHЫX CЛУЧAЙHЫX
C     ЧИCEЛ, OCHOBAHHЫЙ HA TEOPИИ И ПPEДЛOЖEHИЯX, COДEPЖAЩИXCЯ
C     B KHИГE KHУT (1969),TOM 2.
C       ПEPEД ПEPBЫM OБPAЩEHИEM K URAND ЦEЛOЙ ПEPEMEHHOЙ IY
C     CЛEДУET ПPИCBOИTЬ ПPOИЗBOЛЬHOE ЦEЛOЧИCЛEHHOE HAЧAЛЬ-
C     HOE ЗHAЧEHИE. BЫЗЫBAЮЩAЯ ПPOГPAMMA HE ДOЛЖHA ИЗMEHЯTЬ
C     ЗHAЧEHИE IY MEЖДУ ПOCЛEДOBATEЛЬHЫMИ BЫЗOBAMИ. ЗHAЧEHИЯ
C     ФУHKЦИИ URAND ЯBЛЯЮTCЯ ЧИCЛAMИ ИЗ ИHTEPBAЛA (0,1).
C
      INTEGER IA,IC,ITWO,M2,M,MIC
      DOUBLE PRECISION HALFM
      REAL S
      DOUBLE PRECISION DATAN,DSQRT
      DATA M2/0/,ITWO/2/
      IF(M2.NE.0) GOTO 20
C
C     ECЛИ ЭTO ПEPBЫЙ BXOД, TO BЫЧИCЛИTЬ ДЛИHУ
C     ЦEЛOЧИCЛEHHOГO MAШИHHOГO CЛOBA
C
      M=1
   10 M2=M
      M=ITWO*M2
      IF(M .GT. M2) GO TO 10
      HALFM=M2
C
C     BЫЧИCЛИTЬ MHOЖИTEЛЬ И ПPИPAЩEHИE ЛИHEЙHOГO
C     KOHГPУЭHTHOГO METOДA
C
      IA=8*IDINT(HALFM*DATAN(1.D0)/8.D0)+5
      IC=2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0))+1
      MIC=(M2-IC)+M2
C
C     S-MACШTAБИPУЮЩИЙ MHOЖИTEЛЬ ДЛЯ ПPEOБPAЗOBAHИЯ B ЧИCЛO
C     C ПЛABAЮЩEЙ TOЧKOЙ
C
      S=0.5/HALFM
C
C     BЫЧИCЛИTЬ CЛEДУЮЩEE CЛУЧAЙHOE ЧИCЛO
C
   20 IY=IY*IA
C
C     CЛEДУЮЩИЙ OПEPATOP-ДЛЯ MAШИH,KOTOPЫE HE ДOПУCKAЮT
C     ПEPEПOЛHEHИЯ ЦEЛЫX ЧИCEЛ ПPИ CЛOЖEHИИ
C
      IF(IY.GT.MIC) IY=(IY-M2)-M2
      IY=IY+IC
C
C     CЛEДУЮЩИЙ OПEPATOP - ДЛЯ MAШИH, У KOTOPЫX ДЛИHA CЛOBA
C     ДЛЯ CЛOЖEHИЯ БOЛЬШE, ЧEM ДЛЯ УMHOЖEHИЯ
C
      IF(IY/2.GT.M2) IY=(IY-M2)-M2
C
C     CЛEДУЮЩИЙ OПEPATOP - ДЛЯ MAШИH, У KOTOPЫX ПEPEПOЛHEHИE
C     ЦEЛOГO ЧИCЛA BЛИЯET HA ЗHAKOBЫЙ PAЗPЯД
C
      IF(IY.LT.0) IY=(IY+M2)+M2
      URAND=FLOAT(IY)*S
      RETURN
      END


Я её, конечно, переписал под F95:
Код

function URand( IY ) result(result)
    integer, parameter :: R_ = 4
    integer, parameter :: I_ = R_
    integer, parameter :: R2_ = R_ * 2

    real(R_) result
    integer(I_), intent(inout) :: IY

    integer(I_), save :: Ia, Ic, M2, MIc
    real(R2_), save :: halfM
    real(R_), save :: s
    data M2 /0/

    if( M2 == 0 ) then
        ! Если это первый вход, то вычисление длины целочисленного машинного
        ! слова.

        M2 = Huge( M2 ) / 2 + 1
        halfM = real( M2, R2_ )

        ! Вычисление множителя и приращения линейного конгруэнтного метода.

        Ia = 8*Int( halfM * Atan(1._R2_) / 8._R2_, I_ ) + 5
        Ic = 2*Int( halfM * ( 0.5_R2_ - SqRt(3._R2_)/6._R2_ ), I_ ) + 1
        MIc = ( M2 - Ic ) + M2

        ! S - масштабирующий множитель для преобразования в число с плавающей
        ! точкой.
        s = 0.5_R_ / real( halfM, R_ )
    end if

    !Вычисление следующего случайного числа.
    IY = IY * Ia

    ! Следующий оператор - для машин, которые не допускают переполнения целых
    ! чисел при сложении.
    if( IY > MIc )&
        IY = ( IY - M2 ) - M2

    IY = IY + Ic

    ! Следующий оператор - для машин, у которых длина слова для сложения больше,
    ! чем для умножения.
    if( IY / 2 > M2 ) then
        IY = ( IY - M2 ) - M2
    ! Следующий оператор - для машин, у которых переполнение целого числа
    ! влияет на знаковый разряд.
    else if( IY < 0 ) then
        IY = ( IY + M2 ) + M2
    end if

    result = real( IY, R_ ) * s
end function URand


Примерчик по использованию последней:
Код

program ExURand
    implicit none

    integer IY, i

    IY = 0
    do i = 1, 10
        if( URand( IY ) < 0.5 ) then
            print *, "tails"
        else
            print *, "heads"
        end if
    end do

    read *

end program ExURand

Процедура заточена на генерацию числа за числом. Но её можно модифицировать и для возвращения случайного массива. Наверное, эта процедура тоже очень надуманная, но, тем не менее, обладает хорошими статистическими характеристиками.
Могу ещё написать, если что-то не понравится. smile

Автор: valvliv 2.3.2006, 19:55
Мда... Бум изучать. Я вообще-то хотела если не как в матлабе, то хоть типа itime, но у меня проект крутит функции циклически-регулярно, и такая генерилка дает аккуратную синусоиду, что не комильфо с точки зрения случайных чисел.

Простоты хочется, простоты...
smile

Автор: Cr@$h 2.3.2006, 20:13
В чем нужна простота? В интерфейсе, в объеме кода, генерим одно число или массив. С синусоидой не совсем понял... При очередном входе можно задавать разное начальное число, например, через Cpu_Time, тогда выборка не будет повторяться, если это про это.
Добавлено @ 20:15
Али целые числа нужны...

Автор: valvliv 2.3.2006, 20:52
> В чем нужна простота? В интерфейсе, в объеме кода, генерим одно число или массив.

В объеме кода, точно. Я вот сейчас размышляю - а не сгенерить ли ряд в матлабе и сунуть его кушать фортрану... Там будет красивый гауссов белый шум, без обмана...

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

Я брала секунды момента обращения к процедуре и рескейлила в интервал [0,1). Ожидала типа равномерно распределенный шум, но код работает очень методично smile

Л.

Автор: Cr@$h 2.3.2006, 22:05
Странно, я URand проверял... Быть может, не так сильно на неповторяемость.
Не понял делему: либо брать из MATLAB, либо небольшой код на Fortran. Почему так? И что еще больше будет занимать...
Тем не менее, убежден, что хороший генератор будет комбинированным, а поэтому немаленьким по коду. А код нужен из-за производительности или эстетического удовольствия?
Буду исследовать...

Автор: valvliv 2.3.2006, 22:08
Думаю, с урандом все в порядке smile Это я балда ленивая smile

Автор: Cr@$h 2.3.2006, 23:51
Ну, слава Форсайту, в смысле Фортрану. Если она все-таки не подходит или её нужно будет для массивов переделать -- тему продолжить мона.
Не знаю для каких целей и платформ тебе, но URand писалась робасто. Совет по упрощению её кода при некоторых условиях:
  • Система представления чисел, на которой будет запускаться (компилиться) прога одинаковая. Например, проект будет компилиться и под винды, и под линуху и где бы то ни было ещё, но представление чисел будет то же: integer и в Африке -- integer, а real будет имееть такие же границы, скажем. В этом случае весь первый if-блок выполняем заранее для машины/платформы, сосчитав эти статические параметры датчика.
  • Генерить нужно будет только real(4) или только real(8). В этом случае все константы и параметры типов задаем явно, а не через константы (R_, I_). Одно условие: integer должен занимать столько же байт, сколько генерируемый real, чтобы уместить в себе все её значащие цифры.

Например, если я хочу иметь поменьше код генерации для определенной (моей) платформы, то делаю из процедуры:
Код

function URand( IY ) result(result)
    integer, parameter :: R_ = 4
    integer, parameter :: I_ = R_
    integer, parameter :: R2_ = R_ * 2

    real(R_) result
    integer(I_), intent(inout) :: IY

    integer(I_), parameter :: Ia = 843314861, Ic = 453816693, M2 = 1073741824,&
        MIc = 1693666955
    real(R_), parameter :: s = 4.6566129E-10
    
    !Вычисление следующего случайного числа.
    IY = IY * Ia

    ! Следующий оператор - для машин, которые не допускают переполнения целых
    ! чисел при сложении.
    if( IY > MIc )&
        IY = ( IY - M2 ) - M2

    IY = IY + Ic

    ! Следующий оператор - для машин, у которых длина слова для сложения больше,
    ! чем для умножения.
    if( IY / 2 > M2 ) then
        IY = ( IY - M2 ) - M2
    ! Следующий оператор - для машин, у которых переполнение целого числа
    ! влияет на знаковый разряд.
    else if( IY < 0 ) then
        IY = ( IY + M2 ) + M2
    end if

    result = real( IY, R_ ) * s
end function URand

Просто я взял да посчитал все константы.
Можно убрать константы типа R_ и I_, если код ориентируется на одну разновидность (типа менять потом его не придется вдруг).
Можно поисследовать вопрос и убрать ненужные из трех if-ов. Здесь я уже не помню всех этих премудростей. Собственно, компактный датчик: убрать вначале константы разновидностей типа, сосчитать параметры датчика заранее, убрать if-ы. Но можем потерять робастость, переносимость, надежность, но полуим точно прирост в скорости, не потеряем точно.

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)