!----------------------------------------------------------------------------
! Version 0.1.0 (000516) Created. Less dynamic array-of-terms version.
! Version 0.1.1 (000825) Generalised Mult_PPP and Power_PPI.
!                        Dropped Mult_PP and Power_PI.
! Version 0.1.2 (001126) Argument unit on I/O-routines.
! PH Lundow (per-hakan.lundow@math.umu.se)
!
! Modules for basic operations on monomials and polynomials.
!
! USE Polynomials   ! The module Monomial_Definitions must be available
! TYPE(Polynomial) :: poly1, poly2, poly3, poly_array( 10 )
! TYPE(Monomial)   :: mono1, mono2, mono_array( 10 )
! INTEGER :: int, deg( variables )
!
! Monomial_Definitions:
!    Choose a module Monomial_Definitions suited to your application, see the
!    files monoreal.f90, monompi.f90, monompr.f90 or write your own module.
!
! Construction and destruction of polynomials:
!    CALL Init_Poly (poly1)        ! Constructs a polynomial.
!    CALL Clear_Poly (poly1)       ! Destroys a polynomial.
!    CALL Clear_Poly ()            ! Empty. For compatibility purposes only.
!    max_length_poly = 99999       ! Redefines the max length of a polynomial.
!    CALL Init_Poly (poly1)        ! poly1 will have max length 99999.
!    CALL Init_Poly (poly2, 99)    ! poly2 will have max length 99.
!    CALL Init_Poly (poly_array, Maxlen=99)
!    Important:
!       Before construction you may wish to set the maximum length. This can
!       be done either by setting the global variable max_length_poly to an
!       appropriate value or through the optional argument maxlen of Init_Poly.
!       If max_length_poly is redefined (1000 by default) then all future
!       polynomials will have max length max_length_poly but using the optional
!       argument maxlen applies only to that particular polynomial (or array
!       of polynomials). Both maxlen and max_length_poly must be at least 2.
!
! Assignments:
!    poly1 = poly2, poly1 = mono1, poly1 = int, mono1 = mono2, mono1 = int,
!    mono1 % coeff = one_coeff, mono1 % coeff = int, mono1 % degree = deg
!
! Arithmetic operators:
!    mono1 * mono2, mono1 / mono2, mono1 * int, int * mono1, mono1 / int,
!    int / mono1, mono1 ** int, -mono1
!
! Binary relations:
!    poly1 == poly2, poly1 /= poly2, mono1 == mono2, mono1 /= mono2
!
! Procedures:
!    Addition:
!    CALL Add_Poly (poly1, poly2, poly3)  ! poly1 <-- poly1 + poly2 * poly3
!    CALL Add_Poly (poly1, poly2, mono1)  ! poly1 <-- poly1 + poly2 * mono1
!    CALL Add_Poly (poly1, poly2, int)    ! poly1 <-- poly1 + poly2 * int
!    CALL Add_Poly (poly1, poly2, deg)    ! poly1 <-- poly1 + poly2 * x^deg
!    CALL Add_Poly (poly1, poly2)         ! poly1 <-- poly1 + poly2
!    CALL Add_Poly (poly1, mono1)         ! poly1 <-- poly1 + mono1
!    CALL Add_Poly (poly1, int)           ! poly1 <-- poly1 + int
!    CALL Add_Poly (poly1, deg)           ! poly1 <-- poly1 + x^deg
!
!    Multiplication:
!    CALL Mult_Poly (poly1, poly2, poly3) ! poly1 <-- poly2 * poly3
!    CALL Mult_Poly (poly1, mono1)        ! poly1 <-- poly1 * mono1
!    CALL Mult_Poly (poly1, int)          ! poly1 <-- poly1 * int
!    CALL Mult_Poly (poly1, deg)          ! poly1 <-- poly1 * x^deg
!
!    Division:
!    CALL Div_Poly (poly1, mono1)         ! poly1 <-- poly1 / mono1
!    CALL Div_Poly (poly1, int)           ! poly1 <-- poly1 / int
!    CALL Div_Poly (poly1, deg)           ! poly1 <-- poly1 / x^deg
!
!    Powers:
!    CALL Power_Poly (poly1, poly2, int)  ! poly1 <-- poly2 ** int
!
!    Sums:
!    CALL Sum_Poly (poly1, poly_array)    ! poly1 <-- SUM (poly_array)
!    CALL Sum_Poly (poly1, mono_array)    ! poly1 <-- SUM (mono_array)
!
! Handy procedures:
!    CALL Diff_Poly (poly1, int)          ! poly1 <-- d(poly1)/dx_int
!    Diff_Mono (mono1, int)               ! d(mono1)/dx_int
!    CALL Swap_Poly (poly1, poly2)        ! poly1 <--> poly2
!    Length_Poly (poly1)                  ! number of terms in poly1
!    ZeroQ (poly1), ZeroQ (mono1)         ! poly1 == 0 ?, mono1 == 0 ?
!
!    Order (poly1, poly2) ! 1 if poly1 precedes poly2 in the canonical order
!                         !-1 if poly2 precedes poly1 in the canonical order
!                         ! 0 if poly1 is equal to poly2
!
!    Order (mono1, mono2) ! Order is determined only by the monomials degrees
!                         ! 1 if mono1 precedes mono2 in the canonical order
!                         !-1 if mono2 precedes mono1 in the canonical order
!                         ! 0 if degrees are equal, coefficients may differ.
!
! Subroutines for I/O:
!    CALL Write_Poly (int, poly1, fmt, iostat)  ! poly1 --> unit int
!    CALL Write_Mono (int, mono1, fmt, iostat)  ! mono1 --> unit int
!    CALL Read_Poly  (int, poly1, iostat)       ! poly1 <-- unit int
!    CALL Read_Mono  (int, mono1, iostat)       ! mono1 <-- unit int
!    ! Arguments fmt (CHARACTER(*)) and iostat (INTEGER) are optional.
!    ! Write_Poly and Read_Poly recognises unformatted I/O.
!
! Enumeration:
!    Type and procedures:
!    TYPE(Enum_Poly), Enumeration, Next, LastQ
!    Example:
!    TYPE(Enum_Poly) :: enum
!    CALL Enumeration (enum, poly1)
!    mono2 = 0
!    DO WHILE ( .NOT. LastQ (enum) )
!       CALL Next (enum, mono1)
!       mono2 % coeff = mono2 % coeff + mono1 % coeff
!    END DO
!    ! mono2 % coeff is the sum of the coefficients.
!
! Array extensions:
!    TYPE(Polynomial) :: matrix1( 10, 10 ), matrix2( 10, 10 )
!    TYPE(Polynomial) :: vector1( 10 ), vector2( 10 )
!
!    CALL Init_Poly (matrix1)
!    CALL Init_Poly (vector1)
!    CALL Clear_Poly (matrix1)
!    CALL Clear_Poly (vector1)
!
!    matrix1 = poly1, matrix1 = mono1, matrix1 = int
!    vector1 = poly1, vector1 = mono1, vector1 = int
!
!    CALL Swap_Poly (matrix1, matrix2)
!    CALL Swap_Poly (vector1, vector2)
!
!    Length_Poly (matrix1)
!    Length_Poly (vector1)
!
!    CALL Read_Poly (5, matrix1)
!    CALL Read_Poly (5, vector1)
!    CALL Write_Poly (6, matrix1)
!    CALL Write_Poly (6, vector1)
!
! Public variables:
!    max_poly_length
!
! Public parameters:
!    zero_coeff, one_coeff, two_coeff, ten_coeff, zero_mono, one_mono,
!    two_mono, ten_mono, x_mono( 1:variables ), variables, signed_coeff
!
! These modules are written in Essential Fortran 90.
!-----------------------------------------------------------------------------
MODULE Monomials
USE Monomial_Definitions
IMPLICIT NONE

INTERFACE ASSIGNMENT(=)
   MODULE PROCEDURE Assign_MI
END INTERFACE

INTERFACE OPERATOR(-)
   MODULE PROCEDURE Minus_M
END INTERFACE

INTERFACE OPERATOR(*)
   MODULE PROCEDURE Times_MM
   MODULE PROCEDURE Times_MI
   MODULE PROCEDURE Times_IM
END INTERFACE

INTERFACE OPERATOR(/)
   MODULE PROCEDURE Divide_MM
   MODULE PROCEDURE Divide_MI
   MODULE PROCEDURE Divide_IM
END INTERFACE

INTERFACE OPERATOR(**)
   MODULE PROCEDURE Power_MI
END INTERFACE

INTERFACE OPERATOR(==)
   MODULE PROCEDURE EQ_MM
END INTERFACE

INTERFACE OPERATOR(/=)
   MODULE PROCEDURE NE_MM
END INTERFACE

INTERFACE ZeroQ
   MODULE PROCEDURE ZeroQ_M
END INTERFACE

INTERFACE Order
   MODULE PROCEDURE Order_MM
END INTERFACE

INTERFACE Diff_Mono
   MODULE PROCEDURE Diff_M
END INTERFACE

PUBLIC  :: &
   ASSIGNMENT(=), OPERATOR(-), OPERATOR(*), OPERATOR(/), OPERATOR(**),     &
   OPERATOR(==), OPERATOR(/=), ZeroQ, Order, Diff_Mono

