SUBROUTINE nied(QUASI)
C
      INTEGER MAXDIM, NBITS
      PARAMETER (MAXDIM=100, NBITS=31)
C
C   The parameter MAXDIM is the maximum dimension that will be used.
C   NBITS is the number of bits (not counting the sign) in a
C   fixed-point integer.
C
      INTEGER CJ(MAXDIM, 0:NBITS-1), DIMEN, COUNT, NEXTQ(MAXDIM)
      COMMON /COMM2/ CJ, DIMEN, COUNT, NEXTQ
      SAVE   /COMM2/
C
C   The common block /COMM2/ :
C     CJ    - Contains the packed values of Niederreiter's C(I,J,R)
C     DIMEN   - The dimension of the sequence to be generated
C     COUNT - The index of the current item in the sequence,
C             expressed as an array of bits.  COUNT(R) is the same
C             as Niederreiter's AR(N) (page 54) except that N is
C             implicit.
C     NEXTQ - The numerators of the next item in the series.  These
C             are like Niederreiter's XI(N) (page 54) except that
C             N is implicit, and the NEXTQ are integers.  To obtain
C             the values of XI(N), multiply by RECIP (see GOLO2).
C
C   Array CJ of the common block is set up by subroutine CALCC2.
C   The other items in the common block are set up by INLO2.
C
C   ------------------------------------------------------------
C
C
C
      double precision  RECIP
      PARAMETER (RECIP=2.0D0**(-NBITS))
C
C   The parameter RECIP is the multiplier which changes the
C   integers in NEXTQ into the required real values in QUASI.
C
      INTEGER I, R
      double precision  QUASI(*)
C
C Multiply the numerators in NEXTQ by RECIP to get the next
C quasi-random vector
C
      DO 5 I = 1, DIMEN
        QUASI(I) = NEXTQ(I) * RECIP
    5 CONTINUE
C
C Find the position of the right-hand zero in COUNT.  This
C is the bit that changes in the Gray-code representation as
C we go from COUNT to COUNT+1.
C
      R = 0
      I = COUNT
   10 IF (MOD(I,2).NE.0) THEN
        R = R + 1
        I = I/2
        GOTO 10
      ENDIF
C
C Check that we have not passed 2**NBITS calls on GOLO2
C
      IF (R .GE. NBITS) THEN
        WRITE (*,*) ' GOLO2 :  Too many calls'
        STOP
      ENDIF
C
C Compute the new numerators in vector NEXTQ
C
      DO 20 I = 1, DIMEN
C       ***** Bitwise exclusive-or is not standard in Fortran
C       ***** This is the Vax version :
       NEXTQ(I) = IEOR(NEXTQ(I), CJ(I,R))
C       ***** This is the Unix version
C      NEXTQ(I) = XOR(NEXTQ(I), CJ(I,R))
   20 CONTINUE
C
      COUNT = COUNT + 1
      RETURN
      END
C
C     *****   end of PROCEDURE GOLO2

C =========================================================

      SUBROUTINE init_nied(DIM, SKIP)
C
C   This version :  12 February 1992
C
C   See the general comments on implementing Niederreiter's
C   low-discrepancy sequences.
C
C   This subroutine calculates the values of Niederreiter's
C   C(I,J,R) and performs other initialisation necessary
C   before calling GOLO2.
C
C INPUT :
C   DIMEN - The dimension of the sequence to be generated.
C        {DIMEN is called DIM in the argument of INLO2,
C        because DIMEN is subsequently passed via COMMON
C        and is called DIMEN there.}
C
C   SKIP  - The number of values to throw away at the beginning
C           of the sequence.
C
C OUTPUT :
C   To GOLO2, labelled common /COMM2/.
C
C USES :
C   Calls CALCC2 to calculate the values of CJ.
C   ***** A non-standard function is used to compute *****
C   ***** bitwise exclusive-ors.                     *****
C
C
C   ------------------------------------------------------------
C
C
C   This file defines the common block /COMM2/ and some
C   associated parameters.  These are for use in the base 2
C   version of the generator.
C
      INTEGER MAXDIM, NBITS
      PARAMETER (MAXDIM=100, NBITS=31)
C
C   The parameter MAXDIM is the maximum dimension that will be used.
C   NBITS is the number of bits (not counting the sign) in a
C   fixed-point integer.
C
      INTEGER CJ(MAXDIM, 0:NBITS - 1), DIMEN, COUNT
      INTEGER NEXTQ(MAXDIM)
      COMMON /COMM2/ CJ, DIMEN, COUNT, NEXTQ
      SAVE   /COMM2/
C
C   The common block /COMM2/ :
C     CJ    - Contains the packed values of Niederreiter's C(I,J,R)
C     DIMEN   - The dimension of the sequence to be generated
C     COUNT - The index of the current item in the sequence,
C             expressed as an array of bits.  COUNT(R) is the same
C             as Niederreiter's AR(N) (page 54) except that N is
C             implicit.
C     NEXTQ - The numerators of the next item in the series.  These
C             are like Niederreiter's XI(N) (page 54) except that
C             N is implicit, and the NEXTQ are integers.  To obtain
C             the values of XI(N), multiply by RECIP (see GOLO2).
C
C   Array CJ of the common block is set up by subroutine CALCC2.
C   The other items in the common block are set up by INLO2.
C
C   ------------------------------------------------------------
C
C
C
      INTEGER I, R, DIM, SKIP, GRAY
