Author 
Message 
Eva Sabrona Simmo #1 / 10

Random Number Generators
I have been trying to use the 'rand()' random number generator, but it doesn't work adequately. Now, I am trying to find other random number generators that might be of use. I have found through the man pages that 'random()', 'ran()', and 'lcrans' are available. It seems that 'lcrans' would be best for my purposes, but I am programming in fortran and this is a C Library Function. Any ideas on how I can access this function  or any other function in the C Library  from a FORTRAN program? Also, does anyone have any ideas for random number generators that might be of use to me? I would prefer a pseudorandom number generator... Thanks, Eva S. Simmons  ***THE SIMMONS FACTOR ****EVA SABRINA SIMMONS******UTAUSTIN GRAD STUDENT*****
WWW Personal Page: http://www.***.com/ ~weevey ****************** WATCH IT, OR IT MIGHT ATTACK!! ****************************

Fri, 29 Aug 1997 09:11:43 GMT 


Dan P #2 / 10

Random Number Generators
Quote: >I have been trying to use the 'rand()' random number generator, >but it doesn't work adequately. Now, I am trying to find other >random number generators that might be of use. I have found >through the man pages that 'random()', 'ran()', and 'lcrans' >are available. It seems that 'lcrans' would be best for my >purposes, but I am programming in FORTRAN and this is a C Library >Function. Any ideas on how I can access this function  or >any other function in the C Library  from a FORTRAN program?
By reading the documentation of your compilers. Since you didn't mention the system you're using, this is all we can provide as help. On VMS, you simply call your function and specify the C library at link time. On Unix, you write a wrapper function in C and don't specify any library at link time (the C library will be searched by default). Quote: >Also, does anyone have any ideas for random number generators that >might be of use to me? I would prefer a pseudorandom number >generator...
Without specialized hardware (e.g. sources and detectors of radiations) all the generators are pseudorandom. One of the volumes of Knuth's "The Art of Computer Programming" is discussing the problem in great detail. Check your library. Dan  Dan Pop CERN, CN Division
Mail: CERN  PPE, Bat. 31 R004, CH1211 Geneve 23, Switzerland

Fri, 29 Aug 1997 19:39:34 GMT 


Byron Bod #3 / 10

