c****Makewhat.f: program to make random number files for DIEHARD
      print 123
123   format(
     &' This program makes a file of random integers ',/,
     &' for tests by DIEHARD.  Select one from this list:',/)
       print 124
124    format(
     &'  1.  A multiply-with-carry (MWC) generator x(n)=a*x(n-1)+carry
     &mod 2^32',/,
     &'  2.  A MWC generator on pairs of 16 bits',/,
     &'  3.  The "Mother of all random number generators"',/,
     &'  4.  The KISS  generator',/,
     &'  5.  The simple but very good generator COMBO',/,
     &'  6.  The lagged Fibonacci-MWC combination ULTRA',/,
     &'  7.  A combination MWC/subtract-with-borrow (SWB) generator, per
     &iod ~ 10^364',/,
     &'  8.  An extended congruential generator',/,
     &'  9.  The Super-Duper generator',/,
     &' 10.  A subtract-with-borrow generator',/,
     &' 11.  Any specified congruential generator',/,
     &' 12.  The 31-bit generator ran2 from Numerical Recipes',/,
     &' 13.  Any specified shift-register generator, 31 or 32 bits',/,
     &' 14.  The system generator in Sun Fortran f77',/,
     &' 15.  Any lagged-Fibonacci generator,  x(n)=x(n-r) op x(n-s)',/,
     &' 16.  An inverse congruential generator',/)
        open(4,file='makef.txt')
        print*,'  Enter your choice, 1 to 16:'
        read*,jch
       if(jch.eq.1 ) call MAKEMWC1
       if(jch.eq.2 ) call MAKE1616
       if(jch.eq.3 ) call MAKEMTHR
       if(jch.eq.4 ) call MAKEKISS
       if(jch.eq.5 ) call MAKECMBO
       if(jch.eq.6 ) call MAKELTRA
       if(jch.eq.7 ) call MAKESBMC
       if(jch.eq.8 ) call MAKEXCNG
       if(jch.eq.9 ) call MAKESUPR
       if(jch.eq.10) call MAKESWB
       if(jch.eq.11) call MAKECONG
       if(jch.eq.12) call MAKERAN2
       if(jch.eq.13) call MAKESHRG
       if(jch.eq.14) call MAKESUNR
       if(jch.eq.15) call MAKEFIBO
       if(jch.eq.16) call MAKEINVC
       end

        subroutine makemwc1
       implicit real*8 (a-h,o-z)
       integer*4 n(4096),x(2),a(2),w(4),c(4),t,b,xx,cc,z(4)
       character*80 text(22)
       character*15 filename
       character*1 op
        t(xx)=rshift(xx,16)
        b(xx)=and(xx,65535)
       read(4,845) text
       print 845,text
845    format(a78)
       print*,'  Pause:  enter any letter or number key to continue:'
679    format(a1)
       read 679,op
       print*,'  Enter filename (<=15 characters) for output:'
       read 2991,filename
2991   format(a15)
       open(2,file=filename,form='unformatted',access='direct',
     & recl=16384)
       print*,'      Select multiplier `a` from this list:'
       print*,'-------------------------------------------------------'
       print*,' 1791398085 1929682203 1683268614 1965537969 1675393560'
       print*,' 1967773755 1517746329 1447497129 1655692410 1606218150'
       print*,' 2051013963 1075433238 1557985959 1781943330 1893513180'
       print*,' 1631296680 2131995753 2083801278 1873196400 1554115554'
       print*,'-------------------------------------------------------'
       print 281
281   format(
     &' (In general, for any choice of `a`, let m=a*2^32-1. If both m',/
     &,' and (m-1)/2 are prime then the period will be (m-1)/2).')
       print*,'   Enter multiplier a:'
       read*,aa
        print*,'  Enter a seed integer x and initial carry c:'
        read*,xx,cc
        print*,'       Please wait..............'
         x(1)=b(xx)
         x(2)=t(xx)
         c(1)=b(cc)
         c(2)=t(cc)
         c(3)=0
         c(4)=0
         a(1)=b(aa)
         a(2)=t(aa)
          jk=0
           do 2 i=1,700
          jk=jk+1
              do 3 j=1,4096
              call prod(x,a,z)
              call sum(z,c,w)
              x(1)=w(1)
              x(2)=w(2)
              c(1)=w(3)
              c(2)=w(4)
              n(j)=lshift(x(2),16)+x(1)
3             continue
2        write(2,rec=jk) n
        print*,'            FINISHED'
      print 602,filename