C
      DIMEN = DIM
C
C       This assignment just relabels the variable for
C       subsequent use.
C
      IF (DIMEN.LE.0 .OR. DIMEN.GT.MAXDIM) THEN
        WRITE (*,*) ' INLO2 :  Bad dimension'
        STOP
      ENDIF
C
      CALL CALCC2
C
C   Translate SKIP into Gray code
C
C   ***** The bitwise exclusive-or is not standard in Fortran
C   ***** This is the Vax version :
      GRAY = IEOR (SKIP, SKIP/2)
C   ***** THIS is the Unix version
C     GRAY = XOR (SKIP, SKIP/2)
C
C   Now set up NEXTQ appropriately for this value of the Gray code
C
      DO 5 I = 1, DIMEN
    5   NEXTQ(I) = 0
C
      R = 0
   10 IF (GRAY .NE. 0) THEN
        IF (MOD(GRAY,2) .NE. 0) THEN
          DO 20 I = 1, DIMEN
C           ***** This is the non-standard exclusive-or again
C           ***** Vax version :
            NEXTQ(I) = IEOR(NEXTQ(I), CJ(I,R))
C           ***** Unix version :
C           NEXTQ(I) = XOR(NEXTQ(I), CJ(I,R))
   20     CONTINUE
        ENDIF
        GRAY = GRAY / 2
        R = R + 1
        GOTO 10
      ENDIF
C
      COUNT = SKIP
      RETURN
      END
C
C     *****  end of SUBROUTINE INLO2

C =======================================================

      SUBROUTINE CALCC2
C
C   This version :  12 February 1992
C
C      *****  For base-2 only.
C
C   See the general comments on implementing Niederreiter's
C   low-discrepancy sequences.
C
C   This program calculates the values of the constants C(I,J,R).
C   As far as possible, we use Niederreiter's notation.
C   For each value of I, we first calculate all the corresponding
C   values of C :  these are held in the array CI.  All these
C   values are either 0 or 1.  Next we pack the values into the
C   array CJ, in such a way that CJ(I,R) holds the values of C
C   for the indicated values of I and R and for every value of
C   J from 1 to NBITS.  The most significant bit of CJ(I,R)
C   (not counting the sign bit) is C(I,1,R) and the least
C   significant bit is C(I,NBITS,R).
C     When all the values of CJ have been calculated, we return
C   this array to the calling program.
C
C  --------------------------------------------------------------
C
C   We define the common block /COMM2/ and some
C   associated parameters.  These are for use in the base 2
C   version of the generator.
C
      INTEGER MAXDIM, NBITS
      PARAMETER (MAXDIM=100, NBITS=31)
C
C   The parameter MAXDIM is the maximum dimension that will be used.
C   NBITS is the number of bits (not counting the sign) in a
C   fixed-point integer.
C
      INTEGER CJ(MAXDIM, 0:NBITS-1), DIMEN, COUNT, NEXTQ(MAXDIM)
      COMMON /COMM2/ CJ, DIMEN, COUNT, NEXTQ
      SAVE   /COMM2/
C
C   The common block /COMM2/ :
C     CJ    - Contains the packed values of Niederreiter's C(I,J,R)
C     DIMEN   - The dimension of the sequence to be generated
C     COUNT - The index of the current item in the sequence,
C             expressed as an array of bits.  COUNT(R) is the same
C             as Niederreiter's AR(N) (page 54) except that N is
C             implicit.
C     NEXTQ - The numerators of the next item in the series.  These
C             are like Niederreiter's XI(N) (page 54) except that
C             N is implicit, and the NEXTQ are integers.  To obtain
C             the values of XI(N), multiply by RECIP (see GOLO2).
C
C   Array CJ of the common block is set up by subroutine CALCC2.
C   The other items in the common block are set up by INLO2.
C
C   --------------------------------------------------------------
C
C   The following COMMON block, used by many subroutines,
C   gives the order Q of a field, its characteristic P, and its
C   addition, multiplication and subtraction tables.
C   The parameter MAXQ gives the order of the largest field to
C   be handled.
C
      INTEGER MAXQ
      PARAMETER (MAXQ=50)

      INTEGER P, Q, ADD(0:MAXQ-1,0:MAXQ-1)
      INTEGER MUL(0:MAXQ-1, 0:MAXQ-1), SUB(0:MAXQ-1, 0:MAXQ-1)
      COMMON /FIELD/ P, Q, ADD, MUL, SUB
      SAVE /FIELD/
C
C   The following definitions concern the representation of
C   polynomials.
C
      INTEGER MAXDEG, DEG
      PARAMETER (MAXDEG=50, DEG=-1)