PRIVATE :: &
   Assign_MI, Minus_M, Times_MM, Times_MI, Times_IM, Divide_MM, Divide_MI, &
   Divide_IM, Power_MI, EQ_MM, NE_MM, ZeroQ_M, Order_MM, Diff_M

CONTAINS

!-----------------------------------------------------------------------------

SUBROUTINE Assign_MI (u, n)
TYPE(Monomial), INTENT(OUT) :: u
INTEGER,        INTENT(IN ) :: n
u % coeff  = n
u % degree = 0
RETURN
END SUBROUTINE Assign_MI

!-----------------------------------------------------------------------------

FUNCTION Minus_M (u) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u
TYPE(Monomial) :: w  ! RESULT
w % coeff  = - u % coeff
w % degree =   u % degree
RETURN
END FUNCTION Minus_M

!-----------------------------------------------------------------------------

FUNCTION Times_MM (u, v) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u, v
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff  * v % coeff
w % degree = u % degree + v % degree
RETURN
END FUNCTION Times_MM

!-----------------------------------------------------------------------------

FUNCTION Times_MI (u, n) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u
INTEGER,        INTENT(IN) :: n
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff * n
w % degree = u % degree
RETURN
END FUNCTION Times_MI

!-----------------------------------------------------------------------------

FUNCTION Times_IM (n, u) RESULT (w)
INTEGER,        INTENT(IN) :: n
TYPE(Monomial), INTENT(IN) :: u
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff * n
w % degree = u % degree
RETURN
END FUNCTION Times_IM

!-----------------------------------------------------------------------------

FUNCTION Divide_MM (u, v) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u, v
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff  / v % coeff
w % degree = u % degree - v % degree
RETURN
END FUNCTION Divide_MM

!-----------------------------------------------------------------------------

FUNCTION Divide_MI (u, n) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u
INTEGER,        INTENT(IN) :: n
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff / n
w % degree = u % degree
RETURN
END FUNCTION Divide_MI

!-----------------------------------------------------------------------------

FUNCTION Divide_IM (n, u) RESULT (w)
INTEGER,        INTENT(IN) :: n
TYPE(Monomial), INTENT(IN) :: u
TYPE(Monomial) :: w  ! RESULT
w % coeff  = n / u % coeff
w % degree = - u % degree
RETURN
END FUNCTION Divide_IM

!-----------------------------------------------------------------------------

FUNCTION Power_MI (u, n) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u
INTEGER,        INTENT(IN) :: n
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff ** n
w % degree = u % degree * n
RETURN
END FUNCTION Power_MI

!-----------------------------------------------------------------------------

FUNCTION EQ_MM (u, v) RESULT (q)
TYPE(Monomial), INTENT(IN) :: u, v
LOGICAL :: q  ! RESULT
INTEGER :: i

i = 0
IF ( ZeroQ (u) ) i = i + 1
IF ( ZeroQ (v) ) i = i + 2

SELECT CASE (i)
CASE (0) ! u /= 0, v /= 0
   IF ( Order (u, v) == 0 ) THEN ! possibly equal
      IF ( u % coeff == v % coeff ) THEN
         q = .TRUE.
      ELSE
         q = .FALSE.
      END IF
   ELSE
      q = .FALSE.
   END IF
CASE (1) ! u == 0, v /= 0
   q = .FALSE.
CASE (2) ! u /= 0, v == 0
   q = .FALSE.
CASE (3) ! u == 0, v == 0
   q = .TRUE.
CASE DEFAULT
   STOP 'ERROR: CASE'
END SELECT

RETURN
END FUNCTION EQ_MM

!-----------------------------------------------------------------------------

FUNCTION NE_MM (u, v) RESULT (q)
TYPE(Monomial), INTENT(IN) :: u, v
LOGICAL :: q  ! RESULT
q = .NOT. ( u == v )
RETURN
END FUNCTION NE_MM

!-----------------------------------------------------------------------------

FUNCTION ZeroQ_M (u) RESULT (q)
! Tests whether u is a zero-monomial.
TYPE(Monomial), INTENT(IN) :: u
LOGICAL :: q  ! RESULT

INTEGER, PARAMETER :: min_degree = -HUGE (0)

IF ( u % coeff == zero_coeff ) THEN
   q = .TRUE.
ELSE IF ( ANY (u % degree <= min_degree) ) THEN
   q = .TRUE.
ELSE
   q = .FALSE.
END IF

RETURN
END FUNCTION ZeroQ_M

!-----------------------------------------------------------------------------

FUNCTION Order_MM (u, v) RESULT (o)
! If u precedes v in a polynomial then 1 else if v precedes u then -1
! else 0 end if.
TYPE(Monomial), INTENT(IN) :: u, v
INTEGER :: o  ! RESULT
INTEGER :: i

!DO i = 1, variables     ! last exponent change most rapidly
DO i = variables, 1, -1  ! first exponent change most rapidly
   IF ( u % degree( i ) /= v % degree( i ) ) THEN
      IF ( u % degree( i ) > v % degree( i ) ) THEN
         o = 1
      ELSE
         o = -1
      END IF
      RETURN
   END IF
END DO
o = 0

RETURN
END FUNCTION Order_MM

!-----------------------------------------------------------------------------

FUNCTION Diff_M (u, n) RESULT (w)
TYPE(Monomial), INTENT(IN) :: u
INTEGER,        INTENT(IN) :: n
TYPE(Monomial) :: w  ! RESULT
w % coeff  = u % coeff * u % degree( n )
w % degree = u % degree
w % degree( n ) = w % degree( n ) - 1
RETURN
END FUNCTION Diff_M

!-----------------------------------------------------------------------------

END MODULE Monomials

!=============================================================================

MODULE Polynomial_Basics
USE Monomials
IMPLICIT NONE

INTEGER, PUBLIC :: max_length_poly = 1000  ! Can be redefined by user.

TYPE, PRIVATE :: Term
   TYPE(Monomial)      :: mono
   TYPE(Term), POINTER :: next
END TYPE Term

TYPE, PUBLIC :: Polynomial
   PRIVATE
   INTEGER             :: length
   TYPE(Term), POINTER :: head
   TYPE(Term), POINTER :: stack
   TYPE(Term), POINTER :: array( : )
END TYPE Polynomial

TYPE, PUBLIC :: Enum_Poly
   PRIVATE
   TYPE(Term), POINTER :: trm
END TYPE Enum_Poly

INTERFACE ASSIGNMENT(=)
   MODULE PROCEDURE Assign_PP
   MODULE PROCEDURE Assign_PM
   MODULE PROCEDURE Assign_PI
   MODULE PROCEDURE Assign_M1P
END INTERFACE

INTERFACE OPERATOR(==)
   MODULE PROCEDURE EQ_PP
END INTERFACE

INTERFACE OPERATOR(/=)
   MODULE PROCEDURE NE_PP
END INTERFACE

INTERFACE Order
   MODULE PROCEDURE Order_PP
END INTERFACE

INTERFACE ZeroQ
   MODULE PROCEDURE ZeroQ_P
END INTERFACE

INTERFACE Init_Poly
   MODULE PROCEDURE Init_P
END INTERFACE

INTERFACE Clear_Poly
   MODULE PROCEDURE Clear_P
   MODULE PROCEDURE Clear_0
END INTERFACE

INTERFACE Length_Poly
   MODULE PROCEDURE Length_P
END INTERFACE

INTERFACE Swap_Poly
   MODULE PROCEDURE Swap_PP
END INTERFACE

INTERFACE Add_Poly
   MODULE PROCEDURE Add_PM
   MODULE PROCEDURE Add_PI
   MODULE PROCEDURE Add_PD
   MODULE PROCEDURE Add_PP
   MODULE PROCEDURE Add_PPI
   MODULE PROCEDURE Add_PPD
   MODULE PROCEDURE Add_PPM
   MODULE PROCEDURE Add_PPP
END INTERFACE

INTERFACE Mult_Poly
   MODULE PROCEDURE Mult_PM
   MODULE PROCEDURE Mult_PI
   MODULE PROCEDURE Mult_PD
   MODULE PROCEDURE Mult_PPP
END INTERFACE

INTERFACE Div_Poly
   MODULE PROCEDURE Div_PM
   MODULE PROCEDURE Div_PI
   MODULE PROCEDURE Div_PD
END INTERFACE

INTERFACE Power_Poly
   MODULE PROCEDURE Power_PPI
END INTERFACE

INTERFACE Sum_Poly
   MODULE PROCEDURE Sum_PP1
   MODULE PROCEDURE Sum_PM1
END INTERFACE

INTERFACE Diff_Poly
   MODULE PROCEDURE Diff_P
END INTERFACE

INTERFACE Write_Poly
   MODULE PROCEDURE Write_P
END INTERFACE

INTERFACE Read_Poly
   MODULE PROCEDURE Read_P
END INTERFACE

INTERFACE Enumeration
   MODULE PROCEDURE Enumeration_EP
END INTERFACE

INTERFACE Next
   MODULE PROCEDURE Next_EM
END INTERFACE

INTERFACE LastQ
   MODULE PROCEDURE LastQ_E
END INTERFACE

