c
c
      subroutine initiran(seqnum, gseed, state)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine initializes the state for the seqnum(th) member        c
c  of the family of pseudorandom sequences.  initiran must             c
c  be called once (and only once) for each sequence from which         c
c  random binary or integer values are wanted.                         c
c                                                                      c
c  L, k, s, and mbits are properties associated with the family        c
c  of generators used by this package.  L and k represent a            c
c  particular primitive trinomial (mod 2), s is the "canonical bit"    c
c  associated with that trinomial, and mbits is the number of          c
c  mantissa bits used by the floating point pseudorandom generator.    c
c  for L = 158, state must have a declared dimension of at least       c
c  4 x 158 = 632 in the calling program.                               c
c                                                                      c
c  After initializing a sequence with initiran, multiple calls         c
c  can be made to iranx or brand64,  and the state vector will         c
c  stay current.   The state vector must not be modified by the user   c
c  code.   If more than one sequence is to be active at any time,      c
c  there must be a searate state vector for each active sequence.      c
c                                                                      c
c  gseed is the integer value of the user-supplied "global seed."      c
c  The global seed is intended to be fixed for the duration of the     c
c  run and should be the same value for every sequence on every        c
c  processor.   It can be used to modify the random data, in a         c
c  predictable way, on a per-run basis.  seqnum must be in the         c
c  range {0,1,2,3,..., 65535 = 2^16 - 1 = z'0000ffff'} and gseed       c
c  must be in the range {0,1,2,3,..., 65534 = 2^16 - 2 = z'0000fffe'}. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer seqnum, gseed, state(4*L)
c----------------------------------------------------------------------c
      integer hi, lo, a, q, r, t31m1, seed, i
      integer treg1, treg2, highbit1, highbit2
      parameter (a = 16807, q = 127773, r = 2836, t31m1 = z'7fffffff')
c
      call tausinit(seqnum, gseed, treg1, treg2)
c
c  initialize first state register:
      seed     = iand(treg1, t31m1)
      highbit1 = ieor(seed, treg1)
      state(L) = 0
      do  i = 1, L-1
         hi        = seed/q
         lo        = seed - (hi*q)
         seed      = a*lo - r*hi
         if (seed.le.0) seed = seed + t31m1
         state(i)  = seed + seed
      enddo
c
c  incorporate the high bit of the taus. reg. in the
c  fibonacci register, in order to maintain uniqueness:
      if (highbit1.eq.0) then
         state(L-1) = iand(state(L-1),t31m1)
      else
         state(L-1) = ior (state(L-1),highbit1)
      endif
      state(s+1)   = ior(state(s+1),1)
c
c  initialize second state register:
      seed       = iand(treg2,t31m1)
      highbit2   = ieor(seed,treg2)
      state(2*L) = 0
      do  i = L+1, 2*L-1
         hi        = seed/q
         lo        = seed - (hi*q)
         seed      = a*lo - r*hi
         if (seed.le.0) seed = seed + t31m1
         state(i)  = seed + seed
      enddo
c
c  incorporate the high bit of the taus. reg. in the
c  fibonacci register, in order to maintain uniqueness:
      if (highbit2.eq.0) then
         state(2*L-1) = iand(state(2*L-1),t31m1)
      else
         state(2*L-1) = ior (state(2*L-1),highbit2)
      endif
      state(L+s+1) = ior(state(L+s+1),1)
      return
      end
c
c
c
c
c
      subroutine initfran(seqnum, gseed, state)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine initializes the state for the seqnum(th) member        c
c  of the family of pseudorandom sequences.  initfran must             c
c  be called once (and only once) for each sequence from which         c
c  random floating point values are wanted.                            c
c                                                                      c
c  L, k, s, and mbits are properties associated with the family        c
c  of generators used by this package.  L and k represent a            c
c  particular primitive trinomial (mod 2), s is the "canonical bit"    c
c  associated with that trinomial, and mbits is the number of          c
c  mantissa bits used by the floating point pseudorandom generator.    c
c  for L = 158, state must have a declared dimension of at least       c
c  2 x 158 = 316 in the calling program.                               c
c                                                                      c
c  After initializing a sequence with initfran, multiple calls         c
c  can be made to frand,  and the state vector will stay current.      c
c  The state vector must not be modified by the user code.   If        c
c  more than one sequence is to be active at any time, there must      c
c  be a searate state vector for each active sequence.                 c
c                                                                      c
c  gseed is the integer value of the user-supplied "global seed."      c
c  The global seed is intended to be fixed for the duration of the     c
c  run and should be the same value for every sequence on every        c
c  processor.   It can be used to modify the random data, in a         c
c  predictable way, on a per-run basis.  seqnum must be in the         c
c  range {0,1,2,3,..., 65535 = 2^16 - 1 = z'0000ffff'} and gseed       c
c  must be in the range {0,1,2,3,..., 65534 = 2^16 - 2 = z'0000fffe'}. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer seqnum, gseed
      real (kind=8) state(2*L)
