A Pseudo-Random Bit Generator

A Simple and Efficient Pseudo-Random Bit Generator

Link to original document
http://www.stud.uni-muenchen.de/%7Emok-kong.shen/#paper1

M. K. Shen

The following is a design sketch of a simple and efficient pseudo-random bit generator that has successfully passed U. M. Maurer’s universal statistical test for random bit generators [1].

(a) Use a good PRNG that produces a real-valued sequence (y_i) (i = 1,2,….) of uniform distribution in [0,1).

(b) Form the 24-bit integer sequence (z_i) with z_i := 2^24 * y_i which provides the blocks (3 bytes) of pseudo-random bits desired.

For the PRNG in (a) we choose to use the compound PRNG designed by the author in [2] which is based on the idea of randomly interlacing the outputs of a number of constituent PRNGs (see below for the algorithm). With L=8, Q=5000, K=1000000 (in Maurer’s notation) we obtained in an experiment the following values of the test parameter ftu, with n denoting the number of constituent generators of our compound PRNG. The max, min and mean values are with respect to 100 test runs done for each value of n (each run with a different compound PRNG).

        n     min ftu     max ftu     mean ftu
        1     6.804041    7.236945    7.178144
        2     7.155620    7.223903    7.184418
        5     7.173252    7.193011    7.183730
       10     7.177052    7.189913    7.183914
       20     7.178847    7.187996    7.183980
       50     7.177911    7.186389    7.183700
      100     7.180321    7.186633    7.183638
      200     7.181037    7.185719    7.183679
      500     7.180933    7.186409    7.183587
     1000     7.180637    7.186052    7.183680

For rejection rate rho=0.01 the threshold values computed from Table
1 of [1] are:

      t1 = 7.180865     t2 = 7.186466

Thus it can be seen that with n of the order of 50 or larger our
scheme of pseudo-random bit generation is indeed very satisfactory.

Discussion and notes:

  1. Our scheme is designed for computers having 24-bit mantissa for real-valued numbers (32-bit word length). Thus it is from the outset not sensible to use in (b) a multiplier of y_i larger than 2^24. Evidently smaller multipliers would be better than larger ones. This has in fact been verified by using sequences (y_i) of inferior quality. That a multiplier smaller than 2^24 is not called for, as indicated by the results of our experiment, is attributable to the superior quality of the sequence (y_i) generated by our compound PRNG.
  2. The algorithm of the compound PRNG we used to generate (y_i) may be very simply described as follows:
            j := 0;
            do  output(G(j));  j := trunc(n*G(j))  od;
    
    

    where G(j) (0 <= j < n) denote calls to n internal constituent PRNGs with output in the real-valued interval [0,1). (As discussed in [2], these n PRNGs can be fairly arbitrarily chosen, subjected only to a few very minor constraints. In particular, the PRNGs may be of any type.) Note that the order of activation of the n constituent PRNGs is determined by pseudo-random numbers produced by these PRNGs themselves and that only one half of the total outputs generated by them appear as the output of the compound PRNG. In the experiment reported above we have, for the sole purpose of convenience, generated the parameters of the n constituent PRNGs using a standard PRNG and a single seed. However, in actual applications the parameters of the n PRNGs may all be individually specified, if desired, thus imposing a large set of unkown values for the analyst to infer in case n is chosen to be correspondingly large. Further note that the computing cost factor of the compound PRNG is independent of n, being a constant 2 with respect to one single constituent PRNG, thus permitting liberal choice of very large values of n with the purpose to defeat analysis and/or to obtain extremely long period length. (For each number output by the compound PRNG exactly 2 calls of one of the n constituent PRNGs are involved.)

  3. While the n constituent PRNGs used in the experiment reported above are based on linear congruential relations (see the program code below) the bits of the integer-valued numbers obtained from the congruential relations are not directly used. Rather, these integers are divided by the respective moduli of the congruential relations to obtain real-valued numbers in the interval [0,1). It is the bits of the 24-bit integers subsequently obtained through multiplication by 2^24 that finally constitute the pseudo-random bits output by our scheme. In other words, the integers output from the different congruential relations are mapped to one and the same common range [0,2^24-1]. The mapping effected is thus different for outputs stemming from different constituent PRNGs because the moduli of the congruential relations are chosen to be different by construction. This normalization process tends to lead to a qualitatively much better bit sequences represented by (z_i) than the n bit sequences contained in the integers directly obtained from the congruential relations underlying the n constituent PRNGs. It constitutes one of the essential factors contributing to the success of our scheme, the other factor being the random interlacing or intermixing of the output streams of the n constituent PRNGs, a fact which is explicit in the algorithm of the compound PRNG presented above.
  4. Using a somewhat optimized program written in Fortran 90, the time for generating 10 million bytes of pseudo-random bits using our scheme on a 200 MHZ PC was found to be 6.4 sec.
  5. The computer program used in the reported experiment is attached below.
