path: root/src/algebra/ffpoly.spad.pamphlet

\documentclass{article}
\usepackage{open-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 not one?(qe rem d) 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
      not one? leadingCoefficient f => 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 :
      --
      not one? lift(x ** qn1)$MM => 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
      not one? leadingCoefficient f => 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 not one?(lcf := leadingCoefficient f) 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
      term : Rec
      for term: free 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 not one?(lcf := leadingCoefficient f) 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
      term : Rec
      for term: free 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 not one?(lcf := leadingCoefficient f) 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
      term : Rec
      for term: free 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 not one?(lcf := leadingCoefficient f) 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
      term : Rec
      for term: free 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(d)$I + 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}