c Dos generadores de números pseudo aleatorios y ejemplo de uso. c Para compilar este programa y ejecutarlo en linux: c gfortran random_main.f ; ./a.out program main c Genera n números independientes distribuidos uniformemente entre 0 c y 1, con dos generadores, random_number y rannr c y calcula el promedio de (k+1)*x**k, para k=1,2,3,4. c El valor esperado es 1 en los cuatro casos. c Lee n y la semilla (idum0). La semilla debe ser un entero no nulo c la primera vez. Las siguientes veces seguirá con la misma serie c hasta que se introduzca un valor no nulo otra vez. c En el caso de rannr, idum0 y -idum0 producen la misma serie. c real mrandom_number dimension a(4),b(4) c 1 write(*,*) ' n,idum0=?' read(*,*) n,idum0 idum1 = idum0 idum2 = idum0 c do k = 1,4 a(k) = 0 b(k) = 0 end do do i = 1,n x = mrandom_number(idum1) y = rannr(idum2) do k = 1,4 a(k) = a(k) + x**k b(k) = b(k) + y**k end do end do do k = 1,4 a(k) = (k+1)*a(k)/n b(k) = (k+1)*b(k)/n end do write(*,*) 'idum0,n,a,b' write(*,*) 'random_number',idum0,n,(a(k),k=1,4) write(*,*) 'rannr ',idum0,n,(b(k),k=1,4) c go to 1 end function mrandom_number(idum) c USO: idum debe ser una variable (no una constante) c La salida mrandom_number es un número aleatorio u.d. entre 0 y 1. c Si idum no es cero al entrar, inicializa la serie aleatoria con c idum como semilla, y produce el primer número de la serie. c Si idum es cero, produce el siguiente número de la serie. c En ambos casos idum es cero al salir. c random_number y random_seed son funciones intrínsecas de FORTRAN90 c pero el algoritmo en sí no es estándar. Aquí usará el de gfortran. parameter (nmx=120) real::mrandom_number integer seed(nmx) if(idum.ne.0) then call random_seed(size = n) if(n.gt.nmx) then write(*,*) 'mrandom_number::n.gt.nmx, n=',n stop end if do i = 1,n seed(i) = idum end do call random_seed(put = seed) idum = 0 end if call random_number(mrandom_number) end FUNCTION rannr(idum) ! Tomado de Numerical Recipes. c Modificado para que inicialice sii idum no es cero (sólo usa abs(idum)). c Al salir idum=0. IMPLICIT NONE INTEGER, PARAMETER :: K4B=selected_int_kind(9) INTEGER(K4B), INTENT(INOUT) :: idum REAL :: rannr ! “Minimal” random number generator of Park and Miller combined with ! a Marsaglia shift sequence. Returns a uniform random deviate ! between 0.0 and 1.0 (exclusive of the endpoint values). This fully ! portable, scalar generator has the “traditional” (not Fortran 90) ! calling sequence with a random deviate as the returned function ! value: call with idum a negative integer to initialize; ! thereafter, do not alter idum except to reinitialize. The period ! of this generator is about 3.1 × 10**18. INTEGER(K4B), PARAMETER :: IA=16807,IM=2147483647,IQ=127773,IR $ =2836 REAL, SAVE :: am INTEGER(K4B), SAVE :: ix=-1,iy=-1,k if (idum .ne. 0 .or. iy < 0) then ! Initialize. am=nearest(1.0,-1.0)/IM iy=ior(ieor(888889999,abs(idum)),1) ix=ieor(777755555,abs(idum)) c idum=abs(idum)+1 ! Set idum positive. idum=0 ! MODIFICADO: Pongo idum a cero. end if ix=ieor(ix,ishft(ix,13)) ! Marsaglia shift sequence with period 232 − 1. ix=ieor(ix,ishft(ix,-17)) ix=ieor(ix,ishft(ix,5)) k=iy/IQ ! Park-Miller sequence by Schrage’s method, iy=IA*(iy-k*IQ)-IR*k ! period 231 − 2. if (iy < 0) iy=iy+IM rannr=am*ior(iand(IM,ieor(ix,iy)),1) ! Combine the two generators with ! masking to ensure nonzero value. END FUNCTION rannr