      SUBROUTINE NBRAN(N,P,NPAR,ISTART,X)
C
C     PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N
C              FROM THE NEGATIVE BINOMIAL DISTRIBUTION
C              WITH SINGLE PRECISION 'BERNOULLI PROBABILITY'
C              PARAMETER = P, 
C              AND INTEGER 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS'
C              PARAMETER = NPAR.
C              THE NEGATIVE BINOMIAL DISTRIBUTION USED
C              HEREIN HAS MEAN = NPAR*(1-P)/P
C              AND STANDARD DEVIATION = SQRT(NPAR*(1-P)/(P*P))).
C              THIS DISTRIBUTION IS DEFINED FOR
C              ALL NON-NEGATIVE INTEGER X--X = 0, 1, 2, ... .
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = C(NPAR+X-1,NPAR) * P**NPAR * (1-P)**X.
C              WHERE C(NPAR+X-1,NPAR) IS THE COMBINATORIAL FUNCTION
C              EQUALING THE NUMBER OF COMBINATIONS OF NPAR+X-1 ITEMS
C              TAKEN NPAR AT A TIME.
C              THE NEGATIVE BINOMIAL DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF FAILURES
C              BEFORE OBTAINING NPAR SUCCESSES IN AN
C              INDEFINITE SEQUENCE OF BERNOULLI (0,1)
C              TRIALS WHERE THE PROBABILITY OF SUCCESS
C              IN A SINGLE TRIAL = P.
C     INPUT  ARGUMENTS--N      = THE DESIRED INTEGER NUMBER 
C                                OF RANDOM NUMBERS TO BE
C                                GENERATED.
C                     --P      = THE SINGLE PRECISION VALUE 
C                                OF THE 'BERNOULLI PROBABILITY'
C                                PARAMETER FOR THE NEGATIVE BINOMIAL
C                                DISTRIBUTION.
C                                P SHOULD BE BETWEEN
C                                0.0 (EXCLUSIVELY) AND
C                                1.0 (EXCLUSIVELY).
C                     --NPAR   = THE INTEGER VALUE
C                                OF THE 'NUMBER OF SUCCESSES
C                                IN BERNOULLI TRIALS' PARAMETER.
C                                NPAR SHOULD BE A POSITIVE INTEGER.
C                     --ISTART = AN INTEGER FLAG CODE WHICH 
C                                (IF SET TO 0) WILL START THE
C                                GENERATOR OVER AND HENCE
C                                PRODUCE THE SAME RANDOM SAMPLE
C                                OVER AND OVER AGAIN
C                                UPON SUCCESSIVE CALLS TO
C                                THIS SUBROUTINE WITHIN A RUN; OR
C                                (IF SET TO SOME INTEGER
C                                VALUE NOT EQUAL TO 0,
C                                LIKE, SAY, 1) WILL ALLOW
C                                THE GENERATOR TO CONTINUE
C                                FROM WHERE IT STOPPED
C                                AND HENCE PRODUCE DIFFERENT
C                                RANDOM SAMPLES UPON
C                                SUCCESSIVE CALLS TO
C                                THIS SUBROUTINE WITHIN A RUN.
C     OUTPUT ARGUMENTS--X      = A SINGLE PRECISION VECTOR
C                                (OF DIMENSION AT LEAST N)
C                                INTO WHICH THE GENERATED
C                                RANDOM SAMPLE WILL BE PLACED.
C     OUTPUT--A RANDOM SAMPLE OF SIZE N 
C             FROM THE NEGATIVE BINOMIAL DISTRIBUTION
C             WITH 'BERNOULLI PROBABILITY' PARAMETER = P
C             AND 'NUMBER OF SUCCESSES IN BERNOULLI TRIALS' 
C             PARAMETER = NPAR.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C                 --NPAR SHOULD BE A POSITIVE INTEGER.
C     OTHER DATAPAC   SUBROUTINES NEEDED--UNIRAN, BINRAN, GEORAN.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN. 
C     COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT
C              FROM THIS DISCRETE RANDOM NUMBER
C              GENERATOR MUST NECESSARILY BE A
C              SEQUENCE OF ***INTEGER*** VALUES,
C              THE OUTPUT VECTOR X IS SINGLE
C              PRECISION IN MODE.
C              X HAS BEEN SPECIFIED AS SINGLE
C              PRECISION SO AS TO CONFORM WITH THE DATAPAC
C              CONVENTION THAT ALL OUTPUT VECTORS FROM ALL
C              DATAPAC SUBROUTINES ARE SINGLE PRECISION.
C              THIS CONVENTION IS BASED ON THE BELIEF THAT
C              1) A MIXTURE OF MODES (FLOATING POINT
C              VERSUS INTEGER) IS INCONSISTENT AND
C              AN UNNECESSARY COMPLICATION
C              IN A DATA ANALYSIS; AND
C              2) FLOATING POINT MACHINE ARITHMETIC
C              (AS OPPOSED TO INTEGER ARITHMETIC) 
C              IS THE MORE NATURAL MODE FOR DOING 
C              DATA ANALYSIS. 
C     REFERENCES--HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGE 95.
C               --JOHNSON AND KOTZ, DISCRETE
C                 DISTRIBUTIONS, 1969, PAGES 122-142.
C               --FELLER, AN INTRODUCTION TO PROBABILITY
C                 THEORY AND ITS APPLICATIONS, VOLUME 1,
C                 EDITION 2, 1957, PAGES 155-157, 210.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 929.
C               --KENDALL AND STUART, THE ADVANCED THEORY OF
C                 STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 130-131.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING LABORATORY (205.03)
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE:  301-921-2315
C     ORIGINAL VERSION--NOVEMBER  1975. 
C
C---------------------------------------------------------------------
C
      DIMENSION X(1)
