random - suma - restas negativas para niños

las recetas numéricas ran3 generan números negativos

Question

Estoy usando un esquema de recetas numéricas para generar números aleatorios ( ran3 , página 7 en este archivo PDF ). No noté nada extraño, pero esta vez obtuve un número negativo en la etapa de "calentamiento" que es más grande que MBIG . El código parece que esto no debería suceder. Puedo solucionarlo fácilmente cambiando la instrucción if por una declaración while en la línea que dice if(mk.lt.MZ)mk=mk+MBIG pero quiero saber cuáles son las consecuencias.

Editar: aquí está la función

 FUNCTION ran3a(idum)
INTEGER idum
INTEGER MBIG,MSEED,MZ
C REAL MBIG,MSEED,MZ
REAL ran3a,FAC
PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
INTEGER i,iff,ii,inext,inextp,k
INTEGER mj,mk,ma(55)
C REAL mj,mk,ma(55)
SAVE iff,inext,inextp,ma
DATA iff /0/
if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=MSEED-iabs(idum)
mj=mod(mj,MBIG)
ma(55)=mj
mk=1
do 11 i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.MZ)mk=mk+MBIG
mj=ma(ii)
11 continue
do 13 k=1,4
do 12 i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
12 continue
13 continue
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.MZ)mj=mj+MBIG
ma(inext)=mj
ran3a=mj*FAC
return
END

Answer 1

Echó un vistazo al pdf. Lo que debes hacer es

1) Semilla: valor = ran3 (-1)

2) Úselo: value = ran3 (0)

Answer 2

Obtenía Seg Faults (usando gfortran 4.8) porque la función intentaba cambiar el valor de entrada idum del número negativo a 1. No hay razón para esa línea (ni nada con iff ), así que la borré e imprimí la array ma en varios lugares diferentes y no encontró números negativos en la matriz.

Una posibilidad, sin embargo, es que si iabs(idum) es más grande que MSEED , podría tener un problema con la línea mj=MSEED - iabs(idum) . Debe protegerse de esto usando mj=abs(MSEED-abs(idum)) como el libro ha escrito.