Topic: Random Numbers
FYI - http://csrc.nist.gov/publications/nistb … 015_08.pdf
Should this affect some Simply Fortran Users?
Bob
Re: Random Numbers
I'm not sure how many users are relying on Simply Fortran for writing encryption software, but some might be. Random number generation is quite a complex, yet fascinating, field. Many users might be interested, but probably more for the science behind the document.
Jeff Armstrong
Approximatrix, LLC
Approximatrix, LLC
Re: Random Numbers
I am interested in some of the details about the intrinsic (pseudo) random number generator. For example, do you know what algorithm is used, and what is period of the pseudo-random sequence?
Thank you.
Re: Random Numbers
From the runtime library source code:
/*
We use the xoshiro256** generator, a fast high-quality generator
that:
- passes TestU1 without any failures
- provides a "jump" function making it easy to provide many
independent parallel streams.
- Long period of 2**256 - 1
A description can be found at
http://prng.di.unimi.it/
or
https://arxiv.org/abs/1805.01407
The paper includes public domain source code which is the basis for
the implementation below.
*/Just to be clear, that is the algorithm the current runtime library is using. This algorithm could always change if something better comes along.
Jeff Armstrong
Approximatrix, LLC
Approximatrix, LLC
Re: Random Numbers
Terrific!
Thank you very much for this information.
Re: Random Numbers
All,
For those that may have an interest in random number generation in SF, I've posted a basic program in SF that uses random numbers.
Even though it doesn't address the original post here, it might be instructive for those interested. Although the code used has been modified by me, it comes from various miscellaneous sources and is not original. I've named it WH_RAND after the WICHMAN-HILL method.
!***********************************************************************
!* WH_RAND: UNIFORM RANDOM NUMBER GENERATOR VERSION 3.0 *
!* *
!* RANF IS A UNIFORM RANDOM NUMBER ROUTINE, ORIGIN IS UNKOWN. *
!* RDUNWH IS USED IN PLACE OF RANF AS A UNIFORM RANDOM NUMBER ROUTINE.*
!* *
!* MODIFIED WITH RANDOM SEEDS RANSEED() TO INITIALIZE THE RANDOM SEED *
!* USING THE SYSTEM CLOCK TO ALLOW FOR NEWLY GENERATED NUMBERS AFTER *
!* EACH EXECUTION OF THE ROUTINE. *
!* *
!* REVISED AND MODIFIED BY F.W. PERRELLA, PH.D. FEBRUARY 29, 2024 *
!***********************************************************************
!
!-----------------------------------------------------------------------
!---------- MODULES ----------------------------------------------------
!-----------------------------------------------------------------------
MODULE NRTYPES
!--- KIND TYPES FOR 64- AND 32-BIT SIGNED INTEGERS
INTEGER, PARAMETER :: I32 = SELECTED_INT_KIND(9)
!INTEGER, PARAMETER :: I64 = SELECTED_INT_KIND(18)
INTEGER, PARAMETER :: SP = SELECTED_REAL_KIND(6)
!INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(15)
END MODULE NRTYPES
MODULE IRANDOM
USE NRTYPES
INTEGER(I32) :: K2
INTEGER(I32) :: IRAND(250)
INTEGER(I32) :: IX,IY,IZ
END MODULE IRANDOM
MODULE WH_RAND
USE NRTYPES
IMPLICIT NONE
PUBLIC :: RANSEED, RANSET, RANF, RDUNWH
INTERFACE RANSEED
MODULE PROCEDURE RANSEED
END INTERFACE RANSEED
INTERFACE RANSET
MODULE PROCEDURE RANSET
END INTERFACE RANSET
INTERFACE RANF
MODULE PROCEDURE RANF
END INTERFACE RANF
INTERFACE RDUNWH
MODULE PROCEDURE RDUNWH
END INTERFACE RDUNWH
PRIVATE
!-------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------
SUBROUTINE RANSEED(ISEED)
IMPLICIT NONE
INTEGER(I32), INTENT(OUT) :: ISEED
INTEGER(I32) :: I
REAL(SP) :: A, AC, AM, R1
PARAMETER(A = 4.5875, AC = 45381.6693, AM = 214748.3648)
CALL INIT_RANDOM_SEED(ISEED)
I = 0
DO WHILE ((ISEED == 0).AND.(I < 10))
I = I + 1
CALL INIT_RANDOM_SEED(ISEED)
END DO
IF(ISEED .EQ. 0) ISEED = 10000
IF(ISEED .LT. 0) ISEED = -ISEED
R1 = 1.0 + MODULO(A*REAL(ISEED) + AC, AM)
ISEED = INT(R1)
RETURN
END SUBROUTINE RANSEED
!--- INITIALIZE A RANDOM SEED FROM THE SYSTEM CLOCK AT EVERY RUN.
SUBROUTINE INIT_RANDOM_SEED(ISEED)
IMPLICIT NONE
INTEGER(I32), INTENT(OUT) :: ISEED
INTEGER(I32) :: I, N, ICLOCK
INTEGER(I32), DIMENSION(:), ALLOCATABLE :: RSEED
CALL RANDOM_SEED(SIZE = N)
ALLOCATE(RSEED(N))
CALL SYSTEM_CLOCK(COUNT = ICLOCK)
IF (ICLOCK < 0) ICLOCK = 30000
RSEED = ICLOCK + 37 * (/ (I - 1, I = 1, N) /)
CALL RANDOM_SEED(PUT = RSEED)
ISEED = INT( RSEED(N) )
DEALLOCATE(RSEED)
RETURN
END SUBROUTINE INIT_RANDOM_SEED
SUBROUTINE RANSET(ISD)
USE IRANDOM
IMPLICIT NONE
INTEGER(I32), INTENT(IN OUT) :: ISD
INTEGER(I32) :: I, J, K, JRAND, J1
ISD = MODULO(ISD, 65536)
DO I = 1, 250
DO J = 1, 5
ISD = MODULO(ISD*8973 + 1, 65536) - 32768
END DO
JRAND = ISD*65536
J1 = MODULO(I, 2)
DO K = 1, 5
ISD = MODULO(ISD*8973 + J1, 65536)
END DO
IRAND(I) = JRAND + ISD
END DO
K2 = 1
RETURN
END SUBROUTINE RANSET
FUNCTION RANF()
USE IRANDOM
IMPLICIT NONE
REAL(SP) :: RANF
INTEGER(I32) :: L1, IRANF
L1 = MODULO(K2 + 146, 250) + 1
IRANF = IEOR(IRAND(K2), IRAND(L1))
RANF = REAL(IRANF)*2.328306E-10 + 0.5
IRAND(K2) = IRANF
K2 = MODULO(K2, 250) + 1
RETURN
END FUNCTION RANF
FUNCTION RDUNWH (IX,IY,IZ)
!--- GENERATES UNIFORM PSEUDO-RANDOM DEVIATES IN THE
!--- INTERVAL [0,1) USING THE WICHMAN-HILL METHOD.
!--- AN EXACT VALUE OF ZERO IS POSSIBLE BUT NOT LIKELY.
!--- THE CYCLE LENGTH EXCEEDS 6.95E12.
!---
!--- CURRENT VERSION COMPLETED MAY 15, 1987
!---
!--- REFERENCE:
!--- GRIFFITHS, P. AND HILL, I.D. (EDS.),
!--- 'APPLIED STATISTICS ALGORITHMS',
!--- THE ROYAL STATISTICAL SOCIETY/ELLIS HARWOOD LIMITED,
!--- ALGORITHM AS 183, PP. 238-242 (1985).
!---
!--- THE INTEGERS IX, IY, AND IZ REQUIRE INTEGER
!--- ARITHMETIC UP TO 5,212,632
!--- BEFORE THE FIRST CALL TO THIS ROUTINE THE INTEGERS
!--- IX, IY, AND IZ SHOULD BE SET TO VALUES WITHIN THE
!--- RANGE [1,30000]. AFTER THAT THEIR VALUES SHOULD NOT
!--- BE CHANGED EXCEPT BY THIS ROUTINE.
IMPLICIT NONE
REAL(SP) :: RDUNWH, SUMX
INTEGER(I32), INTENT(IN OUT) :: IX, IY, IZ
IX = MOD(171*IX,30269)
IY = MOD(172*IY,30307)
IZ = MOD(170*IZ,30323)
SUMX = REAL(IX)/30269.0+REAL(IY)/30307.0+REAL(IZ)/30323.0
RDUNWH = AMOD(SUMX,1.0)
RETURN
END FUNCTION RDUNWH
END MODULE WH_RAND
!-----------------------------------------------------------------------
!---------- END MODULES ------------------------------------------------
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!---------- TEST PROGRAM -----------------------------------------------
!-----------------------------------------------------------------------
PROGRAM RANDOM_TEST
USE NRTYPES
USE WH_RAND
IMPLICIT NONE
CHARACTER(1) :: CR
INTEGER(I32) :: I, J, K, CNT, N, ISEED, IX,IY,IZ, IXYZ
INTEGER(I32) :: NBINS, BIN, YBIN, IFACTOR, ICYCLE, GCOUNT
REAL(SP) :: YMIN, YMAX, YMEAN, YSIGMA
REAL(SP) :: SUMN, SUMSQ, COVCNT, COVSCORE
REAL(SP) :: Y, RANGEX, W, GMEAN, GSIGMA, GSCORE
PARAMETER(IFACTOR = 100, ICYCLE = 5)
ALLOCATABLE :: Y(:), IXYZ(:), BIN(:), YBIN(:), RANGEX(:)
ALLOCATE(Y(1), IXYZ(3))
ALLOCATE(BIN(IFACTOR), YBIN(IFACTOR), RANGEX(IFACTOR))
!--- INITIALIZE SEED VALUE
CALL RANSEED(ISEED)
DO I = 1, 3
!--- RESET SEED VALUE
CALL RANSET(ISEED)
!--- CHECK THAT THE VALUES IN IN RANGE.
IF (ISEED > 30000) THEN
ISEED = 30000*INT(FLOAT(ISEED - 30000)/FLOAT(ISEED))
END IF
IF (ISEED == 0) ISEED = 13849 / I
IXYZ(I) = ISEED
END DO
!--- SEEDS FOR RDUNWH() RANDOM NUMBER GENERATOR
IX = IXYZ(1)
IY = IXYZ(2)
IZ = IXYZ(3)
!--- GENERATE N*N = 10,000 UNIFORM RANDOM NUMBERS
N = 100
1 CONTINUE
GMEAN = 0.0
GSIGMA = 0.0
GSCORE = 0.0
GCOUNT = 0.0
IF (ALLOCATED(Y)) DEALLOCATE (Y)
ALLOCATE(Y(N*N))
WRITE(*,"(/1X,A)") 'WH_RAND UNIFORM RANDOM NUMBER GENERATOR'
WRITE(*,"(/1X,A,3I6)") 'INITIAL RDUNWH SEEDS IX,IY,IZ =', IXYZ(1:3)
WRITE(*,"(A)") REPEAT('~',40)
WRITE(*,"(/A)") REPEAT('-',72)
!--- CALCULATE ICYCLE MEAN VALUES
DO K = 1, ICYCLE
!--- RESET SEED VALUES
J = RDUNWH(IX,IY,IZ)
!--- STORE NEW SEEDS
IXYZ(1) = IX
IXYZ(2) = IY
IXYZ(3) = IZ
CNT = 0
DO J = 1, N
DO I = 1, N
CNT = CNT + 1
!--- GENERATE AND SAVE RANDOM DATA Y
!--- Y(CNT) = RANF()
Y(CNT) = RDUNWH(IX,IY,IZ)
END DO
END DO
!--- CALCULATE THE MINUMUM AND MAXIMUM VALUES
YMIN = MINVAL(Y)
YMAX = MAXVAL(Y)
!--- CALCULATE THE MEAN AND STANDARD DEVIATION
SUMN = SUM( Y(:) )
SUMSQ = SUM( Y(:)**2 )
YMEAN = SUMN / FLOAT(CNT)
!--- YSIGMA = SQRT( SUM( (Y(:) - YMEAN)**2 ) / CNT )
YSIGMA = SQRT( SUMSQ/FLOAT(CNT) - YMEAN**2 )
!--- The covariance score is computed over each pair of the
!--- data set Y, so that Y1 is paired with Y2, Y3 with Y4, etc.
COVCNT = 0.0
DO I = 1, CNT-1
COVCNT = COVCNT + (Y(I+1)-YMEAN)*(Y(I)-YMEAN)
END DO
COVSCORE = COVCNT / FLOAT(CNT-1)
WRITE(*,"(1X,A,1X,I6,1X,A,F10.4,2X,A,F10.4,2X,A,F10.4)") &
"COUNT", CNT, "MEAN=", YMEAN, " SIGMA=", YSIGMA, &
" COV SCORE=", COVSCORE
IF (ICYCLE > 1) THEN
!--- GRAND MEAN, GRAND SIGMA, GRAND SCORE
GMEAN = GMEAN + YMEAN
GSIGMA = GSIGMA + YSIGMA
GSCORE = GSCORE + COVSCORE
GCOUNT = GCOUNT + CNT
END IF
END DO
WRITE(*,"(A)") REPEAT('-',72)
IF (ICYCLE > 1) THEN
!--- GRAND AVERAGES
GMEAN = GMEAN / FLOAT(ICYCLE)
GSIGMA = GSIGMA / FLOAT(ICYCLE)
GSCORE = GSCORE / FLOAT(ICYCLE)
WRITE(*,"(1X,A,1X,I6,1X,A,F10.4,2X,A,F10.4,2X,A,F10.6)") &
"COUNT", GCOUNT, "MEAN=", GMEAN, &
" SIGMA=", GSIGMA, " COV SCORE=", GSCORE
END IF
!-------------------------------------------------------------------
!--- THIS SECTION USES A SCORE ARRAY AND A RANGE ARRAY, COUNTS
!--- THE NUMBER OF EACH SCORES IN EACH RANGE, AND DISPLAYS A
!--- HORIZONTAL BAR HISTOGRAM.
!-------------------------------------------------------------------
!--- CALCULATE THE NUMBER OF BINS FOR Y VALUES
W = 5.0 * (FLOAT(IFACTOR)*YSIGMA + 0.10) * (N*N)**(-1.0/3.0)
!--- ADJUST VALUES BY IFACTOR FOR PLOTTING
Y(:) = FLOAT(IFACTOR)*Y(:)
YMIN = FLOAT(IFACTOR)*YMIN
YMAX = FLOAT(IFACTOR)*YMAX
NBINS = INT((YMAX - YMIN) / W)
WRITE(*,"(/A,F8.3,2X,A,I3)") 'W=', W, ' NBINS=', NBINS
!--- CALCULATE THE RANGE OF BINS STARTING AT Y MINIMUM
RANGEX(:) = YMIN
DO I = 1, NBINS
IF ( I .EQ. 1 ) THEN
RANGEX(I) = RANGEX(I) + W
ELSE
RANGEX(I) = RANGEX(I-1) + W
END IF
END DO
WRITE(*,"(A)") REPEAT('~',40)
!--- STORE THE Y VALUES IN THE APPROPRIATE BIN
WRITE(*,"(/1X,A)") 'NUMERIC HISTOGRAM'
DO I = 1, NBINS
YBIN(I) = COUNT( Y(:) .LE. RANGEX(I) )
IF( I .GT. 1 ) THEN
BIN(I) = YBIN(I) - YBIN(I-1)
ELSE
BIN(I) = YBIN(I)
ENDIF
WRITE(*,'(F10.2,I10)') RANGEX(I), BIN(I)
END DO
WRITE(*,"(A)") REPEAT('~',40)
!--- DISPLAY HORIZONTAL BAR HISTOGRAM
WRITE(*,"(/1X,A)") 'HISTOGRAM PLOT'
DO I = 1, NBINS
J = INT( FLOAT(BIN(I))/ FLOAT(NBINS) / 1.0 )
WRITE(*, "(F4.1, A, A)") &
REAL(I), ": ", REPEAT("=", J)
END DO
WRITE(*,"(A)") REPEAT('~',40)
!...END RANDOM NUMBER ROUTINE
WRITE(*,'(/A$)') "PRESS '1' TO CONTINUE OR 'ENTER' TO QUIT:"
READ(*,'(A)') CR
IF (TRIM(CR) == "1") THEN
GO TO 1
ELSE
STOP 'PROGRAM ENDED...'
END IF
END PROGRAM RANDOM_TEST
!-----------------------------------------------------------------------
!---------- END TEST PROGRAM -------------------------------------------
!-----------------------------------------------------------------------