C
C   The parameter MAXDEG gives the highest degree of polynomial
C   to be handled.  Polynomials stored as arrays have the
C   coefficient of degree n in POLY(N), and the degree of the
C   polynomial in POLY(-1).  The parameter DEG is just to remind
C   us of this last fact.  A polynomial which is identically 0
C   is given degree -1.
C
C   A polynomial can also be stored in an integer I, with
C        I = AN*Q**N + ... + A0.
C   Routines ITOP and PTOI convert between these two formats.
C
C   ---------------------------------------------------------------
C
C
C
C   MAXE   - We need MAXDIM irreducible polynomials over Z2.
C            MAXE is the highest degree among these.
C   MAXV   - The maximum possible index used in V.
C
      INTEGER MAXE, MAXV
      PARAMETER (MAXE=9, MAXV=NBITS+MAXE)
C
C INPUT :
C   The array CJ to be initialised, and DIMEN the number of
C   dimensions we are using, are transmitted through /COMM2/.
C
C OUTPUT :
C   The array CJ is returned to the calling program.
C
C USES :
C   Subroutine SETFLD is used to set up field arithmetic tables.
C   (Although this is a little heavy-handed for the field of
C   order 2, it makes for uniformity with the general program.)
C   Subroutine CALCV is used to for the auxiliary calculation
C   of the values V needed to get the Cs.
C
      INTEGER PX(-1:MAXDEG), B(-1:MAXDEG)
      INTEGER V(0:MAXV), CI(NBITS, 0:NBITS-1)
      INTEGER E, I, J, R, U, TERM
C
      INTEGER IRRED(MAXDIM, -1:MAXE)
      SAVE IRRED
