Date: Mon, 24 Nov 86 18:06:30 est From: David Krowitz Below the dashed line you will find the bundled files for some random number generators for the Apollos that I wrote along with some test programs. Run /bin/sh to unbundle the files. TEST1, TEST2, TEST3, and TEST4 test the distribution of the generators for uniform, exponential, Lorentzian, and Gaussian numbers against the expected distribution. TEST5 tests the uniform random number generator in triplets (x,y, and color) for visual patterns (I once saw a demo of a supposedly good generator which made a pattern of marching squares when this test was applied!) This package is still under development (on the back burner, really) so there is no documentation files yet. However it's relatively simple to understand from the source files. -- David Krowitz ---------------------------------------------------------- # To unbundle, sh this file echo rand.ftn 1>&2 cat >rand.ftn <<'End of rand.ftn' C***************************************************************************** C***** ***** C***** RAND.FTN ***** C***** ***** C***** Pseudo Random Number Generator Package ***** C***** Version 3 ***** C***** David M. Krowitz August 25, 1986. ***** C***** ***** C***** Copyright (c) 1986 ***** C***** David M. Krowitz ***** C***** Massachusetts Institute of Technology ***** C***** Department of Earth, Atmosheric, and Planetary Sciences ***** C***************************************************************************** C**************************************************************************** C C Routine to seed the random number generator using C the system's time of day. C C**************************************************************************** SUBROUTINE RAND_$INIT INTEGER*4 I,J,K INTEGER*2 CLOCK(3) COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3 LOGICAL INITFLAG INTEGER*4 RANTAB(0:54) INTEGER*4 INDEX1,INDEX2,INDEX3 DATA INITFLAG /.FALSE./ C C Get the 48-bit long time of day from the system and turn it C into a 32-bit long integer. C CALL TIME_$CLOCK (CLOCK) I = INT4(CLOCK(1)) J = INT4(CLOCK(2))*2**4 K = INT4(CLOCK(3))*2**8 RANTAB(0) = I+J+K C C Use a linear congruential psuedo random number generator to C initialize the table for the RAND_$UNIF function. C DO 100 I = 1,54 RANTAB(I) = RANTAB(I-1)*314159812+1 100 CONTINUE INITFLAG = .TRUE. INDEX1 = 0 INDEX2 = 23 INDEX3 = 54 RETURN END C**************************************************************************** C C Uniform Random Number Generator. C Generates a psuedo-random number in the range 0.0 to 1.0 with C a uniform distribution. The first time RAND_$UNIF is called the C argument is used to seed the random number generator if the C random number generator has not already been seeded by a call C to the RAND_$INIT routine. If the seed passed on the first call C is greater than 0, then the seed is used to initialize the C random number generator, otherwise a built in seed number is C used to set up the initial table. Since IEEE standard floating C point numbers have a 24 bit mantissa (one bit being the hidden bit) C we generate random integers in the range 0 to 2**24-1 and then C divide them by 2**24-1 to get a real number in the range 0 to C 1.0, inclusive. C C**************************************************************************** REAL*4 FUNCTION RAND_$UNIF (ISEED) COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3 LOGICAL INITFLAG INTEGER*4 RANTAB(0:54) INTEGER*4 INDEX1,INDEX2,INDEX3 INTEGER*4 ISEED C C If random number generator has not already been seeded C by calling the RAND_$INIT routine, then set up the initial C table for the additive congruential method used here. If C the seed number supplied by the caller is 0, then use our C own seed number. C IF (INITFLAG) GOTO 100 IF (ISEED.GT.0) THEN RANTAB(0) = ISEED ELSE RANTAB(0) = 31415926 ENDIF DO 10 I = 1,54 RANTAB(I) = RANTAB(I-1)*314159812+1 10 CONTINUE INITFLAG = .TRUE. INDEX1 = 0 INDEX2 = 23 INDEX3 = 54 C C Generate the next random integer in the range 0 to 2**24-1 C and then convert it into a real number in the range 0.0 to 1.0 C 100 INDEX1 = MOD(INDEX1+1,55) INDEX2 = MOD(INDEX2+1,55) INDEX3 = MOD(INDEX3+1,55) RANTAB(INDEX1) = RANTAB(INDEX2)+RANTAB(INDEX3) RAND_$UNIF = FLOAT(AND(RANTAB(INDEX1),2**24-1))/FLOAT(2**24-1) RETURN END C**************************************************************************** C C Vector Uniform Random Number Generator. C Same as RAND_$UNIF function except it fills a complete table C full of uniformly distributed pseudo-random numbers in the C range 0.0 to 1.0. This routine avoids the overhead of making C repeated functions calls when you need a lot of random numbers C at once. C C**************************************************************************** SUBROUTINE RAND_$VEC_UNIF (ISEED,TABLE,NUM) COMMON /RANDCOM/ INITFLAG,RANTAB,INDEX1,INDEX2,INDEX3 LOGICAL INITFLAG INTEGER*4 RANTAB(0:54) INTEGER*4 INDEX1,INDEX2,INDEX3 INTEGER*4 ISEED,NUM REAL*4 TABLE(NUM) C C If random number generator has not already been seeded C by calling the RAND_$INIT routine, then set up the initial C table for the additive congruential method used here. If C the seed number supplied by the caller is 0, then use our C own seed number. C IF (INITFLAG) GOTO 100 IF (ISEED.GT.0) THEN RANTAB(0) = ISEED ELSE RANTAB(0) = 31415926 ENDIF DO 10 I = 1,54 RANTAB(I) = RANTAB(I-1)*314159812+1 10 CONTINUE INITFLAG = .TRUE. INDEX1 = 0 INDEX2 = 23 INDEX3 = 54 C C Generate the next random integer in the range 0 to 2**24-1 C and then convert it into a real number in the range 0.0 to 1.0 C 100 DO 200 I = 1,NUM INDEX1 = MOD(INDEX1+1,55) INDEX2 = MOD(INDEX2+1,55) INDEX3 = MOD(INDEX3+1,55) RANTAB(INDEX1) = RANTAB(INDEX2)+RANTAB(INDEX3) TABLE(I) = FLOAT(AND(RANTAB(INDEX1),2**24-1))/FLOAT(2**24-1) 200 CONTINUE RETURN END C**************************************************************************** C C Decaying Exponential Random Number Generator. C Generate random numbers in the range 0.0 to 1.0E+38 (largest C positive single precision real number) with a distribution of C the form: C 1 -X/TAU C --- E , TAU > 0 C TAU C C**************************************************************************** REAL FUNCTION RAND_$EXPONENTIAL(TAU) REAL TAU REAL RAND_$UNIF RAND_$EXPONENTIAL = -TAU*ALOG(RAND_$UNIF(0)) RETURN END C**************************************************************************** C C Vector Decaying Exponential Random Number Generator. C Same as the RAND_$EXPONENTIAL function except that it fills a C complete table of pseudo-random numbers with a distribution C of a decaying exponential. This routine avoids the overhead of C making repeated function calls when you need a lot of values C all at once. C C**************************************************************************** SUBROUTINE RAND_$VEC_EXPONENTIAL(TAU,TABLE,NUM) REAL TAU REAL RAND_$UNIF INTEGER*4 NUM REAL TABLE(NUM) DO 100 I = 1,NUM TABLE(I) = -TAU*ALOG(RAND_$UNIF(0)) 100 CONTINUE RETURN END C**************************************************************************** C C Lorentzian Random Number Generator. C Generate random numbers in the range -1.0E+38 to 1.0E+38 (largest C negative single precision real number to largest positive single C precision real number) with a distribution of the form: C C 1 GAMMA/2 C ---- --------------------- , GAMMA > 0 == WIDTH AT HALF MAX C 2 2 C PI (X-MU) + (GAMMA/2) , MU == MEAN VALUE C C**************************************************************************** REAL FUNCTION RAND_$LORENTZIAN (GAMMA,MU) REAL PI PARAMETER (PI = 3.1415926) REAL GAMMA,MU REAL RAND_$UNIF REAL OLD_GAMMA,OLD_MU REAL CONST1,CONST2 SAVE OLD_GAMMA,OLD_MU SAVE CONST1,CONST2 DATA OLD_GAMMA,OLD_MU/-1.0E38,-1.0E38/ C C Avoid overhead of recalculating the ATAN function C if GAMMA and MU haven't change since the last call. C IF ((GAMMA.EQ.OLD_GAMMA).AND.(MU.EQ.OLD_MU)) GOTO 100 CONST1 = GAMMA/2.0 CONST2 = ATAN(-2.0*MU/GAMMA) OLD_MU = MU OLD_GAMMA = GAMMA 100 RAND_$LORENTZIAN = CONST1*TAN(PI*RAND_$UNIF(0)+CONST2)+MU RETURN END C**************************************************************************** C C Vector Lorentzian Random Number Generator. C Same as the RAND_$LORENTZIAN function except that it fills a C complete table of pseudo-random numbers with a distribution C of a Lorentzian function. This routine avoids the overhead of C making repeated function calls when you need a lot of values C all at once. C C**************************************************************************** SUBROUTINE RAND_$VEC_LORENTZIAN (GAMMA,MU,TABLE,NUM) REAL PI PARAMETER (PI = 3.1415926) REAL GAMMA,MU INTEGER*4 NUM REAL TABLE(NUM) REAL RAND_$UNIF REAL OLD_GAMMA,OLD_MU REAL CONST1,CONST2 SAVE OLD_GAMMA,OLD_MU SAVE CONST1,CONST2 DATA OLD_GAMMA,OLD_MU/-1.0E38,-1.0E38/ C C Avoid overhead of recalculating the ATAN function C if GAMMA and MU haven't change since the last call. C IF ((GAMMA.EQ.OLD_GAMMA).AND.(MU.EQ.OLD_MU)) GOTO 100 CONST1 = GAMMA/2.0 CONST2 = ATAN(-2.0*MU/GAMMA) OLD_MU = MU OLD_GAMMA = GAMMA 100 DO 200 I = 1,NUM TABLE(I) = CONST1*TAN(PI*RAND_$UNIF(0)+CONST2)+MU 200 CONTINUE RETURN END C**************************************************************************** C C Gaussian Random Number Generator. C Generate random numbers in the range -1.0E+38 to 1.0E+38 (largest C negative single precision real number to largest positive single C precision real number) with a distribution of the form: C C 2 C 1 -1 ( X - MU ) C _______________ --- (--------) , SIGMA > 0 == STD. DEVIATION C _____ 2 ( SIGMA ) C / E , MU == MEAN VALUE C \/ 2 PI SIGMA C C**************************************************************************** REAL FUNCTION RAND_$GAUSSIAN (SIGMA,MU) REAL PI REAL DELTA INTEGER SAMPLES PARAMETER (PI = 3.1415926) PARAMETER (DELTA = 6.0E0) PARAMETER (SAMPLES = 200) REAL SIGMA,MU LOGICAL INIT INTEGER*2 INTV(SAMPLES) REAL RAND_$UNIF DOUBLE PRECISION X(SAMPLES),ERF(SAMPLES) DOUBLE PRECISION IERF,DX SAVE INIT SAVE INTV SAVE X,ERF DATA INIT/.FALSE./ C C Define the unit Gaussian function. C GAUSS(X) = 1.0E0/(SQRT(2.0E0*PI))*DEXP(-1.0E0/2.0E0*X**2) C C If not already initialized, tabulate the integral of the C unit gaussian funtion using Simpson's method. C IF (INIT) GOTO 1000 DX = 2.0E0*DELTA/FLOAT(SAMPLES-1) X(1) = -DELTA ERF(1) = 0.0E0 DO 10 I = 2,SAMPLES X(I) = X(I-1)+DX ERF(I) = ERF(I-1)+(GAUSS(X(I-1))+4.0E0*GAUSS(X(I)-DX/2.0E0) & +GAUSS(X(I)))*DX/6.0E0 10 CONTINUE C C Set up a table for looking up what interval of C the inverse ERF function that a particular value C would fall into (since the ERF function was C evaluated at equally spaced intervals, its C inverse does not have equally space points). C INTERVAL =1 DO 30 I = 1,SAMPLES IERF = (I-1)*(1.0E0/FLOAT(SAMPLES-1)) 15 IF (IERF.LT.ERF(INTERVAL)) GOTO 20 INTERVAL = INTERVAL+1 GOTO 15 20 INTV(I) = INTERVAL 30 CONTINUE C C Initialization done. C INIT = .TRUE. C C Look up the inverse ERF funtion for a given random number C between 0.0 and 1.0 and interpolate result. C 1000 R = RAND_$UNIF(0) INTERVAL = INTV(INT2(R*(SAMPLES-1)+1)) 1010 IF (R.LT.ERF(INTERVAL)) GOTO 1100 INTERVAL = INTERVAL+1 GOTO 1010 1100 XX = X(INTERVAL-1)+(X(INTERVAL)-X(INTERVAL-1))/ & (ERF(INTERVAL)-ERF(INTERVAL-1))*(R-ERF(INTERVAL-1)) RAND_$GAUSSIAN = SIGMA*XX+MU RETURN END C**************************************************************************** C C Vector Gaussian Random Number Generator. C Same as the RAND_$GAUSSIAN function except that it fills a C complete table of pseudo-random numbers with a distribution C of a Gaussian function. This routine avoids the overhead of C making repeated function calls when you need a lot of values C all at once. C C**************************************************************************** SUBROUTINE RAND_$VEC_GAUSSIAN (SIGMA,MU,TABLE,NUM) REAL PI REAL DELTA INTEGER SAMPLES PARAMETER (PI = 3.1415926) PARAMETER (DELTA = 6.0E0) PARAMETER (SAMPLES = 200) REAL SIGMA,MU INTEGER*4 NUM REAL TABLE(NUM) LOGICAL INIT INTEGER*2 INTV(SAMPLES) REAL RAND_$UNIF DOUBLE PRECISION X(SAMPLES),ERF(SAMPLES) DOUBLE PRECISION IERF,DX SAVE INIT SAVE INTV SAVE X,ERF DATA INIT/.FALSE./ C C Define the unit Gaussian function. C GAUSS(X) = 1.0E0/(SQRT(2.0E0*PI))*DEXP(-1.0E0/2.0E0*X**2) C C If not already initialized, tabulate the integral of the C unit gaussian funtion using Simpson's method. C IF (INIT) GOTO 1000 DX = 2.0E0*DELTA/FLOAT(SAMPLES-1) X(1) = -DELTA ERF(1) = 0.0E0 DO 10 I = 2,SAMPLES X(I) = X(I-1)+DX ERF(I) = ERF(I-1)+(GAUSS(X(I-1))+4.0E0*GAUSS(X(I)-DX/2.0E0) & +GAUSS(X(I)))*DX/6.0E0 10 CONTINUE C C Set up a table for looking up what interval of C the inverse ERF function that a particular value C would fall into (since the ERF function was C evaluated at equally spaced intervals, its C inverse does not have equally space points). C INTERVAL =1 DO 30 I = 1,SAMPLES IERF = (I-1)*(1.0E0/FLOAT(SAMPLES-1)) 15 IF (IERF.LT.ERF(INTERVAL)) GOTO 20 INTERVAL = INTERVAL+1 GOTO 15 20 INTV(I) = INTERVAL 30 CONTINUE C C Initialization done. C INIT = .TRUE. C C Look up the inverse ERF funtion for a given random number C between 0.0 and 1.0 and interpolate result. C 1000 DO 2000 I = 1,NUM R = RAND_$UNIF(0) INTERVAL = INTV(INT2(R*(SAMPLES-1)+1)) 1010 IF (R.LT.ERF(INTERVAL)) GOTO 1100 INTERVAL = INTERVAL+1 GOTO 1010 1100 XX = X(INTERVAL-1)+(X(INTERVAL)-X(INTERVAL-1))/ & (ERF(INTERVAL)-ERF(INTERVAL-1))*(R-ERF(INTERVAL-1)) TABLE(I) = SIGMA*XX+MU 2000 CONTINUE RETURN END End of rand.ftn echo test1.bld 1>&2 cat >test1.bld <<'End of test1.bld' von ftn test1 -cpu ^1 ftn rand -cpu ^1 bind -b test1 test1.bin rand.bin voff End of test1.bld echo test1.ftn 1>&2 cat >test1.ftn <<'End of test1.ftn' PROGRAM TEST1 C C**************************************************************** C Program to test Steve Ketchum's (spelling?) Random C Number Generator Routine for Uniform Distribution. C**************************************************************** C IMPLICIT INTEGER*2 (A-Z) INTEGER*4 NUMSAMPLES PARAMETER (NUMSAMPLES = 500000) C C Required insert files C %INCLUDE '/SYS/INS/BASE.INS.FTN' %INCLUDE '/SYS/INS/GMR.INS.FTN' C C INTEGER*2 BITMAP_SIZE(2) INTEGER*4 SEGMENT_ID INTEGER*4 ST INTEGER*2 POINT1(2),POINT2(2) INTEGER*4 I,J INTEGER*2 DISTRIBUTION(2,1000) LOGICAL COLOR_NODE REAL RAND_$UNIF C DATA BITMAP_SIZE / 1024, 1024/ C C Start GMR package, create file, and create segment. C CALL GM_$INIT(GM_$DIRECT,INT2(1),BITMAP_SIZE,INT2(8),ST) CALL GM_$FILE_CREATE('test1.gmr',INT2(9),GM_$OVERWRITE, &GM_$1W,FILE_ID,ST) CALL GM_$SEGMENT_CREATE('easy',INT2(4),SEGMENT_ID,ST) C C Find out if we have a color node. C CALL GM_$INQ_CONFIG(CONFIG,ST) IF ((CONFIG.EQ.GM_$BW_800x1024).OR. &(CONFIG.EQ.GM_$BW_1024x800)) THEN COLOR_NODE=.FALSE. ELSE COLOR_NODE=.TRUE. ENDIF C C Init the random number generator and the distribution function. C CALL RAND_$INIT DO 5 J = 1,1000 DISTRIBUTION(1,J)=J DISTRIBUTION(2,J)=0 5 CONTINUE C C Generate some random numbers and check their C distribution function. C DO 10 I = 1,NUMSAMPLES J=1000*RAND_$UNIF(0)+1.0 IF ((J.LT.1).OR.(J.GT.1000)) THEN WRITE (*,11) I,J 11 FORMAT (1X,'ON LOOP:',I6,' VALUE OF J IS:',I6) ELSE DISTRIBUTION(2,J)=DISTRIBUTION(2,J)+1 ENDIF 10 CONTINUE C C Plot the distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(3,ST) CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Plot an axis with tick marks. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =0 DISTRIBUTION (2,2) =1000 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =1000 DISTRIBUTION (2,2) =0 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DO 15 I=0,10 DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =I*100 DISTRIBUTION (1,2) =-5 DISTRIBUTION (2,2) =I*100 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 15 CONTINUE DO 20 I=0,10 DISTRIBUTION (1,1) =I*100 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =I*100 DISTRIBUTION (2,2) =-10 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 20 CONTINUE C C Draw the ideal distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) CALL GM_$DRAW_STYLE (GM_$DOTTED,INT2(10),CMPLX(0.0,0.0), &INT2(2),ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =NUMSAMPLES/1000 DISTRIBUTION (1,2) =1000 DISTRIBUTION (2,2) =NUMSAMPLES/1000 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Display the file. C CALL GM_$DISPLAY_FILE(ST) C C Leave it on screen long enough to see it. C DO 100 I = 1, 50000000 100 CONTINUE C C Clean up. C CALL GM_$SEGMENT_CLOSE(.TRUE.,ST) CALL GM_$FILE_CLOSE(.TRUE.,ST) CALL GM_$TERMINATE(ST) END End of test1.ftn echo test2.bld 1>&2 cat >test2.bld <<'End of test2.bld' von ftn test2 -cpu ^1 ftn rand -cpu ^1 bind -b test2 test2.bin rand.bin voff End of test2.bld echo test2.ftn 1>&2 cat >test2.ftn <<'End of test2.ftn' PROGRAM TEST2 C C**************************************************************** C Program to test Random Number Generator C Routine for Exponential Distribution. C**************************************************************** C IMPLICIT INTEGER*2 (A-Z) C C Required insert files C %INCLUDE '/SYS/INS/BASE.INS.FTN' %INCLUDE '/SYS/INS/GMR.INS.FTN' C C INTEGER*2 BITMAP_SIZE(2) INTEGER*4 SEGMENT_ID INTEGER*4 ST INTEGER*2 POINT1(2),POINT2(2) INTEGER*4 I,J INTEGER*2 DISTRIBUTION(2,1000) LOGICAL COLOR_NODE REAL RAND_$UNIF,RAND_$EXPONENTIAL REAL TAU C PARAMETER (NUMSAMPLES = 100000) PARAMETER (TAU = 200.0) C DATA BITMAP_SIZE / 1024, 1024/ C C Start GMR package, create file, and create segment. C CALL GM_$INIT(GM_$DIRECT,INT2(1),BITMAP_SIZE,INT2(8),ST) CALL GM_$FILE_CREATE('test2.gmr',INT2(9),GM_$OVERWRITE, &GM_$1W,FILE_ID,ST) CALL GM_$SEGMENT_CREATE('easy',INT2(4),SEGMENT_ID,ST) C C Find out if we have a color node. C CALL GM_$INQ_CONFIG(CONFIG,ST) IF ((CONFIG.EQ.GM_$BW_800x1024).OR. &(CONFIG.EQ.GM_$BW_1024x800)) THEN COLOR_NODE=.FALSE. ELSE COLOR_NODE=.TRUE. ENDIF C C Init the random number generator and the distribution function. C CALL RAND_$INIT DO 5 J = 1,1000 DISTRIBUTION(1,J)=J DISTRIBUTION(2,J)=0 5 CONTINUE C C Generate some random numbers and check their C distribution function. C DO 10 I = 1,NUMSAMPLES J=RAND_$EXPONENTIAL(TAU)+1.0 IF ((J.LT.1).OR.(J.GT.1000)) THEN WRITE (*,11) I,J 11 FORMAT (1X,'ON LOOP:',I6,' VALUE OF J IS:',I6) ELSE DISTRIBUTION(2,J)=DISTRIBUTION(2,J)+1 ENDIF 10 CONTINUE C C Plot the distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(3,ST) CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Plot an axis with tick marks. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =0 DISTRIBUTION (2,2) =1000 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =1000 DISTRIBUTION (2,2) =0 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DO 15 I=0,10 DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =100*I DISTRIBUTION (1,2) =-5 DISTRIBUTION (2,2) =100*I CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 15 CONTINUE DO 20 I=0,10 DISTRIBUTION (1,1) =I*100 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =I*100 DISTRIBUTION (2,2) =-10 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 20 CONTINUE C C Draw the ideal distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) CALL GM_$DRAW_STYLE (GM_$DOTTED,INT2(10),CMPLX(0.0,0.0), &INT2(2),ST) DO 30 I = 1,1000 DISTRIBUTION (1,I) =I DISTRIBUTION (2,I) =NUMSAMPLES/TAU*EXP(-I/TAU) 30 CONTINUE CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Display the file. C CALL GM_$DISPLAY_FILE(ST) C C Leave it on screen long enough to see it. C DO 100 I = 1, 50000000 100 CONTINUE C C Clean up. C CALL GM_$SEGMENT_CLOSE(.TRUE.,ST) CALL GM_$FILE_CLOSE(.TRUE.,ST) CALL GM_$TERMINATE(ST) END End of test2.ftn echo test3.bld 1>&2 cat >test3.bld <<'End of test3.bld' von ftn test3 -cpu ^1 ftn rand -cpu ^1 bind -b test3 test3.bin rand.bin voff End of test3.bld echo test3.ftn 1>&2 cat >test3.ftn <<'End of test3.ftn' PROGRAM TEST3 C C**************************************************************** C Program to test Random Number Generator C Routine for Lorentzian Distribution. C**************************************************************** C IMPLICIT INTEGER*2 (A-Z) C C Required insert files C %INCLUDE '/SYS/INS/BASE.INS.FTN' %INCLUDE '/SYS/INS/GMR.INS.FTN' C C INTEGER*2 BITMAP_SIZE(2) INTEGER*4 SEGMENT_ID INTEGER*4 ST INTEGER*2 POINT1(2),POINT2(2) INTEGER*4 I,J INTEGER*2 DISTRIBUTION(2,1000) LOGICAL COLOR_NODE REAL RAND_$UNIF,RAND_$LORENTZIAN REAL PI REAL MEAN,HALFWIDTH C PARAMETER (PI = 3.1415926) PARAMETER (NUMSAMPLES = 100000) PARAMETER (MEAN = 500.0) PARAMETER (HALFWIDTH = 100.0) C DATA BITMAP_SIZE / 1024, 1024/ C C Start GMR package, create file, and create segment. C CALL GM_$INIT(GM_$DIRECT,INT2(1),BITMAP_SIZE,INT2(8),ST) CALL GM_$FILE_CREATE('test3.gmr',INT2(9),GM_$OVERWRITE, &GM_$1W,FILE_ID,ST) CALL GM_$SEGMENT_CREATE('easy',INT2(4),SEGMENT_ID,ST) C C Find out if we have a color node. C CALL GM_$INQ_CONFIG(CONFIG,ST) IF ((CONFIG.EQ.GM_$BW_800x1024).OR. &(CONFIG.EQ.GM_$BW_1024x800)) THEN COLOR_NODE=.FALSE. ELSE COLOR_NODE=.TRUE. ENDIF C C Init the random number generator and the distribution function. C CALL RAND_$INIT DO 5 J = 1,1000 DISTRIBUTION(1,J)=J DISTRIBUTION(2,J)=0 5 CONTINUE C C Generate some random numbers and check their C distribution function. C DO 10 I = 1,NUMSAMPLES J=RAND_$LORENTZIAN(HALFWIDTH,MEAN) IF ((J.LT.1).OR.(J.GT.1000)) THEN WRITE (*,11) I,J 11 FORMAT (1X,'ON LOOP:',I6,' VALUE OF J IS:',I6) ELSE DISTRIBUTION(2,J)=DISTRIBUTION(2,J)+1 ENDIF 10 CONTINUE C C Plot the distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(3,ST) CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Plot an axis with tick marks. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =0 DISTRIBUTION (2,2) =1000 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =1000 DISTRIBUTION (2,2) =0 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DO 15 I=0,10 DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =100*I DISTRIBUTION (1,2) =-5 DISTRIBUTION (2,2) =100*I CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 15 CONTINUE DO 20 I=0,10 DISTRIBUTION (1,1) =I*100 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =I*100 DISTRIBUTION (2,2) =-10 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 20 CONTINUE C C Draw the ideal distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) CALL GM_$DRAW_STYLE (GM_$DOTTED,INT2(10),CMPLX(0.0,0.0), &INT2(2),ST) DO 30 I = 1,1000 DISTRIBUTION (1,I) =I DISTRIBUTION (2,I) =NUMSAMPLES*HALFWIDTH/(2*PI*((I-MEAN)**2 & +(HALFWIDTH/2.0)**2)) 30 CONTINUE CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Display the file. C CALL GM_$DISPLAY_FILE(ST) C C Leave it on screen long enough to see it. C DO 100 I = 1, 50000000 100 CONTINUE C C Clean up. C CALL GM_$SEGMENT_CLOSE(.TRUE.,ST) CALL GM_$FILE_CLOSE(.TRUE.,ST) CALL GM_$TERMINATE(ST) END End of test3.ftn echo test4.bld 1>&2 cat >test4.bld <<'End of test4.bld' von ftn test4 -cpu ^1 ftn rand -cpu ^1 bind -b test4 test4.bin rand.bin voff End of test4.bld echo test4.ftn 1>&2 cat >test4.ftn <<'End of test4.ftn' PROGRAM TEST4 C C**************************************************************** C Program to test Random Number Generator C Routine for Gaussian Distribution. C**************************************************************** C IMPLICIT INTEGER*2 (A-Z) C C Required insert files C %INCLUDE '/SYS/INS/BASE.INS.FTN' %INCLUDE '/SYS/INS/GMR.INS.FTN' C C INTEGER*2 BITMAP_SIZE(2) INTEGER*4 SEGMENT_ID INTEGER*4 ST INTEGER*2 POINT1(2),POINT2(2) INTEGER*4 I,J INTEGER*2 DISTRIBUTION(2,1000) LOGICAL COLOR_NODE REAL RAND_$UNIF,RAND_$GAUSSIAN REAL PI REAL MEAN,SIGMA C PARAMETER (PI = 3.1415926) PARAMETER (NUMSAMPLES = 300000) PARAMETER (MEAN = 500.0) PARAMETER (SIGMA = 100.0) C DATA BITMAP_SIZE / 1024, 1024/ C C Start GMR package, create file, and create segment. C CALL GM_$INIT(GM_$DIRECT,INT2(1),BITMAP_SIZE,INT2(8),ST) CALL GM_$FILE_CREATE('test4.gmr',INT2(9),GM_$OVERWRITE, &GM_$1W,FILE_ID,ST) CALL GM_$SEGMENT_CREATE('easy',INT2(4),SEGMENT_ID,ST) C C Find out if we have a color node. C CALL GM_$INQ_CONFIG(CONFIG,ST) IF ((CONFIG.EQ.GM_$BW_800x1024).OR. &(CONFIG.EQ.GM_$BW_1024x800)) THEN COLOR_NODE=.FALSE. ELSE COLOR_NODE=.TRUE. ENDIF C C Init the random number generator and the distribution function. C CALL RAND_$INIT DO 5 J = 1,1000 DISTRIBUTION(1,J)=J DISTRIBUTION(2,J)=0 5 CONTINUE C C Generate some random numbers and check their C distribution function. C DO 10 I = 1,NUMSAMPLES J=RAND_$GAUSSIAN(SIGMA,MEAN) IF ((J.LT.1).OR.(J.GT.1000)) THEN WRITE (*,11) I,J 11 FORMAT (1X,'ON LOOP:',I6,' VALUE OF J IS:',I6) ELSE DISTRIBUTION(2,J)=DISTRIBUTION(2,J)+1 ENDIF 10 CONTINUE C C Plot the distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(3,ST) CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Plot an axis with tick marks. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =0 DISTRIBUTION (2,2) =1000 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =1000 DISTRIBUTION (2,2) =0 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) DO 15 I=0,10 DISTRIBUTION (1,1) =0 DISTRIBUTION (2,1) =100*I DISTRIBUTION (1,2) =-5 DISTRIBUTION (2,2) =100*I CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 15 CONTINUE DO 20 I=0,10 DISTRIBUTION (1,1) =I*100 DISTRIBUTION (2,1) =0 DISTRIBUTION (1,2) =I*100 DISTRIBUTION (2,2) =-10 CALL GM_$POLYLINE_2D16(INT2(2),DISTRIBUTION,.FALSE.,.FALSE.,ST) 20 CONTINUE C C Draw the ideal distribution. C IF (COLOR_NODE) CALL GM_$DRAW_VALUE(1,ST) CALL GM_$DRAW_STYLE (GM_$DOTTED,INT2(10),CMPLX(0.0,0.0), &INT2(2),ST) DO 30 I = 1,1000 DISTRIBUTION (1,I) =I DISTRIBUTION (2,I) =NUMSAMPLES/(SIGMA*SQRT(2.0*PI))* &EXP(-1.0/2.0*(((I-MEAN)/SIGMA)**2)) 30 CONTINUE CALL GM_$POLYLINE_2D16(INT2(1000),DISTRIBUTION,.FALSE.,.FALSE.,ST) C C Display the file. C CALL GM_$DISPLAY_FILE(ST) C C Leave it on screen long enough to see it. C DO 100 I = 1, 50000000 100 CONTINUE C C Clean up. C CALL GM_$SEGMENT_CLOSE(.TRUE.,ST) CALL GM_$FILE_CLOSE(.TRUE.,ST) CALL GM_$TERMINATE(ST) END End of test4.ftn echo test5.bld 1>&2 cat >test5.bld <<'End of test5.bld' von pas test5.pas -cpu ^1 ftn rand.ftn -cpu ^1 bind -b test5 test5.bin rand.bin voff End of test5.bld echo test5.pas 1>&2 cat >test5.pas <<'End of test5.pas' PROGRAM TEST5; %NOLIST; %INCLUDE '/sys/ins/base.ins.pas'; %INCLUDE '/sys/ins/gpr.ins.pas'; %LIST; VAR screen_window: GPR_$WINDOW_T; {Color bitmap GPR window size} dest_window: GPR_$WINDOW_T; {GPR window for writing pixels to grey_scale bitmap} screen_bitmap: GPR_$BITMAP_DESC_T; {Copy of color bitmap to display on screen while working} hi_plane: GPR_$PLANE_T; {Highest plane used in bitmap} screen_plane: GPR_$PLANE_T; {Highest plane used in screen's bitmap} pixel: GPR_$PIXEL_VALUE_T; status: status_$t; {Status returned by GPR calls} i,j: linteger; {Counters for bitmap coordinates} ii,jj: pinteger; {Counters for copying grey-scale pixels into the bitmap} FUNCTION rand_$unif (IN seed: integer32): REAL;EXTERN; BEGIN {Set up initial screen bitmap size and origin for bitmap displayed on screen and then init GPR.} screen_window.window_base.x_coord := 0; screen_window.window_base.y_coord := 0; screen_window.window_size.x_size := 1024; screen_window.window_size.y_size := 1024; screen_plane := 7; GPR_$INIT (GPR_$BORROW,1,screen_window.window_size,screen_plane,screen_bitmap,status); dest_window.window_size.x_size := 1; dest_window.window_size.y_size := 1; FOR i := 0 TO 10000000 DO BEGIN dest_window.window_base.x_coord := ROUND(screen_window.window_size.x_size *rand_$unif(0)); dest_window.window_base.y_coord := ROUND(screen_window.window_size.y_size *rand_$unif(0)); pixel := ROUND(15*rand_$unif(0)); GPR_$WRITE_PIXELS (pixel,dest_window,status); END; GPR_$TERMINATE (FALSE,status); end. End of test5.pas