Date: Mon, 24 Nov 86 18:06:30 est
From: David Krowitz <mit-erl!mit-kermit!tigger!krowitz@EDDIE.MIT.EDU>
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