C ALGORITHM 738, COLLECTED ALGORITHMS FROM ACM.
C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C VOL. 20, NO. 4, DECEMBER, 1994, P. 494-495.
GUIDE TO GENERAL-BASE PROGRAMS
------------------------------
These routines generate well-spaced sets of points
for use in multidimensional integration and in
global optimization.
The programs are written in Fortran 77.
There are three sets, with main routines
GFARIT
GFPLYS
GENIN
respectively. The subroutines associated with these main
routines are listed in the respective comments to those routines
and in the the program directory.
For the first time these routines are used --
1. Compile all main routines and subroutines.
2. Link each main routine with its respective subroutines.
3. Run the main routines GFARIT, GFPLYS, GENIN in that order.
4. Save the files gftabs.dat and irrtabs.dat produced
by GFARIT and GFPLYS.
*****
***** The procedure above works on the Vax.
***** For your computer, it may have to be modified slightly.
*****
***** The Unix version has common blocks FIELD.inc and
***** COMM.inc which are incorporated by an extended FORTRAN
***** compiler automatically. For the Vax versions, these
***** blocks have already been included in the respective
***** routines; the user need do nothing about them.
*****
Subsequently --
It is necessary only to run GENIN.
&&&&&& GENIN is a driver skeleton. It offers the
&&&&&& user a menu, asking him to choose the
&&&&&& number of an integrand [assuming that TESTF
&&&&&& has an integrand with that number], the
&&&&&& dimension, the base, and the skip value.
&&&&&& It tells the user what the "optimal" base is,
&&&&&& without requiring him to choose it. It also
&&&&&& suggests possible skip values. Various
&&&&&& combinations of integrals, dimensions, bases,
&&&&&& and skip values can be tried before choosing
&&&&&& to quit GENIN.
With the exception of note 3 below,
the general-base programs do not use extended FORTRAN functions
(in contrast to the programs tailored for base 2). However, if
you wish to use base 2, the base-2 programs will run much faster
than the general-base programs.
The general-base programs are portable with the following
possible exceptions:
1. OPEN statements are system dependent.
2. The parameter NBITS is these programs is
the number of bits per word, excluding sign.
It is currently set to 31. If not appropriate for
your computer, modify the PARAMETER statements giving
the value of NBITS throughout.
3. The Unix version contains system-dependent timing
statements, which can be deleted or modified. The
Vax version has no timing statements, but these can be
added, in a system-dependent way, if desired.
GUIDE TO BASE-2 PROGRAMS
------------------------
These programs generate well-spaced points for use
in multidimensional integration and in global optimization.
The programs are written in Fortran 77.
***** To run this set of programs on the VAX,
***** first -
***** compile GENIN2, INLO2, GOLO2, CALCC2, CALCV,
SETFLD, PLYMUL, CHARAC, TESTF
&&&&& Our VAX program does not measure elapsed
&&&&& time. It can be modified in a system
&&&&& - dependent way to do so.
&&&&& On the other hand, our Unix version does
&&&&& measure elapsed time, in a system-dependent
&&&&& way. See the comments on the Unix version
&&&&& of GENIN2.
&&&&& In our VAX programs, the files FIELD.INC
&&&&& and COMM2.INC have already been incorporated
&&&&& where needed. The user needs do to nothing
&&&&& about them.
&&&&& On the other hand, our Unix version processes
&&&&& these two files directly, using extended
&&&&& FORTRAN.
&&&&&& The user can replace our TESTF,
&&&&&& which contains four test integrals,
&&&&&& by his own TESTF with up to four
&&&&&& integrals of interest to him.
***** second -
***** link these compiled routines
***** third -
***** run GENIN2
&&&&&& GENIN2 is a driver skeleton. It gives
&&&&&& the user a menu, asking him to choose
&&&&&& the number of the integral to be estimated
&&&&&& [assuming that TESTF has an integral with
&&&&&& the corresponding number], the dimension,
&&&&&& and the skip value -- offering
&&&&&& a suggested value. Various combinations
&&&&&& of integrals, dimensions, and skip values
&&&&&& can be tried before choosing to quit GENIN2.
The procedure for other computers is similar.
***** For the VAX, we use the extended FORTRAN
***** IEOR [for bitwise exclusive-or]. Its name
***** and format may differ on other computers.
***** IEOR is used in INLO2 and GOLO2.
***** Under Unix, we use the extended Fortran function
***** XOR for IEOR.
***** OPEN statements may well have to be modified
***** to conform to local computer-center requirements.
***** For the base-2 programs, the only OPEN statement
***** is in GENIN2. Interactive input-output is normally
***** from and to the user's terminal. Batch output
***** goes to the file called OUTFIL.DAT.
===========================================================
There are several auxiliary subroutines.
Calling structure :
GENIN2 calls INLO2 and GOLO2.
INLO2 calls CALCC2.
CALCC2 calls CALCV and SETFLD.
CALCV calls PLYMUL.
SETFLD calls CHARAC.
GOLO2, PLYMUL, and CHARAC are self-contained.
The subroutines share a number of parameters and two common
blocks, which are described below.
The principal parameters :
MAXDIM - The maximum dimension that will be used. An
irreducible polynomial must be available for each
dimension : see subroutine CALCC2.
NBITS - The number of bits (not counting the sign) in
a fixed-point integer.
*****
***** Currently, NBITS is set to 31. If your computer
***** does not have 31 bits per word, excluding sign,
***** modify the value of NBITS in each corrsponding
***** PARAMETER statement.
*****
DEG - -1 : see the comments on polynomials below.
Common block /COMM2/
CJ - The packed values of Niedderreiter's C(I,J,R).
DIM - The dimension of the sequence to be generated.
COUNT - The index of the current item in the sequence.
NEXTQ - The numerators of the next item in the sequence.
These are like Niederreiter's XI(N) (page 54),
except that N is implicit, afd the NEXTQ are
integers. To obtain the values of XI(N),
multiply by RECIP (see GOLO2).
Array CJ of common block /COMM2/ is set up by subroutine
CALCC2. The other members of this common block are set
up by subroutine INLO2.
Common block /FIELD/
P - The characteristic of the field in use.
Q - The order of the field in use.
ADD, MUL, and SUB - Arithmetic tables for the field in use.
The values of the common block FIELD are set by SETFLD.
Polynomials :
A polynomial POLY (say) is held in an array POLY(-1:MAXDEG).
The coefficient of x**n is held in POLY(N), and the degree
of the polynomial is held in POLY(-1). The parameter
DEG = -1 reminds us of this when necessary ; that is,
the degree of POLY is given by POLY(DEG).
Further comments accompany each subroutine.
c*** genin.f
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.dat) and
C by GFPLYS (to set up the irreducible polynomials in
C irrtabs.dat).
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.dat and irrtab.dat '
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.dat and irrtab.dat 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.dat'
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
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
PARAMETER (MAXE=5, 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.dat', 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.dat' 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.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
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
c*** gfarit.f
PROGRAM GFARIT
C
C This version : 12 December 1991
C
C The program calculates addition and multiplication tables
C for arithmetic in finite fields, and writes them out to
C the file 'gftabs.dat'. Tables are only calculated for fields
C of prime-power order Q, the other cases being trivial.
C For each value of Q, the file contains first Q, then the
C addition table, and lastly the multiplication table.
C
C
C GFARIT and its associated subroutines below must be
C run to set up the file gftabs.dat [on unit 1].
C After gftabs.dat has been set up, run GFPLYS with its
C associated subroutine to set up the file irrtabs.dat
C [on unit 2]; this requires
C reading gftabs.dat from unit 1. The files gftabs.dat
C and irrtabs.dat can be saved for future use. Thus, each
C user [or set of users with access to these files] needs to
C run the respective sets of programs associated with GFARIT
C and GFPLYS just once. This must be done before running the
C set of programs associated with GENIN. The set of programs
C tailored for base 2, using the driver GENIN2, requires
C neither gftabs.dat nor irrtabs.dat, hence neither GFARIT
C nor GFPLYS.
C
C Below we list [the main program] GFARIT and its associated
C subroutines. An asterisk indicates a subroutine also used
C by GFPLYS and by GENIN. An ampersand denotes a subroutine
C also used by GFPLYS but not by GENIN.
C
C GFARIT
C GFTAB (writes unit 1)
C CHARAC *
C SETFLD *
C ITOP &
C PTOI &
C PLYADD
C PLYMUL *
C PLYDIV
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
C
OPEN (UNIT=1, FILE='gftabs.dat', STATUS = 'unknown')
C
C ****** OPEN statements are system dependent
C
WRITE (*,*) ' This is GFARIT'
DO 10 I = 2, MAXQ
CALL GFTAB(I)
WRITE (*,*) ' GFARIT : Case ', I, ' done'
10 CONTINUE
C
WRITE (*,*) ' GFARIT ends'
STOP
END
C
C ***** end of PROGRAM GFARIT
SUBROUTINE GFTAB (QIN)
C
C This version : 12 December 1991
C
C
C This routine writes gftabs.dat onto unit 1.
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 The common /FIELD/ will be set up to work modulo P, the
C characteristic of field QIN. We construct the tables for
C the field of order QIN in GFADD and GFMUL.
C
INTEGER QIN, I, J, PTOI, CHARAC
INTEGER PI(-1:MAXDEG), PJ(-1:MAXDEG), PK(-1:MAXDEG)
INTEGER GFADD(0:MAXQ-1, 0:MAXQ-1), GFMUL(0:MAXQ-1, 0:MAXQ-1)
C
C IRRPLY holds irreducible polynomials for constructing
C prime-power fields. IRRPLY(-2,I) says which field this
C row is used for, and then the rest of the row is a
C polynomial (with the degree in IRRPLY(-1,I) as usual).
C The chosen irreducible poly is copied into MODPLY for use.
C
INTEGER IRRPLY(8, -2:MAXDEG), MODPLY(-1:MAXDEG)
SAVE IRRPLY
C
DATA (IRRPLY(1,J), J=-2,2) /4, 2, 1, 1, 1/
DATA (IRRPLY(2,J), J=-2,3) /8, 3, 1, 1, 0, 1/
DATA (IRRPLY(3,J), J=-2,2) /9, 2, 1, 0, 1/
DATA (IRRPLY(4,J), J=-2,4) /16, 4, 1, 1, 0, 0, 1/
DATA (IRRPLY(5,J), J=-2,2) /25, 2, 2, 0, 1/
DATA (IRRPLY(6,J), J=-2,3) /27, 3, 1, 2, 0, 1/
DATA (IRRPLY(7,J), J=-2,5) /32, 5, 1, 0, 1, 0, 0, 1/
DATA (IRRPLY(8,J), J=-2,2) /49, 2, 1, 0, 1/
C
IF (QIN .LE. 1 .OR. QIN .GT. MAXQ) THEN
WRITE (*,*) ' ARITH : Bad value of Q'
RETURN
ENDIF
C
P = CHARAC(QIN)
C
C If QIN is not a prime-power, we are not interested.
C
IF (P .EQ. 0 .OR. P .EQ. QIN) RETURN
C
C Otherwise, we set up the elements of the common /FIELD/
C ready to do arithmetic mod P, the characteristic of QIN.
C
CALL SETFLD (P)
C
C Next find a suitable irreducible polynomial and copy it
C to array MODPLY.
C
I = 1
20 IF (IRRPLY(I,-2) .NE. QIN) THEN
I = I + 1
GOTO 20
ENDIF
DO 30 J = -1, IRRPLY(I, DEG)
30 MODPLY(J) = IRRPLY(I, J)
DO 40 J = IRRPLY(I,DEG)+1, MAXDEG
40 MODPLY(J) = 0
C
C Deal with the trivial cases ...
C
DO 60 I = 0, QIN-1
GFADD(I,0) = I
GFADD(0,I) = I
GFMUL(I,0) = 0
60 GFMUL(0,I) = 0
DO 70 I = 1, QIN-1
GFMUL(I,1) = I
70 GFMUL(1,I) = I
C
C ... and now deal with the rest. Each integer from 1 to QIN-1
C is treated as a polynomial with coefficients handled mod P.
C Multiplication of polynomials is mod MODPLY.
C
DO 80 I = 1, QIN-1
CALL ITOP (I, PI)
DO 80 J = 1, I
CALL ITOP (J, PJ)
CALL PLYADD (PI, PJ, PK)
GFADD(I,J) = PTOI (PK)
GFADD(J,I) = GFADD(I,J)
IF (I .GT. 1 .AND. J .GT. 1) THEN
CALL PLYMUL (PI, PJ, PK)
CALL PLYDIV (PK, MODPLY, PJ, PK)
GFMUL(I,J) = PTOI (PK)
GFMUL(J,I) = GFMUL(I,J)
ENDIF
80 CONTINUE
C
C Write the newly-calculated tables out to unit 1
C
WRITE (1, 900) QIN
C
C This is the table gftabs.dat.
C SETFLD opens unit 1, reads gftabs.dat, and then closes unit 1.
C
DO 100 I = 0, QIN-1
100 WRITE (1, 900) (GFADD(I,J), J = 0, QIN-1)
DO 110 I = 0, QIN-1
110 WRITE (1, 900) (GFMUL(I,J), J = 0, QIN-1)
900 FORMAT (20I3)
C
RETURN
END
C
C ***** end of SUBROUTINE GFTAB
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.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
SUBROUTINE ITOP (IN, POLY)
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 IN, I, J, POLY(-1:MAXDEG)
C
C Converts an integer to a polynomial with coefficients in the
C field of order Q.
C
DO 10 J = -1, MAXDEG
10 POLY(J) = 0
C
I = IN
J = -1
20 IF (I .GT. 0) THEN
J = J + 1
IF (J .GT. MAXDEG) THEN
WRITE (*,*) ' ITOP : Polynomial exceeds MAXDEG'
STOP
ENDIF
POLY(J) = MOD (I, Q)
I = I / Q
GOTO 20
ENDIF
POLY(DEG) = J
RETURN
END
C
C ***** end of SUBROUTINE ITOP
INTEGER FUNCTION PTOI (POLY)
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, POLY(-1:MAXDEG)
C
C Converts a polynomial with coefficients in the field of
C order Q to an integer.
C
I = 0
DO 10 J = POLY(DEG), 0, -1
10 I = Q * I + POLY(J)
PTOI = I
RETURN
END
C
C ***** end of INTEGER FUNCTION PTOI
SUBROUTINE PLYADD (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, MAXAB, DEGC
INTEGER PA(-1:MAXDEG), PB(-1:MAXDEG), PC(-1:MAXDEG)
C
C Adds polynomials PA and PB putting result in PC.
C Coefficients are elements of the field of order Q.
C
MAXAB = MAX (PA(DEG), PB(DEG))
DEGC = -1
DO 10 I = 0, MAXAB
PC(I) = ADD (PA(I), PB(I))
IF (PC(I) .NE. 0) DEGC = I
10 CONTINUE
PC(DEG) = DEGC
DO 20 I = MAXAB+1, MAXDEG
20 PC(I) = 0
RETURN
END
C
C ***** end of SUBROUTINE PLYADD
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 PLYDIV (PA, PB, PQ, PR)
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, D, M, BINV, DEGB, DEGR, DEGQ
INTEGER PA(-1:MAXDEG), PB(-1:MAXDEG), PQ(-1:MAXDEG)
INTEGER PR(-1:MAXDEG)
INTEGER PTQ(-1:MAXDEG), PTR(-1:MAXDEG)
C
C Divides polynomial PA by PB, putting the quotient in PQ
C and the remainder in PR.
C Coefficients are elements of the field of order Q.
C
IF (PB(DEG) .EQ. -1) THEN
WRITE (*,*) ' PLYDIV : Divide by zero'
STOP
ENDIF
C
DO 10 I = -1, MAXDEG
PTQ(I) = 0
10 PTR(I) = PA(I)
DEGR = PA(DEG)
DEGB = PB(DEG)
DEGQ = DEGR - DEGB
IF (DEGQ .LT. 0) DEGQ = -1
C
C Find the inverse of the leading coefficient of PB.
C
J = PB(DEGB)
DO 15 I = 1, P-1
IF (MUL(I,J) .EQ. 1) BINV = I
15 CONTINUE
C
DO 30 D = DEGQ, 0, -1
M = MUL (PTR(DEGR), BINV)
DO 20 I = DEGB, 0, -1
20 PTR(DEGR+I-DEGB) = SUB (PTR(DEGR+I-DEGB), MUL (M, PB(I)))
DEGR = DEGR - 1
30 PTQ(D) = M
C
DO 40 I = 0, MAXDEG
PQ(I) = PTQ(I)
40 PR(I) = PTR(I)
PQ(DEG) = DEGQ
50 IF (PR(DEGR) .EQ. 0 .AND. DEGR .GE. 0) THEN
DEGR = DEGR - 1
GOTO 50
ENDIF
PR(DEG) = DEGR
C
RETURN
END
C
C ***** end of SUBROUTINE PLYDIV
c*** gfplys.f
PROGRAM GFPLYS
C
C This version : 12 December 1991
C
C The program calculates irreducible polynomials for various
C finite fields, and writes them out to the file 'irrtabs.dat'.
C Finite field arithmetic is carried out with the help of
C precalculated addition and multiplication tables found on
C the file 'gftabs.dat'. The format of the irreducible polys on
C the output file is
C Q
C d1 a1 a2 ... a(d1)
C d2 b1 b2 ... b(d2)
C etc
C where Q is the order of the field, d1 is the degree of the
C first irreducible poly, a1, a2, ..., a(d1) are its
C coefficients, and so on.
C
C The file 'gftabs.dat' is read on unit 1. GFPLYS and the
C associated subroutines assume that GFARIT has been run previously
C to put the required data in this file. GFPLYS writes its
C output on file 'irrtabs.dat' using unit 2.
C
C GFPLYS and its associated subroutines are listed below.
C An asterisk indicates a subroutine also used by GFARIT
C and by GENIN. An ampersand indicates a subroutine also
C used by GFARIT but not by GENIN.
C
C GFPLYS
C IRRED (reads unit 1, writes unit 2)
C CHARAC *
C SETFLD *
C ITOP &
C PTOI &
C PLYMUL *
C FIND
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
C
OPEN (UNIT=1, FILE='gftabs.dat', STATUS = 'old')
OPEN (UNIT=2, FILE='irrtabs.dat', STATUS = 'unknown')
C
C ***** OPEN statements are system dependent
C
WRITE (*,*) ' This is GFPLYS'
DO 10 I = 2, MAXQ
CALL IRRED(I)
WRITE (*,*) ' GFPLYS : Case ', I, ' done'
10 CONTINUE
C
END
C
C ***** end of PROGRAM GFPLYS
SUBROUTINE IRRED (QIN)
C
C This version : 12 December 1991
C
C
C This routine reads the file gftabs.dat from unit 1
C and writes the file irrtabs.dat onto unit 2.
C GFPLYS opens units 1 and 2.
C SETFLD opens unit 1, reads gftabs.dat, and then closes unit 1
C on each call.
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 NPOLS, MAXS
PARAMETER (NPOLS=25, MAXS=400)
LOGICAL SEIVE(MAXS)
INTEGER MONPOL(MAXS)
C
C The parameter NPOLS says how many irreducible polys are to
C be calculated for each field.
C We find the irreducible polys using a seive. Parameter
C MAXS gives the size of this seive. Array MONPOL holds monic
C polys, array SEIVE says whether the poly is still OK.
C
INTEGER QIN, I, J, K, N, PTOI, FIND, CHARAC
INTEGER PI(-1:MAXDEG), PJ(-1:MAXDEG), PK(-1:MAXDEG)
C
IF (QIN .LE. 1 .OR. QIN .GT. MAXQ) THEN
WRITE (*,*) ' IRRED : Bad value of Q'
RETURN
ENDIF
C
P = CHARAC(QIN)
C
C If no field of order QIN exists, there is nothing to do.
C
IF (P .EQ. 0) RETURN
C
C Otherwise, set up the field arithmetic tables.
C
CALL SETFLD (QIN)
C
C Set up seive containing only monic polys
C
I = 0
J = 1
K = Q
DO 50 N = 1, MAXS
I = I + 1
IF (I .EQ. J) THEN
I = K
J = 2 * K
K = Q * K
ENDIF
MONPOL(N) = I
SEIVE(N) = .TRUE.
50 CONTINUE
C
C Write out irreducible polys as they are found
C
N = 0
WRITE (2, 900) QIN
900 FORMAT (20I3)
DO 200 I = 1, MAXS
IF (SEIVE(I)) THEN
CALL ITOP (MONPOL(I), PI)
K = PI(DEG)
WRITE (2, 900) K, (PI(J), J = 0, K)
N = N + 1
IF (N .EQ. NPOLS) RETURN
C
DO 100 J = I, MAXS
CALL ITOP (MONPOL(J), PJ)
CALL PLYMUL (PI, PJ, PK)
K = FIND (PTOI (PK), MONPOL, J, MAXS)
IF (K .NE. 0) SEIVE(K) = .FALSE.
100 CONTINUE
ENDIF
200 CONTINUE
C
WRITE (*,*) ' IRRED : Seive too small'
WRITE (*,*) ' Only', N, ' irreducible polys were found'
RETURN
C
END
C
C ***** end of SUBROUTINE IRRED
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.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
SUBROUTINE ITOP (IN, POLY)
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 IN, I, J, POLY(-1:MAXDEG)
C
C Converts an integer to a polynomial with coefficients in the
C field of order Q.
C
DO 10 J = -1, MAXDEG
10 POLY(J) = 0
C
I = IN
J = -1
20 IF (I .GT. 0) THEN
J = J + 1
IF (J .GT. MAXDEG) THEN
WRITE (*,*) ' ITOP : Polynomial exceeds MAXDEG'
STOP
ENDIF
POLY(J) = MOD (I, Q)
I = I / Q
GOTO 20
ENDIF
POLY(DEG) = J
RETURN
END
C
C ***** end of SUBROUTINE ITOP
INTEGER FUNCTION PTOI (POLY)
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, POLY(-1:MAXDEG)
C
C Converts a polynomial with coefficients in the field of
C order Q to an integer.
C
I = 0
DO 10 J = POLY(DEG), 0, -1
10 I = Q * I + POLY(J)
PTOI = I
RETURN
END
C
C ***** end of INTEGER FUNCTION PTOI
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
INTEGER FUNCTION FIND (N, TAB, I, MAXTAB)
C
C This version : 12 December 1991
C
INTEGER N, I, MAXTAB, TAB(MAXTAB), J
C
C Look up N in ordered TAB(I) to TAB(MAXTAB)
C
FIND = 0
IF (N .GT. TAB(MAXTAB)) RETURN
DO 10 J = I, MAXTAB
IF (TAB(J) .EQ. N) THEN
FIND = J
RETURN
ENDIF
10 CONTINUE
RETURN
END
C
C ***** end of INTEGER FUNCTION FIND
c*** genin2.f
PROGRAM GENIN2
C
C ***** This is the driver for the base-2 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 : 2 Feb 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 INLO2, GOLO2, and
C related programs. Various possible test integrals are
C provided by the function TESTF. GENIN2 generates only
C sequences with base 2.
C
C Interactive input and output go through the Fortran units
C denoted by * ; at most installations, * corresponds to
C the user's terminal. Fortran unit 2 is used to save the output.
C
C These programs assume that your computer's word length
C is 31 bits, excluding sign. If this is not the case,
C modify the parameter NBITS throughout accordingly.
C
C
C
C GENIN2 and its associated subroutines are listed below.
C
C GENIN2
C INLO2
C GOLO2
C CALCC2
C CALCV
C CHARAC
C SETFLD
C PLYMUL
C TESTF
C
C The suffix 2 indicates routines for use only by
C the set of programs tailored for base 2.
C
C The other routines are also used by the general-base programs.
C
INTEGER MAXDIM, OUTUNT, POWER
PARAMETER (MAXDIM=12, OUTUNT=2, POWER = 12)
C
C The parameter MAXDIM gives the maximum dimension that will
C be used. OUTUNT is the unit to save the output.
C POWER is used in a possible warm-up formula.
C
INTEGER I, NUM, DIMEN, SEQLEN, SKIP, STEP, PBASE
REAL QUASI(MAXDIM), TESTF, EXACTF
DOUBLE PRECISION SUM
C
WRITE (*,*) ' This is program GENIN2'
WRITE (*,*) ' Output file name is OUTFIL.DAT '
OPEN (unit = OUTUNT, file = 'OUTFIL.DAT', status = 'UNKNOWN')
C
C ***** OPEN statements may well have to be modified
C ***** to conform to local computer-center requirements.
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 GENIN2'
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
C ***** The sequence length is the number of
C ***** quasi-random points used to estimate the
C ***** integral, excluding warm-up.
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 :'
C
C ***** Except when comparing results with other
C ***** bases, we suggest taking SEQLEN to be a power
C ***** of 2. Examples:
C
WRITE (*,*) ' 2 ** 10 = ', 2 ** 10
WRITE (*,*) ' 2 ** 15 = ', 2 ** 15
WRITE (*,*) ' 2 ** 20 = ', 2 ** 20
WRITE (*,*) ' Enter SEQLEN (possibly a power of 2 above) '
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 (*,*) ' There is reason to believe that BASE * e, '
WRITE (*,*) ' where e is defined for example in '
WRITE (*,*) ' Bratley, Fox, and Niederreiter [1991], '
WRITE (*,*) ' would be a good choice. Our formula has '
WRITE (*,*) ' has the form SKIP = 2 ** POWER, where '
WRITE (*,*) ' POWER is chosen so that SKIP comes nowhere '
WRITE (*,*) ' near the largest possible machine-representable'
WRITE (*,*) ' integer. It does not hurt to choose '
WRITE (*,*) ' POWER larger than e, because skipping is '
WRITE (*,*) ' done implicitly in O(1) time. '
WRITE (*,*) ' The maximum value of e for a fixed dimension '
WRITE (*,*) ' s grows like log s. We allow some "fat" for '
WRITE (*,*) ' the implicit constant. '
WRITE (*,*) ' Numerically, 2 ** POWER = ', 2 ** POWER
WRITE (*,*) ' Enter SKIP (not necessarily the value above)'
READ (*,*) SKIP
IF (SKIP.LT.0) THEN
WRITE (*,*) ' Number must be nonnegative'
GOTO 20
ENDIF
C
C
CALL INLO2 (DIMEN, SKIP)
WRITE (*,*) ' GENIN2 : Initialization complete'
C
C Write titles and the exact value of the integral
C
WRITE (OUTUNT,27) NUM
27 FORMAT (/,' Test integral ',I2)
WRITE (OUTUNT,28) DIMEN, SEQLEN, SKIP
28 FORMAT (/,' Dimension ',I6,', Base 2 (GENIN2)',
1 /,' Sequence ',I7,', Skipped ',I4)
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
SUM = 0
PBASE = 2 ** 6
WRITE (*,*) ' Odd-looking iteration numbers are powers of 2 '
STEP = 500
DO 100 I = 1, SEQLEN
CALL GOLO2 (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 = 2 * PBASE
C
C There is reason to believe that the subsequence
C of estimates along powers of the base [here 2]
C converges faster than the original sequence or
C the subsequence corresponding to STEP.
C
ENDIF
C
99 FORMAT (I12,G24.7)
IF (I .EQ. 5000) STEP = 1000
IF (I .EQ. 10000) STEP = 5000
100 CONTINUE
C
WRITE (*,*) ' GENIN2 : iteration ends'
GOTO 5
C
END
C
C ***** end of PROGRAM GENIN2
SUBROUTINE INLO2 (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=12, 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
SUBROUTINE GOLO2 (QUASI)
C
C This version : 21 February 1992
C
C This routine is for base 2 only. The driver, GENIN2,
C calls it after proper set-up.
C
C See the general comments on implementing Niederreiter's
C low-discrepancy sequences.
C
C This subroutine generates a new quasi-random vector
C on each call.
C
C INPUT
C From INLO2, labelled common /COMM2/, properly initialized.
C
C OUTPUT
C To the caller, the next vector in the sequence in the
C array QUASI.
C
C USES
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=12, 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
REAL RECIP
PARAMETER (RECIP=2.0**(-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
REAL 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
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=12, 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=5, 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 12 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,1) / 1,0,1 /
DATA (IRRED(2,I), I=-1,1) / 1,1,1 /
DATA (IRRED(3,I), I=-1,2) / 2,1,1,1 /
DATA (IRRED(4,I), I=-1,3) / 3,1,1,0,1 /
DATA (IRRED(5,I), I=-1,3) / 3,1,0,1,1 /
DATA (IRRED(6,I), I=-1,4) / 4,1,1,0,0,1 /
DATA (IRRED(7,I), I=-1,4) / 4,1,0,0,1,1 /
DATA (IRRED(8,I), I=-1,4) / 4,1,1,1,1,1 /
DATA (IRRED(9,I), I=-1,5) / 5,1,0,1,0,0,1 /
DATA (IRRED(10,I), I=-1,5) / 5,1,0,0,1,0,1 /
DATA (IRRED(11,I), I=-1,5) / 5,1,1,1,1,0,1 /
DATA (IRRED(12,I), I=-1,5) / 5,1,1,1,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
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.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
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
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