c----------------------------------------------------------------------c
      integer hi, lo, a, q, r, t31m1, seed, i, mask
      integer treg, tregx, highbit, t30m1
      parameter (a = 16807, q = 127773, r = 2836 , t31m1 = z'7fffffff')
      parameter (t30m1 = z'3fffffff')
      real (kind=8) scale
c
      call tausinit(seqnum, gseed, treg, tregx)
      scale    = 1.
      do  i = 1, mbits
         scale = scale * 0.5
      enddo
      mask     = ishft(1, mbits) - 2
      seed     = iand(treg, t31m1)
      highbit  = ieor(seed, treg)
c
c  initialize register:
      state(1) = scale*float(seqnum + seqnum)
      state(L) = 0.
      do  i = 2, L-2
         hi    = seed/q
         lo    = seed - (hi*q)
         seed  = a*lo - r*hi
         if (seed.le.0) seed = seed + t31m1
         state(i) = scale*float(iand(ishft(seed,-(31-mbits)), mask))
      enddo
c
      hi    = seed/q
      lo    = seed - (hi*q)
      seed  = a*lo - r*hi
      if (seed.le.0) seed = seed + t31m1
      if (highbit.eq.0) then
         seed = iand(seed, t30m1)
      else
         seed = ior (seed, ishft(highbit,-1))
      endif
      state(L-1) = scale*float(iand(ishft(seed,-(31-mbits)), mask))
c
      state(s+1)  = state(s+1) + scale
      return
      end
c
c
c
c
c
      subroutine irand(n, state, x)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine returns n 32-bit random integers in the array, x.      c
c                                                                      c
c  The state vector is updated to continue producing random numbers    c
c  from the current sequence by subsequent calls to this routine.      c
c  The state vector should not be modified by the calling routine,     c
c  unless no more values are required from the current sequence.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer (kind=8) n
      integer state(L, 2, 0:1)
      integer x(n)
c----------------------------------------------------------------------c
      integer i, old, new, kk, t31m1, t32m2
      integer (kind=8) xptr, k8, ii
      parameter (t31m1 = z'7fffffff', t32m2 = z'fffffffe')
c
      old  = 1
      new  = 0
      k8   = k  ! so min works
      xptr = 0
   10 if (xptr.ge.n) go to 99
         old  = new
         new  = ieor(new,1)
         kk = min(k8, n-xptr)
         do  i = 1, kk
            state(i,1,new) = state(k-kk+i,1,old) + state(L-kk+i,1,old)
            state(i,2,new) = state(k-kk+i,2,old) + state(L-kk+i,2,old)
         enddo
         do  i = kk+1, L
            state(i,1,new) = state(i-kk,1,old)
            state(i,2,new) = state(i-kk,2,old)
         enddo
         do  ii = xptr+1, min(xptr+kk, n)
            x(ii) = iand(ishft(state(kk+1-(ii-xptr),1,new),   -1),t31m1)
            x(ii) = ieor( iand(state(kk+1-(ii-xptr),2,new),t32m2),x(ii))
         enddo
         xptr = xptr + kk
         go to 10
   99 continue
c
c  update "old" state for subsequent calls:
      if (new.eq.1) then
         do  i = 1, L
            state(i,1,old) = state(i,1,new)
            state(i,2,old) = state(i,2,new)
         enddo
      endif
      return
      end
c
c
c
c
c
      subroutine frand(n, state, x)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine returns n random real numbers in the array, x.         c
c                                                                      c
c  The state vector is updated to continue producing random numbers    c
c  from the current sequence by subsequent calls to this routine.      c
c  The state vector should not be modified by the calling routine,     c
c  unless no more values are required from the current sequence.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer (kind=8) n
      real (kind=8) state(L, 0:1)
      real (kind=8) x(n)
