Initial population.

1 files changed, 1036 insertions, 0 deletions

diff --git a/src/algebra/ffpoly.spad.pamphlet b/src/algebra/ffpoly.spad.pamphlet
new file mode 100644
index 00000000..eba95386
--- /dev/null
+++ b/src/algebra/ffpoly.spad.pamphlet
@@ -0,0 +1,1036 @@
+\documentclass{article}
+\usepackage{axiom}
+\begin{document}
+\title{\$SPAD/src/algebra ffpoly.spad}
+\author{Alexandre Bouyer, Johannes Grabmeier, Alfred Scheerhorn, Robert Sutor, Barry Trager}
+\maketitle
+\begin{abstract}
+\end{abstract}
+\eject
+\tableofcontents
+\eject
+\section{package FFPOLY FiniteFieldPolynomialPackage}
+<<package FFPOLY FiniteFieldPolynomialPackage>>=
+)abbrev package FFPOLY FiniteFieldPolynomialPackage
+++ Author: A. Bouyer, J. Grabmeier, A. Scheerhorn, R. Sutor, B. Trager
+++ Date Created: January 1991
+++ Date Last Updated: 1 June 1994
+++ Basic Operations:
+++ Related Constructors:
+++ Also See:
+++ AMS Classifications:
+++ Keywords: finite field, polynomial, irreducible polynomial, normal
+++   polynomial, primitive polynomial, random polynomials
+++ References:
+++   [LS] Lenstra, H. W. & Schoof, R. J., "Primitivive Normal Bases
+++        for Finite Fields", Math. Comp. 48, 1987, pp. 217-231
+++   [LN] Lidl, R. & Niederreiter, H., "Finite Fields",
+++        Encycl. of Math. 20, Addison-Wesley, 1983
+++  J. Grabmeier, A. Scheerhorn: Finite Fields in Axiom.
+++   Axiom Technical Report Series, to appear.
+++ Description:
+++   This package provides a number of functions for generating, counting
+++   and testing irreducible, normal, primitive, random polynomials
+++   over finite fields.
+
+FiniteFieldPolynomialPackage GF : Exports == Implementation where
+
+  GF : FiniteFieldCategory
+
+  I    ==> Integer
+  L    ==> List
+  NNI  ==> NonNegativeInteger
+  PI   ==> PositiveInteger
+  Rec  ==> Record(expnt:NNI, coeff:GF)
+  Repr ==> L Rec
+  SUP  ==> SparseUnivariatePolynomial GF
+
+  Exports ==> with
+ --    qEulerPhiCyclotomic : PI -> PI
+--      ++ qEulerPhiCyclotomic(n)$FFPOLY(GF) yields the q-Euler's function
+--      ++ of the n-th cyclotomic polynomial over the field {\em GF} of
+--      ++ order q (cf. [LN] p.122);
+--      ++ error if n is a multiple of the field characteristic.
+    primitive? : SUP -> Boolean
+      ++ primitive?(f) tests whether the polynomial f over a finite
+      ++ field is primitive, i.e. all its roots are primitive.
+    normal? : SUP -> Boolean
+      ++ normal?(f) tests whether the polynomial f over a finite field is
+      ++ normal, i.e. its roots are linearly independent over the field.
+    numberOfIrreduciblePoly : PI -> PI
+      ++ numberOfIrreduciblePoly(n)$FFPOLY(GF) yields the number of
+      ++ monic irreducible univariate polynomials of degree n
+      ++ over the finite field {\em GF}.
+    numberOfPrimitivePoly : PI -> PI
+      ++ numberOfPrimitivePoly(n)$FFPOLY(GF) yields the number of
+      ++ primitive polynomials of degree n over the finite field {\em GF}.
+    numberOfNormalPoly : PI -> PI
+      ++ numberOfNormalPoly(n)$FFPOLY(GF) yields the number of
+      ++ normal polynomials of degree n over the finite field {\em GF}.
+    createIrreduciblePoly : PI -> SUP
+      ++ createIrreduciblePoly(n)$FFPOLY(GF) generates a monic irreducible
+      ++ univariate polynomial of degree n over the finite field {\em GF}.
+    createPrimitivePoly : PI -> SUP
+      ++ createPrimitivePoly(n)$FFPOLY(GF) generates a primitive polynomial
+      ++ of degree n over the finite field {\em GF}.
+    createNormalPoly : PI -> SUP
+      ++ createNormalPoly(n)$FFPOLY(GF) generates a normal polynomial
+      ++ of degree n over the finite field {\em GF}.
+    createNormalPrimitivePoly : PI -> SUP
+      ++ createNormalPrimitivePoly(n)$FFPOLY(GF) generates a normal and
+      ++ primitive polynomial of degree n over the field {\em GF}.
+      ++ Note: this function is equivalent to createPrimitiveNormalPoly(n)
+    createPrimitiveNormalPoly : PI -> SUP
+      ++ createPrimitiveNormalPoly(n)$FFPOLY(GF) generates a normal and
+      ++ primitive polynomial of degree n over the field {\em GF}.
+      ++ polynomial of degree n over the field {\em GF}.
+    nextIrreduciblePoly : SUP -> Union(SUP, "failed")
+      ++ nextIrreduciblePoly(f) yields the next monic irreducible polynomial
+      ++ over a finite field {\em GF} of the same degree as f in the following
+      ++ order, or "failed" if there are no greater ones.
+      ++ Error: if f has degree 0.
+      ++ Note: the input polynomial f is made monic.
+      ++ Also, \spad{f < g} if
+      ++ the number of monomials of f is less
+      ++ than this number for g.
+      ++ If f and g have the same number of monomials,
+      ++ the lists of exponents are compared lexicographically.
+      ++ If these lists are also equal, the lists of coefficients
+      ++ are compared according to the lexicographic ordering induced by
+      ++ the ordering of the elements of {\em GF} given by {\em lookup}.
+    nextPrimitivePoly : SUP -> Union(SUP, "failed")
+      ++ nextPrimitivePoly(f) yields the next primitive polynomial over
+      ++ a finite field {\em GF} of the same degree as f in the following
+      ++ order, or "failed" if there are no greater ones.
+      ++ Error: if f has degree 0.
+      ++ Note: the input polynomial f is made monic.
+      ++ Also, \spad{f < g} if the {\em lookup} of the constant term
+      ++ of f is less than
+      ++ this number for g.
+      ++ If these values are equal, then \spad{f < g} if
+      ++ if the number of monomials of f is less than that for g or if
+      ++ the lists of exponents of f are lexicographically less than the
+      ++ corresponding list for g.
+      ++ If these lists are also equal, the lists of coefficients are
+      ++ compared according to the lexicographic ordering induced by
+      ++ the ordering of the elements of {\em GF} given by {\em lookup}.
+    nextNormalPoly : SUP -> Union(SUP, "failed")
+      ++ nextNormalPoly(f) yields the next normal polynomial over
+      ++ a finite field {\em GF} of the same degree as f in the following
+      ++ order, or "failed" if there are no greater ones.
+      ++ Error: if f has degree 0.
+      ++ Note: the input polynomial f is made monic.
+      ++ Also, \spad{f < g} if the {\em lookup} of the coefficient
+      ++ of the term of degree
+      ++ {\em n-1} of f is less than that for g.
+      ++ In case these numbers are equal, \spad{f < g} if
+      ++ if the number of monomials of f is less that for g or if
+      ++ the list of exponents of f are lexicographically less than the
+      ++ corresponding list for g.
+      ++ If these lists are also equal, the lists of coefficients are
+      ++ compared according to the lexicographic ordering induced by
+      ++ the ordering of the elements of {\em GF} given by {\em lookup}.
+    nextNormalPrimitivePoly : SUP -> Union(SUP, "failed")
+      ++ nextNormalPrimitivePoly(f) yields the next normal primitive polynomial
+      ++ over a finite field {\em GF} of the same degree as f in the following
+      ++ order, or "failed" if there are no greater ones.
+      ++ Error: if f has degree 0.
+      ++ Note: the input polynomial f is made monic.
+      ++ Also, \spad{f < g} if the {\em lookup} of the constant
+      ++ term of f is less than
+      ++ this number for g or if
+      ++ {\em lookup} of the coefficient of the term of degree {\em n-1}
+      ++ of f is less than this number for g.
+      ++ Otherwise, \spad{f < g}
+      ++ if the number of monomials of f is less than
+      ++ that for g or if the lists of exponents for f are
+      ++ lexicographically less than those for g.
+      ++ If these lists are also equal, the lists of coefficients are
+      ++ compared according to the lexicographic ordering induced by
+      ++ the ordering of the elements of {\em GF} given by {\em lookup}.
+      ++ This operation is equivalent to nextPrimitiveNormalPoly(f).
+    nextPrimitiveNormalPoly : SUP -> Union(SUP, "failed")
+      ++ nextPrimitiveNormalPoly(f) yields the next primitive normal polynomial
+      ++ over a finite field {\em GF} of the same degree as f in the following
+      ++ order, or "failed" if there are no greater ones.
+      ++ Error: if f has degree 0.
+      ++ Note: the input polynomial f is made monic.
+      ++ Also, \spad{f < g} if the {\em lookup} of the
+      ++ constant term of f is less than
+      ++ this number for g or, in case these numbers are equal, if the
+      ++ {\em lookup} of the coefficient of the term of degree {\em n-1}
+      ++ of f is less than this number for g.
+      ++ If these numbers are equals, \spad{f < g}
+      ++ if the number of monomials of f is less than
+      ++ that for g, or if the lists of exponents for f are lexicographically
+      ++ less than those for g.
+      ++ If these lists are also equal, the lists of coefficients are
+      ++ coefficients according to the lexicographic ordering induced by
+      ++ the ordering of the elements of {\em GF} given by {\em lookup}.
+      ++ This operation is equivalent to nextNormalPrimitivePoly(f).
+--    random : () -> SUP
+--      ++ random()$FFPOLY(GF) generates a random monic polynomial
+--      ++ of random degree over the field {\em GF}
+    random : PI -> SUP
+      ++ random(n)$FFPOLY(GF) generates a random monic polynomial
+      ++ of degree n over the finite field {\em GF}.
+    random : (PI, PI) -> SUP
+      ++ random(m,n)$FFPOLY(GF) generates a random monic polynomial
+      ++ of degree d over the finite field {\em GF}, d between m and n.
+    leastAffineMultiple: SUP  -> SUP
+      ++ leastAffineMultiple(f) computes the least affine polynomial which
+      ++ is divisible by the polynomial f over the finite field {\em GF},
+      ++ i.e. a polynomial whose exponents are 0 or a power of q, the
+      ++ size of {\em GF}.
+    reducedQPowers: SUP  -> PrimitiveArray SUP
+      ++ reducedQPowers(f)
+      ++ generates \spad{[x,x**q,x**(q**2),...,x**(q**(n-1))]}
+      ++ reduced modulo f where \spad{q = size()$GF} and \spad{n = degree f}.
+    --
+    -- we intend to implement also the functions
+    -- cyclotomicPoly: PI -> SUP, order: SUP -> PI,
+    -- and maybe a new version of irreducible?
+
+
+  Implementation ==> add
+
+    import IntegerNumberTheoryFunctions
+    import DistinctDegreeFactorize(GF, SUP)
+
+
+    MM := ModMonic(GF, SUP)
+
+    sizeGF : PI := size()$GF :: PI
+
+    revListToSUP(l:Repr):SUP ==
+        newl:Repr := empty()
+        -- cannot use map since copy for Record is an XLAM
+        for t in l repeat newl := cons(copy t, newl)
+        newl pretend SUP
+
+    listToSUP(l:Repr):SUP ==
+        newl:Repr := [copy t for t in l]
+        newl pretend SUP
+
+    nextSubset : (L NNI, NNI) -> Union(L NNI, "failed")
+      -- for a list s of length m with 1 <= s.1 < ... < s.m <= bound,
+      -- nextSubset(s, bound) yields the immediate successor of s
+      -- (resp. "failed" if s = [1,...,bound])
+      -- where s < t if and only if:
+      -- (i)  #s < #t; or
+      -- (ii) #s = #t and s < t in the lexicographical order;
+      -- (we have chosen to fix the signature with NNI instead of PI
+      --  to avoid coercions in the main functions)
+
+    reducedQPowers(f) ==
+      m:PI:=degree(f)$SUP pretend PI
+      m1:I:=m-1
+      setPoly(f)$MM
+      e:=reduce(monomial(1,1)$SUP)$MM ** sizeGF
+      w:=1$MM
+      qpow:PrimitiveArray SUP:=new(m,0)
+      qpow.0:=1$SUP
+      for i in 1..m1 repeat  qpow.i:=lift(w:=w*e)$MM
+      qexp:PrimitiveArray SUP:=new(m,0)
+      m = 1 =>
+        qexp.(0$I):= (-coefficient(f,0$NNI)$SUP)::SUP
+        qexp
+      qexp.0$I:=monomial(1,1)$SUP
+      h:=qpow.1
+      qexp.1:=h
+      for i in 2..m1 repeat
+        g:=0$SUP
+        while h ^= 0 repeat
+          g:=g + leadingCoefficient(h) * qpow.degree(h)
+          h:=reductum(h)
+        qexp.i:=(h:=g)
+      qexp
+
+    leastAffineMultiple(f) ==
+    -- [LS] p.112
+      qexp:=reducedQPowers(f)
+      n:=degree(f)$SUP
+      b:Matrix GF:= transpose matrix [entries vectorise
+           (qexp.i,n) for i in 0..n-1]
+      col1:Matrix GF:= new(n,1,0)
+      col1(1,1)  := 1
+      ns : List Vector GF := nullSpace (horizConcat(col1,b) )
+      ----------------------------------------------------------------
+      -- perhaps one should use that the first vector in ns is already
+      -- the right one
+      ----------------------------------------------------------------
+      dim:=n+2
+      coeffVector : Vector GF
+      until empty? ns repeat
+        newCoeffVector := ns.1
+        i : PI :=(n+1) pretend PI
+        while newCoeffVector(i) = 0 repeat
+          i := (i - 1) pretend PI
+        if i < dim then
+          dim := i
+          coeffVector := newCoeffVector
+        ns := rest ns
+      (coeffVector(1)::SUP) +(+/[monomial(coeffVector.k, _
+               sizeGF**((k-2)::NNI))$SUP for k in 2..dim])
+
+--    qEulerPhiCyclotomic n ==
+--      n = 1 => (sizeGF - 1) pretend PI
+--      p : PI := characteristic()$GF :: PI
+--      (n rem p) = 0 => error
+--        "cyclotomic polynomial not defined for this argument value"
+--      q  : PI := sizeGF
+--      -- determine the multiplicative order of q modulo n
+--      e  : PI := 1
+--      qe : PI := q
+--      while (qe rem n) ^= 1 repeat
+--        e  := e + 1
+--        qe := qe * q
+--      ((qe - 1) ** ((eulerPhi(n) quo e) pretend PI) ) pretend PI
+
+    numberOfIrreduciblePoly n ==
+      -- we compute the number Nq(n) of monic irreducible polynomials
+      -- of degree n over the field GF of order q by the formula
+      -- Nq(n) = (1/n)* sum(moebiusMu(n/d)*q**d) where the sum extends
+      -- over all divisors d of n (cf. [LN] p.93, Th. 3.25)
+      n = 1 => sizeGF
+      -- the contribution of d = 1 :
+      lastd : PI  := 1
+      qd    : PI  := sizeGF
+      sum   :  I  := moebiusMu(n) * qd
+      -- the divisors d > 1 of n :
+      divisorsOfn : L PI := rest(divisors n) pretend L PI
+      for d in divisorsOfn repeat
+        qd := qd * (sizeGF) ** ((d - lastd) pretend PI)
+        sum := sum + moebiusMu(n quo d) * qd
+        lastd := d
+      (sum quo n) :: PI
+
+    numberOfPrimitivePoly n == (eulerPhi((sizeGF ** n) - 1) quo n) :: PI
+      -- [each root of a primitive polynomial of degree n over a field
+      --  with q elements is a generator of the multiplicative group
+      --  of a field of order q**n (definition), and the number of such
+      --  generators is precisely eulerPhi(q**n - 1)]
+
+    numberOfNormalPoly n ==
+      -- we compute the number Nq(n) of normal polynomials of degree n
+      -- in GF[X], with GF of order q, by the formula
+      -- Nq(n) = (1/n) * qPhi(X**n - 1) (cf. [LN] p.124) where,
+      -- for any polynomial f in GF[X] of positive degree n,
+      -- qPhi(f) = q**n * (1 - q**(-n1)) *...* (1 - q**(-nr)) =
+      -- q**n * ((q**(n1)-1) / q**(n1)) *...* ((q**(nr)-1) / q**(n_r)),
+      -- the ni being the degrees of the distinct irreducible factors
+      -- of f in its canonical factorization over GF
+      -- ([LN] p.122, Lemma 3.69).
+      -- hence, if n = m * p**r where p is the characteristic of GF
+      -- and gcd(m,p) = 1, we get
+      -- Nq(n) = (1/n)* q**(n-m) * qPhi(X**m - 1)
+      -- now X**m - 1 is the product of the (pairwise relatively prime)
+      -- cyclotomic polynomials Qd(X) for which d divides m
+      -- ([LN] p.64, Th. 2.45), and each Qd(X) factors into
+      -- eulerPhi(d)/e (distinct) monic irreducible polynomials in GF[X]
+      -- of the same degree e, where e is the least positive integer k
+      -- such that d divides q**k - 1 ([LN] p.65, Th. 2.47)
+      n = 1 => (sizeGF - 1) :: NNI :: PI
+      m : PI := n
+      p : PI := characteristic()$GF :: PI
+      q : PI := sizeGF
+      while (m rem p) = 0 repeat   -- find m such that
+        m := (m quo p) :: PI       -- n = m * p**r and gcd(m,p) = 1
+      m = 1 =>
+         -- know that n is a power of p
+        (((q ** ((n-1)::NNI) )  * (q - 1) ) quo n) :: PI
+      prod : I := q - 1
+      divisorsOfm : L PI := rest(divisors m) pretend L PI
+      for d in divisorsOfm repeat
+        -- determine the multiplicative order of q modulo d
+        e  : PI := 1
+        qe : PI := q
+        while (qe rem d) ^= 1 repeat
+          e  := e + 1
+          qe := qe * q
+        prod := prod * _
+          ((qe - 1) ** ((eulerPhi(d) quo e) pretend PI) ) pretend PI
+      (q**((n-m) pretend PI) * prod quo n) pretend PI
+
+    primitive? f ==
+      -- let GF be a field of order q; a monic polynomial f in GF[X]
+      -- of degree n is primitive over GF if and only if its constant
+      -- term is non-zero, f divides X**(q**n - 1) - 1 and,
+      -- for each prime divisor d of q**n - 1,
+      -- f does not divide X**((q**n - 1) / d) - 1
+      -- (cf. [LN] p.89, Th. 3.16, and p.87, following Th. 3.11)
+      n : NNI := degree f
+      n = 0 => false
+      leadingCoefficient f ^= 1 => false
+      coefficient(f, 0) = 0 => false
+      q  : PI := sizeGF
+      qn1: PI := (q**n - 1) :: NNI :: PI
+      setPoly f
+      x := reduce(monomial(1,1)$SUP)$MM -- X rem f represented in MM
+      --
+      -- may be improved by tabulating the residues x**(i*q)
+      -- for i = 0,...,n-1 :
+      --
+      lift(x ** qn1)$MM ^= 1 => false -- X**(q**n - 1) rem f in GF[X]
+      lrec  : L Record(factor:I, exponent:I) := factors(factor qn1)
+      lfact : L PI := []              -- collect the prime factors
+      for rec in lrec repeat          -- of q**n - 1
+        lfact := cons((rec.factor) :: PI, lfact)
+      for d in lfact repeat
+        if (expt := (qn1 quo d)) >= n then
+          lift(x ** expt)$MM = 1 => return false
+      true
+
+    normal? f ==
+      -- let GF be a field with q elements; a monic irreducible
+      -- polynomial f in GF[X] of degree n is normal if its roots
+      -- x, x**q, ... , x**(q**(n-1)) are linearly independent over GF
+      n : NNI := degree f
+      n = 0 => false
+      leadingCoefficient f ^= 1 => false
+      coefficient(f, 0) = 0 => false
+      n = 1 => true
+      not irreducible? f => false
+      g:=reducedQPowers(f)
+      l:=[entries vectorise(g.i,n)$SUP for i in 0..(n-1)::NNI]
+      rank(matrix(l)$Matrix(GF)) = n => true
+      false
+
+    nextSubset(s, bound) ==
+      m : NNI := #(s)
+      m = 0 => [1]
+      -- find the first element s(i) of s such that s(i) + 1 < s(i+1) :
+      noGap : Boolean := true
+      i : NNI := 0
+      restOfs : L NNI
+      while noGap and not empty?(restOfs := rest s) repeat
+      -- after i steps (0 <= i <= m-1) we have s = [s(i), ... , s(m)]
+      -- and restOfs = [s(i+1), ... , s(m)]
+        secondOfs := first restOfs    -- s(i+1)
+        firstOfsPlus1 := first s + 1  -- s(i) + 1
+        secondOfs = firstOfsPlus1 =>
+          s := restOfs
+          i := i + 1
+        setfirst_!(s, firstOfsPlus1)  -- s := [s(i)+1, s(i+1),..., s(m)]
+        noGap := false
+      if noGap then                   -- here s = [s(m)]
+        firstOfs := first s
+        firstOfs < bound => setfirst_!(s, firstOfs + 1) -- s := [s(m)+1]
+        m < bound =>
+            setfirst_!(s, m + 1)      -- s := [m+1]
+            i := m
+        return "failed"               -- (here m = s(m) = bound)
+      for j in i..1 by -1 repeat  -- reconstruct the destroyed
+        s := cons(j, s)           -- initial part of s
+      s
+
+    nextIrreduciblePoly f ==
+      n : NNI := degree f
+      n = 0 => error "polynomial must have positive degree"
+      -- make f monic
+      if (lcf := leadingCoefficient f) ^= 1 then f := (inv lcf) * f
+      -- if f = fn*X**n + ... + f{i0}*X**{i0} with the fi non-zero
+      -- then fRepr := [[n,fn], ... , [i0,f{i0}]]
+      fRepr : Repr := f pretend Repr
+      fcopy : Repr := []
+      -- we can not simply write fcopy := copy fRepr because
+      -- the input(!) f would be modified by assigning
+      -- a new value to one of its records
+      for term in fRepr repeat
+        fcopy := cons(copy term, fcopy)
+      if term.expnt ^= 0 then
+        fcopy := cons([0,0]$Rec, fcopy)
+      tailpol : Repr := []
+      headpol : Repr := fcopy  -- [[0,f0], ... , [n,fn]] where
+                               -- fi is non-zero for i > 0
+      fcopy   := reverse fcopy
+      weight  : NNI := (#(fcopy) - 1) :: NNI -- #s(f) as explained above
+      taillookuplist : L NNI := []
+      -- the zeroes in the headlookuplist stand for the fi
+      -- whose lookup's were not yet computed :
+      headlookuplist : L NNI := new(weight, 0)
+      s  : L NNI := [] -- we will compute s(f) only if necessary
+      n1 : NNI := (n - 1) :: NNI
+      repeat
+        -- (run through the possible weights)
+        while not empty? headlookuplist repeat
+          -- find next polynomial in the above order with fixed weight;
+          -- assume at this point we have
+          -- headpol = [[i1,f{i1}], [i2,f{i2}], ... , [n,1]]
+          -- and tailpol = [[k,fk], ... , [0,f0]] (with k < i1)
+          term := first headpol
+          j := first headlookuplist
+          if j = 0 then j := lookup(term.coeff)$GF
+          j := j + 1 -- lookup(f{i1})$GF + 1
+          j rem sizeGF = 0 =>
+            -- in this case one has to increase f{i2}
+            tailpol := cons(term, tailpol) -- [[i1,f{i1}],...,[0,f0]]
+            headpol := rest headpol        -- [[i2,f{i2}],...,[n,1]]
+            taillookuplist := cons(j, taillookuplist)
+            headlookuplist := rest headlookuplist
+          -- otherwise set f{i1} := index(j)$GF
+          setelt(first headpol, coeff, index(j :: PI)$GF)
+          setfirst_!(headlookuplist, j)
+          if empty? taillookuplist then
+            pol := revListToSUP(headpol)
+            --
+            -- may be improved by excluding reciprocal polynomials
+            --
+            irreducible? pol => return pol
+          else
+            -- go back to fk
+            headpol := cons(first tailpol, headpol) -- [[k,fk],...,[n,1]]
+            tailpol := rest tailpol
+            headlookuplist := cons(first taillookuplist, headlookuplist)
+            taillookuplist := rest taillookuplist
+        -- must search for polynomial with greater weight
+        if empty? s then -- compute s(f)
+          restfcopy := rest fcopy
+          for entry in restfcopy repeat s := cons(entry.expnt, s)
+        weight = n => return "failed"
+        s1 := nextSubset(rest s, n1) :: L NNI
+        s := cons(0, s1)
+        weight := #s
+        taillookuplist := []
+        headlookuplist := cons(sizeGF, new((weight-1) :: NNI, 1))
+        tailpol := []
+        headpol := [] -- [[0,0], [s.2,1], ... , [s.weight,1], [n,1]] :
+        s1 := cons(n, reverse s1)
+        while not empty? s1 repeat
+          headpol := cons([first s1, 1]$Rec, headpol)
+          s1 := rest s1
+        headpol := cons([0, 0]$Rec, headpol)
+
+    nextPrimitivePoly f ==
+      n : NNI := degree f
+      n = 0 => error "polynomial must have positive degree"
+      -- make f monic
+      if (lcf := leadingCoefficient f) ^= 1 then f := (inv lcf) * f
+      -- if f = fn*X**n + ... + f{i0}*X**{i0} with the fi non-zero
+      -- then fRepr := [[n,fn], ... , [i0,f{i0}]]
+      fRepr : Repr := f pretend Repr
+      fcopy : Repr := []
+      -- we can not simply write fcopy := copy fRepr because
+      -- the input(!) f would be modified by assigning
+      -- a new value to one of its records
+      for term in fRepr repeat
+        fcopy := cons(copy term, fcopy)
+      if term.expnt ^= 0 then
+        term  := [0,0]$Rec
+        fcopy := cons(term, fcopy)
+      fcopy   := reverse fcopy
+      xn : Rec := first fcopy
+      c0 : GF  := term.coeff
+      l  : NNI := lookup(c0)$GF rem sizeGF
+      n = 1 =>
+        -- the polynomial X + c is primitive if and only if -c
+        -- is a primitive element of GF
+        q1 : NNI  := (sizeGF - 1) :: NNI
+        while l < q1 repeat -- find next c such that -c is primitive
+          l := l + 1
+          c := index(l :: PI)$GF
+          primitive?(-c)$GF =>
+            return [xn, [0,c]$Rec] pretend SUP
+        "failed"
+      weight : NNI := (#(fcopy) - 1) :: NNI -- #s(f)+1 as explained above
+      s  : L NNI := [] -- we will compute s(f) only if necessary
+      n1 : NNI := (n - 1) :: NNI
+      -- a necessary condition for a monic polynomial f of degree n
+      -- over GF to be primitive is that (-1)**n * f(0) be a
+      -- primitive element of GF (cf. [LN] p.90, Th. 3.18)
+      c  : GF  := c0
+      while l < sizeGF repeat
+        -- (run through the possible values of the constant term)
+        noGenerator : Boolean := true
+        while noGenerator and l < sizeGF repeat
+          -- find least c >= c0 such that (-1)^n c0 is primitive
+          primitive?((-1)**n * c)$GF => noGenerator := false
+          l := l + 1
+          c := index(l :: PI)$GF
+        noGenerator => return "failed"
+        constterm : Rec := [0, c]$Rec
+        if c = c0 and weight > 1 then
+          headpol : Repr := rest reverse fcopy -- [[i0,f{i0}],...,[n,1]]
+                                               -- fi is non-zero for i>0
+          -- the zeroes in the headlookuplist stand for the fi
+          -- whose lookup's were not yet computed :
+          headlookuplist : L NNI := new(weight, 0)
+        else
+          -- X**n + c can not be primitive for n > 1 (cf. [LN] p.90,
+          -- Th. 3.18); next possible polynomial is X**n + X + c
+          headpol : Repr := [[1,0]$Rec, xn] -- 0*X + X**n
+          headlookuplist : L NNI := [sizeGF]
+          s := [0,1]
+          weight := 2
+        tailpol : Repr := []
+        taillookuplist : L NNI := []
+        notReady : Boolean := true
+        while notReady repeat
+          -- (run through the possible weights)
+          while not empty? headlookuplist repeat
+            -- find next polynomial in the above order with fixed
+            -- constant term and weight; assume at this point we have
+            -- headpol = [[i1,f{i1}], [i2,f{i2}], ... , [n,1]] and
+            -- tailpol = [[k,fk],...,[k0,fk0]] (k0<...<k<i1<i2<...<n)
+            term := first headpol
+            j := first headlookuplist
+            if j = 0 then j := lookup(term.coeff)$GF
+            j := j + 1 -- lookup(f{i1})$GF + 1
+            j rem sizeGF = 0 =>
+              -- in this case one has to increase f{i2}
+              tailpol := cons(term, tailpol) -- [[i1,f{i1}],...,[k0,f{k0}]]
+              headpol := rest headpol        -- [[i2,f{i2}],...,[n,1]]
+              taillookuplist := cons(j, taillookuplist)
+              headlookuplist := rest headlookuplist
+            -- otherwise set f{i1} := index(j)$GF
+            setelt(first headpol, coeff, index(j :: PI)$GF)
+            setfirst_!(headlookuplist, j)
+            if empty? taillookuplist then
+              pol := revListToSUP cons(constterm, headpol)
+              --
+              -- may be improved by excluding reciprocal polynomials
+              --
+              primitive? pol => return pol
+            else
+              -- go back to fk
+              headpol := cons(first tailpol, headpol) -- [[k,fk],...,[n,1]]
+              tailpol := rest tailpol
+              headlookuplist := cons(first taillookuplist,
+                                              headlookuplist)
+              taillookuplist := rest taillookuplist
+          if weight = n then notReady := false
+          else
+            -- must search for polynomial with greater weight
+            if empty? s then -- compute s(f)
+              restfcopy := rest fcopy
+              for entry in restfcopy repeat s := cons(entry.expnt, s)
+            s1 := nextSubset(rest s, n1) :: L NNI
+            s  := cons(0, s1)
+            weight := #s
+            taillookuplist := []
+            headlookuplist := cons(sizeGF, new((weight-2) :: NNI, 1))
+            tailpol := []
+            -- headpol = [[s.2,0], [s.3,1], ... , [s.weight,1], [n,1]] :
+            headpol := [[first s1, 0]$Rec]
+            while not empty? (s1 := rest s1) repeat
+              headpol := cons([first s1, 1]$Rec, headpol)
+            headpol := reverse cons([n, 1]$Rec, headpol)
+        -- next polynomial must have greater constant term
+        l := l + 1
+        c := index(l :: PI)$GF
+      "failed"
+
+    nextNormalPoly f ==
+      n : NNI := degree f
+      n = 0 => error "polynomial must have positive degree"
+      -- make f monic
+      if (lcf := leadingCoefficient f) ^= 1 then f := (inv lcf) * f
+      -- if f = fn*X**n + ... + f{i0}*X**{i0} with the fi non-zero
+      -- then fRepr := [[n,fn], ... , [i0,f{i0}]]
+      fRepr : Repr := f pretend Repr
+      fcopy : Repr := []
+      -- we can not simply write fcopy := copy fRepr because
+      -- the input(!) f would be modified by assigning
+      -- a new value to one of its records
+      for term in fRepr repeat
+        fcopy := cons(copy term, fcopy)
+      if term.expnt ^= 0 then
+        term  := [0,0]$Rec
+        fcopy := cons(term, fcopy)
+      fcopy     := reverse fcopy -- [[n,1], [r,fr], ... , [0,f0]]
+      xn : Rec  := first fcopy
+      middlepol : Repr := rest fcopy -- [[r,fr], ... , [0,f0]]
+      a0 : GF  := (first middlepol).coeff -- fr
+      l  : NNI := lookup(a0)$GF rem sizeGF
+      n = 1 =>
+        -- the polynomial X + a is normal if and only if a is not zero
+        l = sizeGF - 1 => "failed"
+        [xn, [0, index((l+1) :: PI)$GF]$Rec] pretend SUP
+      n1 : NNI := (n  - 1) :: NNI
+      n2 : NNI := (n1 - 1) :: NNI
+      -- if the polynomial X**n + a * X**(n-1) + ... is normal then
+      -- a = -(x + x**q +...+ x**(q**n)) can not be zero (where q = #GF)
+      a  : GF  := a0
+      -- if a = 0 then set a := 1
+      if l = 0 then
+        l := 1
+        a := 1$GF
+      while l < sizeGF repeat
+        -- (run through the possible values of a)
+        if a = a0 then
+          -- middlepol = [[0,f0], ... , [m,fm]] with m < n-1
+          middlepol := reverse rest middlepol
+          weight : NNI := #middlepol -- #s(f) as explained above
+          -- the zeroes in the middlelookuplist stand for the fi
+          -- whose lookup's were not yet computed :
+          middlelookuplist : L NNI := new(weight, 0)
+          s : L NNI := [] -- we will compute s(f) only if necessary
+        else
+          middlepol := [[0,0]$Rec]
+          middlelookuplist : L NNI := [sizeGF]
+          s : L NNI := [0]
+          weight : NNI := 1
+        headpol : Repr := [xn, [n1, a]$Rec] -- X**n + a * X**(n-1)
+        tailpol : Repr := []
+        taillookuplist : L NNI := []
+        notReady : Boolean := true
+        while notReady repeat
+          -- (run through the possible weights)
+          while not empty? middlelookuplist repeat
+            -- find next polynomial in the above order with fixed
+            -- a and weight; assume at this point we have
+            -- middlepol = [[i1,f{i1}], [i2,f{i2}], ... , [m,fm]] and
+            -- tailpol = [[k,fk],...,[0,f0]] ( with k<i1<i2<...<m)
+            term := first middlepol
+            j := first middlelookuplist
+            if j = 0 then j := lookup(term.coeff)$GF
+            j := j + 1 -- lookup(f{i1})$GF + 1
+            j rem sizeGF = 0 =>
+              -- in this case one has to increase f{i2}
+              -- tailpol = [[i1,f{i1}],...,[0,f0]]
+              tailpol   := cons(term, tailpol)
+              middlepol := rest middlepol -- [[i2,f{i2}],...,[m,fm]]
+              taillookuplist   := cons(j, taillookuplist)
+              middlelookuplist := rest middlelookuplist
+            -- otherwise set f{i1} := index(j)$GF
+            setelt(first middlepol, coeff, index(j :: PI)$GF)
+            setfirst_!(middlelookuplist, j)
+            if empty? taillookuplist then
+              pol := listToSUP append(headpol, reverse middlepol)
+              --
+              -- may be improved by excluding reciprocal polynomials
+              --
+              normal? pol => return pol
+            else
+              -- go back to fk
+              -- middlepol = [[k,fk],...,[m,fm]]
+              middlepol := cons(first tailpol, middlepol)
+              tailpol := rest tailpol
+              middlelookuplist := cons(first taillookuplist,
+                                               middlelookuplist)
+              taillookuplist := rest taillookuplist
+          if weight = n1 then notReady := false
+          else
+            -- must search for polynomial with greater weight
+            if empty? s then -- compute s(f)
+              restfcopy := rest rest fcopy
+              for entry in restfcopy repeat s := cons(entry.expnt, s)
+            s1 := nextSubset(rest s, n2) :: L NNI
+            s  := cons(0, s1)
+            weight := #s
+            taillookuplist := []
+            middlelookuplist := cons(sizeGF, new((weight-1) :: NNI, 1))
+            tailpol   := []
+            -- middlepol = [[0,0], [s.2,1], ... , [s.weight,1]] :
+            middlepol := []
+            s1 := reverse s1
+            while not empty? s1 repeat
+              middlepol := cons([first s1, 1]$Rec, middlepol)
+              s1 := rest s1
+            middlepol := cons([0,0]$Rec, middlepol)
+        -- next polynomial must have greater a
+        l := l + 1
+        a := index(l :: PI)$GF
+      "failed"
+
+    nextNormalPrimitivePoly f ==
+      n : NNI := degree f
+      n = 0 => error "polynomial must have positive degree"
+      -- make f monic
+      if (lcf := leadingCoefficient f) ^= 1 then f := (inv lcf) * f
+      -- if f = fn*X**n + ... + f{i0}*X**{i0} with the fi non-zero
+      -- then fRepr := [[n,fn], ... , [i0,f{i0}]]
+      fRepr : Repr := f pretend Repr
+      fcopy : Repr := []
+      -- we can not simply write fcopy := copy fRepr because
+      -- the input(!) f would be modified by assigning
+      -- a new value to one of its records
+      for term in fRepr repeat
+        fcopy := cons(copy term, fcopy)
+      if term.expnt ^= 0 then
+        term  := [0,0]$Rec
+        fcopy := cons(term, fcopy)
+      fcopy   := reverse fcopy -- [[n,1], [r,fr], ... , [0,f0]]
+      xn : Rec := first fcopy
+      c0 : GF  := term.coeff
+      lc : NNI := lookup(c0)$GF rem sizeGF
+      n = 1 =>
+        -- the polynomial X + c is primitive if and only if -c
+        -- is a primitive element of GF
+        q1 : NNI  := (sizeGF - 1) :: NNI
+        while lc < q1 repeat -- find next c such that -c is primitive
+          lc := lc + 1
+          c  := index(lc :: PI)$GF
+          primitive?(-c)$GF =>
+            return [xn, [0,c]$Rec] pretend SUP
+        "failed"
+      n1 : NNI := (n  - 1) :: NNI
+      n2 : NNI := (n1 - 1) :: NNI
+      middlepol : Repr := rest fcopy -- [[r,fr],...,[i0,f{i0}],[0,f0]]
+      a0 : GF  := (first middlepol).coeff
+      la : NNI := lookup(a0)$GF rem sizeGF
+      -- if the polynomial X**n + a * X**(n-1) +...+ c is primitive and
+      -- normal over GF then (-1)**n * c is a primitive element of GF
+      -- (cf. [LN] p.90, Th. 3.18), and a = -(x + x**q +...+ x**(q**n))
+      -- is not zero (where q = #GF)
+      c : GF  := c0
+      a : GF  := a0
+      -- if a = 0 then set a := 1
+      if la = 0 then
+        la := 1
+        a  := 1$GF
+      while lc < sizeGF repeat
+        -- (run through the possible values of the constant term)
+        noGenerator : Boolean := true
+        while noGenerator and lc < sizeGF repeat
+          -- find least c >= c0 such that (-1)**n * c0 is primitive
+          primitive?((-1)**n * c)$GF => noGenerator := false
+          lc := lc + 1
+          c  := index(lc :: PI)$GF
+        noGenerator => return "failed"
+        constterm : Rec := [0, c]$Rec
+        while la < sizeGF repeat
+        -- (run through the possible values of a)
+          headpol : Repr := [xn, [n1, a]$Rec] -- X**n + a X**(n-1)
+          if c = c0 and a = a0 then
+            -- middlepol = [[i0,f{i0}], ... , [m,fm]] with m < n-1
+            middlepol := rest reverse rest middlepol
+            weight : NNI := #middlepol + 1 -- #s(f)+1 as explained above
+            -- the zeroes in the middlelookuplist stand for the fi
+            -- whose lookup's were not yet computed :
+            middlelookuplist : L NNI := new((weight-1) :: NNI, 0)
+            s : L NNI := [] -- we will compute s(f) only if necessary
+          else
+            pol := listToSUP append(headpol, [constterm])
+            normal? pol and primitive? pol => return pol
+            middlepol := [[1,0]$Rec]
+            middlelookuplist : L NNI := [sizeGF]
+            s : L NNI := [0,1]
+            weight : NNI := 2
+          tailpol : Repr := []
+          taillookuplist : L NNI := []
+          notReady : Boolean := true
+          while notReady repeat
+          -- (run through the possible weights)
+            while not empty? middlelookuplist repeat
+              -- find next polynomial in the above order with fixed
+              -- c, a and weight; assume at this point we have
+              -- middlepol = [[i1,f{i1}], [i2,f{i2}], ... , [m,fm]]
+              -- tailpol = [[k,fk],...,[k0,fk0]] (k0<...<k<i1<...<m)
+              term := first middlepol
+              j := first middlelookuplist
+              if j = 0 then j := lookup(term.coeff)$GF
+              j := j + 1 -- lookup(f{i1})$GF + 1
+              j rem sizeGF = 0 =>
+                -- in this case one has to increase f{i2}
+                -- tailpol = [[i1,f{i1}],...,[k0,f{k0}]]
+                tailpol   := cons(term, tailpol)
+                middlepol := rest middlepol -- [[i2,f{i2}],...,[m,fm]]
+                taillookuplist   := cons(j, taillookuplist)
+                middlelookuplist := rest middlelookuplist
+              -- otherwise set f{i1} := index(j)$GF
+              setelt(first middlepol, coeff, index(j :: PI)$GF)
+              setfirst_!(middlelookuplist, j)
+              if empty? taillookuplist then
+                pol := listToSUP append(headpol, reverse
+                                cons(constterm, middlepol))
+                --
+                -- may be improved by excluding reciprocal polynomials
+                --
+                normal? pol and primitive? pol => return pol
+              else
+                -- go back to fk
+                -- middlepol = [[k,fk],...,[m,fm]]
+                middlepol := cons(first tailpol, middlepol)
+                tailpol := rest tailpol
+                middlelookuplist := cons(first taillookuplist,
+                                                 middlelookuplist)
+                taillookuplist := rest taillookuplist
+            if weight = n1 then notReady := false
+            else
+              -- must search for polynomial with greater weight
+              if empty? s then -- compute s(f)
+                restfcopy := rest rest fcopy
+                for entry in restfcopy repeat s := cons(entry.expnt, s)
+              s1 := nextSubset(rest s, n2) :: L NNI
+              s  := cons(0, s1)
+              weight := #s
+              taillookuplist := []
+              middlelookuplist := cons(sizeGF, new((weight-2)::NNI, 1))
+              tailpol   := []
+              -- middlepol = [[s.2,0], [s.3,1], ... , [s.weight,1] :
+              middlepol := [[first s1, 0]$Rec]
+              while not empty? (s1 := rest s1) repeat
+                middlepol := cons([first s1, 1]$Rec, middlepol)
+              middlepol := reverse middlepol
+          -- next polynomial must have greater a
+          la := la + 1
+          a  := index(la :: PI)$GF
+        -- next polynomial must have greater constant term
+        lc := lc + 1
+        c  := index(lc :: PI)$GF
+        la := 1
+        a  := 1$GF
+      "failed"
+
+    nextPrimitiveNormalPoly f == nextNormalPrimitivePoly f
+
+    createIrreduciblePoly n ==
+      x := monomial(1,1)$SUP
+      n = 1 => x
+      xn := monomial(1,n)$SUP
+      n >= sizeGF => nextIrreduciblePoly(xn + x) :: SUP
+      -- (since in this case there is most no irreducible binomial X+a)
+      odd? n => nextIrreduciblePoly(xn + 1) :: SUP
+      nextIrreduciblePoly(xn) :: SUP
+
+    createPrimitivePoly n ==
+    -- (see also the comments in the code of nextPrimitivePoly)
+      xn := monomial(1,n)$SUP
+      n = 1 => xn + monomial(-primitiveElement()$GF, 0)$SUP
+      c0 : GF := (-1)**n * primitiveElement()$GF
+      constterm : Rec := [0, c0]$Rec
+      -- try first (probably faster) the polynomials
+      -- f = X**n + f{n-1}*X**(n-1) +...+ f1*X + c0 for which
+      -- fi is 0 or 1 for i=1,...,n-1,
+      -- and this in the order used to define nextPrimitivePoly
+      s  : L NNI := [0,1]
+      weight : NNI := 2
+      s1 : L NNI := [1]
+      n1 : NNI := (n - 1) :: NNI
+      notReady : Boolean := true
+      while notReady repeat
+        polRepr : Repr := [constterm]
+        while not empty? s1 repeat
+          polRepr := cons([first s1, 1]$Rec, polRepr)
+          s1 := rest s1
+        polRepr := cons([n, 1]$Rec, polRepr)
+        --
+        -- may be improved by excluding reciprocal polynomials
+        --
+        primitive? (pol := listToSUP polRepr) => return pol
+        if weight = n then notReady := false
+        else
+          s1 := nextSubset(rest s, n1) :: L NNI
+          s  := cons(0, s1)
+          weight := #s
+      -- if there is no primitive f of the above form
+      -- search now from the beginning, allowing arbitrary
+      -- coefficients f_i, i = 1,...,n-1
+      nextPrimitivePoly(xn + monomial(c0, 0)$SUP) :: SUP
+
+    createNormalPoly n  ==
+      n = 1 => monomial(1,1)$SUP + monomial(-1,0)$SUP
+      -- get a normal polynomial f = X**n + a * X**(n-1) + ...
+      -- with a = -1
+      -- [recall that if f is normal over the field GF of order q
+      -- then a = -(x + x**q +...+ x**(q**n)) can not be zero;
+      -- hence the existence of such an f follows from the
+      -- normal basis theorem ([LN] p.60, Th. 2.35) and the
+      -- surjectivity of the trace ([LN] p.55, Th. 2.23 (iii))]
+      nextNormalPoly(monomial(1,n)$SUP
+                       + monomial(-1, (n-1) :: NNI)$SUP) :: SUP
+
+    createNormalPrimitivePoly n ==
+      xn := monomial(1,n)$SUP
+      n = 1 => xn + monomial(-primitiveElement()$GF, 0)$SUP
+      n1  : NNI := (n - 1) :: NNI
+      c0  : GF  := (-1)**n * primitiveElement()$GF
+      constterm  := monomial(c0, 0)$SUP
+      -- try first the polynomials f = X**n + a *  X**(n-1) + ...
+      -- with a = -1
+      pol := xn + monomial(-1, n1)$SUP + constterm
+      normal? pol and primitive? pol => pol
+      res := nextNormalPrimitivePoly(pol)
+      res case SUP => res
+      -- if there is no normal primitive f with a = -1
+      -- get now one with arbitrary (non-zero) a
+      -- (the existence is proved in [LS])
+      pol := xn + monomial(1, n1)$SUP + constterm
+      normal? pol and primitive? pol => pol
+      nextNormalPrimitivePoly(pol) :: SUP
+
+    createPrimitiveNormalPoly n == createNormalPrimitivePoly n
+
+--    qAdicExpansion m ==
+--      ragits : List I := wholeRagits(m :: (RadixExpansion sizeGF))
+--      pol  : SUP := 0
+--      expt : NNI := #ragits
+--      for i in ragits repeat
+--        expt := (expt - 1) :: NNI
+--        if i ^= 0 then pol := pol + monomial(index(i::PI)$GF, expt)
+--      pol
+
+--    random == qAdicExpansion(random()$I)
+
+--    random n ==
+--      pol := monomial(1,n)$SUP
+--      n1 : NNI := (n - 1) :: NNI
+--      for i in 0..n1 repeat
+--        if (c := random()$GF) ^= 0 then
+--          pol := pol + monomial(c, i)$SUP
+--      pol
+
+    random n ==
+      polRepr : Repr := []
+      n1 : NNI := (n - 1) :: NNI
+      for i in 0..n1 repeat
+        if (c := random()$GF) ^= 0 then
+          polRepr := cons([i, c]$Rec, polRepr)
+      cons([n, 1$GF]$Rec, polRepr) pretend SUP
+
+    random(m,n) ==
+      if m > n then (m,n) := (n,m)
+      d : NNI := (n - m) :: NNI
+      if d > 1 then n := ((random()$I rem (d::PI)) + m) :: PI
+      random(n)
+
+@
+\begin{verbatim}
+-- 12.05.92: JG: long lines
+-- 25.02.92: AS: normal? changed. Now using reducedQPowers
+-- 05.04.91: JG: error in createNormalPrimitivePoly was corrected
+\end{verbatim}
+\section{License}
+<<license>>=
+--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
+--All rights reserved.
+--
+--Redistribution and use in source and binary forms, with or without
+--modification, are permitted provided that the following conditions are
+--met:
+--
+--    - Redistributions of source code must retain the above copyright
+--      notice, this list of conditions and the following disclaimer.
+--
+--    - Redistributions in binary form must reproduce the above copyright
+--      notice, this list of conditions and the following disclaimer in
+--      the documentation and/or other materials provided with the
+--      distribution.
+--
+--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
+--      names of its contributors may be used to endorse or promote products
+--      derived from this software without specific prior written permission.
+--
+--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+@
+<<*>>=
+<<license>>
+
+<<package FFPOLY FiniteFieldPolynomialPackage>>
+@
+\eject
+\begin{thebibliography}{99}
+\bibitem{1} nothing
+\end{thebibliography}
+\end{document}