C
      IPR=6
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)GOTO50
      IF(P.LE.0.0.OR.P.GE.1.0)GOTO55
      IF(NPAR.LT.1)GOTO60
      GOTO90
   50 WRITE(IPR, 5) 
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,11) 
      WRITE(IPR,46)P
      RETURN
   60 WRITE(IPR,25) 
      WRITE(IPR,47)NPAR
      RETURN
   90 CONTINUE
    5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 BINRAN SUBROUTINE IS NON-POSITIVE *****)
   11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 BINRAN SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 BINRAN SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C-----START POINT-----------------------------------------------------
C
      CALL UNIRAN(1,ISTART,G) 
C
C     CHECK ON THE MAGNITUDE OF P,
C     AND BRANCH TO THE FASTER
C     GENERATION METHOD ACCORDINGLY.
C
      IF(P.LT.0.1)GOTO450
C
C     IF P IS MODERATE OR LARGE,
C     GENERATE N NEGATIVE BINOMIAL NUMBERS
C     USING THE FACT THAT THE 
C     WAITING TIME FOR NPAR SUCCESSES IN
C     BERNOULLI TRIALS HAS A
C     NEGATIVE BINOMIAL DISTRIBUTION.
C
      DO100I=1,N
      ISUM=0
      J=1 
  150 CALL BINRAN(1,P,1,1,B)
      IB=B+0.5
      ISUM=ISUM+IB
      IF(ISUM.EQ.NPAR)GOTO250 
      J=J+1
      GOTO150
  250 X(I)=J
  100 CONTINUE
      RETURN
C
C     IF P IS SMALL,
C     GENERATE N NEGATIVE BINOMIAL NUMBERS
C     BY USING THE FACT THAT THE SUM
C     OF GEOMETRIC VARIATES IS A
C     NEGATIVE BINOMIAL VARIATE.
C
  450 DO500I=1,N
      ISUM=0
      DO600J=1,NPAR 
      CALL GEORAN(1,P,1,G)
      IG=G+0.5
      ISUM=ISUM+IG
  600 CONTINUE
      X(I)=ISUM
  500 CONTINUE
      RETURN
C
      END 
      SUBROUTINE BINRAN(N,P,NPAR,ISEED,X)