602   format('  2,867,200 32-bit random integers (11,468,800 bytes)',/
     &,'  have been written to the file ',a15)
        print*,'  The period is', aa*2d0**31-1d0
        end

       subroutine make1616
       integer*4 n(4096),a,b,x,y
       character*80 text(27),dum
       open(2,file='mwc1616.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,846) (dum,j=1,22 )
       read(4,846) text
       print 846,(text(j),j=1,21)
846    format(a78)
       pause
       print 846,(text(j),j=22,27)
       print*,'      Select multipliers a and b, a<>b, from this list:'
       read(4,846) (text(j),j=1,10)
       print 846, (text(j),j=1,10)
       print*,'      Enter a and b: (my favorites are 18000 and 30903)'
       read*,a,b
        print*,' Enter two (<=31 bit) seed integers, not zero'
        read*,x,y
        print*,'       Please wait..............'
           jk=0
           do 2 i=1,700
           jk=jk+1
              do 3 j=1,4096
              x=a*and(x,65535)+rshift(x,16)
              y=b*and(y,65535)+rshift(y,16)
3             n(j)=lshift(x,16)+and(y,65535)
2        write(2,rec=jk) n
        print*,'            FINISHED'
        print*,'  2,867,200 32-bit random integers (11,468,800 bytes)'
        print*,'  have been written to the file mwc1616.32          '
        print*,'  They are the concatenated output of two 16-bit'
        print*,'  multiply-with-carry generators.  The period is'
        print 800, (a*2d0**15-1d0)*(b*2d0**15-1d0)
800     format(f40.0)
        print*,'  (the last 4 digits may be lost to roundoff).'
801     format('  The correct last 4 digits are',i8)
        print 801,mod(mod(a*32768-1,10000)*mod(b*32768-1,10000),10000)
        return
        end

       subroutine makemthr
       implicit real*8 (a-h,o-z)
       integer*4 n(4096), x,y,z,w,v,a,b,c,d,carry
       character*80  text(23),dum
       parameter(dm=2d0**32)
       data a,b,c,d/2111111111,1492,1776,5115/
       data da,db,dc,dd/2111111111d0,1492d0,1776d0,5115d0/
       read(4,845) (dum,j=1,59)
       read(4,845) text
       print 845,text
845    format(a78)
745    format(a15)
      open(2,file='mthr.32',form='unformatted',access='direct',
     & recl=16384)
        print*,' Enter four seed integers:'
        read*,x,y,z,w
        print*,'       Please wait..............'
      jk=0
      do 3 ij=1,700
      jk=jk+1
      do 2 i=1,4096
      if(x.lt.0) then
      dx=x+dm
      else
      dx=x
      endif
      if(y.lt.0) then
      dy=y+dm
      else
      dy=y
      endif
      if(z.lt.0) then
      dz=z+dm
      else
      dz=z
      endif
      if(w.lt.0) then
      dw=w+dm
      else
      dw=w
      endif
      v=a*x+b*y+c*z+d*w+carry
      x=y
      y=z
      z=w
      w=v
      dv=da*dx+db*dy+dc*dz+dd*dw+carry
      carry=dint(dv/dm)
2     n(i)=w
3     write(2,rec=jk) n
        print*,'              '
        print*,'            FINISHED'
      print 602,filename
602   format('  2,867,200 32-bit random integers (11,468,800 bytes)',/
     &,'  have been written to the file ',a15)
      print*,'  The period is about 2^158.97'
        return
        end

       subroutine makekiss
       integer*4 b(4096),x,y,z,w,carry
       character*80 text(22),dum
       open(2,file='kiss.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,845) (dum,j=1,82)
       read(4,845) text
       print 845,text
845    format(a78)
        print*,' Enter four seed integers, not zero'
        read*,x,y,z,w
        carry=0
        print*,'       Please wait..............'
        jk=0
           do 2 i=1,700
           jk=jk+1
              do 33 jj=1,4096
              x=69069*x+1
              y=xor(y,lshift(y,13))
              y=xor(y,rshift(y,17))
              y=xor(y,lshift(y,5))
              k=rshift(z,2)+rshift(w,3)+rshift(carry,2)
              m=w+w+z+carry
              z=w
              w=m
              carry=rshift(k,30)
33              b(jj)=x+y+w
2        write(2,rec=jk) b
        print*,'            FINISHED'
        print*,'  2,867,200 32-bit random integers (11,468,800 bytes)'
        print*,'  have been written to the file kiss.32          '
        return
        end
       subroutine makecmbo
c****Simple combo, period> 2^60.5
c***x(n)=x(n-1)*x(n-2) mod 2^32 added to
c*****period of x(n)=x(n-1)*x(n-2) is 3*2^29 if seeds odd, and one is +or-3 mod
c*****easy to ensure: replace seed x with 3*x*x.
c****mwc z=30903*iand(z,65535)+ishft(z,-16)
        implicit integer(a-z)
        integer b(4096)
        character*80 text(16),dum
        read(4,845) (dum,j=1,104)
        read(4,845) text
845     format(a78)
        print 845,text
        print*,'              '
        print*,'  Enter three seed integers:'
        read*,x,y,z
        x=x+x+1
        x=3*x*x
        y=y+y+1
       open(1,file='combo.32',form='unformatted',access='direct',
     & recl=16384)
         print*,'   Please wait..................'
        jk=0
        do 3 i=1,700
        jk=jk+1
        do 2 j=1,4096
        v=x*y
        x=y
        y=v
        z=30903*and(z,65535)+rshift(z,16)
2       b(j)=y+z
3       write(1,rec=jk) b
        print*,'            FINISHED'
        print*,'  2,867,200 32-bit random integers (11,468,800 bytes)'
        print*,'  have been written to the file combo.32          '
        return
        end
       subroutine makeltra
      integer*4 bb(4096), r,s,ju(99)
      real xx(607)
      character*80 text(15),dum
      data r,s/99,33/
      read(4,845) (dum,j=1,120)
      read(4,845) text
      print 845,text
845   format(a78)
       open(1,file='ultra.32',form='unformatted',access='direct',
     & recl=16384)
       print*, '   Enter four positive integers for seeds:'
        read*,i,j,k,l
        mwcs=i+j+k+l
           DO 92 II=1,r
            i=18273*and(i,65535)+rshift(i,16)
            j=23163*and(j,65535)+rshift(j,16)
            k=24984*and(k,65535)+rshift(k,16)
            l=28854*and(l,65535)+rshift(l,16)
           js=lshift(i,16)+and(j,65535)+lshift(k,16)+and(l,65535)
           js=or(js,1)
           xx(ii)=.5+js*.5**32
92    jU(II)=jS
       print*,'    The seed table:'
       print 25, (jU(ig),ig=1,r)
25     format(5i12)
       call kstest(xx,r,px)
       print 271,px
271    format(' p-value for KSTEST on uniformity of seeds: ',f8.6)
        print*,'    Please wait................'
         IP=r
         JP=s
             jk=0
             do 2 i=1,700
             jk=jk+1
             do 3 j=1,4096
                jUNI=jU(IP)*jU(JP)
             jU(IP)=jUNI
            IP=IP-1
           IF(IP.EQ.0) IP=97
            JP=JP-1
            IF(JP.EQ.0) JP=97
            mwcs=30903*and(mwcs,65535)+rshift(mwcs,16)
3           bb(j)=juni+mwcs
2      write(1,rec=jk) bb
       print*,'     '
       print*,'     FINISHED'
      print*,' 2,867,200  32-bit random integers (11,468,800 bytes)'
      print*,'  have been written to the file ultra.32  '
       print*,'::::::::::::::::::::::::::::::::::::::::::::::'
        return
        end
      subroutine makesbmc
      integer*4 bb(4096), r,s,ju(37),x,y
      character*80 text(9),dum
      real xx(37)
      data r,s/37,24/
      read(4,845) (dum,j=1,135)
      read(4,845) text
      print 845,text
845   format(a78)
      print*,'    Enter four positive integers for seeds:'
        read*,i,j,k,l
      open(1,file='swbmwc.32',form='unformatted',access='direct',
     & recl=16384)
        print*,'    Please wait................'
         mwcs=i+j+k+l
           DO 92 II=1,r
           jS=0.
           DO 93 JJ=1,32
           M=MOD(MOD(I*J,179)*K,179)
           I=J
           J=K
           K=M
           L=MOD(53*L+1,169)
           jS=2*jS
           IF(MOD(L*M,64).GE.32) jS=jS+1
93    continue
      xx(ii)=.5+.5**32*jS
 92    jU(II)=jS
       print*,'    The seed table:'
       print 25, (jU(ig),ig=1,r)
25     format(5i12)
       call kstest(xx,r,px)
       print 271,px
271    format(' p-value for KSTEST on uniformity of seeds: ',f8.6)
         IP=r
         JP=s
          jk=0
             do 2 i=1,700
             jk=jk+1
             do 3 j=1,4096
             x=ju(ip)
             y=ju(jp)
             juni=x-y
           if(rshift(x,1).gt.rshift(y,1)) juni=juni-1
             jU(IP)=jUNI
            IP=IP-1
            IF(IP.EQ.0) IP=r
            JP=JP-1
            IF(JP.EQ.0) JP=r
             mwcs=30903*and(mwcs,65535)+rshift(mwcs,16)
 3       bb(j)=mwcs+juni
2      write(1,rec=jk) bb
       print*,'     '
       print*,'     FINISHED'
       print 349
349    format('     The binary file swbmwc.32'
     &/,'     has been created with 32-bit integers'
     &,/,'     from the subtract-with-borrow sequence'
     &,/,'    x(n)=x(n-24)-x(n-37)-borrow mod 2^32'
     &,/,'   combined with the multiply-with-carry sequence'
     &,/,'   y(n)=30903*y(n-1)+carry mod 2^16,'
     &,/,'  the overall sequence having period about 10^364.')
          return
          end

       subroutine makexcng
       integer*4 n(4096)
       real*8 a,b,c,p,x
       character*80 text(17),dum
       parameter(p=2d0**32-5d0,s=2d0**31)
       open(2,file='excong.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,845) (dum,j=1,144)
       read(4,845) text
       print 845,text
845    format(a78)
        print*,' Enter three seed integers, not all zero:'
        read*,a,b,c
        print*,'       Please wait..............'
        r3=2d0**32/34359737519d0
        r4=2d0**32/34359736739d0
       jk=0
      do 3 ij=1,700
      jk=jk+1
      do 2 i=1,4096
         go to(51,52,53,54),jch
51         x=dmod(1776d0*a+1476d0*b+1176d0*c,4294967291d0)
           go to 166
52         x=dmod(8192*(a+b+c),4294967291d0)
           go to 166
53         x=dmod(2001d0*a+1998d0*b+1995d0*c,34359737519d0)
           x=x*r3
           go to 166
54         x=dmod(524288d0*(a+b+c),34359736739d0)
           x=x*r4
166   a=b
      b=c
      c=x
      if(x.lt.s) then
      n(i)=x
      else
      n(i)=x-2d0**32
      endif
2      continue
3      write(2,rec=jk) n
        print*,'              '
        print*,'            FINISHED'
        print*,'  2,867,200 32-bit random integers (11,468,800 bytes)'
        print*,'  have been written to the file excong.32          '
        return
        end
       subroutine makesupr
       integer*4 n(4096),x,y
       character*80 text(20),dum
       character*1 op
       open(2,file='suprdupr.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,845) (dum,j=1,161)
       read(4,845) text
       print 845,text
845    format(a78)
        print*,' Enter two  seed integers, the second not zero:'
        read*,x,y
        print*,'  Choose your method of combination,'
        print*,'  + for addition, x for exclusive-or, in column 1:'
        read 821,op
821     format(a1)
        print*,'       Please wait..............'
        jk=0
        if(op.eq.'+') then
           do 2 i=1,700
           jk=jk+1
              do 3 j=1,4096
              x=69069*x+1
              y=xor(y,lshift(y,13))
              y=xor(y,rshift(y,17))
              y=xor(y,lshift(y,5))
3             n(j)=x+y
2        write(2,rec=jk) n
         end if
         if(op.eq.'x') then
           do 72 i=1,700
           jk=jk+1
              do 73 j=1,4096
              x=69069*x+1
              y=xor(y,lshift(y,13))
              y=xor(y,rshift(y,17))
              y=xor(y,lshift(y,5))
73             n(j)=xor(x,y)
72        write(2,rec=jk) n
         endif
        print*,'            FINISHED'
        print*,'  2,867,200 32-bit random integers (11,468,800 bytes)'
        print*,'  have been written to the file suprdupr.32  '
        return
        end
       subroutine  makeswb
       integer*4 bb(4096), r,s,ju(48),rc(5),sc(5),x,y,jix(5)
      character*15 filename
      character*80 text(17),dum
      character*42 swb(5),period(5)
      data rc,sc/43,37,24,21,48,22,24,19,6,8/
      data jix/-6,0,0,0,2147483647/
       read(4,845) (dum,j=1,181)
       read(4,845) text
       print 845,text
845    format(a78)
       print*,'    Choose your generator from this list:'
       print 945,swb
945    format(10x,a42)
       data swb/
     &'   1:   x(n)=x(n-43)-x(n-22)-c mod 2^32-5',
     &'   2:   x(n)=x(n-37)-x(n-24)-c mod 2^32  ',
     &'   3:   x(n)=x(n-24)-x(n-19)-c mod 2^32  ',
     &'   4:   x(n)=x(n-21)-x(n- 6)-c mod 2^32  ',
     &'   5:   x(n)=x(n-48)-x(n- 8)-c mod 2^31  '/
      print*,'(Choice 5 provides 31-bit integers that are left-justified
     &.)'
       print*,' enter 1,2,3,4 or 5:'
       read*, jchoice
       r=rc(jchoice)
       s=sc(jchoice)
        print*,'     Enter output file name (<=15 characters):'
        read 645,filename
645     format(a15)
      open(1,file=filename,form='unformatted',access='direct',
     & recl=16384)
        print 348
348     format('    Enter four positive integers for seeds:')
        read*,i,j,k,l
        print*,'    Please wait................'
              do 90 ijk=1,i
c 90              jjn=1
90            jjn=irand(1)
           DO 92 II=1,r
           jS=0.
           DO 93 JJ=1,32
           M=MOD(MOD(I*J,179)*K,179)
           I=J
           J=K
           K=M
           L=MOD(53*L+1,169)
           jS=2*jS
           IF(MOD(L*M,64).GE.32) jS=jS+1
93    continue
      if(jchoice.eq.5) js=xor(js,1)
92    jU(II)=jS
         IP=r
         JP=s
         jk=0
             do 2 i=1,700
             jk=jk+1
             do 3 j=1,4096
             x=ju(ip)
             y=ju(jp)
             juni=x-y
           if(rshift(x,1).lt.rshift(y,1)) juni=juni+jix(jchoice)
             jU(IP)=jUNI
            IP=IP-1
           IF(IP.EQ.0) IP=r
            JP=JP-1
            IF(JP.EQ.0) JP=r
            if(jchoice.eq.5) juni=juni+juni
3           bb(j)=juni
2      write(1,rec=jk) bb
       print*,'     '
       print*,'     FINISHED'
       print 349, filename,swb(jchoice)
349    format('     The binary file ',a15,
     &/,'     has been created with 32-bit integers'
     &,/,'     from the subtract-with-borrow sequence',/,a41)
       print 350,period(jchoice)
350    format(15x, 'The period is',/,a41)
       data period/
     &'  (2^32-5)^43 - (2^32-5)^22, about 10^414',
     &'  2^1178 - 2^762, about 10^354           ',
     &'  (2^759 - 2^599)/3, about 10^228        ',
     &'  (2^666 - 2^186)/3, about 10^200        ',
     &'  (2^1478 - 2^247)/105, about 10^445     '/
       print*,'+++++++++++++++++++++++++++++++++++++++++++++++++++++++'
       return
       end
      subroutine makecong
      implicit real*8(a-h,o-z)
      integer*4 bb(4096), a,b,r,s
      real*4 x(2000),y(2000)
       character*80 text(13)
       character*1 op
       character*15 filename
       read(4,845) (dum,j=1,198)
       read(4,845) text
       print 845,text
845    format(a78)
        print*,' Enter a,b,r and s for the generator'
        print*,'  x(n)=a*x(n-1)+b mod 2^r+s'
      print*,' To avoid overflow, make sure a*2^r < 2^53) '
        read*,a,b,r,s
        print*,' Enter a seed integer:'
        read*,jseed
        da=a
        db=b
        dm=2d0**r+s
        print*,' If you want the 2-lattice without generating'
      print*,' the binary file, enter 0, else enter 1:'
      print*,'   Create binary file? 0 for NO, 1 for YES:'
      read*, jlat
            if (jlat.eq.1) then
            print*,'     Enter output file name (<=15 characters):'
            read 645,filename
645         format(a15)
      open(1,file=filename,form='unformatted',access='direct',
     & recl=16384)
         print*,'   Please wait................'
            dj=jseed
            dr=2d0**(32-r)
             jk=0
            do 2 i=1,700
             jk=jk+1
               do 3 j=1,4096
               dj=dmod(da*dj+db,dm)
               dx=dr*dj
                   if(dx.lt.2d0**31) then
                   bb(j)=dx
                   else
                   bb(j)=dx-2d0**32
                   endif
3             continue
2           write(1,rec=jk) bb
            endif
        print 221,filename
221   format('  2,867,200 32-bit random integers (11,468,800 bytes)',/,
     &'  have been written to the file ',a15)
         print*,' To display the 2-lattice of this generator, '
         print*,'  hit any letter or number key:'
          read 983,op
 983     format(a1)
         dL=2d0**(20+r-32)
         kount=0
         klim=da*dL/dm
            do 29 k=1,klim
            j1=max(0,(k*dm-b)/da)
            j2=(k*dm-b+dL)/da
                 do 39 j=j1,j2
                 yy=dmod(da*j+db,dm)
                 if(yy .lt. dL)  then
                 kount=kount+1
                 x(kount)=j/dm
                 y(kount)=yy/dm
                 endif
39                continue
29            continue
          call plot1(x,y,kount,54)
      return
      end
       subroutine makeran2
       integer*4 b(4096)
       character*80 text(16),dum
       open(2,file='ran2.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,845) (dum,j=1,211)
       read(4,845) text
       print 845,text
845    format(a78)
        print*,' enter a seed integer'
        read*,jseed
          idum=-iabs(jseed)
          jk=0
           do 3 ij=1,700
          jk=jk+1
       do 2  i=1,4096
2       b(i)=lshift(iran2(idum),1)
3        write(2,rec=jk) b
        print*,'  Now 2,867,200 left-adjusted integers from ran2'
        print*,'  (11,486,800 bytes) have been written to ran2.32'
        return
        end

      FUNCTION iran2(idum)
      INTEGER idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV
      REAL AM,EPS,RNMX
      PARAMETER (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1,
     *IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791,
     *NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER idum2,j,k,iv(NTAB),iy
      SAVE iv,iy,idum2
      DATA idum2/123456789/, iv/NTAB*0/, iy/0/
      if (idum.le.0) then
        idum=max(-idum,1)
        idum2=idum
        do 11 j=NTAB+8,1,-1
          k=idum/IQ1
          idum=IA1*(idum-k*IQ1)-k*IR1
          if (idum.lt.0) idum=idum+IM1
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ1
      idum=IA1*(idum-k*IQ1)-k*IR1
      if (idum.lt.0) idum=idum+IM1
      k=idum2/IQ2
      idum2=IA2*(idum2-k*IQ2)-k*IR2
      if (idum2.lt.0) idum2=idum2+IM2
      j=1+iy/NDIV
      iy=iv(j)-idum2
      iv(j)=idum
      if(iy.lt.1)iy=iy+IMM1
       iran2=iy
c      ran2=min(AM*iy,RNMX)
      return
      END
        subroutine makeshrg
c*** creates shift register random number file
c  use 13,18 or 7,24 or 6,25 or 3,28 or reverse for 31 bits
c use 15,17 or 13,17,5 for 32 bits
      integer b(4096),L,R
      character*80 text(12),dum
      character*15 filename
      m(k,n)=xor(k,lshift(k,n))
      mr(k,n)=xor(k,rshift(k,n))
       read(4,845) (dum,j=1,227)
       read(4,845) text
       print 845,text
845    format(a78)
      print*,'  Enter output file name (<= 15 characters):'
      read 645,filename
645   format(a15)
      open(1,file=filename,form='unformatted',access='direct',
     & recl=16384)
      print*,'   Choose number of bits, 31 or 32:'
      read*,nbits
      jk=0
      if(nbits.eq.31) then
      print*,'  Choose Left,Right shifts from these choices:'
      print*,'    13,18  18,13  24,7  7,24 6,25  25,6  28,3  3,28'
      print*,'  enter two integers, L and R, in free format:'
      read*,L,R
      mask=2**31-1
      print*,' Enter seed integer, not zero:'
      read*,j
        print*,'  Please wait..............'
          do 2 i1=1,700
           jk=jk+1
             do 3 i=1,4096
              j=mr(and(m(j,L),mask),R)
3              b(i)=j+j
2         write(1,rec=jk) b
      endif
      if(nbits.eq.32) then
      print*,' For 32 bit integers, you have three choices.'
      print*,'  two shifts: (L,R)=(17,15) or (15,17) and'
      print*,'  three shifts (L1,R,L2)=(13,17,5).'
      print*,' Enter three integers, 17 15 0 or 15 17 0 or 13 17 5:'
      read*,L,R,L2
      print*,' Enter seed integer, not zero:'
      read*,j
        print*,'  Please wait..............'
          if(L2.eq.0) then
          do 4 i1=1,700
           jk=jk+1
               do 5 i=1,4096
               j=mr(m(j,L),R)
5              b(i)=j
4         write(1,rec=jk) b
          elseif (L2.eq.5)then
          do 6 i1=1,700
           jk=jk+1
               do 7 i=1,4096
               j=m(mr(m(j,L),R),L2)
7              b(i)=j
6         write(1,rec=jk) b
          endif
        endif
        print*,'      FINISHED'
        print 221,filename
221   format('  2,867,200 32-bit random integers (11,468,800 bytes)',/,
     &'  have been written to the file ',a15)
       if(nbits.eq.31 .or. L2.eq.5) print 225,nbits
225   format(' The period of your generator is 2^',i2,'-1.')
         if(L.eq.15 .or. L.eq.17) print 226
226   format(' The period of your generator is 2^32-2^21-2^11+1.')
          return
          end
       subroutine makesunr
       integer*4 b(4096)
       character*80 text(18),dum
       open(2,file='sunran.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,845) (dum,j=1,239)
       read(4,845) text
       print 845,text
845    format(a78)
        print*,' Enter one  seed integer:'
        read*,jseed
        print*,'       Please wait..............'
       ijk=irand(jseed)
         ijk=1
        jk=0
       do 3 j=1,700
       jk=jk+1
       do 2 i=1,4096
c 2           b(i)=1
2      b(i)=irand(1)
3      write(2,rec=jk) b
        print*,'            FINISHED'
        print*,' 2,867,200  32-bit random integers (11,468,800 bytes)'
        print*,'  have been written to the file sunran.32  '
        print*,'  Note: the rightmost bit is always 0'
        return
        end
        subroutine makefibo
        integer*4 bb(4096), r,s,op,ju(607)
        character*80 text(9),dum
      real*8 period
      character*15 filename
       read(4,845) (dum,j=1,257)
       read(4,845) text
       print 845,text
845    format(a78)
          print 346
346   format(' Enter lags r and s from this list:',/,
     &'    17,5  33,13  39,14  55,24  63,31  73,25  97,33 607,273',/,
     &'     and op code: 1 for + , 2 for -, 3 for *, 4 for xor:')
      read*,r,s,op
            print*,'     Enter output file name (<=15 characters):'
            read 645,filename
645         format(a15)
      open(1,file=filename,form='unformatted',access='direct',
     & recl=16384)
        print 348
348     format('    Enter four positive integers for seeds:')
        read*,i,j,k,l
        print*,'    Please wait................'
              do 90 ijk=1,i
90            jjn=irand(1)
           DO 92 II=1,97
           jS=0.
           DO 93 JJ=1,32
           M=MOD(MOD(I*J,179)*K,179)
           I=J
           J=K
           K=M
           L=MOD(53*L+1,169)
           jS=2*jS
           IF(MOD(L*M,64).GE.32) jS=jS+1
93    continue
      if(op.eq.3) js=xor(js,1)
92    jU(II)=jS
         IP=r
         JP=s
         jk=0
             do 2 i=1,700
             jk=jk+1
             do 3 j=1,4096
                if(op.eq.1) then
                jUNI=jU(IP)+jU(JP)
                else if(op.eq.2) then
                jUNI=jU(IP)-jU(JP)
                else if(op.eq.3) then
                jUNI=jU(IP)*jU(JP)
                else
                jUNI=xor(jU(IP),jU(JP))
                endif
             jU(IP)=jUNI
            IP=IP-1
           IF(IP.EQ.0) IP=97
            JP=JP-1
            IF(JP.EQ.0) JP=97
3           bb(j)=juni
2      write(1,rec=jk) bb
       period=(2d0**r-1)*2d0**32
       if(op.eq.3) period=.25*period
       if(op.eq.4) period=(2d0**r-1)
       print*,'     '
       print*,'     FINISHED'
       print 349, filename
349    format('     The binary file ',a15
     &,/,'     has been created with 32-bit integers'
     &,/,'     from the lagged-Fibonacci sequence')
       if(op.eq.1) then
       print*,'+++++++++++++++++++++++++++++++++++++++++++++++++++++++'
       print 351,r,s
351    format(15x,'x(n) = x(n-',i3,') + x(n-',i3,') mod 2^32')
       else if(op.eq.2) then
       print 352,r,s
352    format(15x,'x(n) = x(n-',i3,') - x(n-',i3,') mod 2^32')
       else if(op.eq.3) then
       print 353,r,s
353    format(15x,'x(n) = x(n-',i3,') * x(n-',i3,') mod 2^32')
       else
       print 354,r,s
354    format(15x,'x(n) = x(n-',i3,') xor x(n-',i3,')')
       endif
       print 355
355    format('   The period of the sequence is ')
       print*,'     ',period
       return
       end
        subroutine prod(x,y,z)
        integer*4 x(2),y(2),z(4),t,b,r(4),s(4)
        integer*4 d
        t(d)=rshift(d,16)
        b(d)=and(d,65535)
        d=x(1)*y(1)
        z(1)=b(d)
        d=t(d)+x(1)*y(2)
        r(2)=b(d)
        r(3)=t(d)
        d=x(2)*y(1)
        s(2)=b(d)
        d=t(d)+x(2)*y(2)
        s(3)=b(d)
        s(4)=t(d)
        d=r(2)+s(2)
        z(2)=b(d)
        d=t(d)+r(3)+s(3)
        z(3)=b(d)
        z(4)=t(d)+s(4)
        return
        end
        subroutine sum(x,y,z)
        integer*4 x(4),y(4),z(4),t,b
        integer*4 d
        t(d)=rshift(d,16)
        b(d)=and(d,65535)
        d=x(1)+y(1)
        z(1)=b(d)
        d=t(d)+x(2)+y(2)
        z(2)=b(d)
        d=t(d)+x(3)+y(3)
        z(3)=b(d)
        d=t(d)+x(4)+y(4)
        z(4)=b(d)
        return
        end

      subroutine makeinvc
      integer*4 n(4096),a,b,icong
      character*80 text(15),dum
      real*8 b0,b1,b2,a0,a1,a2,dz
      data a,b,dz/69069,362436069,123456789d0/
      b0=4294967296d0
       open(1,file='invcong.32',form='unformatted',access='direct',
     & recl=16384)
       read(4,845) (dum,j=1,266)
       read(4,845) text
       print 845,text
845    format(a78)
      print*,'  Enter a and b and seed integer, free format:'
      read*, a,b,dz
      print*,'    This may take a long time, as this is a very'
      print*,'    very slow generator.  Please wait ...........'
      jk=0
      do 3 ij=1,700
       jk=jk+1
       do 2 i=1,4096
      b1=dz
      a0=0
      a1=1
      b0=4294967296d0
      do while(b1.gt.0)
      b2=dmod(b0,b1)
      a2=a0-a1*((b0-b2)/b1)
      b0=b1
      b1=b2
      a0=a1
      a1=a2
      end do
      dz=dmod(a*b0*a0+b,4294967296d0)
      if(dz.lt.0) dz=dz+4294967296d0
      if(dz.gt.2d0**31) then
      icong=dz-4294967296d0
      else
      icong=dz
      endif
2     n(i)=icong
3     write(1,rec=jk) n
      print 471
 471  format(
     &' The file invcong.32 now contains 11,428,800 bytes of 32-bit',/,
     &' integers from the inverse congruential generator.')
      end

       subroutine asort(list, n)
       real list(n)
       dimension iu(33), il(33)
       m=1
       i=1
       j=n
5       if(i.ge.j) goto 70
10       k=i
       ij= (i+j)/2
       t= list(ij)
       if(list(i).le.t) goto 20
       list(ij)= list(i)
       list(i)= t
       t= list(ij)
20       l=j
       if(list(j).ge.t) goto 40
       list(ij)= list(j)
       list(j)= t
       t= list(ij)
       if(list(i).le.t) goto 40
       list(ij)= list(i)
       list(i)= t
       t= list(ij)
       goto 40
30       list(l)= list(k)
       list(k)=tt
40       l=l-1
       if(list(l).gt.t) goto 40
       tt= list(l)
50       k= k+1
       if(list(k).lt.t) goto 50
       if(k.le.l) goto 30
       if(l-i.le.j-k) goto 60
       il(m)=i
       iu(m)=l
       i=k
       m=m+1
       goto 80
60       il(m)=k
       iu(m)=j
       j=l
       m=m+1
       goto 80
70       m=m-1
       if(m.le.0) return
       i=il(m)
       j=iu(m)
80       if(j-i.ge.11) goto 10
       if(i.eq.1) goto 5
       i=i-1
90       i=i+1
       if(i.eq.j) goto 70
       t= list(i+1)
       if(list(i).le.t) goto 90
       k=i
100       list(k+1)=list(k)
       k=k-1
       if(t.lt.list(k)) goto 100
       list(k+1)=t
       goto 90
       end
      SUBROUTINE KSTEST(Y,N,P)
C      TO TEST WHETHER A SET OF N REAL NUMBERS IS DRAWN
C      FROM A UNIFORM DISTRIBUTION (KOLMOROGOV-SMIRNOV METHOD)
C      THE TEST IS BASED ON THE DISTANCE BETWEEN THE EMPIRICAL
C      AND THEORETICAL DISTRIBUTION FUNCTIONS
C       USAGE: CALL KSTEST(Y,N,P)
C      Y ...   ARRAY OF REAL NUMBERS HYPOTHETICALLY DRAWN
C              FROM A UNIFORM DISTRIBUTION ON (0,1)
C      N ...   NUMBER OF ELEMENTS IN 'Y'
C      P IS THE PROBABILITY ASSOCIATED WITH THE OBSERVED VALUE
C      OF THE ANDERSON-DARLING STATISTIC: N TIMES THE INTEGRAL
C      OF (FN(X)-X)**2/(X*(1-X))
      DIMENSION Y(N)
      integer l(8,10)
      DATA L/40,46,37,34,27,24,20,20,88,59,43,37,29,27,20,
     &22,92,63,48,41,30,30,25,24,82,59,42,37,26,28,26,22,62,48,33,30
     &,23,23,22,18,49,34,22,20,16,17,17,12,17,17,7,8,4,7,5,1,40,18,
     &19,14,16,13,10,9,59,20,10,4,1,1,0,-1,41,43,36,112,15,95,32,58/
      CALL ASORT(Y,N)
      Z=-N*(N+0.)
      DO 2 I=1,N
      T=Y(I)*(1.-Y(N+1-I))
      IF(T.LT.1.E-20)T=1.E-20
  2   Z=Z-(I+I-1)*ALOG(T)
      Z=Z/N
      P=0.
      IF(Z .LT. .01) GO TO 5
      IF(Z .GT. 2.) GO TO 3
      P=2.*EXP(-1.2337/Z)*(1.+Z/8.-.04958*Z*Z/(1.325+Z))/SQRT(Z)
      GO TO 5
  3   IF(Z .GT. 4.) GO TO 4
      P=1.-.6621361*EXP(-1.091638*Z)-.95059*EXP(-2.005138*Z)
      GO TO 5
  4   P=1.-.4938691*EXP(-1.050321*Z)-.5946335*EXP(-1.527198*Z)
  5   M=MIN0(N-2,8)
      E=0.
      DO 6 J=1,10
  6   E=E+L(M,J)*SP(P,J)*.0001
      IF(N .GT. 10) E=10.*E/N
      A=P+E
      RETURN
      END
      FUNCTION SP(X,I)
      SP=0.
      GO TO (7,7,7,7,7,7,7,8,9,10),I
  7   T=ABS(10.*X-.5-I)
      IF(T.GT.1.5) RETURN
      IF(T.LE. .5) THEN
                   SP=1.5-2.*T*T
      ELSE
                   SP=2.25-T*(3.-T)
      ENDIF
      RETURN
  8   IF(X .LE. .8 .OR. X .GE. 1.) RETURN
      SP=100.*(X-.9)**2 -1.
      RETURN
  9   IF(X .LE. 0. .OR. X .GE. .05) RETURN
      IF(X .LE. .01) THEN
                     SP=-100.*X
      ELSE
                     SP=25.*(X-.05)
      ENDIF
      RETURN
  10  IF(X .LE. .98 .OR. X .GE. 1.) RETURN
      SP=.1-10.*ABS(X-.99)
      RETURN
      END

      SUBROUTINE PLOT1(X,Y,N,NC)
      REAL X(1),Y(1),Z(1)
      CHARACTER*1 C(90,54),STAR,BLANK,PLUS,ZERO,F(4860)
      CHARACTER*4860 D,E
      LOGICAL    P2
      EQUIVALENCE (C,D),(E,F)
      PARAMETER(STAR='*',BLANK=' ',PLUS='+',ZERO='0')
      DATA F/4860*' '/
      P2=.FALSE.
      GO TO 1
      ENTRY PLOT2(X,Y,Z,N,NC)
      P2=.TRUE.
  1   NR=48*NC/100
      IF(NC.GT.60) NR=6*NC/10
      XL=X(1)
      XR=XL
      YB=Y(1)
      YT=YB
      DO 2 I=1,N
      XL=AMIN1(XL,X(I))
      XR=AMAX1(XR,X(I))
      IF(P2) YB=AMIN1(YB,Z(I))
      IF(P2) YT=AMAX1(YT,Z(I))
      YB=AMIN1(YB,Y(I))
  2   YT=AMAX1(YT,Y(I))
      PRINT 22,XL,XR,YB,YT
 22   FORMAT('   X RANGE:',G12.3,' TO ',G12.3,'  Y,Z RANGE:',G12.3,'  TO
     + ',G12.3)
      XD=1.
      IF(XL.NE.XR) XD=NC/(XR-XL)
      YD=1.
      IF(YB.NE.YT) YD=NR/(YT-YB)
      D=E
      DO 4 M=1,N
      I=1.+(X(M)-XL)*XD
      J=1.+(Y(M)-YB)*YD
      C(I,J)=STAR
      IF(.NOT. P2) GO TO 4
      K=1.+(Z(M)-YB)*YD
      C(I,K)=PLUS
      IF(J.EQ.K) C(I,K)=ZERO
  4   CONTINUE
      DO 5 J=NR,1,-1
  5   PRINT 21,(C(I,J),I=1,NC)
  21  FORMAT(10X,90A1)
      RETURN
      END