PROGRAM GENIN
C
C ***** This is the driver for the general base programs.
C ***** More accurately, it is a driver skeleton.
C ***** There is a default set of test integrals
C ***** in TESTF. The user can replace TESTF with
C ***** another subroutine called TESTF containing
C ***** integrals of interest to him.
C
C This version : 14 Aug 1992
C
C This program tests the accuracy of numerical integration
C using the low-discrepancy binary sequences of
C H. Niederreiter (1988) as implemented in INLO, GOLO, and
C related programs. Various possible test integrals are
C provided by the function TESTF.
C
C Interactive input and output go through the Fortran units
C denoted by * ; at most installations, * corresponds to
C the user's terminal. For a prime-power base, Fortran unit 1
C is used to read field arithmetic tables ; for any base, it
C is used to read irreducible polynomials.
C Fortran unit 2 is used to save the output.
C
C Our VAX implementation does not measure elapsed time.
C It can be modified, in a system-dependent way, to do so.
C
C GENIN and its associated subroutines are listed below.
C An asterisk denotes a subroutine also used by GFARIT
C (to set up the field-arithmetic tables gftabs.txt) and
C by GFPLYS (to set up the irreducible polynomials in
C irrtabs.txt).
C
C GENIN
C INLO
C CALCC
C CALCV %
C CHARAC * %
C SETFLD * %
C PLYMUL * %
C GOLO
C TESTF %
C
C A percent sign above denotes a routine also used
C in the set of programs tailored for base 2.
C
C
C Both the general-base and base-2 programs assume that
C your computer's word length is 31 bits, excluding sign.
C If this is not the case, modify the parameter NBITS
C throughout the PARAMETER statements in this set of
C programs accordingly.
C
INTEGER MAXDIM, OUTUNT, MAXBAS, READY
PARAMETER (MAXDIM=12, OUTUNT=2, MAXBAS = 13)
C
C The parameter MAXDIM gives the maximum dimension that will
C be used. OUTUNT is the unit to save the output.
C MAXBAS gives the maximum asymptotically-optimal base
C used up to MAXDIM. The user enters an appropriate value
C of READY depending on whether or not the required files
C indicated above have been set up.
C
INTEGER I, NUM, DIMEN, SEQLEN, SKIP, STEP, BASE, CHARAC
INTEGER OPTBAS(2:MAXDIM), PBASE, POWER(2:MAXBAS)
REAL QUASI(MAXDIM), TESTF, EXACTF
DOUBLE PRECISION SUM
C
C The DATA statement below gives the asymptotically-optimal
C base for the respective dimension.
C
DATA (OPTBAS(I), I = 2,MAXDIM) / 2,3,3,5,7,7,9,9,11,11,13 /
C
C
C The data statement below gives values used in a possible
C warm-up calculation.
C
DATA (POWER(I), I = 2,MAXBAS) / 12,8,8,6,6,6,4,4,4,4,4,4 /
C
C There are theoretical reasons to believe that BASE ** e,
C where e is defined for example in Bratley, Fox, and
C Niederreiter (1991), would be a good choice. However,
C we don't want to come anywhere near the largest possible
C machine-representable integer; hence, the compromise
C exponents above. Note: subject to this conditon,
C it can't hurt to take an exponent greater than e, because
C warm-up skipping of initial values is done implicitly
C in O(1) time. The maximum value of e for a fixed dimension
C s grows like log s. We allow some "fat" for the implicit
C constant in our choice of POWER.
C
WRITE (*,*) ' This is program GENIN'
C
WRITE(*,*) ' If you wish to use the base 2, the '
WRITE(*,*) ' alternative set of programs tailored for '
WRITE(*,*) ' base 2 will run much faster. '
C
WRITE (*,*) ' If the files gftabs.txt and irrtab.txt '
WRITE (*,*) ' have not already been set up, then '
WRITE (*,*) ' consult the guide to the programs for '
WRITE (*,*) ' how to do so, set up those files per '
WRITE (*,*) ' the guide, and [temporarily] quit this '
WRITE (*,*) ' set of programs by entering the integer 0 '
WRITE (*,*) ' below; otherwise, enter any other integer. '
WRITE (*,*) ' ENTER an appropriate integer. '
READ (*,*) READY
IF (READY .EQ. 0) THEN
WRITE (*,*) ' Set up gftabs.txt and irrtab.txt now. '
WRITE (*,*) ' Exit GENIN. '
STOP
ENDIF
C
WRITE (*,*) 'If the number of bit per word, excluding sign'
WRITE (*,*) 'is not 31, enter the integer 0 below '
WRITE (*,*) 'and fix the parameter NBITS everywhere; '
WRITE (*,*) 'otherwise, enter a positive integer below.'
WRITE (*,*) 'ENTER an appropriate integer. '
READ (*,*) READY
IF (READY .EQ. 0) THEN
WRITE(*,*) 'Fix NBITS.'
WRITE(*,*) 'Exit GENIN.'
STOP
ENDIF
C
WRITE (*,*) ' Output file name is OUTFIL.DAT'
OPEN (unit = OUTUNT, file = 'OUTFIL.DAT', status = 'UNKNOWN')
C
C ***** OPEN statements are system-dependent.
C ***** Therefore the statement above may have to be
C ***** modified for use at your computer center.
C
C
5 WRITE (*,*) ' Choose a test integral (1 to 4) or 0 to quit :'
READ (*,*) NUM
IF (NUM.LE.0) THEN
WRITE (*,*) ' End of program GENIN'
CLOSE (OUTUNT)
STOP
ENDIF
IF (NUM.GT.4) THEN
WRITE (*,*) ' No such test integral'
GOTO 5
ENDIF
C
C ***** Each test integral is parameterized by
C ***** its dimension.
C
10 WRITE (*,*) ' Enter dimension :'
READ (*,*) DIMEN
IF (DIMEN.GT.MAXDIM) THEN
WRITE (*,*) ' Dimension may not exceed', MAXDIM
GOTO 10
ENDIF
C
12 WRITE (*,*) ' Choose a prime or prime-power base .'
WRITE (*,*) ' The asymptotically-optimal base '
WRITE (*,*) ' for this dimension is OPTBAS(DIMEN) = '
WRITE (*,*) OPTBAS(DIMEN)
WRITE (*,*) ' This base may not be empirically optimal. '
WRITE (*,*) ' Enter BASE: '
C
READ (*,*) BASE
IF (CHARAC(BASE) .EQ. 0) THEN
WRITE (*,*) ' Out of range or bad value : try again'
GOTO 12
ENDIF
C
C
C ***** The sequence length is the number of
C ***** quasi-random points used to estimate the
C ***** integral, excluding warm-up.
C ***** The number of initial quasi-random points
C ***** deleted during warm-up is given by SKIP,
C ***** chosen below.
C
C ***** Some users may wish to rewrite
C ***** the driver to test a [heuristic] "convergence"
C ***** criterion, stopping the generation of points
C ***** when it is passed or when a specified number of
C ***** points have been generated -- whichever occurs
C ***** first.
C
15 WRITE (*,*) ' Choose sequence length :'
WRITE (*,*) ' A power of the base is recommended; e.g., '
WRITE (*,*) BASE ** POWER(BASE)
WRITE (*,*) BASE ** ((POWER(BASE) + 1))
WRITE (*,*) BASE ** ((POWER(BASE) + 2))
WRITE (*,*) BASE ** ((POWER(BASE) + 3))
WRITE (*,*) ' Enter SEQLEN '
READ (*,*) SEQLEN
IF (SEQLEN.LE.0) THEN
WRITE (*,*) ' Length must be strictly positive'
GOTO 15
ENDIF
C
20 WRITE (*,*) ' Choose the number of values to skip.'
WRITE (*,*) ' One possibility is given by the heuristic '
WRITE (*,*) ' formula SKIP = BASE ** POWER(BASE) '
WRITE (*,*) ' when BASE <= MAXBAS, = 10000 otherwise '
IF (BASE .LE. MAXBAS) THEN
SKIP = BASE ** POWER(BASE)
ELSE
SKIP = 10000
ENDIF
WRITE (*,*) ' Numerically, this equals ', SKIP
WRITE (8,*) ' Enter SKIP (not necessarily the value above) : ' C
READ (*,*) SKIP
IF (SKIP.LT.0) THEN
WRITE (*,*) ' Number must be nonnegative'
GOTO 20
ENDIF
C
C
C
CALL INLO (DIMEN, BASE, SKIP)
WRITE (*,*) ' GENIN : Initialization complete'
C
C Write title and the exact value of the integral
C
WRITE (OUTUNT,27) NUM
27 FORMAT (/,' Test integral ',I2)
WRITE (OUTUNT,28) DIMEN, BASE, SEQLEN, SKIP
28 FORMAT (/,' Dimension ',I6,', Base ', I9,
1 /,' Sequence ',I8,', Skipped ',I7)
WRITE (OUTUNT,29) EXACTF(NUM, DIMEN)
29 FORMAT (/,' Correct value is ',G16.7)
WRITE (OUTUNT,30)
30 FORMAT(/,' Iteration Estimate of integral',/)
C
C Now estimate the integral
C
WRITE (*,*) ' The odd-looking iteration numbers '
WRITE (*,*) ' in the output are powers of the base. '
C
PBASE = BASE
SUM = 0
STEP = 500
DO 100 I = 1, SEQLEN
CALL GOLO (QUASI)
SUM = SUM + TESTF(NUM, DIMEN, QUASI)
IF (MOD(I,STEP).EQ.0) THEN
WRITE (OUTUNT,99) I, SUM/I
ENDIF
IF (MOD(I,PBASE) .EQ. 0) THEN
WRITE (OUTUNT,99) I, SUM/I
PBASE = PBASE * BASE
C
C This finds the next power of the base.
C There is reason to believe that convergenence
C properties of the sequence of estimates is
C better along the subsequence corrsponding to
C powers of the base.
C
ENDIF
99 FORMAT (I12,G24.7)
IF (I .EQ. 5000) STEP = 1000
IF (I .EQ. 10000) STEP = 5000
100 CONTINUE
C
WRITE (*,*) ' GENIN : iteration ends'
GOTO 5
C
END
C
C ***** end of PROGRAM GENIN
SUBROUTINE INLO (DIM, BASE, 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 initialization necessary
C before calling GOLO.
C
C INPUT :
C DIMEN - The dimension of the sequence to be generated.
C {In the argument of INLO, DIMEN is called DIM because
C DIMEN is subsequently passed via COMMON and is called
C DIMEN there.}
C
C BASE - The prime or prime-power base to be used.
C SKIP - The number of values to throw away at the beginning
C of the sequence.
C
C OUTPUT :
C To GOLO, labelled common /COMM/.
C
C USES :
C Calls CALCC to calculate the values of CJ.
C Calls SETFLD to set up field arithmetic tables.
C Calls CHARAC to check that base is a prime or prime-power
C in the range we can handle.
C
C -------------------------------------------------------------
C
C
C This segment defines the common block /COMM/ and some
C associated parameters. These are for use in the general base
C version of the generator.
C
INTEGER MAXDIM, MAXFIG, NBITS
PARAMETER (MAXDIM=12, MAXFIG=20, NBITS=31)
C
C The parameter MAXDIM is the maximum dimension that will be used.
C MAXFIG is the maximum number of base-Q digits we can handle.
C MAXINT is the largest fixed point integer we can represent.
C NBITS is the number of bits in a fixed-point integer, not
C counting the sign.
C ***** NBITS is machine dependent *****
C
INTEGER C(MAXDIM, MAXFIG, 0:MAXFIG-1)
INTEGER COUNT(0:MAXFIG-1), D(MAXDIM, MAXFIG)
INTEGER NEXTQ(MAXDIM), QPOW(MAXFIG)
INTEGER DIMEN, NFIGS
REAL RECIP
COMMON /COMM/ C, COUNT, D, NEXTQ, QPOW, DIMEN, NFIGS, RECIP
SAVE /COMM/
C
C The common block /COMM/ :
C C - Contains the values of Niederreiter's C(I,J,R)
C COUNT - The index of the current item in the sequence,
C expressed as an array of base-Q digits. COUNT(R)
C is the same as Niederreiter's AR(N) (page 54)
C except that N is implicit.
C D - The values of D(I,J).
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.
C QPOW - To speed things up a bit. QPOW(I) = Q ** (NFIGS-I).
C DIMEN - The dimension of the sequence to be generated
C NFIGS - The number of base-Q digits we are using.
C RECIP - 1.0 / (Q ** NFIGS)
C
C Array C of the common block is set up by subroutine CALCC.
C The other items in the common block are set up by INLO.
C
C ------------------------------------------------------------
C
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, NQ, R, DIM, SKIP , BASE, CHARAC
REAL TEMP
C
DIMEN = DIM
C
C This assignment just relabels the variable
C for subsequent use.
C
IF (DIMEN.LE.0 .OR. DIMEN.GT.MAXDIM) THEN
WRITE (*,*) ' INLO : Bad dimension'
STOP
ENDIF
C
C
IF (CHARAC(BASE) .EQ. 0) THEN
WRITE (*,*) ' INLO : Base not prime power or out of range'
STOP
ENDIF
C
CALL SETFLD (BASE)
C
C Calculate how many figures to use in base Q = BASE
C
TEMP = LOG(2.0 ** NBITS - 1)/LOG(REAL(Q))
NFIGS = MIN(MAXFIG, INT(TEMP))
C
CALL CALCC
C
C Set up RECIP and QPOW(I) = Q ** (NFIGS-I)
C
RECIP = 1.0 / (Q ** NFIGS)
QPOW(NFIGS) = 1
DO 5 I = NFIGS-1, 1, -1
5 QPOW(I) = Q * QPOW(I+1)
C
C Initialize COUNT
C
I = SKIP
DO 10 R = 0, NFIGS-1
COUNT(R) = MOD(I, Q)
I = I / Q
10 CONTINUE
IF (I .NE. 0) THEN
WRITE (*,*) ' INLO : SKIP too long'
STOP
ENDIF
C
C Initialize D
C
DO 20 I = 1, DIMEN
DO 20 J = 1, NFIGS
20 D(I,J) = 0
C
DO 50 R = 0, NFIGS-1
IF (COUNT(R) .NE. 0) THEN
DO 40 I = 1, DIMEN
DO 30 J = 1, NFIGS
D(I,J) = ADD (D(I,J), MUL (C(I,J,R), COUNT(R)))
30 CONTINUE
40 CONTINUE
ENDIF
50 CONTINUE
C
C Initialize NEXTQ
C
DO 70 I = 1, DIMEN
NQ = 0
DO 60 J = 1, NFIGS
NQ = NQ + D(I,J) * QPOW(J)
60 CONTINUE
NEXTQ(I) = NQ
70 CONTINUE
C
RETURN
END
C
C ***** end of SUBROUTINE INLO
SUBROUTINE CALCC
C
C This version : 12 February 1992
C
C See the general comments on implementing Niederreiter's
C low-discrepancy sequences.
C
C This routine calculates the values of the constants C(I,J,R).
C As far as possible, we use Niederreiter's notation.
C We calculate the values of C for each I in turn.
C When all the values of C have been calculated, we return
C this array to the calling program.
C
C Irreducible polynomials are read from file 'irrtabs.txt'
C through Fortran unit 1. These polys must have been put on
C the file beforehand by GFPOLYS. Unit 1 is closed before
C entry to CALCC and after returning from the subroutine.
C
C Thanks to Michael Baudin for pointing out that MAXE should
C be increased from 5 to 7, since one of the irreducible polynomials
C that must be stored in PX has degree 7, 07 June 2010.
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)
C
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 This segment defines the common block /COMM/ and some
C associated parameters. These are for use in the general base
C version of the generator.
C
INTEGER MAXDIM, MAXFIG, NBITS
PARAMETER (MAXDIM=12, MAXFIG=20, NBITS=31)
C
C The parameter MAXDIM is the maximum dimension that will be used.
C MAXFIG is the maximum number of base-Q digits we can handle.
C MAXINT is the largest fixed point integer we can represent.
C NBITS is the number of bits in a fixed-point integer, not
C counting the sign.
C ***** NBITS is machine dependent *****
C
INTEGER C(MAXDIM, MAXFIG, 0:MAXFIG-1)
INTEGER COUNT(0:MAXFIG-1), D(MAXDIM, MAXFIG)
INTEGER NEXTQ(MAXDIM), QPOW(MAXFIG)
INTEGER DIMEN, NFIGS
REAL RECIP
COMMON /COMM/ C, COUNT, D, NEXTQ, QPOW, DIMEN, NFIGS, RECIP
SAVE /COMM/
C
C The common block /COMM/ :
C C - Contains the values of Niederreiter's C(I,J,R)
C COUNT - The index of the current item in the sequence,
C expressed as an array of base-Q digits. COUNT(R)
C is the same as Niederreiter's AR(N) (page 54)
C except that N is implicit.
C D - The values of D(I,J).
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.
C QPOW - To speed things up a bit. QPOW(I) = Q ** (NFIGS-I).
C DIMEN - The dimension of the sequence to be generated
C NFIGS - The number of base-Q digits we are using.
C RECIP - 1.0 / (Q ** NFIGS)
C
C Array C of the common block is set up by subroutine CALCC.
C The other items in the common block are set up by INLO.
C
C ------------------------------------------------------------
C
C
C
C
C
C MAXE - We need MAXDIM irreducible polynomials over GF(Q).
C MAXE is the highest degree among these.
C MAXV - The maximum index used in V.
C GFUNIT - The unit number to read field tables.
C NPOLS - The number of precalculated irreducible polys.
C
INTEGER MAXE, MAXV, GFUNIT, NPOLS
C PARAMETER (MAXE=5, GFUNIT=1, NPOLS=25)
PARAMETER (MAXE=7, GFUNIT=1, NPOLS=25)
PARAMETER (MAXV = MAXFIG + MAXE)
C
C INPUT :
C DIMEN, the number of dimensions in use, and NFIGS, the number
C of base-Q figures to be used, are passed in through common COMM.
C Necessary field arithmetic tables are passed through common
C FIELD.
C
C OUTPUT
C The array C is returned to the calling program.
C
C USES
C Subroutine CALCV is used for the auxiliary calculation
C of the values V needed to get the Cs.
C
INTEGER PX(-1:MAXE), B(-1:MAXDEG)
INTEGER V(0:MAXV)
INTEGER E, I, J, R, U
C
C Prepare to read irreducible polynomials on unit 1.
C
OPEN (UNIT=GFUNIT, FILE='irrtabs.txt', STATUS='old')
C
C ***** OPEN statements are system dependent
C
10 READ (GFUNIT, 900, END=500) I
900 FORMAT (20I3)
IF (I .NE. Q) THEN
DO 20 J = 1, NPOLS
20 READ (GFUNIT, 900)
GOTO 10
ENDIF
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 Read 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 The polynomial PX is stored in file 'irrtabs.txt' in the
C format
C n a0 a1 a2 ... an
C where n is the degree of the polynomial and the ai are
C its coefficients.
C
READ (GFUNIT, 900) E, (PX(J), J = 0,E)
PX(DEG) = E
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, NFIGS
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. Neiderreiter
C obtains A from V (page 65, near the bottom), and
C then gets C from A (page 56, equation (7)).
C However this can be done in one step.
C
DO 50 R = 0, NFIGS-1
C(I,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
1000 CONTINUE
C
CLOSE (GFUNIT)
RETURN
C
500 WRITE (*,*) ' CALCC : Tables for q =', Q, ' not found'
STOP
END
C
C ***** end of SUBROUTINE CALCC
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
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
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.txt' 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.txt'.
C
ELSE
OPEN (UNIT=1, FILE='gftabs.txt', 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
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
SUBROUTINE GOLO (QUASI)
REAL QUASI(*)
C
C This version : 21 February 1992
C
C See the general comments on implementing Niederreiter's
C low-discrepancy sequences.
C
C Call subroutine INLO before calling GOLO. Thereafter
C GOLO generates a new quasi-random vector on each call,
C
C INPUT
C From INLO, labelled common /COMM/ and labelled common
C /FIELD/, properly initialized. (INLO calls SETFLD
C to initialize /FIELD/).
C
C OUTPUT
C To the user's program, the next vector in the sequence in
C array QUASI.
C
C -------------------------------------------------------------
C
C
C This segment defines the common block /COMM/ and some
C associated parameters. These are for use in the general base
C version of the generator.
C
INTEGER MAXDIM, MAXFIG, NBITS
PARAMETER (MAXDIM=12, MAXFIG=20, NBITS=31)
C
C The parameter MAXDIM is the maximum dimension that will be used.
C MAXFIG is the maximum number of base-Q digits we can handle.
C MAXINT is the largest fixed point integer we can represent.
C NBITS is the number of bits in a fixed-point integer, not
C counting the sign.
C ***** NBITS is machine dependent *****
C
INTEGER C(MAXDIM, MAXFIG, 0:MAXFIG-1)
INTEGER COUNT(0:MAXFIG-1), D(MAXDIM, MAXFIG)
INTEGER NEXTQ(MAXDIM), QPOW(MAXFIG)
INTEGER DIMEN, NFIGS
REAL RECIP
COMMON /COMM/ C, COUNT, D, NEXTQ, QPOW, DIMEN, NFIGS, RECIP
SAVE /COMM/
C
C The common block /COMM/ :
C C - Contains the values of Niederreiter's C(I,J,R)
C COUNT - The index of the current item in the sequence,
C expressed as an array of base-Q digits. COUNT(R)
C is the same as Niederreiter's AR(N) (page 54)
C except that N is implicit.
C D - The values of D(I,J).
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.
C QPOW - To speed things up a bit. QPOW(I) = Q ** (NFIGS-I).
C DIMEN - The dimension of the sequence to be generated
C NFIGS - The number of base-Q digits we are using.
C RECIP - 1.0 / (Q ** NFIGS)
C
C Array C of the common block is set up by subroutine CALCC.
C The other items in the common block are set up by INLO.
C
C ------------------------------------------------------------
C
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, R, OLDCNT, DIFF, NQ
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 Update COUNT, treated as a base-Q integer. Instead of
C recalculating the values of D from scratch, we update
C them for each digit of COUNT which changes. In terms of
C Niederreiter page 54, NEXTQ(I) corresponds to XI(N), with
C N being implicit, and D(I,J) corresponds to XI(N,J), again
C with N implicit. Finally COUNT(R) corresponds to AR(N).
C
R = 0
10 IF (R .EQ. NFIGS) THEN
WRITE (*,*) ' Too many calls on subroutine GOLO'
STOP
ENDIF
OLDCNT = COUNT(R)
IF (COUNT(R) .LT. Q-1) THEN
COUNT(R) = COUNT(R) + 1
ELSE
COUNT(R) = 0
ENDIF
DIFF = SUB(COUNT(R), OLDCNT)
C
C Digit R has just changed. DIFF says how much it changed
C by. We use this to update the values of array D.
C
DO 40 I = 1, DIMEN
DO 30 J = 1, NFIGS
D(I,J) = ADD(D(I,J), MUL(C(I,J,R), DIFF))
30 CONTINUE
40 CONTINUE
C
C If COUNT(R) is now zero, we must propagate the carry
C
IF (COUNT(R).EQ.0) THEN
R = R + 1
GOTO 10
ENDIF
C
C Now use the updated values of D to calculate NEXTQ.
C Array QPOW helps to speed things up a little :
C QPOW(J) is Q ** (NFIGS-J).
C
DO 60 I = 1, DIMEN
NQ = 0
DO 50 J = 1, NFIGS
NQ = NQ + D(I,J) * QPOW(J)
50 CONTINUE
NEXTQ(I) = NQ
60 CONTINUE
C
RETURN
END
C
C ***** end of SUBROUTINE GOLO
REAL FUNCTION TESTF (N, DIMEN, QUASI)
INTEGER I, N, DIMEN
REAL X, EXACTF, QUASI(*)
C
C This version : 4 Mar 1992
C
C Provides a variety of test integrals for quasi-random
C sequences. A call on TESTF computes an estimate of the
C integral ; a call on EXACTF computes the exact value.
C
GOTO (100, 200, 300, 400) N
C
ENTRY EXACTF (N, DIMEN)
GOTO (1100, 1200, 1300, 1400) N
C
C Test integral 1
C
100 TESTF = 1.0
DO 110 I = 1, DIMEN
TESTF = TESTF * ABS(4 * QUASI(I) - 2)
110 CONTINUE
RETURN
C
1100 EXACTF = 1.0
RETURN
C
C Test integral 2
C
200 TESTF = 1.0
DO 210 I = 1, DIMEN
TESTF = TESTF * I * COS(I * QUASI(I))
210 CONTINUE
RETURN
C
1200 EXACTF = 1.0
DO 1210 I = 1, DIMEN
EXACTF = EXACTF * SIN(FLOAT(I))
1210 CONTINUE
RETURN
C
C Test integral 3
C
300 TESTF = 1.0
DO 350 I = 1, DIMEN
X = 2 * QUASI(I) - 1
GOTO (310, 320, 330, 340) MOD(I, 4)
310 TESTF = TESTF * X
GOTO 350
320 TESTF = TESTF * (2*X*X - 1)
GOTO 350
330 TESTF = TESTF * (4*X*X - 3) * X
GOTO 350
340 X = X * X
TESTF = TESTF * (8*X*X - 8*X + 1)
350 CONTINUE
RETURN
C
1300 EXACTF = 0.0
RETURN
C
C Test integral 4
C
400 TESTF = 0
X = 1
DO 410 I = 1, DIMEN
X = - X * QUASI(I)
TESTF = TESTF + X
410 CONTINUE
RETURN
C
C
1400 X = 1.0 / (2 ** (DIMEN ))
IF (MOD(DIMEN, 2) .EQ. 0) THEN
EXACTF = (X - 1) / 3
ELSE
EXACTF = (X + 1) / 3
ENDIF
RETURN
C
END