C
C     PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N
C              FROM THE BINOMIAL DISTRIBUTION
C              WITH SINGLE PRECISION 'BERNOULLI PROBABILITY'
C              PARAMETER = P,
C              AND INTEGER 'NUMBER OF BERNOULLI TRIALS'
C              PARAMETER = NPAR.
C              THE BINOMIAL DISTRIBUTION USED
C              HEREIN HAS MEAN = NPAR*P
C              AND STANDARD DEVIATION = SQRT(NPAR*P*(1-P)).
C              THIS DISTRIBUTION IS DEFINED FOR ALL
C              DISCRETE INTEGER X BETWEEN 0 (INCLUSIVELY)
C              AND NPAR (INCLUSIVELY).
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = C(NPAR,X) * P**X * (1-P)**(NPAR-X).
C              WHERE C(NPAR,X) IS THE COMBINATORIAL FUNCTION
C              EQUALING THE NUMBER OF COMBINATIONS OF NPAR ITEMS
C              TAKEN X AT A TIME.
C              THE BINOMIAL DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF
C              SUCCESSES IN NPAR BERNOULLI (0,1)
C              TRIALS WHERE THE PROBABILITY OF SUCCESS
C              IN A SINGLE TRIAL = P.
C     INPUT  ARGUMENTS--N      = THE DESIRED INTEGER NUMBER
C                                OF RANDOM NUMBERS TO BE
C                                GENERATED.
C                     --P      = THE SINGLE PRECISION VALUE
C                                OF THE 'BERNOULLI PROBABILITY'
C                                PARAMETER FOR THE BINOMIAL
C                                DISTRIBUTION.
C                                P SHOULD BE BETWEEN
C                                0.0 (EXCLUSIVELY) AND
C                                1.0 (EXCLUSIVELY).
C                     --NPAR   = THE INTEGER VALUE
C                                OF THE 'NUMBER OF BERNOULLI TRIALS'
C                                PARAMETER.
C                                NPAR SHOULD BE A POSITIVE INTEGER.
C     OUTPUT ARGUMENTS--X      = A SINGLE PRECISION VECTOR
C                                (OF DIMENSION AT LEAST N)
C                                INTO WHICH THE GENERATED
C                                RANDOM SAMPLE WILL BE PLACED.
C     OUTPUT--A RANDOM SAMPLE OF SIZE N
C             FROM THE BINOMIAL DISTRIBUTION
C             WITH 'BERNOULLI PROBABILITY' PARAMETER = P
C             AND 'NUMBER OF BERNOULLI TRIALS' PARAMETER = NPAR.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C                 --NPAR SHOULD BE A POSITIVE INTEGER.
C     OTHER DATAPAC   SUBROUTINES NEEDED--UNIRAN, GEORAN.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT
C              FROM THIS DISCRETE RANDOM NUMBER
C              GENERATOR MUST NECESSARILY BE A
C              SEQUENCE OF ***INTEGER*** VALUES,
C              THE OUTPUT VECTOR X IS SINGLE
C              PRECISION IN MODE.
C              X HAS BEEN SPECIFIED AS SINGLE
C              PRECISION SO AS TO CONFORM WITH THE DATAPAC
C              CONVENTION THAT ALL OUTPUT VECTORS FROM ALL
C              DATAPAC SUBROUTINES ARE SINGLE PRECISION.
C              THIS CONVENTION IS BASED ON THE BELIEF THAT
C              1) A MIXTURE OF MODES (FLOATING POINT
C              VERSUS INTEGER) IS INCONSISTENT AND
C              AN UNNECESSARY COMPLICATION
C              IN A DATA ANALYSIS; AND
C              2) FLOATING POINT MACHINE ARITHMETIC
C              (AS OPPOSED TO INTEGER ARITHMETIC)
C              IS THE MORE NATURAL MODE FOR DOING
C              DATA ANALYSIS.
C     REFERENCES--JOHNSON AND KOTZ, DISCRETE
C                 DISTRIBUTIONS, 1969, PAGES 50-86.
C               --HASTINGS AND PEACOCK, STATISTICAL
C                 DISTRIBUTIONS--A HANDBOOK FOR
C                 STUDENTS AND PRACTITIONERS, 1975,
C                 PAGE 41.
C               --FELLER, AN INTRODUCTION TO PROBABILITY
C                 THEORY AND ITS APPLICATIONS, VOLUME 1,
C                 EDITION 2, 1957, PAGES 135-142.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 929.
C               --KENDALL AND STUART, THE ADVANCED THEORY OF
C                 STATISTICS, VOLUME 1, EDITION 2, 1963, PAGES 120-125.
C               --MOOD AND GRABLE, INTRODUCTION TO THE THEORY
C                 OF STATISTICS, EDITION 2, 1963, PAGES 64-69.
C               --TOCHER, THE ART OF SIMULATION,
C                 1963, PAGES 39-40.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE--301-921-3651
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL BUREAU OF STANDARDS.
C           THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED,
C           MODIFIED, OR OTHERWISE USED IN A CONTEXT
C           OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM.
C     LANGUAGE--ANSI FORTRAN (1966)
C               EXCEPTION--HOLLERITH STRINGS IN FORMAT STATEMENTS
C                          DENOTED BY QUOTES RATHER THAN NH.
C     VERSION NUMBER--82/7
C     ORIGINAL VERSION--NOVEMBER  1975.
C     UPDATED         --DECEMBER  1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
C
C---------------------------------------------------------------------
C
CCCCC CHARACTER*4 IFEEDB
CCCCC CHARACTER*4 IPRINT
C
CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW
CCCCC COMMON /PRINT/IFEEDB,IPRINT
C
      IPR=6