C
C   This DATA statement supplies the coefficients and the
C   degrees of the first 100 irreducible polynomials over Z2.
C   They are taken from Lidl and Niederreiter, FINITE FIELDS,
C   Cambridge University Press (1984), page 553.
C   The order of these polynomials is the same as the order in
C   file 'irrtabs.dat' used by the general program.
C
C   In this block PX(I, -1) is the degree of the Ith polynomial,
C   and PX(I, N) is the coefficient of x**n in the Ith polynomial.
C
      DATA (IRRED(  1,I),I=-1,9) /  1,0,1,0,0,0,0,0,0,0,0 /
      DATA (IRRED(  2,I),I=-1,9) /  1,1,1,0,0,0,0,0,0,0,0 /
      DATA (IRRED(  3,I),I=-1,9) /  2,1,1,1,0,0,0,0,0,0,0 /
      DATA (IRRED(  4,I),I=-1,9) /  3,1,1,0,1,0,0,0,0,0,0 /
      DATA (IRRED(  5,I),I=-1,9) /  3,1,0,1,1,0,0,0,0,0,0 /
      DATA (IRRED(  6,I),I=-1,9) /  4,1,1,0,0,1,0,0,0,0,0 /
      DATA (IRRED(  7,I),I=-1,9) /  4,1,0,0,1,1,0,0,0,0,0 /
      DATA (IRRED(  8,I),I=-1,9) /  4,1,1,1,1,1,0,0,0,0,0 /
      DATA (IRRED(  9,I),I=-1,9) /  5,1,0,1,0,0,1,0,0,0,0 /
      DATA (IRRED( 10,I),I=-1,9) /  5,1,0,0,1,0,1,0,0,0,0 /
      DATA (IRRED( 11,I),I=-1,9) /  5,1,1,1,1,0,1,0,0,0,0 /
      DATA (IRRED( 12,I),I=-1,9) /  5,1,1,1,0,1,1,0,0,0,0 /
      DATA (IRRED( 13,I),I=-1,9) /  5,1,1,0,1,1,1,0,0,0,0 /
      DATA (IRRED( 14,I),I=-1,9) /  5,1,0,1,1,1,1,0,0,0,0 /
      DATA (IRRED( 15,I),I=-1,9) /  6,1,1,0,0,0,0,1,0,0,0 /
      DATA (IRRED( 16,I),I=-1,9) /  6,1,0,0,1,0,0,1,0,0,0 /
      DATA (IRRED( 17,I),I=-1,9) /  6,1,1,1,0,1,0,1,0,0,0 /
      DATA (IRRED( 18,I),I=-1,9) /  6,1,1,0,1,1,0,1,0,0,0 /
      DATA (IRRED( 19,I),I=-1,9) /  6,1,0,0,0,0,1,1,0,0,0 /
      DATA (IRRED( 20,I),I=-1,9) /  6,1,1,1,0,0,1,1,0,0,0 /
      DATA (IRRED( 21,I),I=-1,9) /  6,1,0,1,1,0,1,1,0,0,0 /
      DATA (IRRED( 22,I),I=-1,9) /  6,1,1,0,0,1,1,1,0,0,0 /
      DATA (IRRED( 23,I),I=-1,9) /  6,1,0,1,0,1,1,1,0,0,0 /
      DATA (IRRED( 24,I),I=-1,9) /  7,1,1,0,0,0,0,0,1,0,0 /
      DATA (IRRED( 25,I),I=-1,9) /  7,1,0,0,1,0,0,0,1,0,0 /
      DATA (IRRED( 26,I),I=-1,9) /  7,1,1,1,1,0,0,0,1,0,0 /
      DATA (IRRED( 27,I),I=-1,9) /  7,1,0,0,0,1,0,0,1,0,0 /
      DATA (IRRED( 28,I),I=-1,9) /  7,1,0,1,1,1,0,0,1,0,0 /
      DATA (IRRED( 29,I),I=-1,9) /  7,1,1,1,0,0,1,0,1,0,0 /
      DATA (IRRED( 30,I),I=-1,9) /  7,1,1,0,1,0,1,0,1,0,0 /
      DATA (IRRED( 31,I),I=-1,9) /  7,1,0,0,1,1,1,0,1,0,0 /
      DATA (IRRED( 32,I),I=-1,9) /  7,1,1,1,1,1,1,0,1,0,0 /
      DATA (IRRED( 33,I),I=-1,9) /  7,1,0,0,0,0,0,1,1,0,0 /
      DATA (IRRED( 34,I),I=-1,9) /  7,1,1,0,1,0,0,1,1,0,0 /
      DATA (IRRED( 35,I),I=-1,9) /  7,1,1,0,0,1,0,1,1,0,0 /
      DATA (IRRED( 36,I),I=-1,9) /  7,1,0,1,0,1,0,1,1,0,0 /
      DATA (IRRED( 37,I),I=-1,9) /  7,1,0,1,0,0,1,1,1,0,0 /
      DATA (IRRED( 38,I),I=-1,9) /  7,1,1,1,1,0,1,1,1,0,0 /
      DATA (IRRED( 39,I),I=-1,9) /  7,1,0,0,0,1,1,1,1,0,0 /
      DATA (IRRED( 40,I),I=-1,9) /  7,1,1,1,0,1,1,1,1,0,0 /
      DATA (IRRED( 41,I),I=-1,9) /  7,1,0,1,1,1,1,1,1,0,0 /
      DATA (IRRED( 42,I),I=-1,9) /  8,1,1,0,1,1,0,0,0,1,0 /
      DATA (IRRED( 43,I),I=-1,9) /  8,1,0,1,1,1,0,0,0,1,0 /
      DATA (IRRED( 44,I),I=-1,9) /  8,1,1,0,1,0,1,0,0,1,0 /
      DATA (IRRED( 45,I),I=-1,9) /  8,1,0,1,1,0,1,0,0,1,0 /
      DATA (IRRED( 46,I),I=-1,9) /  8,1,0,0,1,1,1,0,0,1,0 /
      DATA (IRRED( 47,I),I=-1,9) /  8,1,1,1,1,1,1,0,0,1,0 /
      DATA (IRRED( 48,I),I=-1,9) /  8,1,0,1,1,0,0,1,0,1,0 /
      DATA (IRRED( 49,I),I=-1,9) /  8,1,1,1,1,1,0,1,0,1,0 /
      DATA (IRRED( 50,I),I=-1,9) /  8,1,1,0,0,0,1,1,0,1,0 /
      DATA (IRRED( 51,I),I=-1,9) /  8,1,0,1,0,0,1,1,0,1,0 /
      DATA (IRRED( 52,I),I=-1,9) /  8,1,0,0,1,0,1,1,0,1,0 /
      DATA (IRRED( 53,I),I=-1,9) /  8,1,0,0,0,1,1,1,0,1,0 /
      DATA (IRRED( 54,I),I=-1,9) /  8,1,1,1,0,1,1,1,0,1,0 /
      DATA (IRRED( 55,I),I=-1,9) /  8,1,1,0,1,1,1,1,0,1,0 /
      DATA (IRRED( 56,I),I=-1,9) /  8,1,1,1,0,0,0,0,1,1,0 /
      DATA (IRRED( 57,I),I=-1,9) /  8,1,1,0,1,0,0,0,1,1,0 /
      DATA (IRRED( 58,I),I=-1,9) /  8,1,0,1,1,0,0,0,1,1,0 /
      DATA (IRRED( 59,I),I=-1,9) /  8,1,1,1,1,1,0,0,1,1,0 /
      DATA (IRRED( 60,I),I=-1,9) /  8,1,1,0,0,0,1,0,1,1,0 /
      DATA (IRRED( 61,I),I=-1,9) /  8,1,0,0,1,0,1,0,1,1,0 /
      DATA (IRRED( 62,I),I=-1,9) /  8,1,0,0,0,1,1,0,1,1,0 /
      DATA (IRRED( 63,I),I=-1,9) /  8,1,0,1,1,1,1,0,1,1,0 /
      DATA (IRRED( 64,I),I=-1,9) /  8,1,1,0,0,0,0,1,1,1,0 /
      DATA (IRRED( 65,I),I=-1,9) /  8,1,1,1,1,0,0,1,1,1,0 /
      DATA (IRRED( 66,I),I=-1,9) /  8,1,1,1,0,1,0,1,1,1,0 /
      DATA (IRRED( 67,I),I=-1,9) /  8,1,0,1,1,1,0,1,1,1,0 /
      DATA (IRRED( 68,I),I=-1,9) /  8,1,1,1,0,0,1,1,1,1,0 /
      DATA (IRRED( 69,I),I=-1,9) /  8,1,1,0,0,1,1,1,1,1,0 /
      DATA (IRRED( 70,I),I=-1,9) /  8,1,0,1,0,1,1,1,1,1,0 /
      DATA (IRRED( 71,I),I=-1,9) /  8,1,0,0,1,1,1,1,1,1,0 /
      DATA (IRRED( 72,I),I=-1,9) /  9,1,1,0,0,0,0,0,0,0,1 /
      DATA (IRRED( 73,I),I=-1,9) /  9,1,0,0,0,1,0,0,0,0,1 /
      DATA (IRRED( 74,I),I=-1,9) /  9,1,1,1,0,1,0,0,0,0,1 /
      DATA (IRRED( 75,I),I=-1,9) /  9,1,1,0,1,1,0,0,0,0,1 /
      DATA (IRRED( 76,I),I=-1,9) /  9,1,0,0,0,0,1,0,0,0,1 /
      DATA (IRRED( 77,I),I=-1,9) /  9,1,0,1,1,0,1,0,0,0,1 /
      DATA (IRRED( 78,I),I=-1,9) /  9,1,1,0,0,1,1,0,0,0,1 /
      DATA (IRRED( 79,I),I=-1,9) /  9,1,1,0,1,0,0,1,0,0,1 /
      DATA (IRRED( 80,I),I=-1,9) /  9,1,0,0,1,1,0,1,0,0,1 /
      DATA (IRRED( 81,I),I=-1,9) /  9,1,1,1,1,1,0,1,0,0,1 /
      DATA (IRRED( 82,I),I=-1,9) /  9,1,0,1,0,0,1,1,0,0,1 /
      DATA (IRRED( 83,I),I=-1,9) /  9,1,0,0,1,0,1,1,0,0,1 /
      DATA (IRRED( 84,I),I=-1,9) /  9,1,1,1,1,0,1,1,0,0,1 /
      DATA (IRRED( 85,I),I=-1,9) /  9,1,1,1,0,1,1,1,0,0,1 /
      DATA (IRRED( 86,I),I=-1,9) /  9,1,0,1,1,1,1,1,0,0,1 /
      DATA (IRRED( 87,I),I=-1,9) /  9,1,1,1,0,0,0,0,1,0,1 /
      DATA (IRRED( 88,I),I=-1,9) /  9,1,0,1,0,1,0,0,1,0,1 /
      DATA (IRRED( 89,I),I=-1,9) /  9,1,0,0,1,1,0,0,1,0,1 /
      DATA (IRRED( 90,I),I=-1,9) /  9,1,1,0,0,0,1,0,1,0,1 /
      DATA (IRRED( 91,I),I=-1,9) /  9,1,0,1,0,0,1,0,1,0,1 /
      DATA (IRRED( 92,I),I=-1,9) /  9,1,1,1,1,0,1,0,1,0,1 /
      DATA (IRRED( 93,I),I=-1,9) /  9,1,1,1,0,1,1,0,1,0,1 /
      DATA (IRRED( 94,I),I=-1,9) /  9,1,0,1,1,1,1,0,1,0,1 /
      DATA (IRRED( 95,I),I=-1,9) /  9,1,1,1,1,0,0,1,1,0,1 /
      DATA (IRRED( 96,I),I=-1,9) /  9,1,0,0,0,1,0,1,1,0,1 /
      DATA (IRRED( 97,I),I=-1,9) /  9,1,1,0,1,1,0,1,1,0,1 /
      DATA (IRRED( 98,I),I=-1,9) /  9,1,0,1,0,1,1,1,1,0,1 /
      DATA (IRRED( 99,I),I=-1,9) /  9,1,0,0,1,1,1,1,1,0,1 /
      DATA (IRRED(100,I),I=-1,9) /  9,1,0,0,0,0,0,0,0,1,1 /