c----------------------------------------------------------------------c
      integer i, old, new, kk
      integer (kind=8) xptr, k8, ii
c
      old  = 1
      new  = 0
      k8   = k  ! so min works
      xptr = 0
   10 if (xptr.ge.n) go to 99
         old  = new
         new  = ieor(new,1)
         kk   = min(k8, n-xptr)
         do  i = 1, kk
            state(i, new) = state(k-kk+i, old) + state(L-kk+i, old)
         enddo
         do  i = 1, kk
            if (state(i,new).ge.1.0) state(i,new) = state(i,new) - 1.0
         enddo
         do  i = kk+1, L
            state(i,new) = state(i-kk, old)
         enddo
         do  ii = xptr+1,  min(xptr+kk, n)
            x(ii) = state(kk+1 - (ii-xptr), new)
         enddo
         xptr = xptr + kk
         go to 10
   99 continue
c
c  update "old" state for subsequent calls:
      if (new.eq.1) then
         do  i = 1, L
            state(i,old) = state(i,new)
         enddo
      endif
      return
      end
c
c
c
c
c
      subroutine brand64(n, state, x)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Used for producing random bit streams on machines with 64-bit       c
c  integers, this routine returns n words of random 64-bit integers    c
c  in the array, x.  This is the same string of bits that would        c
c  be produced on a 32-bit machine by calling irand with a             c
c  (2 x n)-long x vector.  Thus, for creating binary streams, use      c
c  irand on 32-bit machines, and brand64 on 64-bit machines.           c
c                                                                      c
c  The state vector is updated to continue producing random numbers    c
c  from the current sequence by subsequent calls to this routine.      c
c  The state vector should not be modified by the calling routine,     c
c  unless no more values are required from the current sequence.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer wrklen
c  wrklen must be even!!
      parameter (wrklen = 8192)
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer (kind=8) n
      integer state(4*L)
      integer (kind=8) x(n)
c----------------------------------------------------------------------c
      integer (kind=8) wrk8(wrklen), low32
      integer (kind=8) xptr, nleft, wrklen8, ii
      integer wrk(wrklen), i
      data low32/ z'00000000FFFFFFFF'/
c
      xptr  = 1
      nleft = n
      wrklen8 = wrklen  ! for call to irand
c
   10 if ( nleft.le.0) go to 99
c
      if (2*nleft.gt.wrklen) then
c
         call irand(wrklen8, state, wrk)
c
         do  i = 1, wrklen
            wrk8(i) = wrk(i)
         enddo

         do  i = 1, wrklen
            wrk8(i) = iand(wrk8(i), low32)
         enddo

         do  ii = xptr, xptr+(wrklen/2)-1
            x(ii) = ior(       wrk8(2*(ii+1-xptr) - 1),
     $                  ishft( wrk8(2*(ii+1-xptr)    ), 32) )
         enddo

         xptr  = xptr  + (wrklen/2)
         nleft = nleft - (wrklen/2)

      else

         call irand ( 2*nleft, state, wrk )

         do  i = 1, 2*nleft
            wrk8(i) = wrk(i)
         enddo

         do  i = 1, 2*nleft
            wrk8(i) = iand(wrk8(i), low32)
         enddo

         do  ii = xptr, xptr+nleft-1
            x(ii) = ior(       wrk8(2*(ii+1-xptr) - 1),
     $                  ishft( wrk8(2*(ii+1-xptr)    ), 32) )
         enddo

         xptr  = xptr + nleft
         nleft = 0

      endif
      go to 10

   99 continue
      return
      end
c
c
c
c
c
      subroutine get_penum(seqnum, npes, nseqs, penum)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  given  a sequence number (seqnum), the number of pes (npes), and    c
c  the number of sequences (nseqs) this routine determines the number  c
c  of the pe on which sequence number seqnum is located.               c
c  Assumes that penum is in the range {0,1,2,...,npes-1} and that      c
c  sequences are numbered starting from 0.                             c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer seqnum, npes, nseqs, penum
c----------------------------------------------------------------------c
      integer q, r
c
      q     = nseqs/npes
      r     = nseqs - npes*q
      penum = seqnum/(q+1)
      if (seqnum.ge. ((q+1)*r) ) penum = r + (seqnum - ((q+1)*r))/q
      return
      end