References

[1] U. M. Maurer, A Universal Statistical Test for Random Bit
    Generators, J. Cryptology (1992) 5:89-105.

[2] http://www.stud.uni-muenchen.de/~mok-kong.shen/#problem2

[3] http://www.stud.uni-muenchen.de/~mok-kong.shen/#paper1

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       program ranbitgen1
c pseudo-random bit generator
c Maurer's test with L=8, V=256, rho=0.01  
c symbol M is used for L in this program
c for details of the compound PRNG being used see
c http://www.stud.uni-muenchen.de/~mok-kong.shen/#problem2
       implicit integer(a-z)
       parameter(maxint=2147483647)
       parameter(maxn=1000)
       parameter(m=8,v=256,mqk=2500000)
       integer a(0:maxn),b(0:maxn),c(0:maxn),x(0:maxn)
       common /cc/a,b,c,x
       real initgen,g
       common /cs/seed
       integer tab(0:v-1),block(1:mqk),bymask(3),bysh(3)
       double precision sum,ftu,zz,maxftu,minftu,sftu,meanftu
       real t1,t2
       common /cr/t1,t2
c
c building masks and numbers of shifts
       h=255
       sh=0
       do i=3,1,-1
         bymask(i)=h
         h=ishft(h,8)
         bysh(i)=sh
         sh=sh-8
       end do
c
c
 100   print *,'input seed [1, 2147483646] (for the standard PRNG)'
       read *,seed
       if (seed .lt. 1 .or. seed .gt. maxint-1) then
         print *,'input values invalid'
         goto 100
       end if
c
c
 200   print *,'input n [1, 1000] (no. of constituent generators)'
       read *,n
       if (n .lt. 1 .or. n .gt. maxn) then
         print *,'input values invalid'
         goto 200
       end if
c
c
 300   print *,'input q,k'
       read *,q,k
       qk=q+k
       if (q .lt. 10*v .or. qk .gt. mqk) then
         print *,'invalid input'
         goto 300
       end if
c compute threshold values t1, t2 of the test parameter ftu
       call crit(m,k)
c
c
 400   print *,'input number of runs'
       read *,nftu
       base=2**24
       minftu=1.0d100    
       maxftu=0.0d0
       sftu=0.0d0
       do ff=1,nftu
c
c taken from program listing of problem 2
c setting up the n constituent generators of the compound PRNG
c initgen delivers a pseudo-random number in [0,1).
       minc=1000000
       maxc=10000000
       diffc=maxc-minc+1
       do j=0,n-1
 20      c(j)=minc+initgen()*diffc
         cj1=c(j)-1
         b(j)=1+initgen()*cj1
         maxa=(maxint-b(j))/cj1
         a(j)=initgen()*(maxa+1)
c redo if a(j) is too small
         if (a(j) .lt. 10) goto 20
         x(j)=1+initgen()*cj1
       end do
       j=0
