Schrage's Method

      SUBROUTINE sub2(npts,seed,output)

      INTEGER npts,seed,output(npts)
      INTEGER magic1,magic2,j,hi,lo,test

      magic1 = 7**5
      magic2 = 2*(2**30 - 1) + 1
      quot = magic2/magic1
      rem = MOD(magic2,magic1)
      output(1) = seed
      DO j = 2,npts
         hi = output(j - 1)/quot
         lo = MOD(output(j - 1),quot)
         test = magic1*lo - rem*hi
         IF (test.GT.0) THEN
            output(j) = test
         ELSE
            output(j) = test + magic2
         END IF
      END DO
      hi = output(npts)/quot
      lo = MOD(output(npts),quot)
      test = magic1*lo - rem*hi
      IF (test.GT.0) THEN
         seed = test
      ELSE
         seed = test + magic2
      END IF
 
      RETURN
      END