c
c
c
c
c
      subroutine get_seq_nums(penum, npes, nseqs, seqfirst, seqlast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine computes the indices of the first and last sequences   c
c  that will be associated with pe number penum, when there are npes   c
c  pes, and nseqs sequences.  Assumes that penum is in the range       c
c  {0,1,2,...,npes-1} and that sequences are numbered starting from 0. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer penum, npes, nseqs, seqfirst, seqlast
c----------------------------------------------------------------------c
      integer q, r
c
      q        = nseqs/npes
      r        = nseqs - npes*q
      seqfirst = penum*q + min(penum, r)
      seqlast  = seqfirst + q - 1
      if (penum.lt.r) seqlast = seqlast + 1
      return
      end
c
c
c
c
c
      subroutine get_seq_bounds(seqnum, nseqs, lendata, chnksize,
     $                          datafirst, datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Given a particular sequence number (seqnum), the number of          c
c  sequences (nseqs), and the total length of the data (lendata) ,     c
c  this will compute datafirst and datalast, the positions in the      c
c  overall data set that sequence number seqnum occupies.   This       c
c  routine will force each sequence to have an integer number of data  c
c  values of length chnksize (with the exception of the last           c
c  sequence, which will have any partial chunk left over from a total  c
c  data length of lendata).  Assumes that data arrays are numbered     c
c  starting from 1, and that sequences are numbered starting from 0.   c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer seqnum, nseqs, chnksize
      integer (kind=8) lendata, datafirst, datalast
c----------------------------------------------------------------------c
      integer (kind=8) nchunks, q, r, seqnum8
c
      nchunks   = lendata/chnksize
      q         = nchunks/nseqs
      r         = nchunks - nseqs*q
      seqnum8   = seqnum  ! so min works
      datafirst = chnksize * (seqnum*q + min(seqnum8, r)) + 1
      datalast  = datafirst + q*chnksize - 1
      if (seqnum.lt.r        ) datalast = datalast + chnksize
      if (seqnum.eq.(nseqs-1)) datalast = lendata
      return
      end
c
c
c
c
c
      subroutine get_my_idata(penum, npes, nseqs, lendata, chnksize,
     $                        gseed, x, datafirst, datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  returns the integer data belonging to the PE numbered penum,        c
c  given that there are npes PEs and given that the user has           c
c  specified nseqs sequences, a total data length of lendata           c
c  integers, a chunk size of chnksize integers, and a global seed      c
c  value of gseed.  The data belonging to PE penum will be returned in c
c  the integer array, x, starting with x(1),  but will be used by the  c
c  calling program as data items datafirst through datalast in the     c
c  overall lendata-long set that resides across all PEs.  For this     c
c  routine to work correctly, npes must be less than the total number  c
c  of chunks in the data set.  I.e. npes must be less than             c
c  lendata/chnksize.  The size of x must be declared to be at least as c
c  large as the value of maxlen returned by the routine check_maxlen.  c
c  Of course, this will be about lendata/npes, but the value of maxlen c
c  returned by check_maxlen is the exact value needed by the the       c
c  worst-case PE.                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer penum, npes, nseqs, chnksize, gseed
      integer (kind=8) lendata, datafirst, datalast
      integer x(lendata)
c----------------------------------------------------------------------c
      integer state(4*L)
      integer seq_frst
      integer (kind=8) nchnks, q, lendatac, r, penum8, ch_frst, ch_last,
     $        mydata_len, dfrst, dlast, dlen, nskip, ix, mydata_left,
     $        len_data1, nseq_len2, nseq_len1, seq_frst2, datafirst2,
     $        len_data2, seq_len1, seq_len2c, seq_len2, rc, seq_len1c
c
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      penum8      = penum  ! so min works
      ch_frst     =  penum   *q + min( penum8   , r) + 1
      ch_last     = (penum+1)*q + min((penum8+1), r)
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if (penum.eq.npes-1) datalast = lendata
      mydata_len  = datalast - datafirst + 1
c
      lendatac = lendata/chnksize
      seq_len2c = lendatac/nseqs
      seq_len1c = seq_len2c + 1
      rc        = lendatac - seq_len2c*nseqs
c
      seq_len2  = chnksize*seq_len2c
      seq_len1  = chnksize*seq_len1c
      nseq_len1 = rc
      nseq_len2 = nseqs - rc
c
      len_data1 = nseq_len1*seq_len1
      len_data2 = nseq_len2*seq_len2
c
      if (datafirst.gt.len_data1) then
         datafirst2 = datafirst - len_data1
         seq_frst2  = datafirst2/seq_len2
         nskip      = datafirst2 - seq_frst2*seq_len2 - 1
         seq_frst   = seq_frst2 + nseq_len1
      else
         seq_frst   = datafirst/seq_len1
         nskip      = datafirst - seq_frst*seq_len1 - 1
      endif
c
      call get_seq_bounds(seq_frst,nseqs,lendata,chnksize,dfrst,dlast)
      dlen        = dlast - dfrst + 1
c
      call initiran(seq_frst,gseed,state)
      call irand(nskip,state,x)
c
      dlen        = dlen - nskip
      mydata_left = mydata_len
      ix          = 1
c
      call irand(dlen,state,x(ix))
c
      mydata_left = mydata_left - dlen
      ix          = ix + dlen
c
   10 if (mydata_left.le.0) go to 99
         seq_frst = seq_frst + 1
         call get_seq_bounds(seq_frst,nseqs,lendata,
     $                       chnksize,dfrst,dlast)
         dlen     = min(dlast-dfrst+1, mydata_left)
c
         call initiran(seq_frst,gseed,state)
         call irand(dlen,state,x(ix))
c
         mydata_left = mydata_left - dlen
         ix          = ix + dlen
         go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine get_my_fdata(penum, npes, nseqs, lendata, chnksize,
     $                        gseed, x, datafirst, datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  returns the floating point data belonging to the PE numbered penum, c
c  given that there are npes PEs and given that the user has           c
c  specified nseqs sequences, a total data length of lendata           c
c  integers, a chunk size of chnksize reals, and a global seed         c
c  value of gseed.  The data belonging to PE penum will be returned in c
c  the real array, x, starting with x(1),  but will be used by the     c
c  calling program as data items datafirst through datalast in the     c
c  overall lendata-long set that resides across all PEs.  For this     c
c  routine to work correctly, npes must be less than the total number  c
c  of chunks in the data set.  I.e. npes must be less than             c
c  lendata/chnksize.  The size of x must be declared to be at least as c
c  large as the value of maxlen returned by the routine check_maxlen.  c
c  Of course, this will be about lendata/npes, but the value of maxlen c
c  returned by check_maxlen is the exact value needed by the the       c
c  worst-case PE.                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer penum, npes, nseqs, chnksize, gseed
      integer (kind=8) lendata, datafirst, datalast
      real (kind=8) x(lendata)
c----------------------------------------------------------------------c
      real (kind=8) state(2*L)
      integer seq_frst
      integer (kind=8) nchnks, q, lendatac, r, penum8, ch_frst, ch_last,
     $        mydata_len, dfrst, dlast, dlen, nskip, ix, mydata_left,
     $        len_data1, nseq_len2, nseq_len1, seq_frst2, datafirst2,
     $        len_data2, seq_len1, seq_len2c, seq_len2, rc, seq_len1c
c
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      penum8      = penum  ! so min works
      ch_frst     =  penum   *q + min( penum8   , r) + 1
      ch_last     = (penum+1)*q + min((penum8+1), r)
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if (penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
c
      lendatac = lendata/chnksize
      seq_len2c = lendatac/nseqs
      seq_len1c = seq_len2c + 1
      rc        = lendatac - seq_len2c*nseqs
c
      seq_len2  = chnksize*seq_len2c
      seq_len1  = chnksize*seq_len1c
      nseq_len1 = rc
      nseq_len2 = nseqs - rc
c
      len_data1 = nseq_len1*seq_len1
      len_data2 = nseq_len2*seq_len2
c
      if (datafirst.gt.len_data1) then
         datafirst2 = datafirst - len_data1
         seq_frst2  = datafirst2/seq_len2
         nskip      = datafirst2 - seq_frst2*seq_len2 - 1
         seq_frst   = seq_frst2 + nseq_len1
      else
         seq_frst   = datafirst/seq_len1
         nskip      = datafirst - seq_frst*seq_len1 - 1
      endif
c
      call get_seq_bounds(seq_frst,nseqs,lendata,chnksize,dfrst,dlast)
      dlen        = dlast - dfrst + 1
c
      call initfran(seq_frst,gseed,state)
      call frand(nskip,state,x)
c
      dlen        = dlen - nskip
      mydata_left = mydata_len
      ix          = 1
c
      call frand(dlen,state,x(ix))
c
      mydata_left = mydata_left - dlen
      ix          = ix + dlen
c
   10 if (mydata_left.le.0) go to 99
         seq_frst = seq_frst + 1
         call get_seq_bounds(seq_frst,nseqs,lendata,
     $                       chnksize,dfrst,dlast)
         dlen     = min(dlast-dfrst+1, mydata_left)
c
         call initfran(seq_frst,gseed,state)
         call frand(dlen,state,x(ix))
c
         mydata_left = mydata_left - dlen
         ix          = ix + dlen
         go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine get_my_bdata64(penum, npes, nseqs, lendata, chnksize,
     $                        gseed, x, datafirst, datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  returns the 64-bit integer data belonging to the PE numbered penum, c
c  given that there are npes PEs and given that the user has           c
c  specified nseqs sequences, a total data length of lendata 64-bit    c
c  integers, a chunk size of chnksize 64-bit integers, and a global    c
c  seed value of gseed.  The data belonging to PE penum will be        c
c  returned in the 64-bit integer array, x, starting with x(1),  but   c
c  will be used by the calling program as data items datafirst         c
c  through datalast in the overall lendata-long set that resides       c
c  across all PEs.  For this routine to work correctly, npes must be   c
c  less than the total number of chunks in the data set.  I.e. npes    c
c  must be less than lendata/chnksize.  The size of x must be          c
c  declared to be at least as large as the value of maxlen returned    c
c  by the routine check_maxlen.  Of course, this will be about         c
c  lendata/npes, but the value of maxlen returned by check_maxlen is   c
c  the exact value needed by the the worst-case PE.                    c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer L, k, s, mbits
      parameter (L = 158, k = 128, s =  63, mbits = 23)
      integer   penum, npes, nseqs, chnksize, gseed
      integer (kind=8) lendata, datafirst, datalast
      integer (kind=8) x(lendata)
c----------------------------------------------------------------------c
      integer state(4*L)
      integer seq_frst
      integer (kind=8) nchnks, q, lendatac, r, penum8, ch_frst, ch_last,
     $        mydata_len, dlen, dfrst, dlast, nskip, ix, mydata_left,
     $        len_data1, nseq_len2, nseq_len1, seq_frst2, datafirst2,
     $        len_data2, seq_len1, seq_len2c, seq_len2, rc, seq_len1c
c
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      penum8      = penum  ! so min works
      ch_frst     =  penum   *q + min( penum8   , r) + 1
      ch_last     = (penum+1)*q + min((penum8+1), r)
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if (penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
c
      lendatac = lendata/chnksize
      seq_len2c = lendatac/nseqs
      seq_len1c = seq_len2c + 1
      rc        = lendatac - seq_len2c*nseqs
c
      seq_len2  = chnksize*seq_len2c
      seq_len1  = chnksize*seq_len1c
      nseq_len1 = rc
      nseq_len2 = nseqs - rc
c
      len_data1 = nseq_len1*seq_len1
      len_data2 = nseq_len2*seq_len2
c
      if (datafirst.gt.len_data1) then
         datafirst2 = datafirst - len_data1
         seq_frst2  = datafirst2/seq_len2
         nskip      = datafirst2 - seq_frst2*seq_len2 - 1
         seq_frst   = seq_frst2 + nseq_len1
      else
         seq_frst   = datafirst/seq_len1
         nskip      = datafirst - seq_frst*seq_len1 - 1
      endif
c
      call get_seq_bounds(seq_frst,nseqs,lendata,chnksize,dfrst,dlast)
      dlen        = dlast - dfrst + 1
c
      call initiran(seq_frst,gseed,state)
      call brand64(nskip,state,x)
c
      dlen        = dlen - nskip
      mydata_left = mydata_len
      ix          = 1
c
      call brand64(dlen,state,x(ix))
c
      mydata_left = mydata_left - dlen
      ix          = ix + dlen
c
   10 if (mydata_left.le.0) go to 99
         seq_frst = seq_frst + 1
         call get_seq_bounds(seq_frst,nseqs,lendata,
     $                       chnksize,dfrst,dlast)
         dlen     = min(dlast-dfrst+1, mydata_left)
c
         call initiran(seq_frst,gseed,state)
         call brand64(dlen,state,x(ix))
c
         mydata_left = mydata_left - dlen
         ix          = ix + dlen
         go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine check_maxlen(npes, nseqs, lendata, chnksize, maxlen)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer npes, nseqs, chnksize
      integer (kind=8) lendata, maxlen
c----------------------------------------------------------------------c
      integer (kind=8) ch_frst, ch_last, datafirst, datalast,
     $        mydata_len, nchnks, penum, q, r
c
      penum       = 0
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      r           = nchnks - npes*q
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      ch_last     = (penum+1)*q + min((penum+1),r)
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if (npes.eq.0) datalast = lendata
      mydata_len  = datalast - datafirst + 1
      maxlen      = mydata_len
c
      penum       = npes-1
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      r           = nchnks - npes*q
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      ch_last     = (penum+1)*q + min((penum+1),r)
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if (penum.eq.npes-1) datalast = lendata
      mydata_len  = datalast - datafirst + 1
      maxlen      = max(maxlen,mydata_len)
c
      return
      end
c
c
c
c
c
      subroutine tausinit(seqnum, gseed, treg64, treg96)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  uses a 32-bit tausworthe generator to map the pair (seqnum,gseed)   c
c  into a 32-bit pseudorandom number, for subsequent use in filling    c
c  out the fibonacci register.  This routine produces two values,      c
c  treg64 and treg96, which represent 64 and 96 steps of the           c
c  initial fill.  Since the calling routine will only be using         c
c  the low 31 bits (mod 2^31 - 1), there are 6 cases to avoid.         c
c  The bad values for the tausworthe register are (in order):          c
c    T^(-64)(80000000),T^(-64)(7fffffff),T^(-64)(ffffffff),            c
c    T^(-96)(80000000),T^(-96)(7fffffff),T^(-96)(ffffffff).            c
c  these cases correspond to (gseed,seqnum) pairs (in decimal):        c
c  (43553,65535), (26141,65535), (52283,65535),                        c
c  (50425,37887), (17320,35839), (34642, 6143).                        c
c  these are mapped to multiples of the prime number 9349, which       c
c  is approximately 65536/7 and which do not correspond to valid       c
c  initial values of seqnum and gseed that would be passed into        c
c  this routine.                                                       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer seqnum, gseed, treg64, treg96
c----------------------------------------------------------------------c
      integer seqmask, gshft, i, treg
      parameter (seqmask = z'0000ffff')
      parameter (gshft   = 16)
      integer tbad(6), t_ok(6)
      data tbad/ z'AA220000', z'661E0000', z'CC3C0000',
     $           z'C4FA9400', z'43A98C00', z'87531800'/
      data t_ok/ z'00002485', z'0000490a', z'00006d8f',
     $           z'00009214', z'0000b699', z'0000db1e'/
c
      treg = ishft ( gseed + 1, gshft )
      treg = ior ( treg, iand(seqnum+1, seqmask) )
      do  i = 1, 6
         if (treg.eq.tbad(i)) treg = t_ok(i)
      enddo
c
      call taus32(treg,64)
      treg64 = treg
      call taus32(treg,32)
      treg96 = treg
      return
      end
c
c
c
c
c
      subroutine taus32(treg, nsteps)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  steps the 32-bit tausworthe generator nsteps times, using the       c
c  mod 2 primitive polynomial (32,7,5,3,2,1,0).  The table entries     c
c  are used to look up the parity of the least significant 7 bits.     c
c  Tausworthe generators are discussed in Numerical Recipes by         c
c  Press, et al.  Cambridge University Press, c1986, pp 209-213.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer treg, nsteps
c----------------------------------------------------------------------c
      integer i, t31, t32m1, table0(0:127), table1(0:127), bit31, low7
      parameter (t31     = z'80000000')
      parameter (t32m1   = z'FFFFFFFF')
c
      data table0 /
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1/
      data table1 /
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0/
c
      do  i = 1, nsteps
         bit31 = iand(treg,t31)
         low7  = iand(treg,127)
         treg  = iand(ishft(treg,1),t32m1)
         if (bit31.eq.0) then
            treg = ior(treg,table0(low7))
         else
            treg = ior(treg,table1(low7))
         endif
      enddo
      return
      end