c
c
       kj=0
       byi=3
       do
         kj=kj+1
         if (kj .gt. qk) exit
         byi=byi+1
         if (byi .gt. 3) then
c 
c the value in [0,1) obtained from g(j) is mutliplied by 2^24
           bu=base*g(j)
c update the index j. See the algorithm in problem 2 
           j=g(j)*n
c
           byi=1
         end if
         block(kj)=ishft(iand(bu,bymask(byi)),bysh(byi))
       end do
c
c begin of Maurer's algorithm for computing the test parameter ftu
c we use double precision computations to reduce rounding errors
       do i=0,v-1
         tab(i)=0
         end do
       do p=1,q
         tab(block(p))=p
       end do
       sum=0.0d0
       do p=q+1,qk
         zz=p-tab(block(p))
         sum=sum+log(zz)
         tab(block(p))=p
       end do
       ftu=(sum/k)/log(2.0d0)
c
c end of Maurer's algorithm
c
       if (ftu .lt. minftu) minftu=ftu
       if (ftu .gt. maxftu) maxftu=ftu
       sftu=sftu+ftu
c
       end do
c
       meanftu=sftu/nftu
c t1 and t2 are for rho=0.01
       print *,'meanftu : ',real(meanftu),'      t1, t2 : ',t1,t2
       print *,'minftu  : ',real(minftu),'     maxftu  : ',real(maxftu)
c
 500   print *,'input 1 for new seed'
       print *,'      2 for new n'
       print *,'      3 for new q, k'
       print *,'      4 for new number of runs'
       print *,'      5 for termination'
       read *,kk
       if (kk .lt. 1 .or. kk .gt. 5) goto 500
       goto (100,200,300,400,900) kk
 900   end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c computation of threshold values of ftu for rejection rate rho=0.01
c see Maurer's paper in J. Cryptology (1992) 5:89-105
       subroutine crit(m,k)
       real t1,t2
       common /cr/t1,t2
       real ef(16),vr(16)
       data ef/0.7326495,1.5374383,2.4016068,3.3112247,4.2534266,
     1         5.2177052,6.1962507,7.1836656,8.1764248,9.1723243,
     2         10.170032,11.168765,12.168070,13.167693,14.167488,
     3         15.167379/
       data vr/0.690,1.338,1.901,2.358,2.705,2.954,3.125,3.238,
     1         3.311,3.356,3.384,3.401,3.410,3.416,3.419,3.421/
c rho=0.01
       y=2.58
       c=0.7-0.8/m+(1.6+12.8/m)*k**(-4.0/m)
       s=c*sqrt(vr(m)/k)
       t1=ef(m)-y*s
       t2=ef(m)+y*s
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c taken from program listing of problem 2
       real function initgen()
       implicit integer(a-z)
       common /cs/seed
       real random
       hi=seed/127773
       lo=mod(seed,127773)
       test=16807*lo-2836*hi
       if (test .gt. 0) then
         seed=test
       else
         seed=test+2147483647
       end if
       random=real(seed)/2147483647
c Avoid rounding to 1.0.
       if (random .ge. 0.999999) random=0.999999
c Return random to caller.
       initgen=random
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c taken from program listing of problem 2. maxn is here 1000
       real function g(j)
       implicit integer(a-z)
       parameter(maxn=1000)
       integer a(0:maxn),b(0:maxn),c(0:maxn),x(0:maxn)
       common /cc/a,b,c,x
       real rr
c
c Linear congruential generator of the mixed type, i.e.
c          x[new] = a * x[old] + b   (mod c)
c
       x(j)=mod(a(j)*x(j)+b(j),c(j))
c Obtaining a real number in [0,1).
       rr=real(x(j))/c(j)
c Avoid rounding to 1.0.
       if (rr .gt. 0.999999) rr=0.999999
c Return rr to caller.
       g=rr
       end

The World Blackjack Portal!