PUBLIC :: &
   ASSIGNMENT(=), OPERATOR(==), OPERATOR(/=), Order, ZeroQ, Init_Poly,       &
   Clear_Poly, Length_Poly, Swap_Poly, Add_Poly, Mult_Poly, Div_Poly,        &
   Power_Poly, Sum_Poly, Diff_Poly, Write_Poly, Read_Poly, Enumeration,      &
   Next, LastQ

PRIVATE :: &
   Assign_PP, Assign_PM, Assign_PI, Assign_M1P, EQ_PP, NE_PP, Order_PP,      &
   ZeroQ_P, Init_P, Clear_P, Clear_0, Length_P, Swap_PP, Add_PM, Add_PI,     &
   Add_PD,Add_PP, Add_PPM, Add_PPM_1, Add_PPM_2, Add_PPM_3, Add_PPI, Add_PPD,&
   Add_PPP, Mult_PM, Mult_PI, Mult_PD, Mult_PPP, Div_PM, Div_PI, Div_PD,     &
   Power_PPI, Power_1, Power_2, Sum_PP1, Sum_PM1, Diff_P, Write_P, Read_P,   &
   Sort_Terms, Add_Terms, Delete_Terms, Reverse_Terms,TermQ, Push_Stack,     &
   Pop_Stack, Enumeration_EP, Next_EM, LastQ_E

CONTAINS

!-----------------------------------------------------------------------------

SUBROUTINE Assign_PP (p, q)
! p <-- q
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
TYPE(Term), POINTER :: trm1, trm2

p % length = q % length
trm1 => p % head
trm2 => q % head % next

! Copy while both polynomials have terms

DO WHILE ( TermQ (trm1 % next) .AND. TermQ (trm2) )
   trm1 => trm1 % next
   trm1 % mono = trm2 % mono
   trm2 => trm2 % next
END DO

! trm1 points to the last term in p and/or trm2 points to the head of q.
! Copy to new terms in p while q still has terms.

DO WHILE ( TermQ (trm2) )
   CALL Pop_Stack (p, trm1 % next)
   trm1 => trm1 % next
   trm1 % mono = trm2 % mono
   trm1 % next => p % head
   trm2 => trm2 % next
END DO

! trm1 points to the last term in p.
! Delete remaining terms from p

CALL Delete_Terms (p, trm1)

RETURN
END SUBROUTINE Assign_PP

!-----------------------------------------------------------------------------

SUBROUTINE Assign_PM (p, u)
! p <-- u
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm

p % length = 0
CALL Delete_Terms (p, p % head)

IF ( .NOT. ZeroQ (u) ) THEN
   CALL Pop_Stack (p, trm)
   trm % mono = u
   p % head % next => trm
   trm % next => p % head
   p % length = 1
END IF

RETURN
END SUBROUTINE Assign_PM

!-----------------------------------------------------------------------------

SUBROUTINE Assign_PI (p, n)
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: n
TYPE(Monomial) :: u
u = n
p = u
RETURN
END SUBROUTINE Assign_PI

!-----------------------------------------------------------------------------

SUBROUTINE Assign_M1P (u, p)
! u( : ) <-- p
TYPE(Monomial),   INTENT(OUT) :: u( : )
TYPE(Polynomial), INTENT(IN ) :: p
TYPE(Term), POINTER :: trm
INTEGER :: i, j

IF ( p % length > SIZE (u) ) STOP 'ERROR: Assign_M1P'

trm => p % head % next
i = 0

DO WHILE ( TermQ (trm) )
   i = i + 1
   u( i ) = trm % mono
   trm => trm % next
END DO

DO j = i+1, SIZE (u)
   u( j ) = zero_mono
END DO

RETURN
END SUBROUTINE Assign_M1P

!-----------------------------------------------------------------------------

FUNCTION EQ_PP (p, q) RESULT (c)
TYPE(Polynomial), INTENT(IN) :: p, q
LOGICAL :: c  ! RESULT
c = Order (p, q) == 0
RETURN
END FUNCTION EQ_PP

!-----------------------------------------------------------------------------

FUNCTION NE_PP (p, q) RESULT (c)
TYPE(Polynomial), INTENT(IN) :: p, q
LOGICAL :: c  ! RESULT
c = Order (p, q) /= 0
RETURN
END FUNCTION NE_PP

!-----------------------------------------------------------------------------

FUNCTION Order_PP (p, q) RESULT (o)
! If p precedes q then o = 1 else if q precedes p then o = -1
! else o = 0 end if. The 'larger' polynomial comes first.
TYPE(Polynomial), INTENT(IN) :: p, q
INTEGER :: o  ! RESULT
TYPE(Term), POINTER :: trm1, trm2

! First we compare the lengths of the polynomials...

IF ( p % length /= q % length ) THEN
   o = SIGN (1, p % length - q % length)
   RETURN
END IF

! Thus p % length = q % length

trm1 => p % head % next
trm2 => q % head % next

! Compare the degrees instead...

DO WHILE ( TermQ (trm1) )
   o = Order (trm1 % mono, trm2 % mono)
   IF ( o /= 0 ) RETURN ! o has correct value
   trm1 => trm1 % next
   trm2 => trm2 % next
END DO

! The degrees are also equal

trm1 => p % head % next
trm2 => q % head % next

! Compare the coefficients...

DO WHILE ( TermQ (trm1) )
   IF ( trm1 % mono % coeff /= trm2 % mono % coeff ) THEN
      IF ( trm1 % mono % coeff > trm2 % mono % coeff ) THEN
         o = 1
      ELSE
         o = -1
      END IF
      RETURN
   END IF
   trm1 => trm1 % next
   trm2 => trm2 % next
END DO

! The polynomials are indeed equal!

o = 0

RETURN
END FUNCTION Order_PP

!-----------------------------------------------------------------------------

FUNCTION ZeroQ_P (p) RESULT (q)
TYPE(Polynomial), INTENT(IN) :: p
LOGICAL :: q  ! RESULT
q = ASSOCIATED (p % head % next, p % head)
RETURN
END FUNCTION ZeroQ_P

!-----------------------------------------------------------------------------

SUBROUTINE Init_P (p, maxlen)
! p <-- 0
TYPE(Polynomial),  INTENT(OUT) :: p
INTEGER, OPTIONAL, INTENT(IN ) :: maxlen
INTEGER :: i, k

p % length = 0
NULLIFY (p % head)
NULLIFY (p % stack)
NULLIFY (p % array)

IF ( PRESENT (maxlen) ) THEN
   IF ( maxlen < 2 ) STOP 'ERROR: Init_P 1'
   ALLOCATE (p % array( maxlen ), STAT=i)
ELSE
   IF ( max_length_poly < 2 ) STOP 'ERROR: Init_P 2'
   ALLOCATE (p % array( max_length_poly ), STAT=i)
END IF
IF ( i /= 0 ) STOP 'ERROR: STAT'

p % head => p % array( 1 )
p % head % next => p % head
p % head % mono = zero_mono
p % stack => p % array( 2 )

k = SIZE (p % array)

DO i = 2, k-1
   p % array( i ) % next => p % array( i+1 )
END DO

NULLIFY (p % array( k ) % next)

RETURN
END SUBROUTINE Init_P

!-----------------------------------------------------------------------------

SUBROUTINE Clear_P (p)
! p <-- NULL
TYPE(Polynomial), INTENT(IN OUT) :: p

p % length = -1
DEALLOCATE (p % array)
NULLIFY (p % head)
NULLIFY (p % stack)
NULLIFY (p % array)

RETURN
END SUBROUTINE Clear_P

!-----------------------------------------------------------------------------

SUBROUTINE Clear_0 ()
RETURN
END SUBROUTINE Clear_0

!-----------------------------------------------------------------------------

FUNCTION Length_P (p) RESULT (l)
TYPE(Polynomial), INTENT(IN) :: p
INTEGER :: l  ! RESULT
l = p % length
RETURN
END FUNCTION Length_P

!-----------------------------------------------------------------------------

SUBROUTINE Swap_PP (p, q)
TYPE(Polynomial), INTENT(IN OUT) :: p, q
TYPE(Term), POINTER :: trm, trm2( : )
INTEGER :: k

k = p % length
p % length = q % length
q % length = k

trm => p % head
p % head => q % head
q % head => trm

trm => p % stack
p % stack => q % stack
q % stack => trm

trm2 => p % array
p % array => q % array
q % array => trm2

RETURN
END SUBROUTINE Swap_PP

!-----------------------------------------------------------------------------

SUBROUTINE Add_PM (p, u)
! p <-- p + u
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm1, trm2
INTEGER :: o

IF ( ZeroQ (u) ) RETURN

trm1 => p % head
trm2 => p % head % next

DO
   o = Order (u, trm2 % mono)
   SELECT CASE (o)
   CASE (1)
      CALL Pop_Stack (p, trm1 % next)
      trm1 => trm1 % next
      trm1 % mono = u
      trm1 % next => trm2
      p % length = p % length + 1
      EXIT
   CASE (0)
      trm2 % mono % coeff = trm2 % mono % coeff + u % coeff
      IF ( trm2 % mono % coeff == zero_coeff ) THEN
         trm1 % next => trm2 % next
         CALL Push_Stack (p, trm2)
         p % length = p % length - 1
      END IF
      EXIT
   CASE (-1)
      trm1 => trm1 % next
      trm2 => trm2 % next
   CASE DEFAULT
      STOP 'ERROR: CASE'
   END SELECT