C
C-----START POINT-----------------------------------------------------
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)GOTO50
      IF(P.LE.0.0.OR.P.GE.1.0)GOTO55
      IF(NPAR.LT.1)GOTO60
      GOTO90
   50 WRITE(IPR, 5)
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,11)
      WRITE(IPR,46)P
      RETURN
   60 WRITE(IPR,25)
      WRITE(IPR,47)NPAR
      RETURN
   90 CONTINUE
    5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 BINRAN SUBROUTINE IS NON-POSITIVE *****)
   11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 BINRAN SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   25 FORMAT(1H , 91H***** FATAL ERROR--THE THIRD  INPUT ARGUMENT TO THE
     1 BINRAN SUBROUTINE IS NON-POSITIVE *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C     CHECK ON THE MAGNITUDE OF P,
C     AND BRANCH TO THE FASTER
C     GENERATION METHOD ACCORDINGLY.
C
      IF(P.LT.0.1)GOTO450
C
C     IF P IS MODERATE OR LARGE,
C     GENERATE N BINOMIAL RANDOM NUMBERS
C     USING THE REJECTION METHOD.
C
      DO100I=1,N
      ISUM=0
      DO200J=1,NPAR
      CALL UNIRAN(1,ISEED,U)
      IF(U.LE.P)ISUM=ISUM+1
  200 CONTINUE
      X(I)=ISUM
  100 CONTINUE
      RETURN
C
C     IF P IS SMALL,
C     GENERATE N BINOMIAL NUMBERS
C     USING THE FACT THAT THE
C     WAITING TIME FOR 1 SUCCESS IN
C     BERNOULLI TRIALS HAS A
C     GEOMETRIC DISTRIBUTION.
C
  450 DO500I=1,N
      ISUM=0
      J=1
  550 CALL GEORAN(1,P,ISEED,G)
      IG=G+0.5
      ISUM=ISUM+IG+1
      IF(ISUM.GT.NPAR)GOTO650
      J=J+1
      GOTO550
  650 X(I)=J-1
  500 CONTINUE
      RETURN
C
      END
      SUBROUTINE UNIRAN(N,ISEED,X)
C
C     PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N
C              FROM THE UNIFORM (RECTANGULAR)
C              DISTRIBUTION ON THE UNIT INTERVAL (0,1).
C              THIS DISTRIBUTION HAS MEAN = 0.5
C              AND STANDARD DEVIATION = SQRT(1/12) = 0.28867513.
C              THIS DISTRIBUTION HAS THE PROBABILITY
C              DENSITY FUNCTION F(X) = 1.
C
C     INPUT  ARGUMENTS--N      = THE DESIRED INTEGER NUMBER 
C                                OF RANDOM NUMBERS TO BE
C                                GENERATED.
C                     --ISEED  = AN INTEGER ISEED VALUE
C     OUTPUT ARGUMENTS--X      = A SINGLE PRECISION VECTOR
C                                (OF DIMENSION AT LEAST N)
C                                INTO WHICH THE GENERATED
C                                RANDOM SAMPLE WILL BE PLACED.
C     OUTPUT--A RANDOM SAMPLE OF SIZE N 
C             FROM THE RECTANGULAR DISTRIBUTION ON (0,1).
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS. 
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C     OTHER           SUBROUTINES NEEDED--NONE.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--NONE.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C
C     ALGORITHM--FIBONACCI GENERATOR
C                AS DEFINED BY GEORGE MARSAGLIA.
C
C     NOTE--THIS GENERATOR IS TRANSPORTABLE.
C           IT IS NOT MACHINE-INDEPENDENT
C           IN THE SENSE THAT FOR A GIVEN VALUE
C           OF THE INPUT SEED ISEED AND FOR A GIVEN VALUE
C           OF MDIG (TO BE DEFINED BELOW),
C           THE SAME SEQUENCE OF UNIRFORM RANDOM
C           NUMBERS WILL RESULT ON DIFFERENT COMPUTERS
C           (VAX, PRIME, PERKIN-ELMER, IBM, UNIVAC, HONEYWELL, ETC.)
C
C     NOTE--IF MDIG = 32 AND IF ISEED = 305,
C           THEN THE OUTPUT FROM THIS GENERATOR SHOULD BE AS FOLLOWS--
C           THE FIRST      NUMBER TO RESULT IS .4771580...
C           THE SECOND     NUMBER TO RESULT IS .4219293...
C           THE THIRD      NUMBER TO RESULT IS .6646181...
C           ...
C           THE THOUSANDTH NUMBER TO RESULT IS .2036834...
C
C     NOTE--IF MDIG = 16 AND IF ISEED = 305,
C           THEN THE OUTPUT FROM THIS GENERATOR SHOULD BE AS FOLLOWS--
C           THE FIRST      NUMBER TO RESULT IS .027832881...
C           THE SECOND     NUMBER TO RESULT IS .56102176... 
C           THE THIRD      NUMBER TO RESULT IS .41456343... 
C           ...
C           THE THOUSANDTH NUMBER TO RESULT IS .19797357... 
C
C     NOTE--IT IS RECOMMENDED THAT UPON 
C           IMPLEMENTATION OF DATAPLOT, THE OUTPUT
C           FROM UNIRAN BE CHECKED FOR AGREEMENT
C           WITH THE ABOVE SAMPLE OUTPUT.
C           ALSO, THERE ARE MANY ANALYSIS AND DIAGNOSTIC
C           TOOLS IN DATAPLOT THAT WILL ALLOW THE 
C           TESTING OF THE RANDOMNESS AND UNIFORMITY
C           OF THIS GENERATOR.
C           SUCH CHECKING IS ESPECIALLY IMPORTANT 
C           IN LIGHT OF THE FACT THAT OTHER DATAPLOT RANDOM 
C           NUMBER GENERATOR SUBROUTINES (NORRAN--NORMAL,
C           LOGRAN--LOGISTIC, ETC.) ALL MAKE USE OF INTERMEDIATE
C           OUTPUT FROM UNIRAN.
C
C     NOTE--THE OUTPUT FROM THIS SUBROUTINE DEPENDS
C           ON THE INPUT SEED (ISEED) AND ON THE
C           VALUE OF MDIG.
C           MDIG MAY NOT BE SMALLER THAN 16.
C           MDIG MAY NOT BE LARGER THAN MAX INTEGER ON YOUR COMPUTER. 
C
C     NOTE--BECAUSE OF THE PREPONDERANCE OF MAINFRAMES
C           WHICH HAVE WORDS OF 32 BITS AND LARGER
C           (E.G, VAX (= 32 BITS), UNIVAC (= 36 BITS), CDC (= 60 BITS), ETC.)
C           MDIG HAS BEEN SET TO 32.
C           THUS THE SAME SEQUENCE OF RANDOM NUMBERS SHOULD RESULT
C           ON ALL OF THESE COMPUTERS.
C
C     NOTE--FOR SMALLER WORD SIZE COMPUTERS (E.G., 24-BIT AND 16-BIT),
C           THE VALUE OF MDIG SHOULD BE CHANGED TO 24 OR 16.
C           IN SUCH CASE, THE OUTPUT WILL NOT BE IDENTICAL TO
C           THE OUTPUT WHEN MDIG = 32.
C
C     NOTE--THE CYCLE OF THE RANDOM NUMBERS DEPENDS ON MDIG.
C           THE CYCLE FROM MDIG = 32 IS LONG ENOUGH FOR MOST
C           PRACTICAL APPLICATIONS.
C           IF A LONGER CYCLE IS DESIRED, THEN INCREASE MDIG.
C
C     NOTE--THE SEED MAY BE ANY POSITIVE INTEGER. 
C           NO APPRECIABLE DIFFERENCE IN THE QUALITY
C           OF THE RANDOM NUMBERS HAS BEEN NOTED
C           BY THE CHOICE OF THE SEED.  THERE IS NO
C           NEED TO USE PRIMES, NOR TO USE EXCEPTIONALLY
C           LARGE NUMBERS, ETC.
C
C     REFERENCES--MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
C               --JOHNSON AND KOTZ, CONTINUOUS UNIVARIATE
C                 DISTRIBUTIONS--2, 1970, PAGES 57-74.
C     WRITTEN BY--JAMES BLUE
C                 SCIENTIFIC COMPUTING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C               --DAVID KAHANER
C                 SCIENTIFIC COMPUTING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C               --GEORGE MARSAGLIA
C                 COMPUTER SCIENCE DEPARTMENT
C                 WASHINGTON STATE UNIVERSITY
C               --JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C
C     LANGUAGE--ANSI FORTRAN (1977)
C     ORIGINAL VERSION--JUNE      1972. 
C     UPDATED         --AUGUST    1974. 
C     UPDATED         --SEPTEMBER 1975. 
C     UPDATED         --NOVEMBER  1975. 
C     UPDATED         --NOVEMBER  1981. 
C     UPDATED         --MAY       1982. 
C     UPDATED         --MARCH     1984. 
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
C
      DIMENSION M(17)
C
C---------------------------------------------------------------------
C
CCCCC CHARACTER*4 IFEEDB
CCCCC CHARACTER*4 IPRINT
C
CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW
CCCCC COMMON /PRINT/IFEEDB,IPRINT
C
C-----SAVE STATEMENTS-------------------------------------------------
C
      SAVE I,J,M,M1,M2
C
C-----DATA STATEMENTS-------------------------------------------------
C
      DATA M(1),M(2),M(3),M(4),M(5),M(6),M(7),M(8),M(9),M(10),M(11),
     1     M(12),M(13),M(14),M(15),M(16),M(17)
     1/    30788,23052,2053,19346,10646,19427,23975,
     1     19049,10949,19693,29746,26748,2796,23890,
     1     29168,31924,16499/ 
      DATA M1,M2,I,J / 32767,256,5,17 / 
C
      IPR=6
C
C-----START POINT-----------------------------------------------------
C
C               ********************************************
C               **  STEP 1--                              **
C               **  CHECK THE INPUT ARGUMENTS FOR ERRORS  **
C               ********************************************
C
      IF(N.GE.1)GOTO90
      WRITE(IPR,999)
  999 FORMAT(1H )
      WRITE(IPR,51) 
   51 FORMAT(1H ,'***** ERROR IN UNIRAN--')
      WRITE(IPR,52) 
   52 FORMAT(1H ,'      THE INPUT NUMBER OF OBSERVATIONS IS ',
     1'NON-POSITIVE.')
      WRITE(IPR,53)N
   53 FORMAT(1H ,'      N = ',I8)
      GOTO9000
   90 CONTINUE
C
C               *******************************************************
C               **  STEP 2--                                         **
C               **  IF A POSITIVE INPUT SEED HAS BEEN GIVEN,         **
C               **  THEN THIS INDICATES THAT THE GENERATOR           **
C               **  SHOULD HAVE ITS INTERNAL M(.) ARRAY REDEFINED--  **
C               **  DO SO IN THIS SECTION.                           **
C               **  IF A NON-POSITIVE INPUT SEED HAS BEEN GIVEN,     **
C               **  THEN THIS INDICATES THAT THE GENERATOR           **
C               **  SHOULD CONTINUE ON FROM WHERE IT LEFT OFF,       **
C               **  AND THEREFORE THIS SECTION IS SKIPPED.           **
C               *******************************************************
C
      IF(ISEED.LE.0)GOTO290
C
CCCCC MDIG=16
      MDIG=32
C
      M1=2**(MDIG-2)+(2**(MDIG-2)-1)
      M2=2**(MDIG/2)
CCCCC ISEED3=MIN0(IABS(ISEED),M1)
      ISEED3=IABS(ISEED)
      IF(M1.LT.IABS(ISEED))ISEED3=M1
      IF(MOD(ISEED3,2).EQ.0)ISEED3=ISEED3-1
      K0=MOD(9069,M2)
      K1=9069/M2
      J0=MOD(ISEED3,M2)
      J1=ISEED3/M2
C
      DO200I=1,17
      ISEED3=J0*K0
      J1=MOD(ISEED3/M2+J0*K1+J1*K0,M2/2)
      J0=MOD(ISEED3,M2)
      M(I)=J0+M2*J1 
  200 CONTINUE
C
      I=5 
      J=17
C
  290 CONTINUE
C
C               *************************************
C               **  STEP 3--                       **
C               **  GENERATE THE N RANDOM NUMBERS  **
C               *************************************
C
      DO300L=1,N
      K=M(I)-M(J)
      IF(K.LT.0)K=K+M1
      M(J)=K
      I=I-1
      IF(I.EQ.0)I=17
      J=J-1
      IF(J.EQ.0)J=17
      AK=K
      AM1=M1
      X(L)=AK/AM1
  300 CONTINUE
C
C               ***************************************************** 
C               **  STEP 4--                                       ** 
C               **  REGARDLESS OF THE VALUE OF THE INPUT SEED,     ** 
C               **  REDEFINE THE VALUE OF ISEED UPON EXIT HERE     ** 
C               **  TO -1 WITH THE NET EFFECT THAT                 ** 
C               **  IF THE USER DOES NOT REDEFINE THE SEED         ** 
C               **  VALUE BEFORE THE NEXT CALL TO THIS GENERATOR,  ** 
C               **  THEN THIS GENERATOR WILL PICK UP               ** 
C               **  WHERE IT LEFT OFF.                             ** 
C               ***************************************************** 
C
      ISEED=(-1)
C
C               *****************
C               **  STEP 90--  **
C               **  EXIT       **
C               *****************
C
 9000 CONTINUE
      RETURN
CCCCC DEBUG TRACE,INIT
CCCCC AT 90
CCCCC TRACE ON
      END 
      SUBROUTINE GEORAN(N,P,ISEED,X)
C
C     PURPOSE--THIS SUBROUTINE GENERATES A RANDOM SAMPLE OF SIZE N
C              FROM THE GEOMETRIC DISTRIBUTION
C              WITH SINGLE PRECISION 'BERNOULLI PROBABILITY'
C              PARAMETER = P.
C              THE GEOMETRIC DISTRIBUTION USED
C              HEREIN HAS MEAN = (1-P)/P
C              AND STANDARD DEVIATION = SQRT((1-P)/(P*P))).
C              THIS DISTRIBUTION IS DEFINED FOR
C              ALL NON-NEGATIVE INTEGER X--X = 0, 1, 2, ... .
C              THIS DISTRIBUTION HAS THE PROBABILITY FUNCTION
C              F(X) = P * (1-P)**X.
C              THE GEOMETRIC DISTRIBUTION IS THE
C              DISTRIBUTION OF THE NUMBER OF FAILURES
C              BEFORE OBTAINING 1 SUCCESS IN AN
C              INDEFINITE SEQUENCE OF BERNOULLI (0,1)
C              TRIALS WHERE THE PROBABILITY OF SUCCESS
C              IN A SINGLE TRIAL = P.
C     INPUT  ARGUMENTS--N      = THE DESIRED INTEGER NUMBER
C                                OF RANDOM NUMBERS TO BE
C                                GENERATED.
C                     --P      = THE SINGLE PRECISION VALUE
C                                OF THE 'BERNOULLI PROBABILITY'
C                                PARAMETER FOR THE GEOMETRIC
C                                DISTRIBUTION.
C                                P SHOULD BE BETWEEN
C                                0.0 (EXCLUSIVELY) AND
C                                1.0 (EXCLUSIVELY).
C     OUTPUT ARGUMENTS--X      = A SINGLE PRECISION VECTOR
C                                (OF DIMENSION AT LEAST N)
C                                INTO WHICH THE GENERATED
C                                RANDOM SAMPLE WILL BE PLACED.
C     OUTPUT--A RANDOM SAMPLE OF SIZE N
C             FROM THE GEOMETRIC DISTRIBUTION
C             WITH 'BERNOULLI PROBABILITY' PARAMETER = P.
C     PRINTING--NONE UNLESS AN INPUT ARGUMENT ERROR CONDITION EXISTS.
C     RESTRICTIONS--THERE IS NO RESTRICTION ON THE MAXIMUM VALUE
C                   OF N FOR THIS SUBROUTINE.
C                 --P SHOULD BE BETWEEN 0.0 (EXCLUSIVELY)
C                   AND 1.0 (EXCLUSIVELY).
C     OTHER DATAPAC   SUBROUTINES NEEDED--UNIRAN.
C     FORTRAN LIBRARY SUBROUTINES NEEDED--ALOG.
C     MODE OF INTERNAL OPERATIONS--SINGLE PRECISION.
C     LANGUAGE--ANSI FORTRAN (1977)
C     COMMENT--NOTE THAT EVEN THOUGH THE OUTPUT
C              FROM THIS DISCRETE RANDOM NUMBER
C              GENERATOR MUST NECESSARILY BE A
C              SEQUENCE OF ***INTEGER*** VALUES,
C              THE OUTPUT VECTOR X IS SINGLE
C              PRECISION IN MODE.
C              X HAS BEEN SPECIFIED AS SINGLE
C              PRECISION SO AS TO CONFORM WITH THE DATAPAC
C              CONVENTION THAT ALL OUTPUT VECTORS FROM ALL
C              DATAPAC SUBROUTINES ARE SINGLE PRECISION.
C              THIS CONVENTION IS BASED ON THE BELIEF THAT
C              1) A MIXTURE OF MODES (FLOATING POINT
C              VERSUS INTEGER) IS INCONSISTENT AND
C              AN UNNECESSARY COMPLICATION
C              IN A DATA ANALYSIS; AND
C              2) FLOATING POINT MACHINE ARITHMETIC
C              (AS OPPOSED TO INTEGER ARITHMETIC)
C              IS THE MORE NATURAL MODE FOR DOING
C              DATA ANALYSIS.
C     REFERENCES--TOCHER, THE ART OF SIMULATION,
C                 1963, PAGES 14-15.
C               --HAMMERSLEY AND HANDSCOMB, MONTE CARLO METHODS,
C                 1964, PAGE 36.
C               --FELLER, AN INTRODUCTION TO PROBABILITY
C                 THEORY AND ITS APPLICATIONS, VOLUME 1,
C                 EDITION 2, 1957, PAGES 155-157, 210.
C               --NATIONAL BUREAU OF STANDARDS APPLIED MATHEMATICS
C                 SERIES 55, 1964, PAGE 929.
C     WRITTEN BY--JAMES J. FILLIBEN
C                 STATISTICAL ENGINEERING DIVISION
C                 CENTER FOR APPLIED MATHEMATICS
C                 NATIONAL BUREAU OF STANDARDS
C                 WASHINGTON, D. C. 20234
C                 PHONE--301-921-3651
C     NOTE--DATAPLOT IS A REGISTERED TRADEMARK
C           OF THE NATIONAL BUREAU OF STANDARDS.
C           THIS SUBROUTINE MAY NOT BE COPIED, EXTRACTED,
C           MODIFIED, OR OTHERWISE USED IN A CONTEXT
C           OUTSIDE OF THE DATAPLOT LANGUAGE/SYSTEM.
C     LANGUAGE--ANSI FORTRAN (1966)
C               EXCEPTION--HOLLERITH STRINGS IN FORMAT STATEMENTS
C                          DENOTED BY QUOTES RATHER THAN NH.
C     VERSION NUMBER--82/7
C     ORIGINAL VERSION--NOVEMBER  1975.
C     UPDATED         --DECEMBER  1981.
C     UPDATED         --MAY       1982.
C
C-----CHARACTER STATEMENTS FOR NON-COMMON VARIABLES-------------------
C
C---------------------------------------------------------------------
C
      DIMENSION X(*)
