random number generator random() questionable!!!?

[Alex Martelli]

ben at dit.lth.se (Ben Smeets) writes:

>Hi:

>After doing some simulations I got suspicious about the UNIX random
>number generator random() available on the SUN3-xx machines (ver 4.1).
	...
>	2) Has anybody else noted some stange things with
>	   this random number generator.?

Don't know about this particular one, but random number generators are
notoriously easy-to-botch...

>	3) Are the people who can provide me with a good random
>	   number generator?

Knuth, "Art of Computer programming", volume 2, is a great reference.
Here's an implementation which is SLOW (boy, is it ever!) but gives
random uniform variates with GREAT spectral characteristics, and
(very important to ME) the SAME sequence of numbers on ANY platform 
whatsoever (should not be hard to recode in C):

       FUNCTION RAN1(IDUM)
*
* PORTABLE random-number generation routine, computationally costly
* but with nonpareil characteristics.  From "Numerical recipes", p.196;
* identifiers as in the original; added explicit declarations & SAVE.
*
* Call with IDUM < 0 to re-seed (not needed the first time), IDUM == 0
* for normal operation.
*
* Note the SAME sequence is generated on ANY machine (crucial for use
* of random numbers for any intra-machine test!).
*
       REAL*4 RAN1,R(97),RM1,RM2
       INTEGER*4 IDUM,M1,M2,M3,IA1,IA2,IA3,IC1,IC2,IC3,IFF,IX1,IX2,IX3,J
       PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=1.0/M1)
       PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1.0/M2)
       PARAMETER (M3=243000,IA3=4561,IC3=51349)
       SAVE IFF, R, IX1, IX2, IX3
       DATA IFF /0/
*
* Initialize on first call, or whenever the parm is negative
*
      IF(IDUM.lt.0 .or. IFF.eq.0) THEN
          IFF=1
*                                   seed first routine
          IX1=MOD(IC1-IDUM,M1)
          IX1=MOD(IA1*IX1+IC1,M1)
*                                   seed second routine
          IX2=MOD(IX1,M2)
          IX1=MOD(IA1*IX1+IC1,M1)
*                                   seed third routine
          IX3=MOD(IX1,M3)
*                                   fill table w. 1st/2nd routines
          DO 11, J=1,97
                 IX1=MOD(IA1*IX1+IC1,M1)
                 IX2=MOD(IA2*IX2+IC2,M2)
*                                   combine hi with lo
                 R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
11        CONTINUE
          IDUM=1
      ENDIF
	  NCALLS = NCALLS + 1
*
* Starting point, save when (re)initializing: get next number on 
* each of the three sequences
*
       IX1=MOD(IA1*IX1+IC1,M1)
       IX2=MOD(IA2*IX2+IC2,M2)
       IX3=MOD(IA3*IX3+IC3,M3)
*                                   take index from 3rd sequence
       J=1+(97*IX3)/M3
*                                   sanity check, may omit in production
       IF(J.GT.97.OR.J.LT.1) THEN
		   WRITE(*,*) 'RAN1: J=',J,' IX3=',IX3
		   STOP
	   ENDIF
*                                   return this entry of table...
       RAN1=R(J)
*                                   ...and refill it
       R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
       RETURN
       END

-- 
Alex Martelli - CAD.LAB s.p.a., v. Stalingrado 45, Bologna, Italia
Email: (work:) staff at cadlab.sublink.org, (home:) alex at am.sublink.org
Phone: (work:) ++39 (51) 371099, (home:) ++39 (51) 250434; 
Fax: ++39 (51) 366964 (work only; any time of day or night).