Random Number Generators
Quote: > I have been trying to use the 'rand()' random number generator, > but it doesn't work adequately. Now, I am trying to find other > random number generators that might be of use. I have found > through the man pages that 'random()', 'ran()', and 'lcrans' > are available. It seems that 'lcrans' would be best for my > purposes, but I am programming in FORTRAN and this is a C Library > Function. Any ideas on how I can access this function  or > any other function in the C Library  from a FORTRAN program? > Also, does anyone have any ideas for random number generators that > might be of use to me? I would prefer a pseudorandom number > generator... > Thanks, > Eva S. Simmons
Here's an earlier Marsaglia generator that is easier to implement. bb ====================================================== C This random number generator originally appeared in "Toward a Universal C Random Number Generator" by George Marsaglia and Arif Zaman. C Florida State University Report: FSUSCRI8750 (1987) C C It was later modified by F. James and published in "A Review of Pseudo C random Number Generators" C C THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE. C (However, a newly discovered technique can yield C a period of 10^600. But that is still in the development stage.) C C It passes ALL of the tests for random number generators and has a period C of 2^144, is completely portable (gives bit identical results on all C machines with at least 24bit mantissas in the floating point C representation). C C The algorithm is a combination of a Fibonacci sequence (with lags of 97 C and 33, and operation "subtraction plus one, modulo one") and an C "arithmetic sequence" (using subtraction). C C On a Vax 11/780, this random number generator can produce a number in C 13 microseconds. C======================================================================== PROGRAM TstRAN REAL temp(100) INTEGER IJ, KL, len C These are the seeds needed to produce the test case results IJ = 1802 KL = 9373 C Do the initialization call rmarin(ij,kl) C Generate 20000 random numbers len = 100 do 10 i = 1, 200 call RANMAR(temp, len) 10 continue C If the random number generator is working properly, the next six random C numbers should be: C 6533892.0 14220222.0 7275067.0 C 6172232.0 8354498.0 10633180.0 len = 6 call RANMAR(temp, len) write(6,20) (4096.0*4096.0*temp(I), I=1,6) 20 format (3f12.1) end subroutine RMARIN(IJ,KL) C This is the initialization routine for the random number generator RANMAR() C NOTE: The seed variables can have values between: 0 <= IJ <= 31328 C 0 <= KL <= 30081 C The random number sequences created by these two seeds are of sufficient C length to complete an entire calculation with. For example, if sveral C different groups are working on different parts of the same calculation, C each group could be assigned its own IJ seed. This would leave each group C with 30000 choices for the second seed. That is to say, this random C number generator can create 900 million different subsequences  with C each subsequence having a length of approximately 10^30. C C Use IJ = 1802 & KL = 9373 to test the random number generator. The C subroutine RANMAR should be used to generate 20000 random numbers. C Then display the next six random numbers generated multiplied by 4096*4096 C If the random number generator is working properly, the random numbers C should be: C 6533892.0 14220222.0 7275067.0 C 6172232.0 8354498.0 10633180.0 real U(97), C, CD, CM integer I97, J97 logical TEST data TEST /.FALSE./ common /raset1/ U, C, CD, CM, I97, J97, TEST if( IJ .lt. 0 .or. IJ .gt. 31328 .or. * KL .lt. 0 .or. KL .gt. 30081 ) then print '(A)', ' The first random number seed must have a value *between 0 and 31328' print '(A)',' The second seed must have a value between 0 and *30081' stop endif i = mod(IJ/177, 177) + 2 j = mod(IJ , 177) + 2 k = mod(KL/169, 178) + 1 l = mod(KL, 169) do 2 ii = 1, 97 s = 0.0 t = 0.5 do 3 jj = 1, 24 m = mod(mod(i*j, 179)*k, 179) i = j j = k k = m l = mod(53*l+1, 169) if (mod(l*m, 64) .ge. 32) then s = s + t endif t = 0.5 * t 3 continue U(ii) = s 2 continue C = 362436.0 / 16777216.0 CD = 7654321.0 / 16777216.0 CM = 16777213.0 /16777216.0 I97 = 97 J97 = 33 TEST = .TRUE. return end subroutine ranmar(RVEC, LEN) C This is the random number generator proposed by George Marsaglia in C Florida State University Report: FSUSCRI8750 C It was slightly modified by F. James to produce an array of pseudorandom C numbers. REAL RVEC(*) real U(97), C, CD, CM integer I97, J97 logical TEST common /raset1/ U, C, CD, CM, I97, J97, TEST integer ivec if( .NOT. TEST ) then print '(A)',' Call the init routine (RMARIN) before calling RAN *MAR' stop endif do 100 ivec = 1, LEN uni = U(I97)  U(J97) if( uni .lt. 0.0 ) uni = uni + 1.0 U(I97) = uni I97 = I97  1 if(I97 .eq. 0) I97 = 97 J97 = J97  1 if(J97 .eq. 0) J97 = 97 C = C  CD if( C .lt. 0.0 ) C = C + CM uni = uni  C if( uni .lt. 0.0 ) uni = uni + 1.0 RVEC(ivec) = uni 100 continue return end 

Fri, 29 Aug 1997 22:14:44 GMT 


Byron Bod #4 / 10