C
C---------------------------------------------------------------------
C
CCCCC CHARACTER*4 IFEEDB
CCCCC CHARACTER*4 IPRINT
C
CCCCC COMMON /MACH/IRD,IPR,CPUMIN,CPUMAX,NUMBPC,NUMCPW,NUMBPW
CCCCC COMMON /PRINT/IFEEDB,IPRINT
C
      IPR=6
C
C-----START POINT-----------------------------------------------------
C
C     CHECK THE INPUT ARGUMENTS FOR ERRORS
C
      IF(N.LT.1)GOTO50
      IF(P.LE.0.0.OR.P.GE.1.0)GOTO55
      GOTO90
   50 WRITE(IPR, 5)
      WRITE(IPR,47)N
      RETURN
   55 WRITE(IPR,11)
      WRITE(IPR,46)P
      RETURN
   90 CONTINUE
    5 FORMAT(1H , 91H***** FATAL ERROR--THE FIRST  INPUT ARGUMENT TO THE
     1 GEORAN SUBROUTINE IS NON-POSITIVE *****)
   11 FORMAT(1H ,115H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE
     1 GEORAN SUBROUTINE IS OUTSIDE THE ALLOWABLE (0,1) INTERVAL *****)
   46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****)
   47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8   ,6H *****)
C
C     GENERATE N UNIFORM (0,1) RANDOM NUMBERS;
C
      CALL UNIRAN(N,ISEED,X)
C
C     GENERATE N GEOMETRIC RANDOM NUMBERS
C     USING THE PERCENT POINT FUNCTION TRANSFORMATION METHOD.
C
      DO100I=1,N
      IF(X(I).EQ.0.0)GOTO100
      ARG1=1.0-X(I)
      ARG2=1.0-P
      ANUM=ALOG(ARG1)
      ADEN=ALOG(ARG2)
      RATIO=ANUM/ADEN
      IRATIO=RATIO
      X(I)=IRATIO
      ARATIO=IRATIO
      IF(ARATIO.EQ.RATIO)X(I)=IRATIO-1
  100 CONTINUE
C
      RETURN
      END