C
C   Prepare to work in Z2
C
      CALL SETFLD (2)
C
      DO 1000 I = 1, DIMEN
C
C   For each dimension, we need to calculate powers of an
C   appropriate irreducible polynomial :  see Niederreiter
C   page 65, just below equation (19).
C     Copy the appropriate irreducible polynomial into PX,
C   and its degree into E.  Set polynomial B = PX ** 0 = 1.
C   M is the degree of B.  Subsequently B will hold higher
C   powers of PX.
C
        E = IRRED(I, DEG)
        DO 10 J = -1, E
          PX(J) = IRRED(I,J)
   10   CONTINUE
        B(DEG) = 0
        B(0) = 1
C
C   Niederreiter (page 56, after equation (7), defines two
C   variables Q and U.  We do not need Q explicitly, but we
C   do need U.
C
        U = 0
C
        DO 90 J = 1, NBITS
C
C   If U = 0, we need to set B to the next power of PX
C   and recalculate V.  This is done by subroutine CALCV.
C
          IF (U .EQ. 0) CALL CALCV (PX, B, V, MAXV)
C
C Now C is obtained from V.  Niederreiter
C obtains A from V (page 65, near the bottom), and then gets
C C from A (page 56, equation (7)).  However this can be done
C in one step.  Here CI(J,R) corresponds to
C Niederreiter's C(I,J,R).
C
          DO 50 R = 0, NBITS-1
            CI(J,R) = V(R+U)
   50     CONTINUE
C
C Increment U.  If U = E, then U = 0 and in Niederreiter's
C paper Q = Q + 1.  Here, however, Q is not used explicitly.
C
          U = U + 1
          IF (U .EQ. E) U = 0
  90    CONTINUE
C
C  The array CI now holds the values of C(I,J,R) for this value
C  of I.  We pack them into array CJ so that CJ(I,R) holds all
C  the values of C(I,J,R) for J from 1 to NBITS.
C
        DO 120 R = 0, NBITS-1
          TERM = 0
          DO 110 J = 1, NBITS
            TERM = 2 * TERM + CI(J,R)
  110     CONTINUE
          CJ(I,R) = TERM
  120   CONTINUE
C
 1000 CONTINUE
      END
C
C      *****  end of CALCC2

C =======================================================

      SUBROUTINE CALCV (PX, B, V, MAXV)
C
C   This version :  12 February 1991
C
C   See the general comments on implementing Niederreiter's
C   low-discrepancy sequences.
C
C   This program calculates the values of the constants V(J,R) as
C   described in BFN section 3.3.  It is called from either CALCC or
C   CALCC2.  The values transmitted through common /FIELD/ determine
C   which field we are working in.
C
C INPUT :
C   PX is the appropriate irreducible polynomial for the dimension
C     currently being considered.  The degree of PX will be called E.
C   B is the polynomial defined in section 2.3 of BFN.  On entry to
C     the subroutine, the degree of B implicitly defines the parameter
C     J of section 3.3, by degree(B) = E*(J-1).
C   MAXV gives the dimension of the array V.
C   On entry, we suppose that the common block /FIELD/ has been set
C     up correctly (using SETFLD).
C
C OUTPUT :
C   On output B has been multiplied by PX, so its degree is now E*J.
C   V contains the values required.
C
C USES :
C   The subroutine PLYMUL is used to multiply polynomials.
C
C
C   ------------------------------------------------------------
C
C   The following COMMON block, used by many subroutines,
C   gives the order Q of a field, its characteristic P, and its
C   addition, multiplication, and subtraction tables.
C   The parameter MAXQ gives the order of the largest field to
C   be handled.
C
      INTEGER MAXQ
      PARAMETER (MAXQ=50)

      INTEGER P, Q, ADD(0:MAXQ-1,0:MAXQ-1)
      INTEGER MUL(0:MAXQ-1, 0:MAXQ-1), SUB(0:MAXQ-1, 0:MAXQ-1)
      COMMON /FIELD/ P, Q, ADD, MUL, SUB
      SAVE /FIELD/
C
C   The following definitions concern the representation of
C   polynomials.
C
      INTEGER MAXDEG, DEG
      PARAMETER (MAXDEG=50, DEG=-1)
C
C   The parameter MAXDEG gives the highest degree of polynomial
C   to be handled.  Polynomials stored as arrays have the
C   coefficient of degree n in POLY(N), and the degree of the
C   polynomial in POLY(-1).  The parameter DEG is just to remind
C   us of this last fact.  A polynomial which is identically 0
C   is given degree -1.
C
C   A polynomial can also be stored in an integer I, with
C        I = AN*Q**N + ... + A0.
C   Routines ITOP and PTOI convert between these two formats.
C
C   -----------------------------------------------------------
C
C
C
      INTEGER MAXV, E, I, J, KJ, M, BIGM, R, TERM
      INTEGER PX(-1:MAXDEG), B(-1:MAXDEG), V(0:MAXV)
      INTEGER H(-1:MAXDEG)
C
      INTEGER ARBIT, NONZER
      ARBIT() = 1
C
C   We use ARBIT() to indicate where the user can place
C   an arbitrary element of the field of order Q, while NONZER
C   shows where he must put an arbitrary non-zero element
C   of the same field.  For the code,
C   this means 0 <= ARBIT < Q and 0 < NONZER < Q.  Within these
C   limits, the user can do what he likes.  ARBIT is declared as
C   a function as a reminder that a different arbitrary value may
C   be returned each time ARBIT is referenced.
C
C    BIGM is the M used in section 3.3.
C    It differs from the [little] m used in section 2.3,
C    denoted here by M.
C
      NONZER = 1
C
      E = PX(DEG)
C
C   The poly H is PX**(J-1), which is the value of B on arrival.
C   In section 3.3, the values of Hi are defined with a minus sign :
C   don't forget this if you use them later !
C
      DO 10 I = -1, B(DEG)
   10   H(I) = B(I)
      BIGM = H(DEG)
C
C   Now multiply B by PX so B becomes PX**J.
C   In section 2.3, the values of Bi are defined with a minus sign :
C   don't forget this if you use them later !
C
      CALL PLYMUL (PX, B, B)
      M = B(DEG)
C
C   We don't use J explicitly anywhere, but here it is just in case.
C
      J = M / E
C
C   Now choose a value of Kj as defined in section 3.3.
C   We must have 0 <= Kj < E*J = M.
C   The limit condition on Kj does not seem very relevant
C   in this program.
C
      KJ = BIGM
C
C   Now choose values of V in accordance with the conditions in
C   section 3.3
C
      DO 20 R = 0, KJ-1
   20   V(R) = 0
      V(KJ) = 1
C
      IF (KJ .LT. BIGM) THEN
C
        TERM = SUB (0, H(KJ))
C
        DO 30 R = KJ+1, BIGM-1
          V(R) = ARBIT()
C
C         Check the condition of section 3.3,
C         remembering that the H's have the opposite sign.
C
          TERM = SUB (TERM, MUL (H(R), V(R)))
   30   CONTINUE
C
C         Now V(BIGM) is anything but TERM
C
          V(BIGM) = ADD (NONZER, TERM)
C
        DO 40 R = BIGM+1, M-1
   40     V(R) = ARBIT()
C
      ELSE
C       This is the case KJ .GE. BIGM
C
        DO 50 R = KJ+1, M-1
   50     V(R) = ARBIT()
C
      ENDIF
C
C   Calculate the remaining V's using the recursion of section 2.3,
C   remembering that the B's have the opposite sign.
C
      DO 70 R = 0, MAXV-M
        TERM = 0
        DO 60 I = 0, M-1
          TERM = SUB (TERM, MUL (B(I), V(R+I)))
   60   CONTINUE
        V(R+M) = TERM
   70 CONTINUE
C
      RETURN
      END
C
C     ***** end of SUBROUTINE CALCV

C ========================================================

      SUBROUTINE PLYMUL (PA, PB, PC)
C
C   This version :  12 December 1991
C
C
C   ------------------------------------------------------------
C
C   The following COMMON block, used by many subroutines,
C   gives the order Q of a field, its characteristic P, and its
C   addition, multiplication and subtraction tables.
C   The parameter MAXQ gives the order of the largest field to
C   be handled.
C
      INTEGER MAXQ
      PARAMETER (MAXQ=50)

      INTEGER P, Q, ADD(0:MAXQ-1,0:MAXQ-1)
      INTEGER MUL(0:MAXQ-1, 0:MAXQ-1), SUB(0:MAXQ-1, 0:MAXQ-1)
      COMMON /FIELD/ P, Q, ADD, MUL, SUB
      SAVE /FIELD/
C
C   The following definitions concern the representation of
C   polynomials.
C
      INTEGER MAXDEG, DEG
      PARAMETER (MAXDEG=50, DEG=-1)
C
C   The parameter MAXDEG gives the highest degree of polynomial
C   to be handled.  Polynomials stored as arrays have the
C   coefficient of degree n in POLY(N), and the degree of the
C   polynomial in POLY(-1).  The parameter DEG is just to remind
C   us of this last fact.  A polynomial which is identically 0
C   is given degree -1.
C
C   A polynomial can also be stored in an integer I, with
C        I = AN*Q**N + ... + A0.
C   Routines ITOP and PTOI convert between these two formats.
C
C   -----------------------------------------------------------
C
C
C
      INTEGER I, J, DEGA, DEGB, DEGC, TERM
      INTEGER PA(-1:MAXDEG), PB(-1:MAXDEG), PC(-1:MAXDEG)
      INTEGER PT(-1:MAXDEG)
C
C   Multiplies polynomial PA by PB putting the result in PC.
C   Coefficients are elements of the field of order Q.
C
      DEGA = PA(DEG)
      DEGB = PB(DEG)
      IF (DEGA .EQ. -1 .OR. DEGB .EQ. -1) THEN
        DEGC = -1
      ELSE
        DEGC = DEGA + DEGB
      ENDIF
      IF (DEGC .GT. MAXDEG) THEN
        WRITE (*,*) ' PLYMUL :  Degree of product exceeds MAXDEG'
        STOP
      ENDIF
C
      DO 20 I = 0, DEGC
        TERM = 0
        DO 10 J = MAX(0, I-DEGA), MIN(DEGB, I)
   10     TERM = ADD(TERM, MUL(PA(I-J), PB(J)))
   20   PT(I) = TERM
C
      PC(DEG) = DEGC
      DO 30 I = 0, DEGC
   30   PC(I) = PT(I)
      DO 40 I = DEGC+1, MAXDEG
   40   PC(I) = 0
      RETURN
      END
C
C     *****   end of SUBROUTINE PLYMUL

C =======================================================

      SUBROUTINE SETFLD (QIN)
      INTEGER QIN
C
C   This version : 12 December 1991
C
C   This subroutine sets up addition, multiplication, and
C   subtraction tables for the finite field of order QIN.
C   If necessary, it reads precalculated tables from the file
C   'gftabs.dat' using unit 1.  These precalculated tables
C   are supposed to have been put there by GFARIT.
C
C      *****  For the base-2 programs, these precalculated
C      *****  tables are not needed and, therefore, neither
C      *****  is GFARIT.
C
C
C   Unit 1 is closed both before and after the call of SETFLD.
C
C USES
C   Integer function CHARAC gets the characteristic of a field.
C
C
C   ------------------------------------------------------------
C
C   The following COMMON block, used by many subroutines,
C   gives the order Q of a field, its characteristic P, and its
C   addition, multiplication and subtraction tables.
C   The parameter MAXQ gives the order of the largest field to
C   be handled.
C
      INTEGER MAXQ
      PARAMETER (MAXQ=50)

      INTEGER P, Q, ADD(0:MAXQ-1,0:MAXQ-1)
      INTEGER MUL(0:MAXQ-1, 0:MAXQ-1), SUB(0:MAXQ-1, 0:MAXQ-1)
      COMMON /FIELD/ P, Q, ADD, MUL, SUB
      SAVE /FIELD/
C
C   The following definitions concern the representation of
C   polynomials.
C
      INTEGER MAXDEG, DEG
      PARAMETER (MAXDEG=50, DEG=-1)
C
C   The parameter MAXDEG gives the highest degree of polynomial
C   to be handled.  Polynomials stored as arrays have the
C   coefficient of degree n in POLY(N), and the degree of the
C   polynomial in POLY(-1).  The parameter DEG is just to remind
C   us of this last fact.  A polynomial which is identically 0
C   is given degree -1.
C
C   A polynomial can also be stored in an integer I, with
C        I = AN*Q**N + ... + A0.
C   Routines ITOP and PTOI convert between these two formats.
C
C   -----------------------------------------------------------
C
C
C
      INTEGER I, J, N, CHARAC
C
      IF (QIN .LE. 1 .OR. QIN .GT. MAXQ) THEN
        WRITE (*,*) ' SETFLD :  Bad value of Q'
        STOP
      ENDIF
C
      Q = QIN
      P = CHARAC(Q)
C
      IF (P .EQ. 0) THEN
        WRITE (*,*) ' SETFLD :  There is no field of order', Q
        STOP
      ENDIF
C
C Set up to handle a field of prime order :  calculate ADD and MUL.
C
      IF (P .EQ. Q) THEN
        DO 10 I = 0, Q-1
          DO 10 J = 0, Q-1
            ADD(I,J) = MOD(I+J, P)
            MUL(I,J) = MOD(I*J, P)
   10   CONTINUE
C
C Set up to handle a field of prime-power order :  tables for
C ADD and MUL are in the file 'gftabs.dat'.
C
      ELSE
        OPEN (UNIT=1, FILE='gftabs.dat', STATUS='old')
C
C    *****  OPEN statements are system dependent.
C
   20   READ (1, 900, END=500) N
  900   FORMAT (20I3)
        DO 30 I = 0, N-1
          READ (1, 900) (ADD(I,J), J = 0, N-1)
   30   CONTINUE
        DO 40 I = 0, N-1
          READ (1, 900) (MUL(I,J), J = 0, N-1)
   40   CONTINUE
        IF (N .NE. Q) GOTO 20
        CLOSE (1)
      ENDIF
C
C Now use the addition table to set the subtraction table.
C
      DO 60 I = 0, Q-1
        DO 50 J = 0, Q-1
          SUB(ADD(I,J), I) = J
   50   CONTINUE
   60 CONTINUE
      RETURN
C
  500 WRITE (*,*) ' SETFLD :  Tables for q =', Q, ' not found'
      STOP
C
      END
C
C     ***** end of SUBROUTINE SETFLD

C ===========================================================

      INTEGER FUNCTION CHARAC (QIN)
C
C   This version :  12 December 1991
C
C   This function gives the characteristic for a field of
C   order QIN.  If no such field exists, or if QIN is out of
C   the range we can handle, returns 0.
C
C
C   ------------------------------------------------------------
C
C   The following COMMON block, used by many subroutines,
C   gives the order Q of a field, its characteristic P, and its
C   addition, multiplication and subtraction tables.
C   The parameter MAXQ gives the order of the largest field to
C   be handled.
C
      INTEGER MAXQ
      PARAMETER (MAXQ=50)

      INTEGER P, Q, ADD(0:MAXQ-1,0:MAXQ-1)
      INTEGER MUL(0:MAXQ-1, 0:MAXQ-1), SUB(0:MAXQ-1, 0:MAXQ-1)
      COMMON /FIELD/ P, Q, ADD, MUL, SUB
      SAVE /FIELD/
C
C   The following definitions concern the representation of
C   polynomials.
C
      INTEGER MAXDEG, DEG
      PARAMETER (MAXDEG=50, DEG=-1)
C
C   The parameter MAXDEG gives the highest degree of polynomial
C   to be handled.  Polynomials stored as arrays have the
C   coefficient of degree n in POLY(N), and the degree of the
C   polynomial in POLY(-1).  The parameter DEG is just to remind
C   us of this last fact.  A polynomial which is identically 0
C   is given degree -1.
C
C   A polynomial can also be stored in an integer I, with
C        I = AN*Q**N + ... + A0.
C   Routines ITOP and PTOI convert between these two formats.
C
C   -----------------------------------------------------------
C
C
C
      INTEGER QIN, CH(MAXQ)
      SAVE CH
C
      DATA CH  / 0,  2,  3,  2,  5,  0,  7,  2,  3,  0,
     1          11,  0, 13,  0,  0,  2, 17,  0, 19,  0,
     2           0,  0, 23,  0,  5,  0,  3,  0, 29,  0,
     3          31,  2,  0,  0,  0,  0, 37,  0,  0,  0,
     4          41,  0, 43,  0,  0,  0, 47,  0,  7,  0/
C
      IF (QIN .LE. 1 .OR. QIN .GT. MAXQ) THEN
        CHARAC = 0
      ELSE
        CHARAC = CH(QIN)
      ENDIF
C
      END
C
C     ***** end of INTEGER FUNCTION CHARAC