Random Number Generators
Quote: > I have been trying to use the 'rand()' random number generator, > but it doesn't work adequately. Now, I am trying to find other > random number generators that might be of use. I have found > through the man pages that 'random()', 'ran()', and 'lcrans' > are available. It seems that 'lcrans' would be best for my > purposes, but I am programming in FORTRAN and this is a C Library > Function. Any ideas on how I can access this function  or > any other function in the C Library  from a FORTRAN program? > Also, does anyone have any ideas for random number generators that > might be of use to me? I would prefer a pseudorandom number > generator... > Thanks, > Eva S. Simmons
Here's a supposedly superior generator that was posted a while back. I haven't tried it yet. bb ======================================================== Hello everyone! Well, here's what you've all been waiting for. Again, would someone please crosspost this to the relevant groups (numerical analysis and so on)? I am including, separated by a row of 72 stars, my response (received today) from Fred James from CERN, who wrote the code, and my original email from James from about a year ago. What you should cut out and compile starts right after the row of plus signs in this original email. This consists of the SUBROUTINE RANLUX, up until the END statement, and then the PROGRAM RLXTST, which you can use to make sure everything is working correctly on you computer, which includes the output which should be produced in the comment lines at the end. All I did was delete James's opening remarks, separate the subroutine and the program, and compile them separately. Except for the C's at the beginning of the lines, of course, the test output on my computer (VAX 3100) was identical to the test output in the RLXTST comments. All in all, less than 2 minutes work. Which speaks a lot for coding in absolutely standard FORTRAN, as well as for the VAX compiler. (While most compilers, if not all, which offer extensions of course support `only' standard FORTRAN, some, such as IBM AIX on a RISC 6000, have to be given options in order _NOT_ to support extensions which are INCOMPATIBLE with the standard (like folding to lower case). If you start a compiler with the defaults, then my view is that it is OK to support extensions, but it should at least compile standard FORTRAN without having to be explicitly told that this is standard FORTRAN or otherwise produce results incompatible with the standard.) The Martin Luescher mentioned is a professor at DESY/University of Hamburg who does theoretical elementary particle physics and who wrote an article (DESY 93133, September 1993, ISSN 04189833) entitled `A Portable HighQuality Random Number Generator for Lattice Field Theory Simulations' which I accidentally ran into in the Hamburg Observatory library. Luescher directed me to Fred James, and since then we've used the generator quite a bit, tested it extensively just to be safe and have noticed no problems. For most purposes, one could probably use the highest level with no problems as far as CPU time goes. Have fun!! Phillip Helbig Tel. .............. +49 40 7252 4110
Gojenbergsweg 112 Fax ............... +49 40 7252 4198 D21029 Hamburg Telex ............... 217884 hamst d ************************************************************************ Hello, Thank you for your message about RANLUX. I am happy to hear of satisfied users. Yes, I know a little about some of the various rng's that are being proposed these days, in particular by the Numerical Recipes people. I have told Bill Press about RANLUX and he has the references. The area of rng's was always one of the weakest points in Num. Recip., and I guess that will continue to be true, although their algorithms will improve of course. They do a very good job in so many areas, but you can't really expect them to be best in everything. Press seems to be working closely with Marsaglia, which I also did at one time, until I discovered Martin Luescher and Pierre L'Ecuyer who have a more rigorous and more convincing approach to the problem. You may know that Marsaglia once said: "Random numbers are like sex: ... Even when they are bad they are still pretty good." However, we must remember that the actual RANLUX algorithm owes much to several people including Marsaglia: 1. Marsaglia invented the "subtractwithborrow" algorithm which has many practical advantages (fast and portable because all computations are done in portable floatingpoint, very long period, and reasonably easy to initialize, restart, etc.), but is not sufficiently random. 2. L'Ecuyer and Tuzaka (and perhaps others) recognized and proved that Marsaglia's algorithm is in fact equivalent to a linear congruential generator (which is bad) with a very big multiplier (which is good). This gave a first basis for understanding the randomness properties and in particular the flaws. 3. Luescher, using Kolmogorov's chaos theory, showed rigorously how to improve the algorithm by skipping, and how much you have to skip to make it lose all memory of its past. He also can calculate the period exactly, with or without skipping. 4. I wrote the Fortran version of the generator. Feel free to propagate RANLUX, but make sure to give proper credit always to Luescher and the two articles in Computer Physics Communications where it all is published. Thanks again, Fred. ************************************************************************ I guess that theoretical cosmology is easier than experimental cosmology (less timeconsuming) but this random number generator should be good enough for both. The references are Lu:scher, Computer Physics Communications 79 (1994) 100, and James, CPC 79 (1994) 111. +++++++++++++++++++++++++++++++++++ SUBROUTINE RANLUX(RVEC,LENV) C Subtractandborrow random number generator proposed by C Marsaglia and Zaman, implemented by F. James with the name C RCARRY in 1991, and later improved by Martin Luescher C in 1993 to produce "Luxury Pseudorandom Numbers". C Fortran 77 coded by F. James, 1993 C C references: C M. Luscher, Computer Physics Communications 79 (1994) 100 C F. James, Computer Physics Communications 79 (1994) 111 C C LUXURY LEVELS. C   The available luxury levels are: C C level 0 (p=24): equivalent to the original RCARRY of Marsaglia C and Zaman, very long period, but fails many tests. C level 1 (p=48): considerable improvement in quality over level 0, C now passes the gap test, but still fails spectral test. C level 2 (p=97): passes all known tests, but theoretically still C defective. C level 3 (p=223): DEFAULT VALUE. Any theoretically possible C correlations have very small chance of being observed. C level 4 (p=389): highest possible luxury, all 24 bits chaotic. C C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C!!! Calling sequences for RANLUX: ++ C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++ C!!! 32bit random floating point numbers between ++ C!!! zero (not included) and one (also not incl.). ++ C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++ C!!! one 32bit integer INT and sets Luxury Level LUX ++ C!!! which is integer between zero and MAXLEV, or if ++ C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++ C!!! should be set to zero unless restarting at a break++ C!!! point given by output of RLUXAT (see RLUXAT). ++ C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++ C!!! which can be used to restart the RANLUX generator ++ C!!! at the current point by calling RLUXGO. K1 and K2++ C!!! specify how many numbers were generated since the ++ C!!! initialization with LUX and INT. The restarting ++ C!!! skips over K1+K2*E9 numbers, so it can be long.++ C!!! A more efficient but less convenient way of restarting is by: ++ C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++ C!!! ISVEC of 25 32bit integers (see RLUXUT) ++ C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++ C!!! 32bit integer seeds, to be used for restarting ++ C!!! ISVEC must be dimensioned 25 in the calling program ++ C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ DIMENSION RVEC(LENV) DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25) PARAMETER (MAXLEV=4, LXDFLT=3) DIMENSION NDSKIP(0:MAXLEV) DIMENSION NEXT(24) PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265) PARAMETER (ITWO24=2**24, ICONS=2147483563) SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED INTEGER LUXLEV LOGICAL NOTYET DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/ DATA I24,J24,CARRY/24,10,0./ C default C Luxury Level 0 1 2 *3* 4 DATA NDSKIP/0, 24, 73, 199, 365 / Corresponds to p=24 48 97 223 389 C time factor 1 2 3 6 10 on slow workstation C 1 1.5 2 3 5 on fast mainframe C C NOTYET is .TRUE. if no initialization has been performed yet. C Default Initialization by Multiplicative Congruential IF (NOTYET) THEN NOTYET = .FALSE. JSEED = JSDFLT INSEED = JSEED WRITE(6,'(A,I12)') '
... read more »

Fri, 29 Aug 1997 22:05:34 GMT 


che.. #5 / 10

Random Number Generators
Quote: >I have been trying to use the 'rand()' random number generator, >but it doesn't work adequately. Now, I am trying to find other >random number generators that might be of use. I have found >through the man pages that 'random()', 'ran()', and 'lcrans' >are available. It seems that 'lcrans' would be best for my >purposes, but I am programming in FORTRAN and this is a C Library >Function. Any ideas on how I can access this function  or >any other function in the C Library  from a FORTRAN program? >Also, does anyone have any ideas for random number generators that >might be of use to me? I would prefer a pseudorandom number >generator... > Thanks, > Eva S. Simmons
The book "Numerical Recipes in Fortran" by William H. Press et al. contains an entire chapter on random number generators. This will probably have all that you need. Erik Quote: > >***THE SIMMONS FACTOR ****EVA SABRINA SIMMONS******UTAUSTIN GRAD STUDENT*****
> WWW Personal Page: http://ccwf.cc.utexas.edu/~weevey >****************** WATCH IT, OR IT MIGHT ATTACK!! ****************************

Sat, 30 Aug 1997 08:07:58 GMT 


Bob Apthor #6 / 10

Random Number Generators
[snip] w;Also, does anyone have any ideas for random number generators that w;might be of use to me? I would prefer a pseudorandom number w;generator... Eva, Find a copy of Numerical Recipes  it should be available in a good college bookstore or engineering library. It gives some reasonable algorithms including one based on the Data Encryption Standard. Bob  ~ CMPQwk #1.4220~ UNREGISTERED EVALUATION COPY Lousiana  Darwin Takes A Holiday! '[1;35;40m=> Delphi Internet Jet SST v2.009  (C) PBE

Sun, 31 Aug 1997 11:13:46 GMT 


Eva Sabrona Simmo #7 / 10

Random Number Generators
Thanks for the helpful advice about random number generators. I've decided to implement the RAN4 routine from the Numerical Recipes book for FORTRAN. So far, it has worked just fine...:) Enclosed is a copy of this code:  c ********************************************************************** function ran4(idum) c ********************************************************************** c PURPOSE: Returns a uniform random deviate in the range 0.0 to 1.0 c generated by pseudoDES (DESlike) hasing of the 64bit word c (idums,idum), where idums was set by a previous call with negative c idum. Also increments idum. Routine can be used to generate a c random sequence by successive calls, leaving idum unaltered between c calls; or it can randomly access the nth deviate in a sequence by c calling with idum=n. Different sequences are initialized by calls c with differing negative values of idum. c ********************************************************************** integer idum real ran4 integer idums,irword,itemp,jflmsk,jflone,lword real ftemp equivalence (itemp,ftemp) save idums,jflone,jflmsk data idums/0/, jflone/Z'3F800000'/, jflmsk /Z'007FFFFF'/ c ********************************************************************** if (idum.lt.0) then idums=idum idum=1 endif irword=idum lword=idums call psdes(lword,irword) itemp=ior(jflone,iand(jflmsk,irword)) ran4=ftemp1.0 idum=idum+1 c ********************************************************************** return end c ********************************************************************** c ********************************************************************** subroutine psdes(lword,irword) c ********************************************************************** c PURPOSE: "PseudoDES" hashing of the 64bit word (lword,irword). c Both 32bit arguments are turned hashed on all bits. NOTE: This c routine assumes that arbitrary 32bit integers can be added without c overflow. To accomplish this, you may need to compile with a c special directive (e.g., /check=nooverflow for VMS). In other c languages, such as C, one can instead type the integers as "unsigned". c ********************************************************************** integer i,ia,ib,iswap,itmph,itmpl,c1(4),c2(4),irword,lword,NITER parameter (NITER=4) save c1,c2 data c1 /Z'BAA96887',Z'1E17D32C',Z'03BCDC3C',Z'0F33D1B2'/, + c2 /Z'4B0F3B58',Z'E874F0C3',Z'6955C5A6',Z'55A7CA46'/ c ********************************************************************** do 500 i=1,NITER iswap=irword ia=ieor(irword,c1(i)) itmpl=iand(ia,65535) itmph=iand(ishft(ia,16),65535) ib=itmpl**2+not(itmph**2) ia=ior(ishft(ib,i6),iand(ishft(ib,16),65535)) irword=ieor(lword,ieor(c2(i),ia)+itmpl*itmph) lword=iswap 500 continue c ********************************************************************** return end c **********************************************************************  Thanks again for the help!:) Eva S. Simmons  ***THE SIMMONS FACTOR ****EVA SABRINA SIMMONS******UTAUSTIN GRAD STUDENT*****
WWW Personal Page: http://ccwf.cc.utexas.edu/~weevey ****************** WATCH IT, OR IT MIGHT ATTACK!! ****************************

Mon, 01 Sep 1997 06:16:41 GMT 


Rick Venab #8 / 10

Random Number Generators
On 15 Mar 1995 16:16:41 0600 Eva Sabrona Simmons pontificated: Quote: > Thanks for the helpful advice about random number generators. I've > decided to implement the RAN4 routine from the Numerical Recipes book > for FORTRAN. So far, it has worked just fine...:) > Enclosed is a copy of this code:
Oops! Posting Numerical Recipes code is technically a violation of their copyright; they (Flannery, Press, et. al.) do not permit free distribution of the code in any form, especially electronically. The book is quite clear about this; you are expected to pay royalties for each routine you distribute. Besides, better random number routines have been posted to this newsgroup, and/or are available from NETLIB. There is also some work by that venerable sage Knuth, who studied the problem in the same detailed way as he did with sorting.  Rick Venable =====\ = "Eschew Obfuscation" FDA/CBER Biophysics Lab ____/ = Bethesda, MD U.S.A.  \ / = / Not an official statement \

Tue, 02 Sep 1997 12:53:33 GMT 


H. D. Knobl #9 / 10

Random Number Generators
"Numerical Recipes" codes are in many cases good pedagogically, but not the bestimplementations of (semi)numerical computations to be used for for production/research work.In general, copying code from textbooks, even if diskettes are available, should be done with caution.

Tue, 02 Sep 1997 18:18:42 GMT 


Dan P #10 / 10

Random Number Generators
Quote: >"Numerical Recipes" codes are in many cases good pedagogically, but not >the bestimplementations of (semi)numerical computations to be used for >for production/research work.In general, copying code from textbooks, >even if diskettes are available, should be done with caution.
This is particularly true when the book in question is "Numerical Recipes". Those who use or plan to use this book, should have a look at: http://www.lysator.liu.se/c/numericalrecipes.html Dan  Dan Pop CERN, CN Division
Mail: CERN  PPE, Bat. 31 R004, CH1211 Geneve 23, Switzerland

Wed, 03 Sep 1997 06:03:18 GMT 


