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.txt'. 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.txt [on unit 1].
C After gftabs.txt has been set up, run GFPLYS with its
C associated subroutine to set up the file irrtabs.txt
C [on unit 2]; this requires
C reading gftabs.txt from unit 1. The files gftabs.txt
C and irrtabs.txt 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.txt nor irrtabs.txt, 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.txt', 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.txt 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 ( QIN )
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.txt.
C SETFLD opens unit 1, reads gftabs.txt, 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.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
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 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