uniran(3f) - [M_datapac:RANDOM] generate Uniform random numbers
SUBROUTINE UNIRAN(N,Iseed,X)
INTEGER,intent(in) :: N
INTEGER,intent(inout) :: Iseed
REAL(kind=wp),intent(out) :: X(:)
UNIRAN(3f) generates a random sample of size N from the uniform
(rectangular) distribution on the unit interval (0,1).
This distribution has mean = 0.5 and standard deviation = sqrt(1/12)
= 0.28867513. This distribution has the probability density function
f(X) = 1
N The desired integer number of random numbers to be generated.
ISEED An integer iseed value. Should be set to a non-negative value
to start a new sequence of values. Will be set to -1 on return
to indicate the next call should continue the current random
sequence walk.
X A vector (of dimension at least N) into which the generated
random sample of size N from the rectangular distribution on
(0,1) will be placed.
Sample program:
program demo_uniran
use M_datapac, only : uniran, plotxt, sort, label
implicit none
integer,parameter :: n=400
real :: x(n)
integer :: iseed
call label('uniran')
iseed=1234
call UNIRAN(n,Iseed,X)
call plotxt(x,n) ! plot random values
call sort(x,n,x) ! sort values
call plotxt(x,n) ! should display the f(x)=1 nature
! of the distribution
end program demo_uniran
Results:
THE FOLLOWING IS A PLOT OF X(I) (VERTICALLY) VERSUS I (HORIZONTALLY
I-----------I-----------I-----------I-----------I
0.9982013E+00 - X X X X X X X
0.9566447E+00 I X XX X XX XXXX X XX
0.9150882E+00 I X XXXX X XXXX X X X X X XX X X
0.8735316E+00 I XX XXX XX XXXXXX X X X X X
0.8319750E+00 I XXX X X X X X XX X X X
0.7904184E+00 I XXX XX X X X X X XXX
0.7488618E+00 - X XX X X X X X X X X X XXX X
0.7073053E+00 I X X X X X X X X X X X X
0.6657487E+00 I X XXX X XX XXX X X X X X X X
0.6241921E+00 I X X X X XX X X X X
0.5826355E+00 I X X X XX X X X X X
0.5410789E+00 I XX X XX X X XXX XX X X
0.4995224E+00 - X X XX X XXX XX
0.4579658E+00 I XXX X XXX X X X XX XX X X X X
0.4164092E+00 I XX X X X X X X X X X
0.3748527E+00 I XX XX XX X X XX XX XX X X
0.3332961E+00 I X XXX X X X X X XXXX X X
0.2917395E+00 I X X X X X X XX XX XX X
0.2501829E+00 - XX XXX X X X X XX XXXX X
0.2086263E+00 I X X X X X X XX X X XX
0.1670697E+00 I X X X XX X XX XX XXX X X XXX X
0.1255132E+00 I X XXX X X X X XX X
0.8395660E-01 I X X X X XX X X X X X X X XX
0.4240000E-01 I X XX X X X X X X X X X
0.8433913E-03 - X X X X X
I-----------I-----------I-----------I-----------I
0.1000E+01 0.1008E+03 0.2005E+03 0.3002E+03 0.4000E+03
THE FOLLOWING IS A PLOT OF X(I) (VERTICALLY) VERSUS I (HORIZONTALLY
I-----------I-----------I-----------I-----------I
0.9982013E+00 - XX
0.9566447E+00 I XXX
0.9150882E+00 I XXXX
0.8735316E+00 I XXX
0.8319750E+00 I XXX
0.7904184E+00 I XXX
0.7488618E+00 - XXX
0.7073053E+00 I XXX
0.6657487E+00 I XXX
0.6241921E+00 I XX
0.5826355E+00 I XXX
0.5410789E+00 I XXX
0.4995224E+00 - XX
0.4579658E+00 I XXXX
0.4164092E+00 I XX
0.3748527E+00 I XXX
0.3332961E+00 I XXX
0.2917395E+00 I XX
0.2501829E+00 - XXXX
0.2086263E+00 I XX
0.1670697E+00 I XXXX
0.1255132E+00 I XXX
0.8395660E-01 I XXX
0.4240000E-01 I XX
0.8433913E-03 - X
I-----------I-----------I-----------I-----------I
0.1000E+01 0.1008E+03 0.2005E+03 0.3002E+03 0.4000E+03
The original DATAPAC library was written by James Filliben of the
Statistical Engineering Division, National Institute of Standards
and Technology.
* James Blue
Scientific Computing Division
Center for Applied Mathematics
National Bureau of Standards
Washington, D. C. 20234
* David Kahaner
Scientific Computing Division
Center for Applied Mathematics
National Bureau of Standards
* George Marsaglia
Computer Science Department
Washington State University
* James J. Filliben
Statistical Engineering Division
Center for Applied Mathematics
National Bureau of Standards
John Urban, 2022.05.31
CC0-1.0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | N | |||
| integer, | intent(inout) | :: | Iseed | |||
| real(kind=wp), | intent(out) | :: | X(:) |
SUBROUTINE UNIRAN(N,Iseed,X) INTEGER,intent(in) :: N INTEGER,intent(inout) :: Iseed REAL(kind=wp),intent(out) :: X(:) REAL(kind=wp) :: ak , am1 INTEGER i, iseed3, j, j0, j1, k, k0, k1, l, m1, m2, mdig INTEGER m(17) ! MODE OF INTERNAL OPERATIONS--. ! ! ALGORITHM--FIBONACCI GENERATOR AS DEFINED BY GEORGE MARSAGLIA. ! ! NOTE--THIS GENERATOR IS TRANSPORTABLE. ! IT IS NOT MACHINE-INDEPENDENT ! IN THE SENSE THAT FOR A GIVEN VALUE ! OF THE INPUT SEED ISEED AND FOR A GIVEN VALUE ! OF MDIG (TO BE DEFINED BELOW), ! THE SAME SEQUENCE OF UNIFORM RANDOM ! NUMBERS WILL RESULT ON DIFFERENT COMPUTERS ! (VAX, PRIME, PERKIN-ELMER, IBM, UNIVAC, HONEYWELL, ETC.) ! ! NOTE--IF MDIG = 32 AND IF ISEED = 305, ! THEN THE OUTPUT FROM THIS GENERATOR SHOULD BE AS FOLLOWS-- ! THE FIRST NUMBER TO RESULT IS .4771580... ! THE SECOND NUMBER TO RESULT IS .4219293... ! THE THIRD NUMBER TO RESULT IS .6646181... ! ... ! THE THOUSANDTH NUMBER TO RESULT IS .2036834... ! ! NOTE--IF MDIG = 16 AND IF ISEED = 305, ! THEN THE OUTPUT FROM THIS GENERATOR SHOULD BE AS FOLLOWS-- ! THE FIRST NUMBER TO RESULT IS .027832881... ! THE SECOND NUMBER TO RESULT IS .56102176... ! THE THIRD NUMBER TO RESULT IS .41456343... ! ... ! THE THOUSANDTH NUMBER TO RESULT IS .19797357... ! ! NOTE--IT IS RECOMMENDED THAT UPON ! IMPLEMENTATION OF DATAPLOT, THE OUTPUT ! FROM UNIRAN BE CHECKED FOR AGREEMENT ! WITH THE ABOVE SAMPLE OUTPUT. ! ALSO, THERE ARE MANY ANALYSIS AND DIAGNOSTIC ! TOOLS IN DATAPLOT THAT WILL ALLOW THE ! TESTING OF THE RANDOMNESS AND UNIFORMITY ! OF THIS GENERATOR. ! SUCH CHECKING IS ESPECIALLY IMPORTANT ! IN LIGHT OF THE FACT THAT OTHER DATAPLOT RANDOM ! NUMBER GENERATOR SUBROUTINES (NORRAN--NORMAL, ! LOGRAN--LOGISTIC, ETC.) ALL MAKE USE OF INTERMEDIATE ! OUTPUT FROM UNIRAN. ! ! NOTE--THE OUTPUT FROM THIS SUBROUTINE DEPENDS ! ON THE INPUT SEED (ISEED) AND ON THE ! VALUE OF MDIG. ! MDIG MAY NOT BE SMALLER THAN 16. ! MDIG MAY NOT BE LARGER THAN MAX INTEGER ON YOUR COMPUTER. ! ! NOTE--BECAUSE OF THE PREPONDERANCE OF MAINFRAMES ! WHICH HAVE WORDS OF 32 BITS AND LARGER ! (E.G, VAX (= 32 BITS), UNIVAC (= 36 BITS), CDC (= 60 BITS), ETC.) ! MDIG HAS BEEN SET TO 32. ! THUS THE SAME SEQUENCE OF RANDOM NUMBERS SHOULD RESULT ! ON ALL OF THESE COMPUTERS. ! ! NOTE--FOR SMALLER WORD SIZE COMPUTERS (E.G., 24-BIT AND 16-BIT), ! THE VALUE OF MDIG SHOULD BE CHANGED TO 24 OR 16. ! IN SUCH CASE, THE OUTPUT WILL NOT BE IDENTICAL TO ! THE OUTPUT WHEN MDIG = 32. ! ! NOTE--THE CYCLE OF THE RANDOM NUMBERS DEPENDS ON MDIG. ! THE CYCLE FROM MDIG = 32 IS LONG ENOUGH FOR MOST ! PRACTICAL APPLICATIONS. ! IF A LONGER CYCLE IS DESIRED, THEN INCREASE MDIG. ! ! NOTE--THE SEED MAY BE ANY POSITIVE INTEGER. ! NO APPRECIABLE DIFFERENCE IN THE QUALITY ! OF THE RANDOM NUMBERS HAS BEEN NOTED ! BY THE CHOICE OF THE SEED. THERE IS NO ! NEED TO USE PRIMES, NOR TO USE EXCEPTIONALLY ! LARGE NUMBERS, ETC. ! !-----SAVE STATEMENTS------------------------------------------------- ! SAVE i , j , m , m1 , m2 ! !-----DATA STATEMENTS------------------------------------------------- ! DATA m(1) , m(2) , m(3) , m(4) , m(5) , m(6) , m(7) , m(8) , & & m(9) , m(10) , m(11) , m(12) , m(13) , m(14) , m(15) , & & m(16) , m(17)/30788 , 23052 , 2053 , 19346 , 10646 , 19427 , & & 23975 , 19049 , 10949 , 19693 , 29746 , 26748 , 2796 , & & 23890 , 29168 , 31924 , 16499/ DATA m1 , m2 , i , j/32767 , 256 , 5 , 17/ ! !-----START POINT----------------------------------------------------- ! ! ******************************************** ! ** STEP 1-- ** ! ** CHECK THE INPUT ARGUMENTS FOR ERRORS ** ! ******************************************** ! IF ( N>=1 ) THEN ! ! ******************************************************* ! ** STEP 2-- ** ! ** IF A POSITIVE INPUT SEED HAS BEEN GIVEN, ** ! ** THEN THIS INDICATES THAT THE GENERATOR ** ! ** SHOULD HAVE ITS INTERNAL M(.) ARRAY REDEFINED-- ** ! ** DO SO IN THIS SECTION. ** ! ** IF A NON-POSITIVE INPUT SEED HAS BEEN GIVEN, ** ! ** THEN THIS INDICATES THAT THE GENERATOR ** ! ** SHOULD CONTINUE ON FROM WHERE IT LEFT OFF, ** ! ** AND THEREFORE THIS SECTION IS SKIPPED. ** ! ******************************************************* ! IF ( Iseed>0 ) THEN ! mdig = 32 ! m1 = 2**(mdig-2) + (2**(mdig-2)-1) m2 = 2**(mdig/2) iseed3 = IABS(Iseed) IF ( m1<IABS(Iseed) ) iseed3 = m1 IF ( MOD(iseed3,2)==0 ) iseed3 = iseed3 - 1 k0 = MOD(9069,m2) k1 = 9069/m2 j0 = MOD(iseed3,m2) j1 = iseed3/m2 ! DO i = 1 , 17 iseed3 = j0*k0 j1 = MOD(iseed3/m2+j0*k1+j1*k0,m2/2) j0 = MOD(iseed3,m2) m(i) = j0 + m2*j1 ENDDO ! i = 5 j = 17 ENDIF ! ! ! ************************************* ! ** STEP 3-- ** ! ** GENERATE THE N RANDOM NUMBERS ** ! ************************************* ! DO l = 1 , N k = m(i) - m(j) IF ( k<0 ) k = k + m1 m(j) = k i = i - 1 IF ( i==0 ) i = 17 j = j - 1 IF ( j==0 ) j = 17 ak = k am1 = m1 X(l) = ak/am1 ENDDO ! ! ***************************************************** ! ** STEP 4-- ** ! ** REGARDLESS OF THE VALUE OF THE INPUT SEED, ** ! ** REDEFINE THE VALUE OF ISEED UPON EXIT HERE ** ! ** TO -1 WITH THE NET EFFECT THAT ** ! ** IF THE USER DOES NOT REDEFINE THE SEED ** ! ** VALUE BEFORE THE NEXT CALL TO THIS GENERATOR, ** ! ** THEN THIS GENERATOR WILL PICK UP ** ! ** WHERE IT LEFT OFF. ** ! ***************************************************** ! Iseed = (-1) ELSE WRITE (G_IO,99001) 99001 FORMAT (' ') WRITE (G_IO,99002) 99002 FORMAT (' ','***** Error in UNIRAN--') WRITE (G_IO,99003) 99003 FORMAT (' ',' The input number of observations is non-positive.') WRITE (G_IO,99004) N 99004 FORMAT (' ',' N = ',I0) ENDIF ! ! ***************** ! ** STEP 90-- ** ! ** EXIT ** ! ***************** ! END SUBROUTINE UNIRAN