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.txt'.
C   Finite field arithmetic is carried out with the help of
C   precalculated addition and multiplication tables found on
C   the file 'gftabs.txt'.  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.txt' 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.txt' 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.txt', STATUS = 'old')
      OPEN (UNIT=2, FILE='irrtabs.txt', 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.txt from unit 1
C   and writes the file irrtabs.txt onto unit 2.
C   GFPLYS opens units 1 and 2.
C   SETFLD opens unit 1, reads gftabs.txt, 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 sieve.  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 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.
          end do
        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.txt' using unit 1.  These precalculated tables
C   are supposed to have been put there by GFARIT.
C
C      *****  For the base-2 programs, these precalculated
C      *****  tables are not needed and, therefore, neither
C      *****  is GFARIT.
C
C
C   Unit 1 is closed both before and after the call of SETFLD.
C
C USES
C   Integer function CHARAC gets the characteristic of a field.
C
C
C   ------------------------------------------------------------
C
C   The following COMMON block, used by many subroutines,
C   gives the order Q of a field, its characteristic P, and its
C   addition, multiplication and subtraction tables.
C   The parameter MAXQ gives the order of the largest field to
C   be handled.
C
      INTEGER MAXQ
      PARAMETER (MAXQ=50)
 
      INTEGER P, Q, ADD(0:MAXQ-1,0:MAXQ-1)
      INTEGER MUL(0:MAXQ-1, 0:MAXQ-1), SUB(0:MAXQ-1, 0:MAXQ-1)
      COMMON /FIELD/ P, Q, ADD, MUL, SUB
      SAVE /FIELD/
C
C   The following definitions concern the representation of
C   polynomials.
C
      INTEGER MAXDEG, DEG
      PARAMETER (MAXDEG=50, DEG=-1)
C
C   The parameter MAXDEG gives the highest degree of polynomial
C   to be handled.  Polynomials stored as arrays have the
C   coefficient of degree n in POLY(N), and the degree of the
C   polynomial in POLY(-1).  The parameter DEG is just to remind
C   us of this last fact.  A polynomial which is identically 0
C   is given degree -1.
C
C   A polynomial can also be stored in an integer I, with
C        I = AN*Q**N + ... + A0.
C   Routines ITOP and PTOI convert between these two formats.
C
C   -----------------------------------------------------------
C
C
C
      INTEGER I, J, N, CHARAC
C
      IF (QIN .LE. 1 .OR. QIN .GT. MAXQ) THEN
        WRITE (*,*) ' SETFLD :  Bad value of Q'
        STOP
      ENDIF
C
      Q = QIN
      P = CHARAC(Q)
C
      IF (P .EQ. 0) THEN
        WRITE (*,*) ' SETFLD :  There is no field of order', Q
        STOP
      ENDIF
C
C Set up to handle a field of prime order :  calculate ADD and MUL.
C
      IF (P .EQ. Q) THEN
        DO I = 0, Q-1
          DO J = 0, Q-1
            ADD(I,J) = MOD(I+J, P)
            MUL(I,J) = MOD(I*J, P)
          end do
        end do
C
C Set up to handle a field of prime-power order :  tables for
C ADD and MUL are in the file 'gftabs.txt'.
C
      ELSE
        OPEN (UNIT=1, FILE='gftabs.txt', STATUS='old')
C
C    *****  OPEN statements are system dependent.
C
   20   READ (1, 900, END=500) N
  900   FORMAT (20I3)
        DO I = 0, N-1
          READ (1, 900) (ADD(I,J), J = 0, N-1)
        end do
        DO I = 0, N-1
          READ (1, 900) (MUL(I,J), J = 0, N-1)
        end do
        IF (N .NE. Q) GOTO 20
        CLOSE (1)
      ENDIF
C
C Now use the addition table to set the subtraction table.
C
      DO I = 0, Q-1
        DO J = 0, Q-1
          SUB(ADD(I,J), I) = J
        end do
      end do

      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 J = -1, MAXDEG
        POLY(J) = 0
      end do
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 J = POLY(DEG), 0, -1
        I = Q * I + POLY(J)
      end do
      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 I = 0, DEGC
        TERM = 0
        DO J = MAX(0, I-DEGA), MIN(DEGB, I)
          TERM = ADD(TERM, MUL(PA(I-J), PB(J)))
        end do
        PT(I) = TERM
      end do
C
      PC(DEG) = DEGC
      DO I = 0, DEGC
        PC(I) = PT(I)
      end do
      DO I = DEGC+1, MAXDEG
        PC(I) = 0
      end do
      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