END DO

RETURN
END SUBROUTINE Add_PM

!-----------------------------------------------------------------------------

SUBROUTINE Add_PI (p, n)
! p <-- p + n
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,  INTENT(IN    ) :: n
TYPE(Monomial) :: u
u = n
CALL Add_Poly (p, u)
RETURN
END SUBROUTINE Add_PI

!-----------------------------------------------------------------------------

SUBROUTINE Add_PD (p, degree)
! p <-- p + x^degree
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: degree( : )
TYPE(Monomial) :: u
IF ( SIZE (degree) /= variables ) STOP 'ERROR: Add_PD'
u % coeff  = one_coeff
u % degree = degree
CALL Add_PM (p, u)
RETURN
END SUBROUTINE Add_PD

!-----------------------------------------------------------------------------

SUBROUTINE Add_PP (p, q)
! p <-- p + q
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
TYPE(Term), POINTER :: trm1, trm2, trm3
INTEGER :: o

IF ( ASSOCIATED (p % head, q % head) ) STOP 'ERROR: Add_PP'

trm1 => p % head
trm2 => p % head % next
trm3 => q % head % next

IF ( signed_coeff ) THEN ! signed coefficients

   DO WHILE ( TermQ (trm3) )
      o = Order (trm3 % mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono = trm3 % mono
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! trm3 has the same degree as trm2
         trm2 % mono % coeff = trm2 % mono % coeff + trm3 % mono % coeff
         IF ( trm2 % mono % coeff == zero_coeff ) THEN
            trm1 % next => trm2 % next
            CALL Push_Stack (p, trm2)
            trm2 => trm1 % next
            p % length = p % length - 1
         ELSE
            trm1 => trm1 % next
            trm2 => trm2 % next
         END IF
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

ELSE ! unsigned coefficients

   DO WHILE ( TermQ (trm3) )
      o = Order (trm3 % mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono = trm3 % mono
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! trm3 has the same degree as trm2
         trm2 % mono % coeff = trm2 % mono % coeff + trm3 % mono % coeff
         trm1 => trm1 % next
         trm2 => trm2 % next
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

END IF

RETURN
END SUBROUTINE Add_PP

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPM (p, q, u)
! p <-- p + q * u
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
TYPE(Monomial),   INTENT(IN    ) :: u
INTEGER :: i

IF ( ASSOCIATED (p % head, q % head) ) STOP 'ERROR: Add_PPM'

IF ( ZeroQ (u) ) RETURN

i = 0
IF ( u % coeff /= one_coeff ) i = i + 1
IF ( ANY (u % degree /= 0)  ) i = i + 2

SELECT CASE (i)
CASE (0)  ! coeff == 1, degree == 0
   CALL Add_PP (p, q)
CASE (1)  ! coeff /= 1, degree == 0
   CALL Add_PPM_1 (p, q, u)
CASE (2)  ! coeff == 1, degree /= 0
   CALL Add_PPM_2 (p, q, u)
CASE (3)  ! coeff /= 1, degree /= 0
   CALL Add_PPM_3 (p, q, u)
CASE DEFAULT
   STOP 'ERROR: CASE'
END SELECT

RETURN
END SUBROUTINE Add_PPM

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPM_1 (p, q, u)
! coeff /= 1, degree == 0
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm1, trm2, trm3
INTEGER :: o

trm1 => p % head
trm2 => p % head % next
trm3 => q % head % next

IF ( signed_coeff ) THEN ! signed coefficients

   DO WHILE ( TermQ (trm3) )
      o = Order (trm3 % mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! u * trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono % coeff  = trm3 % mono % coeff * u % coeff
         trm1 % mono % degree = trm3 % mono % degree
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! u * trm3 has the same degree as trm2
         trm2%mono%coeff = trm2%mono%coeff + trm3%mono%coeff * u%coeff
         IF ( trm2 % mono % coeff == zero_coeff ) THEN
            trm1 % next => trm2 % next
            CALL Push_Stack (p, trm2)
            trm2 => trm1 % next
            p % length = p % length - 1
         ELSE
            trm1 => trm1 % next
            trm2 => trm2 % next
         END IF
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes u * trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

ELSE ! unsigned coefficients

   DO WHILE ( TermQ (trm3) )
      o = Order (trm3 % mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! u * trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono % coeff  = trm3 % mono % coeff * u % coeff
         trm1 % mono % degree = trm3 % mono % degree
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! u * trm3 has the same degree as trm2
         trm2%mono%coeff = trm2%mono%coeff + trm3%mono%coeff * u%coeff
         trm1 => trm1 % next
         trm2 => trm2 % next
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes u * trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

END IF

RETURN
END SUBROUTINE Add_PPM_1

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPM_2 (p, q, u)
! coeff == 1, degree /= 0
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm1, trm2, trm3
TYPE(Monomial) :: tmp_mono
INTEGER :: o

trm1 => p % head
trm2 => p % head % next
trm3 => q % head % next

IF ( signed_coeff ) THEN ! signed coefficients

   DO WHILE ( TermQ (trm3) )
      tmp_mono % degree = trm3 % mono % degree + u % degree
      o = Order (tmp_mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! u * trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono % coeff  = trm3 % mono % coeff
         trm1 % mono % degree = tmp_mono % degree
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! u * trm3 has the same degree as trm2
         trm2 % mono % coeff = trm2 % mono % coeff + trm3 % mono % coeff
         IF ( trm2 % mono % coeff == zero_coeff ) THEN
            trm1 % next => trm2 % next
            CALL Push_Stack (p, trm2)
            trm2 => trm1 % next
            p % length = p % length - 1
         ELSE
            trm1 => trm1 % next
            trm2 => trm2 % next
         END IF
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes u * trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

ELSE ! unsigned coefficients

   DO WHILE ( TermQ (trm3) )
      tmp_mono % degree = trm3 % mono % degree + u % degree
      o = Order (tmp_mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! u * trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono % coeff  = trm3 % mono % coeff
         trm1 % mono % degree = tmp_mono % degree
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! u * trm3 has the same degree as trm2
         trm2 % mono % coeff = trm2 % mono % coeff + trm3 % mono % coeff
         trm1 => trm1 % next
         trm2 => trm2 % next
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes u * trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

END IF

RETURN
END SUBROUTINE Add_PPM_2

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPM_3 (p, q, u)
! coeff /= 1, degree /= 0
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm1, trm2, trm3
TYPE(Monomial) :: tmp_mono
INTEGER :: o

trm1 => p % head
trm2 => p % head % next
trm3 => q % head % next

IF ( signed_coeff ) THEN ! signed coefficients

   DO WHILE ( TermQ (trm3) )
      tmp_mono % degree = trm3 % mono % degree + u % degree
      o = Order (tmp_mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! u * trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono = trm3 % mono * u
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! u * trm3 has the same degree as trm2
         trm2%mono%coeff = trm2%mono%coeff + u%coeff * trm3%mono%coeff
         IF ( trm2 % mono % coeff == zero_coeff ) THEN
            trm1 % next => trm2 % next
            CALL Push_Stack (p, trm2)
            trm2 => trm1 % next
            p % length = p % length - 1
         ELSE
            trm1 => trm1 % next
            trm2 => trm2 % next
         END IF
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes u * trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

ELSE ! unsigned coefficients

   DO WHILE ( TermQ (trm3) )
      tmp_mono % degree = trm3 % mono % degree + u % degree
      o = Order (tmp_mono, trm2 % mono)
      SELECT CASE (o)
      CASE (1)  ! u * trm3 precedes trm2
         CALL Pop_Stack (p, trm1 % next)
         trm1 => trm1 % next
         trm1 % next => trm2
         trm1 % mono = trm3 % mono * u
         trm3 => trm3 % next
         p % length = p % length + 1
      CASE (0)  ! u * trm3 has the same degree as trm2
         trm2%mono%coeff = trm2%mono%coeff + u%coeff * trm3%mono%coeff
         trm1 => trm1 % next
         trm2 => trm2 % next
         trm3 => trm3 % next
      CASE (-1)  ! trm2 precedes u * trm3
         trm1 => trm1 % next
         trm2 => trm2 % next
      CASE DEFAULT
         STOP 'ERROR: CASE'
      END SELECT
   END DO

END IF

RETURN
END SUBROUTINE Add_PPM_3

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPI (p, q, n)
! p <-- p + q * n
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
INTEGER,  INTENT(IN    ) :: n
TYPE(Monomial) :: u
u = n
CALL Add_Poly (p, q, u)
RETURN
END SUBROUTINE Add_PPI

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPD (p, q, degree)
! p <-- p + q * x^degree
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
INTEGER,          INTENT(IN    ) :: degree( : )
TYPE(Monomial) :: u
IF ( SIZE (degree) /= variables ) STOP 'ERROR: Add_PPD'
u % coeff  = one_coeff
u % degree = degree
CALL Add_PPM (p, q, u)
RETURN
END SUBROUTINE Add_PPD

!-----------------------------------------------------------------------------

SUBROUTINE Add_PPP (p, q, r)
! p <-- p + q * r
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q, r
TYPE(Term), POINTER :: trm

IF ( ASSOCIATED (p % head, q % head) ) STOP 'ERROR: Add_PPP 1'
IF ( ASSOCIATED (p % head, r % head) ) STOP 'ERROR: Add_PPP 2'

IF ( q % length > r % length ) THEN

   CALL Reverse_Terms (r % head)
   trm => r % head % next
   DO WHILE ( TermQ (trm) )
      CALL Add_Poly (p, q, trm % mono)
      trm => trm % next
   END DO
   CALL Reverse_Terms (r % head)

ELSE IF ( q % length < r % length ) THEN

   CALL Reverse_Terms (q % head)
   trm => q % head % next
   DO WHILE ( TermQ (trm) )
      CALL Add_Poly (p, r, trm % mono)
      trm => trm % next
   END DO
   CALL Reverse_Terms (q % head)

ELSE IF ( ASSOCIATED (q % head, r % head) ) THEN ! q === r

   trm => q % head % next
   DO WHILE ( TermQ (trm) )
      CALL Add_Poly (p, r, trm % mono)
      trm => trm % next
   END DO

ELSE ! q =!= r

   CALL Reverse_Terms (r % head)
   trm => r % head % next
   DO WHILE ( TermQ (trm) )
      CALL Add_Poly (p, q, trm % mono)
      trm => trm % next
   END DO
   CALL Reverse_Terms (r % head)

END IF

RETURN
END SUBROUTINE Add_PPP

!-----------------------------------------------------------------------------

SUBROUTINE Mult_PM (p, u)
! p <-- p * u
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm
INTEGER :: i

IF ( ZeroQ (u) ) THEN
   p = zero_mono
   RETURN
END IF

i = 0
IF ( u % coeff /= one_coeff ) i = i + 1
IF ( ANY (u % degree /= 0)  ) i = i + 2
!! 0 <= i <= 3

trm => p % head % next

SELECT CASE (i)
CASE (0)  ! coeff = 1, degree = 0
   RETURN
CASE (1)  ! coeff /= 1, degree = 0
   DO WHILE ( TermQ (trm) )
      trm % mono % coeff  = trm % mono % coeff * u % coeff ! result /= 0
      trm => trm % next
   END DO
CASE (2)  ! coeff = 1, degree /= 0
   DO WHILE ( TermQ (trm) )
      trm % mono % degree = trm % mono % degree + u % degree
      trm => trm % next
   END DO
CASE (3) ! coeff /= 1, degree /= 0
   DO WHILE ( TermQ (trm) )
      trm % mono = trm % mono * u
      trm => trm % next
   END DO
CASE DEFAULT
   STOP 'ERROR: CASE'
END SELECT

RETURN
END SUBROUTINE Mult_PM

!-----------------------------------------------------------------------------

SUBROUTINE Mult_PI (p, n)
! p <-- p * n
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: n
TYPE(Monomial) :: u
u = n
CALL Mult_Poly (p, u)
RETURN
END SUBROUTINE Mult_PI

!-----------------------------------------------------------------------------

SUBROUTINE Mult_PD (p, degree)
! p <-- p * x^degree
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: degree( : )
TYPE(Monomial) :: u
IF ( SIZE (degree) /= variables ) STOP 'ERROR: Mult_PD'
u % coeff  = one_coeff
u % degree = degree
CALL Mult_PM (p, u)
RETURN
END SUBROUTINE Mult_PD

!-----------------------------------------------------------------------------

SUBROUTINE Mult_PPP (p, q, r)
! p <-- q * r  (not efficient)
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q, r
TYPE(Polynomial) :: s
CALL Init_Poly (s)
CALL Add_Poly (s, q, r)
p = s
CALL Clear_Poly (s)
RETURN
END SUBROUTINE Mult_PPP

!-----------------------------------------------------------------------------

SUBROUTINE Div_PM (p, u)
! p <-- p / u
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Monomial),   INTENT(IN    ) :: u
TYPE(Term), POINTER :: trm1, trm2
INTEGER :: i

IF ( ZeroQ (u) ) STOP 'ERROR: Div_PM'

i = 0
IF ( u % coeff /= one_coeff ) i = i + 1
IF ( ANY (u % degree /= 0)  ) i = i + 2
!! 0 <= i <= 3

trm1 => p % head
trm2 => p % head % next

SELECT CASE (i)
CASE (0)  ! coeff = 1, degree = 0

   RETURN

CASE (1)  ! coeff /= 1, degree = 0

   DO WHILE ( TermQ (trm2) )
      trm2 % mono % coeff  = trm2 % mono % coeff / u % coeff
      IF ( trm2 % mono % coeff == zero_coeff ) THEN
         trm1 % next => trm2 % next
         CALL Push_Stack (p, trm2)
         trm2 => trm1 % next
         p % length = p % length - 1
      ELSE
         trm1 => trm1 % next
         trm2 => trm2 % next
      END IF
   END DO

CASE (2)  ! coeff = 1, degree /= 0

   DO WHILE ( TermQ (trm2) )
      trm2 % mono % degree = trm2 % mono % degree - u % degree! result /= 0
      trm2 => trm2 % next
   END DO

CASE (3) ! coeff /= 1, degree /= 0

   DO WHILE ( TermQ (trm2) )
      trm2 % mono = trm2 % mono / u
      IF ( trm2 % mono % coeff == zero_coeff ) THEN
         trm1 % next => trm2 % next
         CALL Push_Stack (p, trm2)
         trm2 => trm1 % next
         p % length = p % length - 1
      ELSE
         trm1 => trm1 % next
         trm2 => trm2 % next
      END IF
   END DO

CASE DEFAULT

   STOP 'ERROR: CASE'

END SELECT

RETURN
END SUBROUTINE Div_PM

!-----------------------------------------------------------------------------

SUBROUTINE Div_PI (p, n)
! p <-- p / n
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: n
TYPE(Monomial) :: u
u = n
CALL Div_Poly (p, u)
RETURN
END SUBROUTINE Div_PI

!-----------------------------------------------------------------------------

SUBROUTINE Div_PD (p, degree)
! p <-- p / x^degree
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: degree( : )
TYPE(Monomial) :: u
IF ( SIZE (degree) /= variables ) STOP 'ERROR: Div_PD'
u % coeff  = one_coeff
u % degree = degree
CALL Div_PM (p, u)
RETURN
END SUBROUTINE Div_PD

!-----------------------------------------------------------------------------

SUBROUTINE Power_PPI (p, q, n)
! p <-- q ** n (n is an integer)
! p may be the same variable as q.
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
INTEGER,          INTENT(IN    ) :: n
TYPE(Polynomial) :: r

IF ( n == 1 ) THEN
   p = q
   RETURN
END IF

CALL Init_Poly (r)

SELECT CASE (q % length)
CASE (0)
   IF ( n == 0 ) STOP 'ERROR: Power_PPI 2'
   r = zero_mono
CASE (1)
   r = q % head % next % mono ** n
CASE (2)
   IF ( n < 0 ) STOP 'ERROR: Power_PPI 3'
   CALL Power_1 (r, q, n)
CASE DEFAULT ! q % length > 2
   IF ( n < 0 ) STOP 'ERROR: Power_PPI 4'
   CALL Power_2 (r, q, n)
END SELECT

p = r
CALL Clear_Poly (r)

RETURN
END SUBROUTINE Power_PPI

!-----------------------------------------------------------------------------

SUBROUTINE Power_1 (p, q, n)
! p <-- q ** n(n is an integer, q has length 2)
! WARNING: p = q not allowed!
! Computes the power of q using the binomial theorem
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
INTEGER,          INTENT(IN    ) :: n
TYPE(Term), POINTER :: trm
TYPE(Monomial) :: c, uu, vv, u, v
INTEGER :: k

p % length = n + 1
trm => q % head % next
v = trm % mono
trm => trm % next
u  = trm % mono
c  = one_mono
uu = one_mono
vv = v ** n
trm => p % head

DO k = 0, n
   IF ( TermQ (trm % next) ) THEN
      trm => trm % next
   ELSE
      CALL Pop_Stack (p, trm % next)
      trm => trm % next
      trm % next => p % head
   END IF
   trm % mono = c * uu * vv
   IF ( k == n ) EXIT
   c % coeff = (c % coeff * (n - k)) / (k + 1)
   uu = uu * u
   vv = vv / v
END DO

CALL Delete_Terms (p, trm)

RETURN
END SUBROUTINE Power_1

!-----------------------------------------------------------------------------

RECURSIVE SUBROUTINE Power_2 (p, q, n)
! p <-- q ** n(n is an integer)
! WARNING: p = q not allowed!
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q
INTEGER,          INTENT(IN    ) :: n
TYPE(Polynomial) :: r, s

p = zero_mono

SELECT CASE (n)
CASE (0)
   p = one_mono
CASE (1)
   p = q
CASE (2)
   CALL Add_Poly (p, q, q)
CASE (3)
   CALL Init_Poly (r)
   CALL Add_Poly (r, q, q)
   CALL Add_Poly (p, q, r)
   CALL Clear_Poly (r)
CASE (4)
   CALL Init_Poly (r)
   CALL Add_Poly (r, q, q)
   CALL Add_Poly (p, r, r)
   CALL Clear_Poly (r)
CASE DEFAULT  ! n > 4

   CALL Init_Poly (r)
   CALL Init_Poly (s)

   CALL Power_2 (r, q, n/2)
   IF ( MOD (n, 2) == 0 ) THEN
      CALL Add_Poly (p, r, r)
   ELSE
      CALL Add_Poly (s, r, r)
      CALL Add_Poly (p, q, s)
   END IF

   CALL Clear_Poly (r)
   CALL Clear_Poly (s)

END SELECT

RETURN
END SUBROUTINE Power_2

!-----------------------------------------------------------------------------

SUBROUTINE Sum_PP1 (p, q)
! p <-- SUM (q)
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Polynomial), INTENT(IN    ) :: q( : )
INTEGER :: i
p = zero_mono
DO i = 1, SIZE (q)
   CALL Add_Poly (p, q( i ))
END DO
RETURN
END SUBROUTINE Sum_PP1

!-----------------------------------------------------------------------------

SUBROUTINE Sum_PM1 (p, u)
! p <-- SUM (u)
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Monomial),   INTENT(IN    ) :: u( : )
TYPE(Monomial) :: old_mono
TYPE(Term), POINTER :: trm
LOGICAL :: sorted
INTEGER :: i

sorted = .TRUE.
old_mono % degree = HUGE (0)
p % length = 0
trm => p % head

DO i = 1, SIZE (u)
   IF ( ZeroQ (u( i )) ) CYCLE
   sorted = sorted .AND. Order (old_mono, u( i )) >= 0
   old_mono % degree = u( i ) % degree
   p % length = p % length + 1
   IF ( TermQ (trm % next) ) THEN
      trm => trm % next
   ELSE
      CALL Pop_Stack (p, trm % next)
      trm => trm % next
      trm % next => p % head
   END IF
   trm % mono = u( i )
END DO

! trm points to the last non-zero term. Delete the rest.

CALL Delete_Terms (p, trm)

! All non-zero terms are in p, but possibly in the wrong order.

IF ( .NOT. sorted ) THEN

   trm => p % head
   DO i = 1, p % length
      trm => trm % next
   END DO

   CALL Sort_Terms (p % head % next, trm)
   trm % next => p % head

END IF

! Now the terms are in correct order, but two consecutive terms
! may have the same degrees. If so we add them.

CALL Add_Terms (p)

RETURN
END SUBROUTINE Sum_PM1

!-----------------------------------------------------------------------------

SUBROUTINE Diff_P (p, x)
! p <-- dp/(dx)
TYPE(Polynomial), INTENT(IN OUT) :: p
INTEGER,          INTENT(IN    ) :: x
TYPE(Term), POINTER :: trm1, trm2

IF ( x < 1 .OR. x > variables ) STOP 'ERROR: Diff_P'

trm1 => p % head
trm2 => p % head % next

DO WHILE ( TermQ (trm2) )
   trm2 % mono = Diff_Mono (trm2 % mono, x)
   IF ( trm2 % mono % coeff == zero_coeff ) THEN
      trm1 % next => trm2 % next
      CALL Push_Stack (p, trm2)
      trm2 => trm1 % next
      p % length = p % length - 1
   ELSE
      trm1 => trm1 % next
      trm2 => trm2 % next
   END IF
END DO

RETURN
END SUBROUTINE Diff_P

!-----------------------------------------------------------------------------

SUBROUTINE Write_P (unit, p, fmt, iostat)
! Writes a polynomial p to unit u. Can control formatting of each
! monomial with optional argument fmt. Error-handling through
! optional argument iostat. Recognises unformatted output with
! INQUIRE-statement.
INTEGER,                INTENT(IN ) :: unit
TYPE(Polynomial),       INTENT(IN ) :: p
CHARACTER(*), OPTIONAL, INTENT(IN ) :: fmt
INTEGER,      OPTIONAL, INTENT(OUT) :: iostat
TYPE(Term), POINTER :: trm
CHARACTER(20) :: charvar
INTEGER :: err

err = 0
trm => p % head % next

INQUIRE (UNIT=unit, FORM=charvar)

IF ( TRIM (charvar) == 'UNDEFINED' ) STOP 'ERROR: Write_P 1'

IF ( TRIM (charvar) == 'UNFORMATTED' ) THEN

   IF ( PRESENT (fmt) ) STOP 'ERROR: Write_P 2'

   DO WHILE ( TermQ (trm) )
      WRITE (UNIT=unit, IOSTAT=err) trm % mono
      IF ( err /= 0 ) EXIT
      trm => trm % next
   END DO

   IF ( err == 0 ) THEN
      WRITE (UNIT=unit, IOSTAT=err) zero_mono
   END IF

ELSE IF ( PRESENT (fmt) ) THEN ! charvar == 'FORMATTED'

   DO WHILE ( TermQ (trm) )
      CALL Write_Mono (unit, trm % mono, fmt, iostat=err)
      IF ( err /= 0 ) EXIT
      trm => trm % next
   END DO

   IF ( err == 0 ) THEN
      WRITE (UNIT=unit, FMT='(A)', IOSTAT=err) &
            TRIM (REPEAT ('0 ', variables+1))
   END IF

ELSE ! charvar == 'FORMATTED'

   DO WHILE ( TermQ (trm) )
      CALL Write_Mono (unit, trm % mono, iostat=err)
      IF ( err /= 0 ) EXIT
      trm => trm % next
   END DO

   IF ( err == 0 ) THEN
      WRITE (UNIT=unit, FMT='(A)', IOSTAT=err) &
            TRIM (REPEAT ('0 ', variables+1))
   END IF

END IF

IF ( PRESENT (iostat) ) THEN
   iostat = err
ELSE IF ( err /= 0 ) THEN
   STOP 'ERROR: Write_P 3'
END IF

RETURN
END SUBROUTINE Write_P

!-----------------------------------------------------------------------------

SUBROUTINE Read_P (unit, p, iostat)
! Reads a polynomial p from unit. Error-handling with
! optional argument iostat. Recongnises unformatted input
! with INQUIRE-statement.
INTEGER,           INTENT(IN    ) :: unit
TYPE(Polynomial),  INTENT(IN OUT) :: p
INTEGER, OPTIONAL, INTENT(   OUT) :: iostat
TYPE(Term), POINTER :: trm
TYPE(Monomial) :: mono, old_mono
CHARACTER(20)  :: charvar
INTEGER :: err
LOGICAL :: sorted

sorted = .TRUE.
old_mono % degree = HUGE (0)
p % length = 0
trm => p % head

! Read monomials until we see a zero-coefficient or something goes wrong.

INQUIRE (UNIT=unit, FORM=charvar)

IF ( TRIM (charvar) == 'UNDEFINED') STOP 'ERROR: Read_P 1'

IF ( TRIM (charvar) == 'UNFORMATTED' ) THEN

   ! WARNING: It is assumed that the monomials are in correct order!
   ! The polynomial must have been written by Write_Poly!

   DO
      READ (UNIT=unit, IOSTAT=err) mono
      IF ( err /= 0 ) EXIT
      IF ( ZeroQ (mono) ) EXIT
      IF ( TermQ (trm % next) ) THEN
         trm => trm % next
      ELSE
         CALL Pop_Stack (p, trm % next)
         trm => trm % next
         trm % next => p % head
      END IF
      trm % mono = mono
      p % length = p % length + 1
   END DO

ELSE ! charvar == 'FORMATTED'

   DO
      CALL Read_Mono (unit, mono, iostat=err)
      IF ( err /= 0 ) EXIT
      IF ( ZeroQ (mono) ) EXIT
      sorted = sorted .AND. Order (old_mono, mono) >= 0
      old_mono % degree = mono % degree
      IF ( TermQ (trm % next) ) THEN
         trm => trm % next
      ELSE
         CALL Pop_Stack (p, trm % next)
         trm => trm % next
         trm % next => p % head
      END IF
      trm % mono = mono
      p % length = p % length + 1
   END DO

END IF

IF ( err == 0 ) THEN

   ! Everything is ok, we encountered a zero-coefficient, ie,
   ! an END-OF-POLYNOMIAL.

   CALL Delete_Terms (p, trm) ! Delete the terms beyond trm.

   IF ( .NOT. sorted ) THEN
      CALL Sort_Terms (p % head % next, trm)
      trm % next => p % head
   END IF

   CALL Add_Terms (p)

END IF

IF ( PRESENT (iostat) ) THEN
   iostat = err
ELSE IF ( err /= 0 ) THEN
   STOP 'ERROR: Read_P 2'
END IF

RETURN
END SUBROUTINE Read_P

!-----------------------------------------------------------------------------

RECURSIVE SUBROUTINE Sort_Terms (trma, trmz)
! Sorts the terms between trma and trmz.
TYPE(Term), POINTER :: trma, trmz  ! ARGUMENTS
TYPE(Term), POINTER :: trm1, trm2, trm3

IF ( ASSOCIATED (trma, trmz) ) RETURN ! Not much to sort

trm1 => trma
trm2 => trmz
trm3 => trma % next

DO WHILE ( .NOT. ASSOCIATED (trm3, trmz) )
   IF ( Order (trma % mono, trm3 % mono) >= 0 ) THEN
      ! trm3 comes after or equals trma
      trm1 % next => trm3
      trm1 => trm3
   ELSE
      ! trm3 comes strictly before trma
      trm2 % next => trm3
      trm2 => trm3
   END IF
   trm3 => trm3 % next
END DO

IF ( Order (trma % mono, trmz % mono) >= 0 ) THEN
   trm1 % next => trmz
   IF ( ASSOCIATED (trmz, trm2) ) THEN
      CALL Sort_Terms (trma % next, trmz)
   ELSE
      trm1 => trmz % next
      CALL Sort_Terms (trm1, trm2)
      CALL Sort_Terms (trma % next, trmz)
      trm2 % next => trma
      trma => trm1
   END IF
ELSE
   CALL Sort_Terms (trmz, trm2)
   IF ( .NOT. ASSOCIATED (trma, trm1) ) THEN
      CALL Sort_Terms (trma % next, trm1)
   END IF
   trm2 % next => trma
   trma => trmz
   trmz => trm1
END IF

RETURN
END SUBROUTINE Sort_Terms

!-----------------------------------------------------------------------------

SUBROUTINE Add_Terms (p)
! Adds the terms of p and deletes zero-terms.
! WARNING: The terms must be sorted, see Sort_Terms.
! WARNING: Terms must be non-zero.
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Term), POINTER :: trm1, trm2
INTEGER :: o

trm1 => p % head % next
trm2 => trm1 % next

DO WHILE ( TermQ (trm2) )
   o = Order (trm1 % mono, trm2 % mono)
   IF ( o == 0 ) THEN
      trm1 % mono % coeff = trm1 % mono % coeff + trm2 % mono % coeff
      trm1 % next => trm2 % next
      CALL Push_Stack (p, trm2)
      trm2 => trm1 % next
      p % length = p % length - 1
   ELSE
      trm1 => trm1 % next
      trm2 => trm2 % next
   END IF
END DO

! The remaining part checks for zero-terms and deletes them.

IF ( .NOT. signed_coeff ) RETURN

trm1 => p % head
trm2 => p % head % next

DO WHILE ( TermQ (trm2) )
   IF ( trm2 % mono % coeff == zero_coeff ) THEN
      trm1 % next => trm2 % next
      CALL Push_Stack (p, trm2)
      trm2 => trm1 % next
      p % length = p % length - 1
   ELSE
      trm1 => trm1 % next
      trm2 => trm2 % next
   END IF
END DO

RETURN
END SUBROUTINE Add_Terms

!-----------------------------------------------------------------------------

SUBROUTINE Delete_Terms (p, trm)
! Deletes all terms (not the head) *after* trm.
! Upon leave trm points to the term it pointed to on entry.
! OK even if trm points to the head.
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Term),       POINTER        :: trm
TYPE(Term), POINTER :: trm2

DO WHILE ( TermQ (trm % next) )
   trm2 => trm % next
   trm % next => trm2 % next
   CALL Push_Stack (p, trm2)
END DO

! trm => |term_0| => |term_1| => ... => |term_n| => |head|
! call delete_terms (trm)
! trm => |term_0| => |head|

RETURN
END SUBROUTINE Delete_Terms

!-----------------------------------------------------------------------------

SUBROUTINE Reverse_Terms (trm)
! Reverse the order of the terms beginning with the head trm.
TYPE(Term), POINTER :: trm
TYPE(Term), POINTER :: trm1, trm2, trm3

trm1 => trm
trm2 => trm1 % next
trm3 => trm2 % next

DO
   trm2 % next => trm1
   IF ( .NOT. TermQ (trm2) ) EXIT
   trm1 => trm2
   trm2 => trm3
   trm3 => trm3 % next
END DO

RETURN
END SUBROUTINE Reverse_Terms

!-----------------------------------------------------------------------------

FUNCTION TermQ (trm) RESULT (t)
! Tests whether trm points to a proper term, ie not the head, of a polynomial.
TYPE(Term), POINTER :: trm
LOGICAL :: t  ! RESULT
INTEGER, PARAMETER :: min_degree = -HUGE (0)
t = trm % mono % degree( 1 ) > min_degree
RETURN
END FUNCTION TermQ

!-----------------------------------------------------------------------------

SUBROUTINE Push_Stack (p, trm)
! Returns to the stack the term that trm points to.
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Term),       POINTER        :: trm
trm % next => p % stack
p % stack  => trm
RETURN
END SUBROUTINE Push_Stack

!-----------------------------------------------------------------------------

SUBROUTINE Pop_Stack (p, trm)
! On return trm points to a new term.
! If necessary a new chunk is allocated.
TYPE(Polynomial), INTENT(IN OUT) :: p
TYPE(Term),       POINTER :: trm  ! INTENT(OUT)

IF ( .NOT. ASSOCIATED (p % stack) ) STOP 'ERROR: Pop_Stack'

trm => p % stack
p % stack => p % stack % next

RETURN
END SUBROUTINE Pop_Stack

!-----------------------------------------------------------------------------
! Example of use of the Enum_Poly-type:
! Compute the sum of the coefficients in p.
! ---------
! TYPE(Polynomial) :: p
! TYPE(Enum_Poly)  :: e
! TYPE(Monomial)   :: u, v
! ---------
! u = zero_mono
! CALL Enumeration (e, p)
! DO WHILE ( .NOT. LastQ (e) )
!    CALL Next (e, v)
!    u % coeff = u % coeff + v % coeff
! END DO
! CALL Write_Mono (6, u)
! ---------
!-----------------------------------------------------------------------------

SUBROUTINE Enumeration_EP (enum, poly)
! Sets enum to the first term in the polynomial poly.
TYPE(Enum_Poly),  INTENT(OUT) :: enum
TYPE(Polynomial), INTENT(IN ) :: poly
enum % trm => poly % head % next
RETURN
END SUBROUTINE Enumeration_EP

!-----------------------------------------------------------------------------

SUBROUTINE Next_EM (enum, mono)
! Updates enum to the next term and mono
! is set to the current monomial.
TYPE(Enum_Poly), INTENT(IN OUT) :: enum
TYPE(Monomial),  INTENT(   OUT) :: mono
mono = enum % trm % mono
enum % trm => enum % trm % next
RETURN
END SUBROUTINE Next_EM

!-----------------------------------------------------------------------------

FUNCTION LastQ_E (enum) RESULT (q)
! Tests whether we have seen all terms, ie whether
! enum points to the head of a polynomial.
TYPE(Enum_Poly), INTENT(IN) :: enum
LOGICAL :: q  ! RESULT
q = .NOT. TermQ (enum % trm)
RETURN
END FUNCTION LastQ_E

!-----------------------------------------------------------------------------

END MODULE Polynomial_Basics

!=============================================================================

MODULE Polynomial_Extensions
! Contains array extensions of some procedures for up to two indices.
USE Polynomial_Basics
IMPLICIT NONE

INTERFACE ASSIGNMENT(=)
   MODULE PROCEDURE Assign_P1M
   MODULE PROCEDURE Assign_P2M
   MODULE PROCEDURE Assign_P1I
   MODULE PROCEDURE Assign_P2I
END INTERFACE

INTERFACE Init_Poly
   MODULE PROCEDURE Init_P1
   MODULE PROCEDURE Init_P2
END INTERFACE

INTERFACE Clear_Poly
   MODULE PROCEDURE Clear_P1
   MODULE PROCEDURE Clear_P2
END INTERFACE

INTERFACE Swap_Poly
   MODULE PROCEDURE Swap_P1P1
   MODULE PROCEDURE Swap_P2P2
END INTERFACE

INTERFACE Length_Poly
   MODULE PROCEDURE Length_P1
   MODULE PROCEDURE Length_P2
END INTERFACE

INTERFACE Write_Poly
   MODULE PROCEDURE Write_P1
   MODULE PROCEDURE Write_P2
END INTERFACE

INTERFACE Read_Poly
   MODULE PROCEDURE Read_P1
   MODULE PROCEDURE Read_P2
END INTERFACE

PUBLIC  :: &
   ASSIGNMENT(=), Init_Poly, Clear_Poly, Swap_Poly, Length_Poly, Write_Poly, &
   Read_Poly

PRIVATE :: &
   Assign_P1M, Assign_P2M, Assign_P1I, Assign_P2I, Init_P1, Init_P2,Clear_P1,&
   Clear_P2, Swap_P1P1, Swap_P2P2, Length_P1, Length_P2, Write_P1, Write_P2, &
   Read_P1, Read_P2

CONTAINS

!-----------------------------------------------------------------------------

SUBROUTINE Assign_P1M (p, u)
TYPE(Polynomial), INTENT(IN OUT) :: p( : )
TYPE(Monomial),   INTENT(IN    ) :: u
INTEGER :: i
DO i = 1, SIZE (p)
   p( i ) = u
END DO
RETURN
END SUBROUTINE Assign_P1M

!-----------------------------------------------------------------------------

SUBROUTINE Assign_P2M (p, u)
TYPE(Polynomial), INTENT(IN OUT) :: p( :, : )
TYPE(Monomial),   INTENT(IN    ) :: u
INTEGER :: i, j
DO j = 1, SIZE (p, 2)
   DO i = 1, SIZE (p, 1)
      p( i, j ) = u
   END DO
END DO
RETURN
END SUBROUTINE Assign_P2M

!-----------------------------------------------------------------------------

SUBROUTINE Assign_P1I (p, n)
TYPE(Polynomial), INTENT(IN OUT) :: p( : )
INTEGER,          INTENT(IN    ) :: n
INTEGER :: i
DO i = 1, SIZE (p)
   p( i ) = n
END DO
RETURN
END SUBROUTINE Assign_P1I

!-----------------------------------------------------------------------------

SUBROUTINE Assign_P2I (p, n)
TYPE(Polynomial), INTENT(IN OUT) :: p( :, : )
INTEGER,          INTENT(IN    ) :: n
INTEGER :: i, j
DO j = 1, SIZE (p, 2)
   DO i = 1, SIZE (p, 1)
      p( i, j ) = n
   END DO
END DO
RETURN
END SUBROUTINE Assign_P2I

!-----------------------------------------------------------------------------

SUBROUTINE Init_P1 (p, maxlen)
TYPE(Polynomial),  INTENT(IN OUT) :: p( : )
INTEGER, OPTIONAL, INTENT(IN    ) :: maxlen
INTEGER :: i
IF ( PRESENT (maxlen) ) THEN
   DO i = 1, SIZE (p)
      CALL Init_Poly (p( i ), maxlen)
   END DO
ELSE
   DO i = 1, SIZE (p)
      CALL Init_Poly (p( i ))
   END DO
END IF
RETURN
END SUBROUTINE Init_P1

!-----------------------------------------------------------------------------

SUBROUTINE Init_P2 (p, maxlen)
TYPE(Polynomial),  INTENT(IN OUT) :: p( :, : )
INTEGER, OPTIONAL, INTENT(IN    ) :: maxlen
INTEGER :: i, j
IF ( PRESENT (maxlen) ) THEN
   DO j = 1, SIZE (p, 2)
      DO i = 1, SIZE (p, 1)
         CALL Init_Poly (p( i, j ), maxlen)
      END DO
   END DO
ELSE
   DO j = 1, SIZE (p, 2)
      DO i = 1, SIZE (p, 1)
         CALL Init_Poly (p( i, j ))
      END DO
   END DO
END IF
RETURN
END SUBROUTINE Init_P2

!-----------------------------------------------------------------------------

SUBROUTINE Clear_P1 (p)
TYPE(Polynomial), INTENT(IN OUT) :: p( : )
INTEGER :: i
DO i = 1, SIZE (p)
   CALL Clear_Poly (p( i ))
END DO
RETURN
END SUBROUTINE Clear_P1

!-----------------------------------------------------------------------------

SUBROUTINE Clear_P2 (p)
TYPE(Polynomial), INTENT(IN OUT) :: p( :, : )
INTEGER :: i, j
DO j = 1, SIZE (p, 2)
   DO i = 1, SIZE (p, 1)
      CALL Clear_Poly (p( i, j ))
   END DO
END DO
RETURN
END SUBROUTINE Clear_P2

!-----------------------------------------------------------------------------

SUBROUTINE Swap_P1P1 (p, q)
TYPE(Polynomial), INTENT(IN OUT) :: p( : ), q( : )
INTEGER :: i
IF ( SIZE (p) /= SIZE (q) ) STOP 'ERROR: Swap_P1P1'
DO i = 1, SIZE (p)
   CALL Swap_Poly (p( i ), q( i ))
END DO
RETURN
END SUBROUTINE Swap_P1P1

!-----------------------------------------------------------------------------

SUBROUTINE Swap_P2P2 (p, q)
TYPE(Polynomial), INTENT(IN OUT) :: p( :, : ), q( :, : )
INTEGER :: i, j
IF ( SIZE (p, 1) /= SIZE (q, 1) ) STOP 'ERROR: Swap_P2P2 1'
IF ( SIZE (p, 2) /= SIZE (q, 2) ) STOP 'ERROR: Swap_P2P2 2'
DO j = 1, SIZE (p, 2)
   DO i = 1, SIZE (p, 1)
      CALL Swap_Poly (p( i, j ), q( i, j ))
   END DO
END DO
RETURN
END SUBROUTINE Swap_P2P2

!-----------------------------------------------------------------------------

FUNCTION Length_P1 (p) RESULT (len)
TYPE(Polynomial), INTENT(IN) :: p( : )
INTEGER :: len( SIZE(p) )  ! RESULT
INTEGER :: i
DO i = 1, SIZE (p)
   len( i ) = Length_Poly (p( i ))
END DO
RETURN
END FUNCTION Length_P1

!-----------------------------------------------------------------------------

FUNCTION Length_P2 (p) RESULT (len)
TYPE(Polynomial), INTENT(IN) :: p( :, : )
INTEGER :: len( SIZE(p,1), SIZE(p,2) )  ! RESULT
INTEGER :: i, j
DO j = 1, SIZE (p, 2)
   DO i = 1, SIZE (p, 1)
      len( i, j ) = Length_Poly (p( i, j ))
   END DO
END DO
RETURN
END FUNCTION Length_P2

!-----------------------------------------------------------------------------

SUBROUTINE Write_P1 (unit, p, fmt, iostat)
INTEGER,                INTENT(IN ) :: unit
TYPE(Polynomial),       INTENT(IN ) :: p( : )
CHARACTER(*), OPTIONAL, INTENT(IN ) :: fmt
INTEGER,      OPTIONAL, INTENT(OUT) :: iostat
INTEGER :: i, err

IF ( PRESENT (fmt) ) THEN
   DO i = 1, SIZE (p)
      CALL Write_Poly (unit, p( i ), fmt, iostat=err)
      IF ( err /= 0 ) EXIT
   END DO
ELSE
   DO i = 1, SIZE (p)
      CALL Write_Poly (unit, p( i ), iostat=err)
      IF ( err /= 0 ) EXIT
   END DO
END IF

IF ( PRESENT (iostat) ) THEN
   iostat = err
ELSE IF ( err /= 0 ) THEN
   STOP 'ERROR: Write_P1'
END IF

RETURN
END SUBROUTINE Write_P1

!-----------------------------------------------------------------------------

SUBROUTINE Write_P2 (unit, p, fmt, iostat)
INTEGER,                INTENT(IN ) :: unit
TYPE(Polynomial),       INTENT(IN ) :: p( :, : )
CHARACTER(*), OPTIONAL, INTENT(IN ) :: fmt
INTEGER,      OPTIONAL, INTENT(OUT) :: iostat
INTEGER :: i, j, err

IF ( PRESENT (fmt) ) THEN
   Loop1 : &
   DO j = 1, SIZE (p, 2)
      DO i = 1, SIZE (p, 1)
         CALL Write_Poly (unit, p( i, j ), fmt, iostat=err)
         IF ( err /= 0 ) EXIT Loop1
      END DO
   END DO Loop1
ELSE
   Loop2 : &
   DO j = 1, SIZE (p, 2)
      DO i = 1, SIZE (p, 1)
         CALL Write_Poly (unit, p( i, j ), iostat=err)
         IF ( err /= 0 ) EXIT Loop2
      END DO
   END DO Loop2
END IF

IF ( PRESENT (iostat) ) THEN
   iostat = err
ELSE IF ( err /= 0 ) THEN
   STOP 'ERROR: Write_P2'
END IF

RETURN
END SUBROUTINE Write_P2

!-----------------------------------------------------------------------------

SUBROUTINE Read_P1 (unit, p, iostat)
INTEGER,           INTENT(IN    ) :: unit
TYPE(Polynomial),  INTENT(IN OUT) :: p( : )
INTEGER, OPTIONAL, INTENT(   OUT) :: iostat
INTEGER :: i, err

DO i = 1, SIZE (p)
   CALL Read_Poly (unit, p( i ), iostat=err)
   IF ( err /= 0 ) EXIT
END DO

IF ( PRESENT (iostat) ) THEN
   iostat = err
ELSE IF ( err /= 0 ) THEN
   STOP 'ERROR: Read_P1'
END IF

RETURN
END SUBROUTINE Read_P1

!-----------------------------------------------------------------------------

SUBROUTINE Read_P2 (unit, p, iostat)
INTEGER,           INTENT(IN    ) :: unit
TYPE(Polynomial),  INTENT(IN OUT) :: p( :, : )
INTEGER, OPTIONAL, INTENT(   OUT) :: iostat
INTEGER :: i, j, err

Loop : &
DO j = 1, SIZE (p, 2)
   DO i = 1, SIZE (p, 1)
      CALL Read_Poly (unit, p( i, j ), iostat=err)
      IF ( err /= 0 ) EXIT Loop
   END DO
END DO Loop

IF ( PRESENT (iostat) ) THEN
   iostat = err
ELSE IF ( err /= 0 ) THEN
   STOP 'ERROR: Read_P2'
END IF

RETURN
END SUBROUTINE Read_P2

!-----------------------------------------------------------------------------

END MODULE Polynomial_Extensions

!=============================================================================

MODULE Polynomials
USE Polynomial_Extensions
IMPLICIT NONE
END MODULE Polynomials

!=============================================================================