! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 2b
!
! numeros aleatorios - subrutina portable utilizando
! el metodo de sustraccion (Knuth, Seminumerical Algorithms)
!
! distribucion uniforme en [0,1)
!
!
program eje2b
implicit none
  integer :: i,j,k,old
  integer, dimension(10) :: arr
  real :: valr
  
! pongo a cero el array
!
  arr=0
  
! inicializo la subrutina que genera numeros aleatorios
!
  write(*,*)'Ingrese semilla (int): '
  read(*,*)j
  valr=-float(j)
  
! itero j veces
!
  write(*,*)'Ingrese nro. iteraciones: '
  read(*,*)j
  
  do i=1,j
  
! obtengo un numero aleatorio
!
    call random1(valr)

! lo guardo en un histograma
!
	k=int(valr*10.)+1
	arr(k)=arr(k)+1
  enddo
  
! encuentro el bin con el mayor valor
!
  old=-12
  k=0
  do i=1,10
    if(arr(i) > old)then
	  old=arr(i)
	  k=i
	endif
  enddo
  
! muestro el resultado
!
  do i=1,10
    write(*,*)i,arr(i),float(arr(i))/float(arr(k))
  enddo
  
end program eje2b


! *******************************************************
! subrutina mejorada para generar numeros
! aleatorios con distribucion uniforme en [0,1)
! y portable (no usa random_number)
!
! la subrutina devuelve en val el numero
! aleatorio. Si en la entrada val es negativo
! la subrutina se inicializa
!
subroutine random1(val)
implicit none
  integer, parameter :: mbig=1000000000, mseed=161803398, mz=0
  real, parameter :: fac=1./mbig
  integer :: i,k,mj,mk,ii
  integer, save :: inext,inextp
  integer, dimension(55), save :: ma
  real :: val
  
! inicializa
!
  if(val < 0.)then
    mj=mseed+int(val)
	mj=mod(mj,mbig)
	ma(55)=mj
	mk=1
	do i=1,54
	  ii=mod(21*i,55)
	  ma(ii)=mk
	  mk=mj-mk
	  if(mk < mz)mk=mk+mbig
	  mj=ma(ii)
	enddo
	do k=1,4
	  do i=1,55
	    ma(i)=ma(i)-ma(1+mod(i+30,55))
		if(ma(i) < mz)ma(i)=ma(i)+mbig
	  enddo
	enddo
	inext=0
	inextp=31
  endif
  
! punto de entrada general
!
  inext=inext+1
  if(inext == 56)inext=1
  inextp=inextp+1
  if(inextp == 56)inextp=1
  mj=ma(inext)-ma(inextp)
  if(mj < mz)mj=mj+mbig
  ma(inext)=mj
  val=float(mj)*fac
  return
end subroutine random1
! *******************************************************