ffdoms.spad

)abbrev domain IPF InnerPrimeField
++ Authors: N.N., J.Grabmeier, A.Scheerhorn
++ Date Created: ?, November 1990, 26.03.1991
++   Revision 05.06.2018 by J. Grabmeier
++ Basic Operations:
++ Related Constructors: PrimeField
++ Also See:
++ AMS Classifications:
++ Keywords: prime characteristic, prime field, finite field
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encycoldia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in Axiom.
++   Axiom Technical Report Series, ATR/5 NP2522.
++  K. Pommerening: lectures on Cryptology,
++   http://www.staff.uni-mainz.de/pommeren/Cryptology
++ Description:
++   InnerPrimeField(p) implements the field with p elements by
++   using IntegerMod p.
++   Note: argument p MUST be a prime (this domain does not check).
++   See \spadtype{PrimeField} for a domain that does check.
++   In addition to the inherited operations of IntegerMod p, the domain
++   provides exploits the structure of the cyclic group of its invertible
++   elements.  It stores a primitive element w, i.a. generator of this
++   group and it stores a logarithm table for w as soon as this is required.
++   sqrt was added in 2018.

InnerPrimeField(p : PositiveInteger) : Exports == Implementation where

  I   ==> Integer
  NNI ==> NonNegativeInteger
  PI  ==> PositiveInteger
  TBL ==> Table(PI, NNI)
  R   ==> Record(key : PI, entry : NNI)
  SUP ==> SparseUnivariatePolynomial
  OF  ==> OutputForm

  Exports ==> Join(FiniteFieldCategory, FiniteAlgebraicExtensionField(%), _
                ConvertibleTo(Integer), Canonical) with
    sqrt : % -> %
      ++ sqrt(x) computes one y such that y^2 = x, error if there is no
      ++ square root, i.e. jacobi(x,p) = -1.
      ++ Implementation according to
      ++ http://www.staff.uni-mainz.de/pommeren/Cryptology/Asymmetric/5_NTh/
    quadraticNonResidue: () ->  %
      ++ quadraticNonResidue() computes the smallest non negative integer,
      ++ which represents a quadratic non residue.

  Implementation ==> IntegerMod p add
-- local variables ==========================================================
    primitiveElementNotPresent? : Boolean := true
    -- gets false after initialization of the primitive Element.
    primitiveElt : PI := 1
    -- primitiveElt stores a primitive element as soon as it is computed by
    -- createPrimitiveElement().
    cyclicGroupSize : NNI := qcoerce(p-1)@NNI
    -- the size of the cyclic group of invertible elements.
    factorsOfCyclicGroupSize: List Record(factor : Integer,
                                exponent : NonNegativeInteger) := []
    -- will hold the factorization of the cyclic group size of the group of
    -- invertible elements, as soon as the logarithm tables are required.
    logarithmTableNotPresent? : Boolean := true
    -- gets false after initialization of the logarithm tables.
    discLogTable : Table(PI, TBL) := table()$Table(PI, TBL)
    -- tables indexed by the prime factors of the size q of the cyclic group
    -- discLogTable.factor is a table with keys
    -- primitiveElement() ^ (i * (q quo factor)) and entries i for
    -- i in 0..n-1, n computed in initialize() in order to use
    -- the minimal size limit 'limit' optimal.

-- local functions ========================================================
    twoPowerOfCyclicGroupSize : NNI :=
      cGS : NNI := cyclicGroupSize
      r : NNI := 0
      while cGS rem 2 = 0 repeat
        cGS := cGS quo 2
        r := r+1
      r

    twoPower: % -> NNI
      ++ twoPower(x: %) computes the highest power of two in the
      ++ factorization of the order of x ~= 0.
    twoPower(x: %): NNI ==
      x = 0 => error "twoPower: argument must not be 0."
      ord : NNI := order x
      r : NNI := 0
      while ord rem 2 = 0 repeat
        ord := ord quo 2
        r := r+1
      r
    initializePrimitiveElement : () -> Void
      ++ initializePrimitiveElement() initializes the computation of a
      ++ primitive Element to be stored.
    initializePrimitiveElement(): Void ==
      factorsOfCyclicGroupSize := factors(factor(cyclicGroupSize)$I
                                         )$(Factored I)
      -- get a primitive element:
      primitiveElt := lookup(createPrimitiveElement())
      -- set initialization flag:
      primitiveElementNotPresent? := false
      void

    initializeLogarithmTable : () -> Void
      ++ initializeLogarithmTable() initializes the computation of a
      ++ logarithm table for the ++ primitive Element to be stored.
    initializeLogarithmTable(): Void ==
      if primitiveElementNotPresent? then initializePrimitiveElement()
      -- now set up tables for discrete logarithm computations
      limit : Integer := 30
      -- the minimum size for the discrete logarithm table
      for rec in factorsOfCyclicGroupSize repeat
        primeDivisor := rec.factor
        base : % := primitiveElement() ^ (cyclicGroupSize quo primeDivisor)
        l : Integer := length primeDivisor
        n : Integer := if odd? l
          then n := shift(primeDivisor, -(l quo 2))
          else n := shift(1, (l quo 2))
        if n < limit then
          d := (primeDivisor - 1) quo limit + 1
          n := (primeDivisor - 1) quo d + 1
        tbl : TBL := table()$TBL
        a : % := 1
        for i in 0..(n-1)::NNI repeat
          insert!([lookup a , i::NNI]$R, tbl)$TBL
          a := a * base
        insert!([primeDivisor::PI, copy(tbl)$TBL]$Record(key: PI, entry: TBL),
                discLogTable)
      -- set initialization flag:
      logarithmTableNotPresent? := false
      void

-- exported functions from % ==============================================
    quadraticNonResidue(): % ==
      found? : Boolean := false
      q : I := 1
      while not found? repeat
        q := q+1
        found? := jacobi(q, p)$IntegerNumberTheoryFunctions = -1
      q :: %
    sqrt(x: %): % ==
      zero? x => x
      jacobi(convert(x)@I, p)$IntegerNumberTheoryFunctions = -1 =>
        error "sqrt: argument does not have a square root by Jacobi symbol."
      3 = (p rem 4) =>
        y : % := x ^ ((p+1) quo 4)
        -- we choose the smaller one of y and -y of their canonical
        -- representatives in 0..p-1
        if convert(y)@I < convert(-y)@I then y else -y
      -- now we must have 1 = (p rem 4) and jacobi(x,p) = 1 (provided p is
      -- really a prime).
      b : % := quadraticNonResidue()
      e : NNI := twoPowerOfCyclicGroupSize
      u : NNI := ((p-1) :: NNI) quo (2^e)
      z : % := x
      r : NNI := twoPower z
      lr : List NNI := [r]
      while r > 0 repeat
        z := z * b ^ (2 ^ ( (e-r) :: NNI))
        r := twoPower z
        lr := cons(r, lr)
      y : % := z ^ ( (u+1) quo 2)
      for r in rest lr  repeat
        y := y / b ^ (2 ^ (e-r-1) :: NNI)
      if convert(y)@I < convert(-y)@I then y else -y

    generator() == 1

    -- This uses Fermat's little theorem x^(p-1)=1 (mod p), so x^(q(p-1)+r)
    -- = x^r (mod p)
    x : % ^ n : Integer ==
      zero?(n) => 1
      zero?(x) => 0
      r := positiveRemainder(n, p-1)::NNI
      ((x pretend IntegerMod p) ^$IntegerMod(p) r) pretend %

    if p <= convert(max()$SingleInteger)@Integer then
      q := p::SingleInteger
      recip x ==
        zero?(y := convert(x)@Integer :: SingleInteger) => "failed"
        invmod(y, q)::Integer::%
    else
      recip x ==
        zero?(y := convert(x)@Integer) => "failed"
        invmod(y, p)::%

    convert(x : %) : I == x pretend I

    normalElement() == 1

    createNormalElement() == 1

    characteristic() == p

    factorsOfCyclicGroupSize() ==
      -- this fixes an infinite loop of function calls, problem was that
      -- factors factor(1) is the empty list:
      p = 2 => factorsOfCyclicGroupSize
      if empty? factorsOfCyclicGroupSize then initializePrimitiveElement()
      factorsOfCyclicGroupSize

    representationType() == "prime"

    tableForDiscreteLogarithm(fac) ==
      if logarithmTableNotPresent? then initializeLogarithmTable()
      tbl := search(fac::PI, discLogTable)$Table(PI, TBL)
      tbl case "failed" =>
        error "tableForDiscreteLogarithm: argument must be prime divisor_
 of the order of the multiplicative group"
      tbl

    primitiveElement() ==
      if primitiveElementNotPresent? then initializePrimitiveElement()
      index(primitiveElt)


    degree(x) : PI == 1::PositiveInteger
    extensionDegree() : PI == 1::PositiveInteger

--    sizeOfGroundField() == p::NonNegativeInteger

    inGroundField?(x)  == true

    coordinates(x : %) : Vector(%) == new(1, x)$(Vector %)

    represents(v)  == v(1)

    retract(x) == x

    retractIfCan(x) == x

    basis() == new(1, 1::%)$(Vector %)
    basis(n : PI) ==
      n = 1 => basis()
      error("basis: argument must divide extension degree")

    definingPolynomial() ==
      monomial(1, 1)$(SUP %) - monomial(1, 0)$(SUP %)


    minimalPolynomial(x) ==
      monomial(1, 1)$(SUP %) - monomial(x, 0)$(SUP %)

    charthRoot(x : %) : % == x

)abbrev domain PF PrimeField
++ Authors: N.N.,
++ Date Created: November 1990, 26.03.1991
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: prime characteristic, prime field, finite field
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encycoldia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++ Description:
++   PrimeField(p) implements the field with p elements if p is a
++   prime number.
++   Error: if p is not prime.
--++   with new compiler, want to put the error check before the add
PrimeField(p : PositiveInteger) : Exp == Impl where
  Exp ==> Join(FiniteFieldCategory, FiniteAlgebraicExtensionField(%), _
      ConvertibleTo(Integer), Canonical) with
    sqrt : % -> %
      ++ sqrt(x) computes one y such that y^2 = x, error if there is no square
      ++ root, i.e. jacobi(x,p) = -1.
  Impl ==>  InnerPrimeField(p) add
    if not prime?(p)$IntegerPrimesPackage(Integer) then
      error "Argument to prime field must be a prime"

)abbrev package FFPOLY FiniteFieldPolynomialPackage
++ Author: A. Bouyer, J. Grabmeier, A. Scheerhorn, R. Sutor, B. Trager
++ Date Created: January 1991
++ 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, ATR/5 NP2522.
++ 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}.
    lexSmaller? : (SUP, SUP) -> Boolean
      ++ lexSmaller?(f, g) compares monic f and g of the same
      ++ degree in the following order.
      ++ Error: if f or g is not monic or if f and g have different
      ++ degrees or if common degree is 0.
      ++ \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}.
    clexSmaller? : (SUP, SUP) -> Boolean
      ++ clexSmaller?(f, g) compares monic f and g of the same
      ++ degree in the following order.
      ++ Error: if f or g is not monic or if f and g have different
      ++ degrees or if common degree is 0.
      ++ \spad{f < g} if the constant term of \spad{f} is zero
      ++ and constant term of \spad{g} is nonzero.
      ++ If both
      ++ constant term of \spad{f} and \spad{g} are nonzero
      ++ then \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{lexSmaller?} is used
      ++ as ordering predicate.
    nlexSmaller? : (SUP, SUP) -> Boolean
      ++ nlexSmaller?(f, g) compares monic f and g of the same
      ++ degree \spad{n} in the following order.
      ++ Error: if f or g is not monic or if f and g have different
      ++ degrees or if common degree is 0.
      ++ \spad{f < g} if the coefficient of the term of degree {\em n-1}
      ++ of \spad{f} is zero and than that for g is nonzero.
      ++ Also, \spad{f < g} if both coefficients are nonzero and
      ++ {\em lookup} of the coefficient of \spad{f} is less than
      ++ that for \spad{g}.
      ++ In case those coefficients are equal, then \spad{lexSmaller?}
      ++ is used as ordering predicate.
    cnlexSmaller? : (SUP, SUP) -> Boolean
      ++ cnlexSmaller?(f, g) compares monic f and g of the same
      ++ degree \spad{n} in the following order.
      ++ Error: if f or g is not monic or if f and g have different
      ++ degrees or if common degree is 0.
      ++ \spad{f < g} if the constant term of \spad{f} is zero
      ++ and constant term of \spad{g} is nonzero.  If both
      ++ constant term of \spad{f} and \spad{g} are nonzero
      ++ then \spad{f < g} if the {\em lookup} of the constant term
      ++ of f is less than
      ++ this number for g.
      ++ If constant terms are equal then \spad{nlexSmaller?}
      ++ is used as ordering predicate.
    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.
      ++ \spad{lexSmaller?} is used as ordering predicate.
    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.
      ++ \spad{clexSmaller?} is used as ordering predicate.
    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.
      ++ \spad{nlexSmaller?} is used as ordering predicate.
    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.
      ++ \spad{cnlexSmaller?} is used as ordering predicate.
      ++ 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.
      ++ \spad{cnlexSmaller?} is used as ordering predicate.
      ++ 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 from IntegerNumberTheoryFunctions

    sizeGF : PI := size()$GF :: PI

    char_GF : PI := characteristic()$GF :: PI

    deg_GF : NNI := 0

    get_deg_GF() : PI ==
        if deg_GF = 0 then
            n : PI := 1$PI;
            ss : NNI := sizeGF;
            while ss > char_GF repeat
                ss := (ss exquo$NNI char_GF)::NNI
                n := n +$PI 1$PI
            deg_GF := n
        deg_GF::PI

    elem1 : GF := index(1)$GF

    poly_or_prime_rep := representationType()$GF case "polynomial"
                         or representationType()$GF case "prime"

    MM := ModMonic(GF, SUP)

    rep_to_SUP(l) ==> l pretend SUP

    listToSUP(l : Repr) : SUP ==
        newl : Repr := [copy t for t in l]
        rep_to_SUP(newl)

    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 := qcoerce(degree(f)$SUP)
      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
      if (lcf := leadingCoefficient f) ~= 1 then f := (inv lcf) * f
      qexp := reducedQPowers(f)
      n := degree(f)$SUP
      b : Matrix GF := transpose(matrix([entries(vectorise(qexp(ii), n))
                                          for ii 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
      for nse in ns repeat
        i : PI := qcoerce(n + 1)
        while nse(i) = 0 repeat
          i := qcoerce(i - 1)
        if i < dim then
          dim := i
          coeffVector := nse
      (coeffVector(1)::SUP) +(+/[monomial(coeffVector.k, _
               sizeGF^((k-2)::NNI))$SUP for k in 2..dim])

--    qEulerPhiCyclotomic n ==
--      n = 1 => qcoerce(sizeGF - 1)
--      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
--      qcoerce((qe - 1) ^ (qcoerce(eulerPhi(n) quo e)) )

    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) ^ (qcoerce(d - lastd))
        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 * _
          qcoerce((qe - 1) ^ (qcoerce(eulerPhi(d) quo e)) )
      qcoerce(q^(qcoerce(n - m)@NNI) * prod quo n)

    -- test if monic irreducible f is primitive
    primitive_i?(f : SUP) : Boolean ==
      -- 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 is irreducible 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
      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
      lrec  := factorList(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

    primitive?(f) ==
        degree(f) = 0 => false
        leadingCoefficient f ~= 1 => false
        coefficient(f, 0) = 0 => false
        -- cheaper and prunes better than computation of x^(q^n - 1)
        not(irreducible?(f)$UnivariateFiniteFieldFactorize(GF, SUP)) => false
        primitive_i?(f)

    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)$UnivariateFiniteFieldFactorize(GF, SUP) => 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

    normal_and_primitive?(f : SUP) : Boolean ==
        normal?(f) and primitive_i?(f)

    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

    lexSmaller?(f : SUP, g : SUP) : Boolean ==
        n := degree(f)
        n ~= degree(g) => error "polynomials must have equal degrees"
        n < 1 => error "polynomials must have positive degree"
        k1 := numberOfMonomials(f)
        k2 := numberOfMonomials(g)
        k1 < k2 => true
        k2 < k1 => false
        f1 := reductum(f)
        g1 := reductum(g)
        while f1 ~= 0 repeat
            k1 := degree(f1)
            k2 := degree(g1)
            k1 < k2 => return true
            k2 < k1 => return false
            f1 := reductum(f1)
            g1 := reductum(g1)
        f1 := reductum(f)
        g1 := reductum(g)
        while f1 ~= 0 repeat
            k1 := lookup(leadingCoefficient(f1))
            k2 := lookup(leadingCoefficient(g1))
            k1 < k2 => return true
            k2 < k1 => return false
            f1 := reductum(f1)
            g1 := reductum(g1)
        false

    ll_cmp(x1 : GF, x2 : GF) : SingleInteger ==
        x1 = 0 and x2 ~= 0 => -1
        x1 ~= 0 and x2 = 0 => 1
        k1 := lookup(x1)
        k2 := lookup(x2)
        k1 < k2 => -1
        k2 < k2 => 1
        0

    clexSmaller?(f : SUP, g : SUP) : Boolean ==
        n := degree(f)
        n ~= degree(g) => error "polynomials must have equal degrees"
        n < 1 => error "polynomials must have positive degree"
        s := ll_cmp(coefficient(f, 0), coefficient(g, 0))
        s < 0 => true
        0 < s => false
        lexSmaller?(f, g)

    nlexSmaller?(f : SUP, g : SUP) : Boolean ==
        n := degree(f)
        n ~= degree(g) => error "polynomials must have equal degrees"
        n < 1 => error "polynomials must have positive degree"
        s := ll_cmp(coefficient(f, (n - 1)::NNI), coefficient(g, (n - 1)::NNI))
        s < 0 => true
        0 < s => false
        lexSmaller?(f, g)

    cnlexSmaller?(f : SUP, g : SUP) : Boolean ==
        n := degree(f)
        n ~= degree(g) => error "polynomials must have equal degrees"
        n < 1 => error "polynomials must have positive degree"
        s := ll_cmp(coefficient(f, 0), coefficient(g, 0))
        s < 0 => true
        0 < s => false
        nlexSmaller?(f, g)

    v_NNI ==> PrimitiveArray(NNI)
    v_GF ==> PrimitiveArray(GF)

    vecs_to_pol(exp_v : v_NNI,  coeff_v : v_GF, w : NNI) : SUP ==
        resl : Repr := []
        for i in 0..w repeat
            resl := cons([exp_v(i), coeff_v(i)], resl)
        rep_to_SUP(resl)

    do_weight(exp_v : v_NNI, ind_v : v_NNI, coeff_v : v_GF,
              min_i : NNI, max_i : NNI, w : NNI,
              tp? : SUP -> Boolean) : Union(SUP, "failed") ==
        -- skip reducible cases
        if sizeGF = 2 and odd?(w) then
            -- fill vectors with expected values
            for i in min_i..(max_i - 1) repeat
                exp_v(i) := i
                coeff_v(i) := elem1
                ind_v(i) := 1
            return "failed"
        repeat
            i := min_i
            while i < max_i repeat
                j := ind_v(i)
                j := j + 1 -- lookup(f_i)$GF + 1
                j = sizeGF =>
                    ind_v(i) := 1
                    coeff_v(i) := elem1
                    i := i + 1
                ind_v(i) := j
                coeff_v(i) := index(j::PI)
                break
            -- if coeffs wrapped around need to take next
            -- exponents tuple
            if i = max_i then
                i := max(min_i, 1)
                while i < max_i repeat
                    j := exp_v(i)
                    j := j + 1
                    j = exp_v(i + 1) =>
                        exp_v(i) := exp_v(i - 1) + 1
                        i := i + 1
                    exp_v(i) := j
                    break
                -- need bigger weight, fail this one
                i = max_i => return "failed"
            pol := vecs_to_pol(exp_v, coeff_v, w)
            -- print(pol::OutputForm)
            tp?(pol) => return pol

    -- run through the possible weights, starting from w
    do_weights(exp_v : v_NNI, ind_v : v_NNI, coeff_v : v_GF,
               min_i : NNI, max_i : NNI, w : NNI,
               tp? : SUP -> Boolean) : Union(SUP, "failed") ==
        n := exp_v(w)
        repeat
            -- skip weight if there nothing else to increase
            if max_i > min_i then
                resu := do_weight(exp_v, ind_v, coeff_v, min_i,
                                  max_i, w, tp?)
                resu case SUP => return resu@SUP
            w = n => return "failed"
            w1 := w
            w := w + 1
            i := w
            for i1 in w1..max_i by -1 repeat
                exp_v(i) := exp_v(i1)
                coeff_v(i) := coeff_v(i1)
                i := i1
            exp_v(max_i) := max_i
            ind_v(max_i) := 1
            coeff_v(max_i) := elem1
            max_i := max_i + 1
            -- artificial, we will bump it in do_weight
            ind_v(min_i) := 0

    get_rep(f) ==> f pretend Repr

    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
        f_rep : Repr := get_rep(f)
        c0 := coefficient(f, 0)
        n = 1 =>
            xn := first(f_rep)
            lc : NNI := (c0 = 0 => 0; lookup(c0))
            lc = sizeGF - 1 => "failed"
            c := index((lc + 1)::PI)
            rep_to_SUP([xn, [0, c]$Rec])
        w : NNI := (#(f_rep) - 1) :: NNI
        if c0 = 0 then w := w + 1
        good_binomials := true
        -- skip binomials if all are reducible
        if n > 1 and w = 1 then
            facs := factors(factor(n))
            good_binomials : Boolean := odd?(char_GF)
            for fac in facs while good_binomials repeat
                fac.factor = 2 and fac.exponent > 1 =>
                    (sizeGF - 1) rem 4 ~= 0 =>
                        good_binomials := false
                if (sizeGF - 1) rem fac.factor ~= 0 then
                    good_binomials := false
            if not(good_binomials) then w := 2
        exp_v : v_NNI := new(n + 1, 0)
        ind_v : v_NNI := new(n + 1, 0)
        coeff_v : v_GF := new(n + 1, 0)
        for term in f_rep for i in 0.. repeat
            exp_v(w - i) := term.expnt
            ci := term.coeff
            coeff_v(w - i) := ci
            ind_v(w - i) := lookup(ci)
        if not(good_binomials) then
            exp_v(1) := 1
            coeff_v(1) := elem1
            ind_v(1) := 1
        do_weights(exp_v, ind_v, coeff_v, 0, w, w,
                   irreducible?$UnivariateFiniteFieldFactorize(GF, SUP))

    lc_Rec ==> Record(nl : PI, nc : GF)

    -- Find nl bigger or equal than l such that m1_to_n*index(nl)
    -- is primitive.  If successful return [nl, index(nl)],
    -- otherwise return "failed"
    get_next_GF_generator(l : NNI, m1_to_n : GF) : Union(lc_Rec, "failed") ==
        l1 := l::PI
        while l1 < sizeGF repeat
            c := index(l1)$GF
            primitive?(m1_to_n*c)$GF => return [l1, c]
            l1 := l1 + 1
        "failed"

    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
        f_rep : Repr := get_rep(f)
        xn : Rec := first(f_rep)
        c0 : GF  := coefficient(f, 0)
        lc : NNI := lookup(c0)$GF
        n = 1 =>
             m1_to_n := -(1$GF)
             -- x + c is primitive iff -c is primitive
             lcu := get_next_GF_generator(lc + 1, m1_to_n)
             lcu case "failed" => "failed"
             c := (lcu@lc_Rec).nc
             rep_to_SUP([xn, [0, c]$Rec])
        w : NNI := (#(f_rep) - 1)::NNI
        if c0 = 0 then w := w + 1
        exp_v : v_NNI := new(n + 1, 0)
        ind_v : v_NNI := new(n + 1, 0)
        coeff_v : v_GF := new(n + 1, 0)
        for term in f_rep for i in 0.. repeat
            exp_v(w - i) := term.expnt
            ci := term.coeff
            coeff_v(w - i) := ci
            ind_v(w - i) := lookup(ci)
        -- 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
        m1_to_n := ((-1$Integer)^n)::GF
        -- if c = 0 then set lc to 1
        if c = 0 then lc := 1
        repeat
            lcu := get_next_GF_generator(lc, m1_to_n)
            lcu case "failed" => return "failed"
            lcr := lcu@lc_Rec
            -- move to next constant term
            lc := lcr.nl + 1
            ind_v(0) := lcr.nl
            coeff_v(0) := c := lcr.nc
            if not(c = c0 and w > 1) then
                -- 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
                w := 2
                exp_v(1) := 1
                ind_v(1) := 0
                coeff_v(1) := 0
                exp_v(w) := n
                ind_v(w) := 1
                coeff_v(w) := 1
            resu := do_weights(exp_v, ind_v, coeff_v, 1, w, w, primitive?)
            resu case SUP => return resu@SUP

    nextNormalPoly f ==
        n : NNI := degree f
        n = 0 => error "polynomial must have positive degree"
        n1 := (n - 1)::NNI
        -- make f monic
        if (lcf := leadingCoefficient(f)) ~= 1 then f := (inv lcf) * f
        f_rep : Repr := get_rep(f)
        a0 : GF  := coefficient(f, n1)
        la  : NNI := lookup(a0)$GF rem sizeGF
        n = 1 =>
            xn : Rec  := first(f_rep)
            -- the polynomial X + a is normal if and only if a is not zero
            la = sizeGF - 1 => "failed"
            rep_to_SUP([xn, [0, index((la + 1)::PI)$GF]$Rec])

        -- 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 la := 1
        if la = 0 then
            la := 1
            a := elem1
        w : NNI := (#(f_rep) - 1) :: NNI
        if coefficient(f, 0) = 0 then w := w + 1
        exp_v : v_NNI := new(n + 1, 0)
        ind_v : v_NNI := new(n + 1, 0)
        coeff_v : v_GF := new(n + 1, 0)
        for term in f_rep for i in 0.. repeat
            exp_v(w - i) := term.expnt
            ci := term.coeff
            coeff_v(w - i) := ci
            ind_v(w - i) := lookup(ci)
        while la < sizeGF repeat
            -- (run through the possible values of a)
            if not(a = a0) then
                w := 2
                ind_v(0) := 0
                coeff_v(0) := 0
                exp_v(w - 1) := n1
                ind_v(w - 1) := la
                coeff_v(w - 1) := a
                exp_v(w) := n
                ind_v(w) := 1
                coeff_v(w) := 1
            resu := do_weights(exp_v, ind_v, coeff_v, 0, (w - 1)::NNI,
                               w, normal?)
            resu case SUP => return resu@SUP
            la := la + 1
            a := index(la::PI)$GF
        "failed"

    nextNormalPrimitivePoly f ==
        n : NNI := degree f
        -- degree 1 primitive polynomial is normal
        n = 1 => nextPrimitivePoly(f)
        n = 0 => error "polynomial must have positive degree"
        n1 := (n - 1)::NNI
        -- make f monic
        if (lcf := leadingCoefficient f) ~= 1 then f := (inv lcf) * f
        f_rep : Repr := get_rep(f)
        c0 : GF  := coefficient(f, 0)
        lc : NNI := lookup(c0)$GF
        a0 : GF  := coefficient(f, n1)
        la : NNI := lookup(a0)$GF
        -- 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 la := 1
        if a = 0 then
            la := 1
            a  := elem1
        m1_to_n := ((-1$Integer)^n)::GF
        w : NNI := (#(f_rep) - 1) :: NNI
        -- if c = 0 then set lc to 1
        if c = 0 then
            lc := 1
            w := w + 1
        exp_v : v_NNI := new(n + 1, 0)
        ind_v : v_NNI := new(n + 1, 0)
        coeff_v : v_GF := new(n + 1, 0)
        for term in f_rep for i in 0.. repeat
            exp_v(w - i) := term.expnt
            ci := term.coeff
            coeff_v(w - i) := ci
            ind_v(w - i) := lookup(ci)
        first_pass := true
        repeat
            lcu := get_next_GF_generator(lc, m1_to_n)
            lcu case "failed" => return "failed"
            lcr := lcu@lc_Rec
            lc := lcr.nl
            c := lcr.nc
            while la < sizeGF repeat
                -- We use original coefficients only on first
                -- pass here and then only when a0 and c0 were OK.
                -- Otherwise we need to reinitialize vectors
                if not (a = a0 and c = c0) then
                    first_pass := false
                    w := 2
                    ind_v(0) := lc
                    coeff_v(0) := c
                    exp_v(w - 1) := n1
                    ind_v(w - 1) := la
                    coeff_v(w - 1) := a
                    exp_v(w) := n
                    ind_v(w) := 1
                    coeff_v(w) := 1
                if w = 2 and not(first_pass) then
                    pol := vecs_to_pol(exp_v, coeff_v, w)
                    normal_and_primitive?(pol) => return pol
                resu := do_weights(exp_v, ind_v, coeff_v, 1,
                                   (w - 1)::NNI, w, normal_and_primitive?)
                resu case SUP => return resu@SUP
                la := la + 1
                a  := index(la :: PI)$GF
            la := 1
            a := 1
            lc := lc + 1

    nextPrimitiveNormalPoly f == nextNormalPrimitivePoly f

    createIrreduciblePoly n ==
        x := monomial(1, 1)$SUP
        n = 1 => x
        xn := monomial(1, n)$SUP
        -- does not make much sense, but we do this to keep output
        -- the same as in earlier versions
        n >= sizeGF => nextIrreduciblePoly(xn + x) :: SUP

        dg := get_deg_GF()
        kg : GF :=
            poly_or_prime_rep =>
                (l := gcd(n, dg)) > 1 =>
                    index((char_GF^((dg exquo l)::NNI) - 1)::PI)
                odd? n => 1
                0
            0
        nextIrreduciblePoly(xn + kg::SUP)::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 (there should be one)
      pol := xn + monomial(-1, n1)$SUP + constterm
      normal? pol and primitive? pol => pol
      res := nextNormalPrimitivePoly(pol)
      res case SUP => res
      print(message("weird, can not find normal and primitive"
                    "poly with a{n-1} = -1")$OutputForm)
      -- if there is no normal primitive f with a = -1
      -- get now one with arbitrary (non-zero) a
      -- (the existence is proved in [LS])
      -- FIXME: elem1 instead of 1 ????
      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)
      rep_to_SUP(cons([n, 1$GF]$Rec, polRepr))

    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)


)abbrev package FFPOLY2 FiniteFieldPolynomialPackage2
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations: rootOfIrreduciblePoly
++ Related Constructors: FiniteFieldCategory
++ Also See:
++ AMS Classifications:
++ Keywords: finite field, zeros of polynomials, Berlekamp's trace algorithm
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in Axiom.
++   Axiom Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldPolynomialPackage2(F, GF) exports some functions concerning
++  finite fields, which depend on a finite field {\em GF} and an
++  algebraic extension F of {\em GF}, e.g. a zero of a polynomial
++  over {\em GF} in F.
FiniteFieldPolynomialPackage2(F, GF) : Exports == Implementation where
  F : FieldOfPrimeCharacteristic with
      coerce : GF -> %
        ++ coerce(x) \undocumented{}
      lookup : % -> PositiveInteger
        ++ lookup(x) \undocumented{}
      basis : PositiveInteger -> Vector %
        ++ basis(n) \undocumented{}
      Frobenius : % -> %
        ++ Frobenius(x) \undocumented{}
  -- F should be a algebraic extension of the finite field GF, either an
  -- algebraic closure of GF or a simple algebraic extension field of GF
  GF : FiniteFieldCategory

  I   ==> Integer
  NNI ==> NonNegativeInteger
  PI  ==> PositiveInteger
  SUP ==> SparseUnivariatePolynomial
  OUT ==> OutputForm
  M   ==> Matrix
  V   ==> Vector
  L   ==> List
  FFPOLY ==> FiniteFieldPolynomialPackage(GF)
  SUPF2 ==> SparseUnivariatePolynomialFunctions2(GF, F)

  Exports ==> with

    rootOfIrreduciblePoly : SUP GF -> F
      ++ rootOfIrreduciblePoly(f) computes one root of the monic,
      ++ irreducible polynomial f, which degree must divide the extension degree
      ++ of {\em F} over {\em GF},
      ++ i.e. f splits into linear factors over {\em F}.


  Implementation ==> add

-- we use berlekamps trace algorithm
-- it is not checked whether the polynomial is irreducible over GF]]
    rootOfIrreduciblePoly(pf) ==
--    not irreducible(pf)$FFPOLY =>
--      error("polynomial has to be irreducible")
      sizeGF := size()$GF
      -- if the polynomial is of degree one, we're ready
      deg := degree(pf)$(SUP GF)::PI
      deg = 0 => error("no roots")
      deg = 1 => -coefficient(pf, 0)$(SUP GF)::F
      p : SUP F := map(coerce, pf)$SUPF2
      -- compute qexp, qexp(i) = x ^(size()GF ^ i) mod p
      -- with this list it's easier to compute the gcd(p(x), trace(x))
      qexp := reducedQPowers(pf)$FFPOLY
      stillToFactor := p
      -- take linear independent elements, the basis of F over GF
      basis : Vector F := basis(deg)$F
      basispointer : I := 1
      -- as p is irreducible over GF, 0 can't be a root of p
      -- therefore we can use the predicate zero?(root) for indicating
      -- whether a root is found
      root := 0$F
      while zero?(root)$F repeat
        beta : F := basis.basispointer
        -- gcd(trace(x)+gf, p(x)) has degree 0, that's why we skip beta=1
        if beta = 1$F then
          basispointer := basispointer + 1
          beta := basis.basispointer
        basispointer := basispointer+1
        -- compute the polynomial trace(beta * x) mod p(x) using explist
        trModp : SUP F := map(coerce, qexp.0)$SUPF2 * beta
        for i in 1..deg-1 repeat
          beta := Frobenius(beta)
          trModp := trModp +$(SUP F) beta *$(SUP F) map(coerce, qexp.i)$SUPF2
        -- if it is of degree 0, it doesn't help us finding a root
        if degree(trModp)$(SUP F) > 0 then
          -- for all elements gf of GF do
          for j in 1..sizeGF repeat
            -- compute gcd(trace(beta * x) + gf, stillToFactor)
            h := gcd(stillToFactor, trModp +$(SUP F) _
             (index(qcoerce(j))$GF::F::(SUP F)))$(SUP F)
            -- make the gcd polynomial monic
            if leadingCoefficient(h)$(SUP F) ~= 1$F then
              h := (inv leadingCoefficient(h)) * h
            degh := degree(h)$(SUP F)
            degSTF := degree(stillToFactor)$(SUP F)
            -- if the gcd has degree one we are ready
            degh = 1 => root := -coefficient(h, 0)$(SUP F)
            -- if the quotient of stillToFactor and the gcd has
            -- degree one, we're also ready
            degSTF - degh = 1 =>
              root := -coefficient(stillToFactor quo h, 0)$(SUP F)
            -- otherwise the gcd helps us finding a root, only if its
            -- degree is between 2 and degree(stillToFactor)-2
            if degh > 1 and degh < degSTF then
              2*degh > degSTF => stillToFactor := stillToFactor quo h
              stillToFactor := h
      root


)abbrev domain FFP FiniteFieldExtensionByPolynomial
++ Authors: R.Sutor, J. Grabmeier, O. Gschnitzer, A. Scheerhorn
++ Basic Operations:
++ Related Constructors:
++ Also See: FiniteFieldCyclicGroupExtensionByPolynomial,
++  FiniteFieldNormalBasisExtensionByPolynomial
++ AMS Classifications:
++ Keywords: field, extension field, algebraic extension,
++   finite extension, finite field, Galois field
++ Reference:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics an
++   Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++   FiniteFieldExtensionByPolynomial(GF, defpol) implements the extension
++   of the finite field {\em GF} generated by the extension polynomial
++   {\em defpol} which MUST be irreducible.
++   Note: the user has the responsibility to ensure that
++   {\em defpol} is irreducible.

FiniteFieldExtensionByPolynomial(GF : FiniteFieldCategory, _
  defpol : SparseUnivariatePolynomial GF) : Exports == Implementation where
--  GF     : FiniteFieldCategory
--  defpol : SparseUnivariatePolynomial GF

  PI   ==> PositiveInteger
  NNI  ==> NonNegativeInteger
  SUP  ==> SparseUnivariatePolynomial
  I    ==> Integer
  R    ==> Record(key : PI, entry : NNI)
  TBL  ==> Table(PI, NNI)
  SAE  ==> SimpleAlgebraicExtension(GF, SUP GF, defpol)
  OUT  ==> OutputForm

  Exports ==> FiniteAlgebraicExtensionField(GF)

  Implementation ==>  add

-- global variables ====================================================

    Rep := SAE

    extdeg : PI        := qcoerce(degree(defpol)$(SUP GF))
    -- the extension degree

    alpha            := new()$Symbol :: OutputForm
    -- a new symbol for the output form of field elements

    sizeCG : Integer := size()$GF^extdeg - 1
    -- the order of the multiplicative group

    facOfGroupSize := []$(List Record(factor : Integer, exponent : NonNegativeInteger))
    -- the factorization of sizeCG

    normalElt : PI := 1
    -- for the lookup of the normal Element computed by
    -- createNormalElement

    primitiveElt : PI := 1
    -- for the lookup of the primitive Element computed by
    -- createPrimitiveElement()

    initlog? : Boolean := true
    -- gets false after initialization of the discrete logarithm table

    initelt? : Boolean := true
    -- gets false after initialization of the primitive and the
    -- normal element


    discLogTable : Table(PI, TBL) := table()$Table(PI, TBL)
    -- tables indexed by the factors of sizeCG,
    -- discLogTable(factor) is a table  with keys
    -- primitiveElement() ^ (i * (sizeCG quo factor)) and entries i for
    -- i in 0..n-1, n computed in initialize() in order to use
    -- the minimal size limit 'limit' optimal.

-- functions ===========================================================

--    createNormalElement() ==
--      a := primitiveElement()
--      nElt := generator()
--      for i in 1.. repeat
--        normal? nElt => return nElt
--        nElt := nElt*a
--      nElt

    generator() == reduce(monomial(1, 1)$SUP(GF))$Rep
    norm x   == resultant(defpol, lift x)

    initializeElt : () -> Void
    initializeLog : () -> Void
    basis(n : PI) ==
      (extdeg rem n) ~= 0 => error "argument must divide extension degree"
      a : % := norm(primitiveElement(), n)
      vector [a^i for i in 0..n-1]

    degree(x) : PI ==
      y : % := 1
      m := zero(extdeg, extdeg+1)$(Matrix GF)
      for i in 1..extdeg+1 repeat
        setColumn!(m, i, coordinates(y))$(Matrix GF)
        y := y*x
      rank(m)::PI

    minimalPolynomial(x : %) ==
      y : % := 1
      m := zero(extdeg, extdeg+1)$(Matrix GF)
      for i in 1..extdeg+1 repeat
        setColumn!(m, i, coordinates(y))$(Matrix GF)
        y := y*x
      v := first nullSpace(m)$(Matrix GF)
      +/[monomial(v.(i+1), i)$(SUP GF) for i in 0..extdeg]


    normal?(x) ==
      l : List List GF := [entries coordinates x]
      a := x
      for i in 2..extdeg repeat
        a := Frobenius(a)
        l := concat(l, entries coordinates a)$(List List GF)
      ((rank matrix(l)$(Matrix GF)) = extdeg::NNI) => true
      false


    a : GF * x : % == a *$Rep x
    n : I * x : % == n *$Rep x
    -x == -$Rep x
    random() == random()$Rep
    coordinates(x : %) == coordinates(x)$Rep
    represents(v) == represents(v)$Rep
    coerce(x : GF) : % == coerce(x)$Rep
    definingPolynomial() == defpol
    retract(x) == retract(x)$Rep
    retractIfCan(x) == retractIfCan(x)$Rep
    index(x) == index(x)$Rep
    lookup(x) == lookup(x)$Rep
    x : %/y : % == x /$Rep y
    x : %/a : GF == x/coerce(a)
--    x: % / a: GF ==
--      a = 0$GF => error "division by zero"
--      x * inv(coerce(a))
    x : % * y : % == x *$Rep y
    x : % + y : % == x +$Rep y
    x : % - y : % == x -$Rep y
    x : % = y : % == x =$Rep y
    basis() == basis()$Rep
    0 == 0$Rep
    1 == 1$Rep

    factorsOfCyclicGroupSize() ==
      if empty? facOfGroupSize then initializeElt()
      facOfGroupSize

    representationType() == "polynomial"

    tableForDiscreteLogarithm(fac) ==
      if initlog? then initializeLog()
      tbl := search(fac::PI, discLogTable)$Table(PI, TBL)
      tbl case "failed" =>
        error "tableForDiscreteLogarithm: argument must be prime divisor_
 of the order of the multiplicative group"
      tbl

    primitiveElement() ==
      if initelt? then initializeElt()
      index(primitiveElt)

    normalElement() ==
      if initelt? then initializeElt()
      index(normalElt)

    initializeElt() ==
      facOfGroupSize := factors(factor(sizeCG)$Integer)
      -- get a primitive element
      pE := createPrimitiveElement()
      primitiveElt := lookup(pE)
      -- create a normal element
      nElt := generator()
      while not normal? nElt repeat
        nElt := nElt*pE
      normalElt := lookup(nElt)
      -- set elements initialization flag
      initelt? := false
      void()$Void

    initializeLog() ==
      if initelt? then initializeElt()
-- set up tables for discrete logarithm
      limit : Integer := 30
    -- the minimum size for the discrete logarithm table
      for f in facOfGroupSize repeat
        fac := f.factor
        base : % := primitiveElement() ^ (sizeCG quo fac)
        l : Integer := length(fac)$Integer
        n : Integer := 0
        if odd?(l)$Integer then n := shift(fac, -(l quo 2))
                           else n := shift(1, (l quo 2))
        if n < limit then
          d := (fac-1) quo limit + 1
          n := (fac-1) quo d + 1
        tbl : TBL := table()$TBL
        a : % := 1
        for i in (0::NNI)..(n-1)::NNI repeat
          insert!([lookup(a), i::NNI]$R, tbl)$TBL
          a := a*base
        insert!([fac::PI, copy(tbl)$TBL]_
               $Record(key : PI, entry : TBL), discLogTable)$Table(PI, TBL)
      -- set logarithm initialization flag
      initlog? := false
      -- tell user about initialization
      --print("discrete logarithm tables initialized"::OUT)
      void()$Void

    coerce(e : %) : OutputForm == outputForm(lift(e), alpha)

    extensionDegree() : PI == extdeg

    size() == qcoerce(sizeCG + 1)

--  sizeOfGroundField() == size()$GF

    inGroundField?(x) ==
      retractIfCan(x) = "failed" => false
      true

    characteristic() == characteristic()$GF

)abbrev domain FFX FiniteFieldExtension
++ Authors: R.Sutor, J. Grabmeier, A. Scheerhorn
++ Basic Operations:
++ Related Constructors: FiniteFieldExtensionByPolynomial,
++  FiniteFieldPolynomialPackage
++ Also See: FiniteFieldCyclicGroupExtension,
++  FiniteFieldNormalBasisExtension
++ AMS Classifications:
++ Keywords: field, extension field, algebraic extension,
++   finite extension, finite field, Galois field
++ Reference:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics an
++   Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++   FiniteFieldExtensionByPolynomial(GF, n) implements an extension
++   of the finite field {\em GF} of degree n generated by the extension
++   polynomial constructed by
++   \spadfunFrom{createIrreduciblePoly}{FiniteFieldPolynomialPackage} from
++   \spadtype{FiniteFieldPolynomialPackage}.
FiniteFieldExtension(GF, n) : Exports == Implementation where
  GF : FiniteFieldCategory
  n : PositiveInteger
  Exports ==> FiniteAlgebraicExtensionField(GF)
  -- MonogenicAlgebra(GF, SUP) with  -- have to check this
  Implementation ==> FiniteFieldExtensionByPolynomial(GF,
       createIrreduciblePoly(n)$FiniteFieldPolynomialPackage(GF))

)abbrev domain IFF InnerFiniteField
++ Author: ???
++ Basic Operations:
++ Related Constructors: FiniteFieldExtensionByPolynomial,
++  FiniteFieldPolynomialPackage
++ Also See: FiniteFieldCyclicGroup, FiniteFieldNormalBasis
++ AMS Classifications:
++ Keywords: field, extension field, algebraic extension,
++   finite extension, finite field, Galois field
++ Reference:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics an
++   Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++   InnerFiniteField(p, n) implements finite fields with \spad{p^n} elements
++   where p is assumed prime but does not check.
++   For a version which checks that p is prime, see \spadtype{FiniteField}.
InnerFiniteField(p : PositiveInteger, n : PositiveInteger) ==
     FiniteFieldExtension(InnerPrimeField p, n)

)abbrev domain FF FiniteField
++ Author: ???
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: field, extension field, algebraic extension,
++  finite extension, finite field, Galois field
++ Reference:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics an
++   Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteField(p, n) implements finite fields with p^n elements.
++  This packages checks that p is prime.
++  For a non-checking version, see \spadtype{InnerFiniteField}.
FiniteField(p : PositiveInteger, n : PositiveInteger) : _
   FiniteAlgebraicExtensionField(PrimeField p) ==_
   FiniteFieldExtensionByPolynomial(PrimeField p, _
       createIrreduciblePoly(n)$FiniteFieldPolynomialPackage(PrimeField p))

)abbrev domain FFCGP FiniteFieldCyclicGroupExtensionByPolynomial
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations:
++ Related Constructors: FiniteFieldFunctions
++ Also See: FiniteFieldExtensionByPolynomial,
++  FiniteFieldNormalBasisExtensionByPolynomial
++ AMS Classifications:
++ Keywords: finite field, primitive elements, cyclic group
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++  AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldCyclicGroupExtensionByPolynomial(GF, defpol)  implements a
++  finite extension field of the ground field {\em GF}. Its elements are
++  represented by powers of a primitive element, i.e. a generator of the
++  multiplicative (cyclic) group. As primitive
++  element we choose the root of the extension polynomial {\em defpol},
++  which MUST be primitive (user responsibility). Zech logarithms are stored
++  in a table of size half of the field size, and use \spadtype{SingleInteger}
++  for representing field elements, hence, there are restrictions
++  on the size of the field.


FiniteFieldCyclicGroupExtensionByPolynomial(GF, defpol) : _
  Exports == Implementation where
  GF    : FiniteFieldCategory                -- the ground field
  defpol : SparseUnivariatePolynomial GF      -- the extension polynomial
  -- the root of defpol is used as the primitive element

  PI    ==> PositiveInteger
  NNI   ==> NonNegativeInteger
  I     ==> Integer
  SI    ==> SingleInteger
  SUP   ==> SparseUnivariatePolynomial
  SAE   ==> SimpleAlgebraicExtension(GF, SUP GF, defpol)
  V     ==> Vector GF
  FFP   ==> FiniteFieldExtensionByPolynomial(GF, defpol)
  FFF   ==> FiniteFieldFunctions(GF)
  OUT   ==> OutputForm
  ARR   ==> PrimitiveArray(SI)
  TBL   ==> Table(PI, NNI)


  Exports ==> FiniteAlgebraicExtensionField(GF)  with

    getZechTable : () -> ARR
      ++ getZechTable() returns the zech logarithm table of the field
      ++ it is used to perform additions in the field quickly.
  Implementation ==>  add

-- global variables ===================================================

    Rep := SI
    -- elements are represented by small integers in the range
    -- (-1)..(size()-2). The (-1) representing the field element zero,
    -- the other small integers representing the corresponding power
    -- of the primitive element, the root of the defining polynomial

    -- it would be very nice if we could use the representation
    -- Rep := Union("zero", IntegerMod(size()$GF ^ degree(defpol) -1)),
    -- why doesn't the compiler like this ?

    qzero ==> (-1)@SI

    extdeg : NNI  := degree(defpol)$(SUP GF)::NNI
    -- the extension degree

    sizeFF : NNI := qcoerce(size()$GF ^ extdeg)
    -- the size of the field

    if sizeFF > 2^20 then
      error "field too large for this representation"

    sizeCG : SI := qconvert(sizeFF - 1)
    -- the order of the cyclic group

    sizeFG : SI := qconvert(sizeCG quo (size()$GF - 1))
    -- the order of the factor group


    zechlog : ARR := new(((sizeFF+1) quo 2)::NNI, -1::SI)$ARR
    -- the table for the zech logarithm

    alpha := new()$Symbol :: OutputForm
    -- get a new symbol for the output representation of
    -- the elements

    primEltGF : GF :=
      odd?(extdeg)$I => -$GF coefficient(defpol, 0)$(SUP GF)
      coefficient(defpol, 0)$(SUP GF)
    -- the corresponding primitive element of the ground field
    -- equals the trace of the primitive element w.r.t. the ground field

    facOfGroupSize := []$(List Record(factor : Integer, exponent : NNI))
    -- the factorization of sizeCG

    initzech? : Boolean := true
    -- gets false after initialization of the zech logarithm array

    initelt? : Boolean := true
    -- gets false after initialization of the normal element

    normalElt : SI := 0
    -- the global variable containing a normal element

-- functions ==========================================================

    -- for completeness we have to give a dummy implementation for
    -- 'tableForDiscreteLogarithm', although this function is not
    -- necessary in the cyclic group representation case

    tableForDiscreteLogarithm(fac) == table()$TBL


    getZechTable() == zechlog
    initializeZech : () -> Void
    initializeElt : () -> Void

    order(x : %) : PI ==
      zero?(x) =>
        error"order: order of zero undefined"
      (sizeCG quo gcd(sizeCG, qcoerce(convert(x)@I)@NNI))::PI

    primitive?(x : %) ==
      zero?(x) or (x = 1) => false
      gcd(x::Rep, sizeCG)$Rep = 1$Rep => true
      false

    coordinates(x : %) ==
      x = 0 => new(extdeg, 0)$(Vector GF)
      primElement : SAE := convert(monomial(1, 1)$(SUP GF))$SAE
-- the primitive element in the corresponding algebraic extension
      coordinates(primElement ^$SAE convert(x)@I)$SAE

    x : % + y : % ==
      if initzech? then initializeZech()
      x =$Rep qzero => y
      y =$Rep qzero => x
      d : Rep := submod(y, x, sizeCG)$Rep
      d <= shift(sizeCG, -$SI (1$SI)) =>
        zechlog.(convert(d)@I) =$SI -1::SI => qzero
        addmod(x, zechlog.(convert(d)@I), sizeCG)$Rep
      --d: Rep := positiveRemainder(x -$Rep y, sizeCG)$Rep
      d : Rep := (sizeCG -$SI d)::Rep
      addmod(y, zechlog.(convert(d)@I), sizeCG)$Rep
      --positiveRemainder(x +$Rep zechlog.(convert(d)@I) -$Rep d, sizeCG)$Rep


    initializeZech() ==
      zechlog := createZechTable(defpol)$FFF
      -- set initialization flag
      initzech? := false
      void()$Void

    basis(n : PI) ==
      extensionDegree() rem n ~= 0 =>
        error("argument must divide extension degree")
      m := sizeCG quo (size()$GF^n-1)
      [index((1+i*m) ::PI) for i in 0..(n-1)]::Vector %

    n : I * x : % == ((n::GF)::%) * x

    minimalPolynomial(a) ==
      f : SUP % := monomial(1, 1)$(SUP %) - monomial(a, 0)$(SUP %)
      u : % := Frobenius(a)
      while not(u = a) repeat
        f := f * (monomial(1, 1)$(SUP %) - monomial(u, 0)$(SUP %))
        u := Frobenius(u)
      p : SUP GF := 0$(SUP GF)
      while not zero?(f)$(SUP %) repeat
        g : GF := retract(leadingCoefficient(f)$(SUP %))
        p := p+monomial(g, _
                      degree(f)$(SUP %))$(SUP GF)
        f := reductum(f)$(SUP %)
      p

    factorsOfCyclicGroupSize() ==
      if empty? facOfGroupSize then initializeElt()
      facOfGroupSize

    representationType() == "cyclic"

    definingPolynomial() == defpol

    random() == random(qconvert(sizeFF))$Rep -$Rep 1$Rep



    represents(v) ==
      u : FFP := represents(v)$FFP
      u =$FFP 0$FFP => 0
      qconvert(discreteLog(u)$FFP)



    coerce(e : GF) : % ==
      zero?(e)$GF => 0
      log : I := discreteLog(primEltGF, e)$GF::NNI *$I sizeFG
      -- now 1$GF is coerced to 0$Rep which is correct.
      qconvert(positiveRemainder(log, sizeCG))


    retractIfCan(x : %) ==
      zero? x => 0$GF
      u := (x::Rep) exquo$Rep qconvert(sizeFG)
      u = "failed" => "failed"
      primEltGF ^$GF (convert(u::%)@I)

    retract(x : %) ==
      a := retractIfCan(x)
      a="failed" => error "element not in ground field"
      a :: GF

    basis() == [index(i :: PI) for i in 1..extdeg]::Vector %


    inGroundField?(x) ==
      zero? x=> true
      positiveRemainder(x::Rep, qconvert(sizeFG))$Rep =$Rep 0$Rep => true
      false

    discreteLog(b : %, x : %) ==
      zero? x => "failed"
      e := extendedEuclidean(b, sizeCG, x)$Rep
      e = "failed" => "failed"
      e1 : Record(coef1 : %, coef2 : %) := e :: Record(coef1 : %, coef2 : %)
      qcoerce(convert(positiveRemainder(e1.coef1, sizeCG)$Rep)@I)

    - x : % ==
        x =$Rep qzero => qzero
        characteristic() =$I 2 => x
        addmod(x, shift(sizeCG, -1)$SI, sizeCG)

    generator() == 1$SI
    createPrimitiveElement() == 1$SI
    primitiveElement() == 1$SI

    discreteLog(x : %) ==
      zero? x => error "discrete logarithm error"
      qcoerce(convert(x)@I)

    normalElement() ==
      if initelt? then initializeElt()
      normalElt::%

    initializeElt() ==
      facOfGroupSize := factors(factor(sizeCG)$Integer)
      normalElt := createNormalElement()::Rep
      initelt? := false
      void()$Void

    extensionDegree() : PI == qcoerce(extdeg)

    characteristic() == characteristic()$GF

    lookup(x : %) ==
      x =$Rep (-$Rep 1$Rep) => qcoerce(sizeFF)
      qcoerce(convert(x +$Rep 1$Rep)@I)

    index(a : PI) ==
      qconvert(positiveRemainder(a, sizeFF)$I) -$Rep 1$Rep

    0 == (-$Rep 1$Rep)

    1 == 0$Rep

-- to get a "exponent like" output form
    coerce(x : %) : OUT ==
        x =$Rep (-$Rep 1$Rep) => message("0")
        x =$Rep 0$Rep => message("1")
        y : I := lookup(x) - 1
        alpha ^$OUT (y::OUT)

    x : % = y : % ==  x =$Rep y

    x : % * y : % ==
      x =$Rep qzero => qzero
      y =$Rep qzero => qzero
      addmod(x, y, sizeCG)$Rep

    a : GF * x : % == coerce(a)@% * x
    x : %/a : GF == x/coerce(a)@%

--    x: % / a: GF ==
--      a = 0$GF => error "division by zero"
--      x * inv(coerce(a))

    inv(x : %)  ==
      zero?(x) => error "inv: not invertible"
      (x = 1) => 1
      sizeCG -$Rep x

    x : % ^ n : PI == x ^ n::I

    x : % ^ n : NNI == x ^ n::I

    x : % ^ n : I ==
      m : Rep := qconvert(positiveRemainder(n, sizeCG)$I)
      m =$Rep 0$Rep => 1
      x = 0 => 0
      mulmod(m, x, sizeCG::Rep)$Rep

)abbrev domain FFCGX FiniteFieldCyclicGroupExtension
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 04.04.1991
++ Basic Operations:
++ Related Constructors:  FiniteFieldCyclicGroupExtensionByPolynomial,
++  FiniteFieldPolynomialPackage
++ Also See: FiniteFieldExtension, FiniteFieldNormalBasisExtension
++ AMS Classifications:
++ Keywords: finite field, primitive elements, cyclic group
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++  AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldCyclicGroupExtension(GF, n)  implements a extension of degree n
++  over the ground field {\em GF}. Its elements are represented by powers of
++  a primitive element, i.e. a generator of the multiplicative (cyclic) group.
++  As primitive element we choose the root of the extension polynomial, which
++  is created by {\em createPrimitivePoly} from
++  \spadtype{FiniteFieldPolynomialPackage}. Zech logarithms are stored
++  in a table of size half of the field size, and use \spadtype{SingleInteger}
++  for representing field elements, hence, there are restrictions
++  on the size of the field.


FiniteFieldCyclicGroupExtension(GF, extdeg) : _
  Exports == Implementation where
  GF       : FiniteFieldCategory
  extdeg   : PositiveInteger
  PI       ==> PositiveInteger
  FFPOLY         ==> FiniteFieldPolynomialPackage(GF)
  SI       ==> SingleInteger
  Exports  ==> FiniteAlgebraicExtensionField(GF) with
    getZechTable : () -> PrimitiveArray(SingleInteger)
      ++ getZechTable() returns the zech logarithm table of the field.
      ++ This table is used to perform additions in the field quickly.
  Implementation ==> FiniteFieldCyclicGroupExtensionByPolynomial(GF, _
                          createPrimitivePoly(extdeg)$FFPOLY)

)abbrev domain FFCG FiniteFieldCyclicGroup
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 04.04.1991
++ Basic Operations:
++ Related Constructors:  FiniteFieldCyclicGroupExtensionByPolynomial,
++  FiniteFieldPolynomialPackage
++ Also See: FiniteField, FiniteFieldNormalBasis
++ AMS Classifications:
++ Keywords: finite field, primitive elements, cyclic group
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++ Description:
++  FiniteFieldCyclicGroup(p, n) implements a finite field extension of degee n
++  over the prime field with p elements. Its elements are represented by
++  powers of a primitive element, i.e. a generator of the multiplicative
++  (cyclic) group. As primitive element we choose the root of the extension
++  polynomial, which is created by {\em createPrimitivePoly} from
++  \spadtype{FiniteFieldPolynomialPackage}. The Zech logarithms are stored
++  in a table of size half of the field size, and use \spadtype{SingleInteger}
++  for representing field elements, hence, there are restrictions
++  on the size of the field.

FiniteFieldCyclicGroup(p, extdeg) : _
  Exports == Implementation where
  p : PositiveInteger
  extdeg   : PositiveInteger
  PI       ==> PositiveInteger
  FFPOLY         ==> FiniteFieldPolynomialPackage(PrimeField(p))
  SI       ==> SingleInteger
  Exports  ==> FiniteAlgebraicExtensionField(PrimeField(p)) with
    getZechTable : () -> PrimitiveArray(SingleInteger)
      ++ getZechTable() returns the zech logarithm table of the field.
      ++ This table is used to perform additions in the field quickly.
  Implementation ==> FiniteFieldCyclicGroupExtensionByPolynomial(PrimeField(p), _
                          createPrimitivePoly(extdeg)$FFPOLY)


)abbrev package INBFF InnerNormalBasisFieldFunctions
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords: finite field, normal basis
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++   Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  D.R.Stinson: Some observations on parallel Algorithms for fast
++   exponentiation in GF(2^n), Siam J. Comp., Vol.19, No.4, pp.711-717,
++   August 1990
++  T.Itoh, S.Tsujii: A fast algorithm for computing multiplicative inverses
++   in GF(2^m) using normal bases, Inf. and Comp. 78, pp.171-177, 1988
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  InnerNormalBasisFieldFunctions(GF) (unexposed):
++  This package has functions used by
++  every normal basis finite field extension domain.

InnerNormalBasisFieldFunctions(GF) : Exports == Implementation where
  GF    : FiniteFieldCategory       -- the ground field

  PI   ==> PositiveInteger
  NNI  ==> NonNegativeInteger
  I    ==> Integer
  SI   ==> SingleInteger
  SUP  ==> SparseUnivariatePolynomial
  VGF  ==> Vector GF
  M    ==> Matrix
  V    ==> Vector
  L    ==> List
  OUT  ==> OutputForm
  TERM ==> Record(value : GF, index : SI)

  Exports ==> with

      setFieldInfo : (V L TERM, GF) -> Void
        ++ setFieldInfo(m, p) initializes the field arithmetic, where m is
        ++ the multiplication table and p is the respective normal element
        ++ of the ground field GF.
      random    : PI           -> VGF
        ++ random(n) creates a vector over the ground field with random entries.
      index     : (PI, PI)      -> VGF
        ++ index(n, m) is a index function for vectors of length n over
        ++ the ground field.
      pol       : VGF          -> SUP GF
        ++ pol(v) turns the vector \spad{[v0, ..., vn]} into the polynomial
        ++ \spad{v0+v1*x+ ... + vn*x^n}.
      xn         : NNI          -> SUP GF
        ++ xn(n) returns the polynomial \spad{x^n-1}.
      dAndcExp  : (VGF, NNI, SI) -> VGF
        ++ dAndcExp(v, n, k) computes \spad{v^e} interpreting v as an element of
        ++ normal basis field. A divide and conquer algorithm similar to the
        ++ one from D.R.Stinson,
        ++ "Some observations on parallel Algorithms for fast exponentiation in
        ++ GF(2^n)", Siam J. Computation, Vol.19, No.4, pp.711-717, August 1990
        ++ is used. Argument k is a parameter of this algorithm.
      repSq     : (VGF, NNI)    -> VGF
        ++ repSq(v, e) computes \spad{v^e} by repeated squaring,
        ++ interpreting v as an element of a normal basis field.
      expPot    : (VGF, SI, SI)  -> VGF
        ++ expPot(v, e, d) returns the sum from \spad{i = 0} to
        ++ \spad{e - 1} of \spad{v^(q^i*d)}, interpreting
        ++ v as an element of a normal basis field and where q is
        ++ the size of the ground field.
        ++ Note: for a description of the algorithm, see T.Itoh and S.Tsujii,
        ++ "A fast algorithm for computing multiplicative inverses in GF(2^m)
        ++ using normal bases",
        ++ Information and Computation 78, pp.171-177, 1988.
      qPot      : (VGF, I)      -> VGF
        ++ qPot(v, e) computes \spad{v^(q^e)}, interpreting v as an element of
        ++ normal basis field, q the size of the ground field.
        ++ This is done by a cyclic e-shift of the vector v.

-- the semantic of the following functions is obvious from the finite field
-- context, for description see category FAXF
      "^"      :(VGF,I)       -> VGF
        ++ x^n \undocumented{}
        ++ See \spadfunFrom{^}{DivisionRing}
      "*"       :(VGF,VGF)     -> VGF
        ++ x*y \undocumented{}
        ++ See \spadfunFrom{*}{SemiGroup}
      "/"       :(VGF,VGF)     -> VGF
        ++ x/y \undocumented{}
        ++ See \spadfunFrom{/}{Field}
      norm      : (VGF, PI)      -> VGF
        ++ norm(x, n) \undocumented{}
        ++ See \spadfunFrom{norm}{FiniteAlgebraicExtensionField}
      trace     : (VGF, PI)      -> VGF
        ++ trace(x, n) \undocumented{}
        ++ See \spadfunFrom{trace}{FiniteAlgebraicExtensionField}
      inv       : VGF          -> VGF
        ++ inv x \undocumented{}
        ++ See \spadfunFrom{inv}{DivisionRing}
      lookup    : VGF          -> PI
        ++ lookup(x) \undocumented{}
        ++ See \spadfunFrom{lookup}{Finite}
      normal?   : VGF          -> Boolean
        ++ normal?(x) \undocumented{}
        ++ See \spadfunFrom{normal?}{FiniteAlgebraicExtensionField}
      basis     : PI           -> V VGF
        ++ basis(n) \undocumented{}
        ++ See \spadfunFrom{basis}{FiniteAlgebraicExtensionField}
      normalElement : PI         -> VGF
        ++ normalElement(n) \undocumented{}
        ++ See \spadfunFrom{normalElement}{FiniteAlgebraicExtensionField}
      minimalPolynomial : VGF   -> SUP GF
        ++ minimalPolynomial(x) \undocumented{}
        ++ See \spadfunFrom{minimalPolynomial}{FiniteAlgebraicExtensionField}

  Implementation ==> add

-- global variables ===================================================

    sizeGF : NNI := size()$GF
    -- the size of the ground field

    multTable : V L TERM := new(1, []$(L TERM))$(V L TERM)
    -- global variable containing the multiplication table

    trGen : GF := 1$GF
    -- controls the embedding of the ground field

    logq : List SI := [0, 10::SI, 16::SI, 20::SI, 23::SI, 0, 28::SI, _
                             30::SI, 32::SI, 0, 35::SI]
    -- logq.i is about 10*log2(i) for the values <12 which
    -- can match sizeGF. It's used by "^"

    expTable : L L SI := [[], _
        [4::SI, 12::SI, 48::SI, 160::SI, 480::SI, 0], _
        [8::SI, 72::SI, 432::SI, 0], _
        [18::SI, 216::SI, 0], _
        [32::SI, 480::SI, 0], [], _
        [72::SI, 0], [98::SI, 0], [128::SI, 0], [], [200::SI, 0]]
    -- expT is used by "^" to optimize the parameter k
    -- before calling dAndcExp(.., .., k)

-- functions ===========================================================

--  computes a^(-1) = a^((q^extDeg)-2)
--  see reference of function expPot
    inv(a) ==
      b : VGF := qPot(expPot(a, (#a-1)::NNI::SI, 1::SI)$%, 1)$%
      erg : VGF := inv((a *$% b).1 *$GF trGen)$GF *$VGF b

-- "^" decides which exponentiation algorithm will be used, in order to
-- get the fastest computation. If dAndcExp is used, it chooses the
-- optimal parameter k for that algorithm.
    a ^ ex  ==
      e : NNI := positiveRemainder(ex, sizeGF^((#a)::PI)-1)$I :: NNI
      zero?(e)$NNI => new(#a, trGen)$VGF
      (e = 1)$NNI  => copy(a)$VGF
--    inGroundField?(a) => new(#a, ((a.1*trGen) ^$GF e))$VGF
      e1 : SI := (length(e)$I)::SI
      sizeGF >$I 11 =>
        q1 : SI := (length(sizeGF)$I)::SI
        logqe : SI := (e1 quo$SI q1) +$SI 1$SI
        10::SI * (logqe + sizeGF-2) > 15::SI * e1 =>
--        print("repeatedSquaring"::OUT)
          repSq(a, e)
--      print("divAndConquer(a,e,1)"::OUT)
        dAndcExp(a, e, 1)
      logqe : SI := ((10::SI *$SI e1) quo$SI (logq.sizeGF)) +$SI 1$SI
      k : SI := 1$SI
      expT : List SI := expTable.sizeGF
      while (logqe >= expT.k) and not zero? expT.k repeat k := k +$SI 1$SI
      mult : I := (sizeGF - 1) *$I sizeGF ^$I (qcoerce(k - 1)) +$I
              ((logqe +$SI k -$SI 1$SI) quo$SI k)::I -$I 2
      (10*mult) >= (15 * (e1::I)) =>
--      print("repeatedSquaring(a,e)"::OUT)
        repSq(a, e)
--    print(hconcat(["divAndConquer(a,e,"::OUT,k::OUT,")"::OUT])$OUT)
      dAndcExp(a, e, k)

-- computes a^e by repeated squaring
    repSq(b, e) ==
      a := copy(b)$VGF
      (e = 1) => a
      odd?(e)$I => a * repSq(a*a, qcoerce(e quo 2))
      repSq(a*a, qcoerce(e quo 2))

-- computes a^e using the divide and conquer algorithm similar to the
-- one from D.R.Stinson,
-- "Some observations on parallel Algorithms for fast exponentiation in
-- GF(2^n)", Siam J. Computation, Vol.19, No.4, pp.711-717, August 1990
    dAndcExp(a, e, k) ==
      plist : List VGF := [copy(a)$VGF]
      qk : I := sizeGF^(qcoerce(k))
      for j in 2..(qk-1) repeat
        if positiveRemainder(j, sizeGF)=0 then b := qPot(plist.(j quo sizeGF), 1)$%
                            else b := a *$% last(plist)$(List VGF)
        plist := concat(plist, b)
      l : List NNI := []
      ex : I := e
      while not(ex = 0) repeat
        l := concat(l, qcoerce(positiveRemainder(ex, qk)))
        ex := ex quo qk
      if first(l)=0 then erg : VGF := new(#a, trGen)$VGF
                    else erg : VGF := plist.(first(l))
      i : SI := k
      for j in rest(l) repeat
        if j ~= 0 then erg := erg *$% qPot(plist.j, i)$%
        i := i+k
      erg

    a * b ==
      e : SI := (#a)::SI
      erg := zero(#a)$VGF
      for t in multTable(1) repeat
        for j in 1..e repeat
          y := t.value  -- didn't work without defining x and y
          x := t.index
          k : SI := addmod(x, j::SI, e)$SI +$SI 1$SI
          erg(k) := erg(k) +$GF a(j) *$GF b(j) *$GF y
      for i in 1..e-1 repeat
        for j in i+1..e repeat
          for t in multTable(j-i+1) repeat
            y := t.value   -- didn't work without defining x and y
            x := t.index
            k : SI := addmod(x, i::SI, e)$SI +$SI 1$SI
            erg(k) := erg(k) +$GF (a(i)*$GF b(j) +$GF a(j)*$GF b(i)) *$GF y
      erg

    lookup(x) ==
      erg : I := 0
      for j in (#x)..1 by -1 repeat
        erg := (erg * sizeGF) + (lookup(x(j))$GF rem sizeGF)
      erg = 0 => (sizeGF^(#x)) :: PI
      erg :: PI

--  computes the norm of a over GF^d, d must divide extdeg
--  see reference of function expPot below
    norm(a, d) ==
      dSI := d::SI
      r := divide((#a)::SI, dSI)
      not(r.remainder = 0) => error "norm: 2.arg must divide extdeg"
      expPot(a, r.quotient, dSI)$%

--  computes expPot(a, e, d) = sum form i=0 to e-1 over a^(q^id))
--  see T.Itoh and S.Tsujii,
--  "A fast algorithm for computing multiplicative inverses in GF(2^m)
--   using normal bases",
--  Information and Computation 78, pp.171-177, 1988
    expPot(a, e, d) ==
      deg : SI := (#a)::SI
      e = 1 => copy(a)$VGF
      k2 : SI := d
      y := copy(a)
      if bit?(e, 0) then
        erg := copy(y)
        qpot : SI := k2
      else
        erg := new(#a, inv(trGen)$GF)$VGF
        qpot : SI := 0
      for k in 1..length(e) repeat
        y := y *$% qPot(y, k2)
        k2 := addmod(k2, k2, deg)$SI
        if bit?(e, k) then
          erg := erg *$% qPot(y, qpot)
          qpot := addmod(qpot, k2, deg)$SI
      erg

-- computes qPot(a, n) = a^(q^n), q=size of GF
    qPot(e, n) ==
      ei := (#e)::SI
      m : SI := positiveRemainder(n::SI, ei)$SI
      zero?(m) => e
      e1 := zero(#e)$VGF
      for i in m+1..ei repeat e1(i) := e(i-m)
      for i in 1..m    repeat e1(i) := e(ei+i-m)
      e1

    trace(a, d) ==
      dSI := d::SI
      r := divide((#a)::SI, dSI)$SI
      not(r.remainder = 0) => error "trace: 2.arg must divide extdeg"
      v := copy(a(1..dSI))$VGF
      sSI : SI := r.quotient
      for i in 1..dSI repeat
        for j in 1..sSI-1 repeat
          v(i) := v(i) + a(i+j::SI*dSI)
      v

    random(n) ==
      v := zero(n)$VGF
      for i in 1..n repeat v(i) := random()$GF
      v


    xn(m) == monomial(1, m)$(SUP GF) - 1$(SUP GF)

    normal?(x) ==
      gcd(xn(#x), pol(x))$(SUP GF) = 1 => true
      false

    x : VGF / y : VGF == x *$% inv(y)$%


    setFieldInfo(m, n) ==
      multTable := m
      trGen := n
      void()$Void

    minimalPolynomial(x) ==
      dx := #x
      y := new(#x, inv(trGen)$GF)$VGF
      m := zero(dx, dx+1)$(M GF)
      for i in 1..dx+1 repeat
        dy := #y
        for j in 1..dy repeat
          for k in 0..((dx quo dy)-1) repeat
            qsetelt!(m, j+k*dy, i, y.j)$(M GF)
        y := y *$% x
      v := first nullSpace(m)$(M GF)
      pol(v)$%

    basis(n) ==
      bas : (V VGF) := new(n, zero(n)$VGF)$(V VGF)
      for i in 1..n repeat
        uniti := zero(n)$VGF
        qsetelt!(uniti, i, 1$GF)$VGF
        qsetelt!(bas, i, uniti)$(V VGF)
      bas

    normalElement(n) ==
       v := zero(n)$VGF
       qsetelt!(v, 1, 1$GF)
       v
--    normalElement(n) == index(n, 1)$%

    index(degm, n) ==
      m : I := n rem$I (sizeGF ^ degm)
      erg := zero(degm)$VGF
      for j in 1..degm repeat
        erg(j) := index(qcoerce(sizeGF + (m rem sizeGF)))$GF
        m := m quo sizeGF
      erg

    pol(x) ==
        +/[monomial(x(i), (i-1)::NNI)$(SUP GF) for i in 1..(#x)::I]

)abbrev domain FFNBP FiniteFieldNormalBasisExtensionByPolynomial
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations:
++ Related Constructors: InnerNormalBasisFieldFunctions, FiniteFieldFunctions,
++ Also See: FiniteFieldExtensionByPolynomial,
++  FiniteFieldCyclicGroupExtensionByPolynomial
++ AMS Classifications:
++ Keywords: finite field, normal basis
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM .
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldNormalBasisExtensionByPolynomial(GF, uni) implements a
++  finite extension of the ground field {\em GF}. The elements are
++  represented by coordinate vectors with respect to. a normal basis,
++  i.e. a basis
++  consisting of the conjugates (q-powers) of an element, in this case
++  called normal element, where q is the size of {\em GF}.
++  The normal element is chosen as a root of the extension
++  polynomial, which MUST be normal over {\em GF}  (user responsibility)
FiniteFieldNormalBasisExtensionByPolynomial(GF, uni) : Exports == _
    Implementation where
  GF    : FiniteFieldCategory            -- the ground field
  uni   : Union(SparseUnivariatePolynomial GF, _
            Vector List Record(value : GF, index : SingleInteger))

  PI   ==> PositiveInteger
  NNI  ==> NonNegativeInteger
  I    ==> Integer
  SI   ==> SingleInteger
  SUP  ==> SparseUnivariatePolynomial
  V    ==> Vector GF
  M    ==> Matrix GF
  OUT  ==> OutputForm
  TERM ==> Record(value : GF, index : SI)
  R    ==> Record(key : PI, entry : NNI)
  TBL  ==> Table(PI, NNI)
  FFF    ==> FiniteFieldFunctions(GF)
  INBFF  ==> InnerNormalBasisFieldFunctions(GF)

  Exports ==> FiniteAlgebraicExtensionField(GF)  with

      getMultiplicationTable : ()   -> Vector List TERM
        ++ getMultiplicationTable() returns the multiplication
        ++ table for the normal basis of the field.
        ++ This table is used to perform multiplications between field elements.
      getMultiplicationMatrix : ()     -> M
        ++ getMultiplicationMatrix() returns the multiplication table in
        ++ form of a matrix.
      sizeMultiplication : ()   -> NNI
        ++ sizeMultiplication() returns the number of entries in the
        ++ multiplication table of the field.
        ++ Note: the time of multiplication
        ++ of field elements depends on this size.
  Implementation ==> add

-- global variables ===================================================

    Rep := V     -- elements are represented by vectors over GF

    alpha       := new()$Symbol :: OutputForm
    -- get a new Symbol for the output representation of the elements

    initlog? : Boolean := true
    -- gets false after initialization of the logarithm table

    initelt? : Boolean := true
    -- gets false after initialization of the primitive element

    initmult? : Boolean := true
    -- gets false after initialization of the multiplication
    -- table or the primitive element

    extdeg : PI   := 1

    defpol : SUP(GF) := 0$SUP(GF)
    -- the defining polynomial

    multTable : Vector List TERM := new(1, []$(List TERM))
    -- global variable containing the multiplication table

    if uni case (Vector List TERM) then
      multTable := uni@(Vector List TERM)
      extdeg := qcoerce(#multTable)
      vv : V := new(extdeg, 0)$V
      vv(1) := 1$GF
      setFieldInfo(multTable, 1$GF)$INBFF
      defpol := minimalPolynomial(vv)$INBFF
      initmult? := false
    else
      defpol := uni :: SUP(GF)
      extdeg := qcoerce(degree(defpol)$(SUP GF))
      multTable : Vector List TERM := new(extdeg, []$(List TERM))

    basisOutput : List OUT :=
      qs : OUT := ('q)::OUT
      append([alpha, alpha ^$OUT qs], _
        [alpha ^$OUT (qs ^$OUT i::OUT) for i in 2..extdeg-1] )


    facOfGroupSize := []$(List Record(factor : Integer, exponent : NNI))
    -- the factorization of the cyclic group size


    traceAlpha : GF := -$GF coefficient(defpol, (degree(defpol)-1)::NNI)
    -- the inverse of the trace of the normalElt
    -- is computed here. It defines the embedding of
    -- GF in the extension field

    primitiveElt : PI := 1
    -- for the lookup the primitive Element computed by createPrimitiveElement()

    discLogTable : Table(PI, TBL) := table()$Table(PI, TBL)
    -- tables indexed by the factors of sizeCG,
    -- discLogTable(factor) is a table with keys
    -- primitiveElement() ^ (i * (sizeCG quo factor)) and entries i for
    -- i in 0..n-1, n computed in initialize() in order to use
    -- the minimal size limit 'limit' optimal.

-- functions ===========================================================

    initializeLog : ()     -> Void
    initializeElt : ()     -> Void
    initializeMult : ()     -> Void


    coerce(v : GF) : %  == new(extdeg, v /$GF traceAlpha)$Rep
    represents(v)   ==  v::%

    degree(a) : PI ==
      d : PI := 1
      b := qPot(a::Rep, 1)$INBFF
      while (b ~= a) repeat
        b := qPot(b::Rep, 1)$INBFF
        d := d+1
      d

    linearAssociatedExp(x, f) ==
      xm : SUP(GF) := monomial(1$GF, extdeg)$(SUP GF) - 1$(SUP GF)
      r := (f * pol(x::Rep)$INBFF) rem xm
      vectorise(r, extdeg)$(SUP GF)
    linearAssociatedLog(x) ==  pol(x::Rep)$INBFF
    linearAssociatedOrder(x) ==
      xm : SUP(GF) := monomial(1$GF, extdeg)$(SUP GF) - 1$(SUP GF)
      xm quo gcd(xm, pol(x::Rep)$INBFF)
    linearAssociatedLog(b, x) ==
      zero? x => 0
      xm : SUP(GF) := monomial(1$GF, extdeg)$(SUP GF) - 1$(SUP GF)
      e := extendedEuclidean(pol(b::Rep)$INBFF, xm, pol(x::Rep)$INBFF)$(SUP GF)
      e = "failed" => "failed"
      e1 := e :: Record(coef1 : (SUP GF), coef2 : (SUP GF))
      e1.coef1

    getMultiplicationTable() ==
      if initmult? then initializeMult()
      multTable
    getMultiplicationMatrix() ==
      if initmult? then initializeMult()
      createMultiplicationMatrix(multTable)$FFF
    sizeMultiplication() ==
      if initmult? then initializeMult()
      sizeMultiplication(multTable)$FFF

    trace(a : %) == retract trace(a, 1)
    norm(a : %) == retract norm(a, 1)
    generator() == normalElement(extdeg)$INBFF
    basis(n : PI) ==
      (extdeg rem n) ~= 0 => error "argument must divide extension degree"
      [Frobenius(trace(normalElement, n), i) for i in 0..(n-1)]::(Vector %)

    a : GF * x : % == a *$Rep x

    x : %/a : GF == x/coerce(a)
--    x: % / a: GF ==
--      a = 0$GF => error "division by zero"
--      x * inv(coerce(a))


    coordinates(x : %)  == x::Rep

    Frobenius(e) == qPot(e::Rep, 1)$INBFF
    Frobenius(e, n) == qPot(e::Rep, n)$INBFF

    retractIfCan(x) ==
      inGroundField?(x) =>
        x(1) *$GF traceAlpha
      "failed"

    retract(x) ==
      inGroundField?(x) =>
        x(1) *$GF traceAlpha
      error("element not in ground field")

-- to get a "normal basis like" output form
    coerce(x : %) : OUT ==
      l : List OUT := []$(List OUT)
      n : PI := extdeg
      (n = 1) => x(1) :: OUT
      for i in 1..n for b in basisOutput repeat
        if not zero?(x(i)) then
          mon : OUT :=
            (x(i) = 1) => b
            (x(i)::OUT) *$OUT b
          l := cons(mon, l)$(List OUT)
      empty?(l)$(List OUT) => (0::OUT)
      r := reduce("+",l)$(List OUT)
      r

    initializeElt() ==
      facOfGroupSize := factors factor(size()$GF^extdeg-1)$I
      -- get a primitive element
      primitiveElt := lookup(createPrimitiveElement())
      initelt? := false
      void()$Void

    initializeMult() ==
      multTable := createMultiplicationTable(defpol)$FFF
      setFieldInfo(multTable, traceAlpha)$INBFF
      -- reset initialize flag
      initmult? := false
      void()$Void

    initializeLog() ==
      if initelt? then initializeElt()
      -- set up tables for discrete logarithm
      limit : Integer := 30
      -- the minimum size for the discrete logarithm table
      for f in facOfGroupSize repeat
        fac := f.factor
        base : % := index(primitiveElt)^((size()$GF^extdeg -$I 1$I) quo$I fac)
        l : Integer := length(fac)$Integer
        n : Integer := 0
        if odd?(l)$I then n := shift(fac, -$I (l quo$I 2))$I
                     else n := shift(1, l quo$I 2)$I
        if n <$I limit then
          d := (fac -$I 1$I) quo$I limit +$I 1$I
          n := (fac -$I 1$I) quo$I d +$I 1$I
        tbl : TBL := table()$TBL
        a : % := 1
        for i in (0::NNI)..(n-1)::NNI repeat
          insert!([lookup(a), i::NNI]$R, tbl)$TBL
          a := a*base
        insert!([fac::PI, copy(tbl)$TBL]_
               $Record(key : PI, entry : TBL), discLogTable)$Table(PI, TBL)
      initlog? := false
      -- tell user about initialization
      --print("discrete logarithm table initialized"::OUT)
      void()$Void

    tableForDiscreteLogarithm(fac) ==
      if initlog? then initializeLog()
      tbl := search(fac::PI, discLogTable)$Table(PI, TBL)
      tbl case "failed" =>
        error "tableForDiscreteLogarithm: argument must be prime _
divisor of the order of the multiplicative group"
      tbl

    primitiveElement() ==
      if initelt? then initializeElt()
      index(primitiveElt)

    factorsOfCyclicGroupSize() ==
      if empty? facOfGroupSize then initializeElt()
      facOfGroupSize

    extensionDegree() : PI == extdeg

    -- sizeOfGroundField() == qcoerce(size()$GF)

    definingPolynomial() == defpol

    trace(a, d) ==
      v := trace(a::Rep, d)$INBFF
      erg := v
      for i in 2..(extdeg quo d) repeat
        erg := concat(erg, v)$Rep
      erg

    characteristic() == characteristic()$GF

    random() == random(extdeg)$INBFF

    x : % * y : % ==
      if initmult? then initializeMult()
      setFieldInfo(multTable, traceAlpha)$INBFF
      x::Rep *$INBFF y::Rep


    1 == new(extdeg, inv(traceAlpha)$GF)$Rep

    0 == zero(extdeg)$Rep

    size() == size()$GF ^ extdeg

    index(n : PI) == index(extdeg, n)$INBFF

    lookup(x : %) == lookup(x::Rep)$INBFF


    basis() ==
      a := basis(extdeg)$INBFF
      vector([e::% for e in entries a])


    x : % ^ e : I ==
      if initmult? then initializeMult()
      setFieldInfo(multTable, traceAlpha)$INBFF
      (x::Rep) ^$INBFF e

    normal?(x) == normal?(x::Rep)$INBFF

    -(x : %) == -$Rep x
    x : % + y : % == x +$Rep y
    x : % - y : % == x -$Rep y
    x : % = y : % == x =$Rep y
    n : I * x : % == x *$Rep (n::GF)




    representationType() == "normal"

    minimalPolynomial(a) ==
      if initmult? then initializeMult()
      setFieldInfo(multTable, traceAlpha)$INBFF
      minimalPolynomial(a::Rep)$INBFF

-- is x an element of the ground field GF ?
    inGroundField?(x) ==
      erg := true
      for i in 2..extdeg repeat
          not(x(i) =$GF x(1)) => erg := false
      erg

    x : % / y : % ==
      if initmult? then initializeMult()
      setFieldInfo(multTable, traceAlpha)$INBFF
      x::Rep /$INBFF y::Rep

    inv(a) ==
      if initmult? then initializeMult()
      setFieldInfo(multTable, traceAlpha)$INBFF
      inv(a::Rep)$INBFF

    norm(a, d) ==
      if initmult? then initializeMult()
      setFieldInfo(multTable, traceAlpha)$INBFF
      norm(a::Rep, d)$INBFF

    normalElement() == normalElement(extdeg)$INBFF

)abbrev domain FFNBX FiniteFieldNormalBasisExtension
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations:
++ Related Constructors: FiniteFieldNormalBasisExtensionByPolynomial,
++  FiniteFieldPolynomialPackage
++ Also See: FiniteFieldExtension, FiniteFieldCyclicGroupExtension
++ AMS Classifications:
++ Keywords: finite field, normal basis
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldNormalBasisExtensionByPolynomial(GF, n)  implements a
++  finite extension field of degree n over the ground field {\em GF}.
++  The elements are represented by coordinate vectors with respect
++  to a normal basis,
++  i.e. a basis consisting of the conjugates (q-powers) of an element, in
++  this case called normal element. This is chosen as a root of the extension
++  polynomial, created by {\em createNormalPoly} from
++  \spadtype{FiniteFieldPolynomialPackage}
FiniteFieldNormalBasisExtension(GF, extdeg) : _
  Exports == Implementation where
  GF    : FiniteFieldCategory                -- the ground field
  extdeg : PositiveInteger                    -- the extension degree
  NNI     ==> NonNegativeInteger
  FFF         ==> FiniteFieldFunctions(GF)
  TERM    ==> Record(value : GF, index : SingleInteger)
  Exports ==> FiniteAlgebraicExtensionField(GF)  with
    getMultiplicationTable : ()   -> Vector List TERM
      ++ getMultiplicationTable() returns the multiplication
      ++ table for the normal basis of the field.
      ++ This table is used to perform multiplications between field elements.
    getMultiplicationMatrix : ()  -> Matrix GF
      ++ getMultiplicationMatrix() returns the multiplication table in
      ++ form of a matrix.
    sizeMultiplication : ()   -> NNI
      ++ sizeMultiplication() returns the number of entries in the
      ++ multiplication table of the field. Note: the time of multiplication
      ++ of field elements depends on this size.

  Implementation ==> FiniteFieldNormalBasisExtensionByPolynomial(GF, _
                    createLowComplexityNormalBasis(extdeg)$FFF)

)abbrev domain FFNB FiniteFieldNormalBasis
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations:
++ Related Constructors: FiniteFieldNormalBasisExtensionByPolynomial,
++  FiniteFieldPolynomialPackage
++ Also See: FiniteField, FiniteFieldCyclicGroup
++ AMS Classifications:
++ Keywords: finite field, normal basis
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++   AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldNormalBasis(p, n) implements a
++  finite extension field of degree n over the prime field with p elements.
++  The elements are represented by coordinate vectors with respect to
++  a normal basis,
++  i.e. a basis consisting of the conjugates (q-powers) of an element, in
++  this case called normal element.
++  This is chosen as a root of the extension polynomial
++  created by \spadfunFrom{createNormalPoly}{FiniteFieldPolynomialPackage}.
FiniteFieldNormalBasis(p, extdeg) : _
  Exports == Implementation where
  p : PositiveInteger
  extdeg : PositiveInteger                    -- the extension degree
  NNI     ==> NonNegativeInteger
  FFF  ==> FiniteFieldFunctions(PrimeField(p))
  TERM    ==> Record(value : PrimeField(p), index : SingleInteger)
  Exports ==> FiniteAlgebraicExtensionField(PrimeField(p))  with
    getMultiplicationTable : ()   -> Vector List TERM
      ++ getMultiplicationTable() returns the multiplication
      ++ table for the normal basis of the field.
      ++ This table is used to perform multiplications between field elements.
    getMultiplicationMatrix : ()  -> Matrix PrimeField(p)
      ++ getMultiplicationMatrix() returns the multiplication table in
      ++ form of a matrix.
    sizeMultiplication : ()   -> NNI
      ++ sizeMultiplication() returns the number of entries in the
      ++ multiplication table of the field. Note: The time of multiplication
      ++ of field elements depends on this size.

  Implementation ==> FiniteFieldNormalBasisExtensionByPolynomial(PrimeField(p), _
                    createLowComplexityNormalBasis(extdeg)$FFF)

)abbrev package FFHOM FiniteFieldHomomorphisms
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Basic Operations:
++ Related Constructors: FiniteFieldCategory, FiniteAlgebraicExtensionField
++ Also See:
++ AMS Classifications:
++ Keywords: finite field, homomorphism, isomorphism
++ References:
++  R.Lidl, H.Niederreiter: Finite Field, Encyclopedia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  J. Grabmeier, A. Scheerhorn: Finite Fields in AXIOM.
++  AXIOM Technical Report Series, ATR/5 NP2522.
++ Description:
++  FiniteFieldHomomorphisms(F1, GF, F2) exports coercion functions of
++  elements between the fields {\em F1} and {\em F2}, which both must be
++  finite simple algebraic extensions of the finite ground field {\em GF}.
FiniteFieldHomomorphisms(F1, GF, F2) : Exports == Implementation where
  F1 : FiniteAlgebraicExtensionField(GF)
  GF : FiniteFieldCategory
  F2 : FiniteAlgebraicExtensionField(GF)
 -- the homomorphism can only convert elements w.r.t. the last extension .
  -- Adding a function 'groundField()' which returns the ground field of GF
  -- as a variable of type FiniteFieldCategory in the new compiler, one
  -- could build up 'convert' recursively to get an homomorphism w.r.t
  -- the whole extension.

  I   ==> Integer
  NNI ==> NonNegativeInteger
  SI  ==> SingleInteger
  PI  ==> PositiveInteger
  SUP ==> SparseUnivariatePolynomial
  M   ==> Matrix GF
  FFP ==> FiniteFieldExtensionByPolynomial
  FFPOL2 ==> FiniteFieldPolynomialPackage2
  FFPOLY ==> FiniteFieldPolynomialPackage
  OUT ==> OutputForm

  Exports ==> with

    coerce : F1  ->  F2
      ++ coerce(x) is the homomorphic image of x from
      ++ {\em F1} in {\em F2}. Thus {\em coerce} is a
      ++ field homomorphism between the fields extensions
      ++ {\em F1} and {\em F2} both over ground field {\em GF}
      ++ (the second argument to the package).
      ++ Error: if the extension degree of {\em F1} doesn't divide
      ++ the extension degree of {\em F2}.
      ++ Note that the other coercion function in the
      ++ \spadtype{FiniteFieldHomomorphisms} is a left inverse.

    coerce : F2  ->  F1
      ++ coerce(x) is the homomorphic image of x from
      ++ {\em F2} in {\em F1}, where {\em coerce} is a
      ++ field homomorphism between the fields extensions
      ++ {\em F2} and {\em F1} both over ground field {\em GF}
      ++ (the second argument to the package).
      ++ Error: if the extension degree of {\em F2} doesn't divide
      ++ the extension degree of {\em F1}.
      ++ Note that the other coercion function in the
      ++ \spadtype{FiniteFieldHomomorphisms} is a left inverse.
    -- coerce(coerce(x: F1)@F2)@F1 = x and coerce(coerce(y: F2)@F1)@F2 = y

  Implementation ==> add

-- global variables ===================================================

    degree1 : NNI := extensionDegree()$F1
    degree2 : NNI := extensionDegree()$F2
    -- the degrees of the last extension

    -- a necessary condition for the one field being an subfield of
    -- the other one is, that the respective extension degrees are
    -- multiples
    if max(degree1, degree2) rem min(degree1, degree2) ~= 0 then
      error "FFHOM: one extension degree must divide the other one"

    conMat1to2 : M := zero(degree2, degree1)$M
    -- conversion matrix for the conversion direction F1 -> F2
    conMat2to1 : M := zero(degree1, degree2)$M
    -- conversion matrix for the conversion direction F2 -> F1

    repType1 := representationType()$F1
    repType2 := representationType()$F2
    -- the representation types of the fields

    init? : Boolean := true
    -- gets false after initialization

    defPol1 := definingPolynomial()$F1
    defPol2 := definingPolynomial()$F2
    -- the defining polynomials of the fields


-- functions ==========================================================


    compare : (SUP GF, SUP GF) -> Boolean
    -- compares two polynomials

    convertWRTsameDefPol12 : F1  ->  F2
    convertWRTsameDefPol21 : F2  ->  F1
    -- homomorphism if the last extension of F1 and F2 was build up
    -- using the same defining polynomials

    convertWRTdifferentDefPol12 : F1  ->  F2
    convertWRTdifferentDefPol21 : F2  ->  F1
    -- homomorphism if the last extension of F1 and F2 was build up
    -- with different defining polynomials

    initialize : () -> Void
    -- computes the conversion matrices

    compare(g : (SUP GF), f : (SUP GF)) ==
      degree(f)$(SUP GF)  >$NNI degree(g)$(SUP GF) => true
      degree(f)$(SUP GF) <$NNI degree(g)$(SUP GF) => false
      equal : Integer := 0
      for i in degree(f)$(SUP GF)..0 by -1 while equal = 0 repeat
        not zero?(coefficient(f, i)$(SUP GF))$GF and _
             zero?(coefficient(g, i)$(SUP GF))$GF => equal := 1
        not zero?(coefficient(g, i)$(SUP GF))$GF and _
             zero?(coefficient(f, i)$(SUP GF))$GF => equal := (-1)
        (f1 := lookup(coefficient(f, i)$(SUP GF))$GF) >$PositiveInteger _
         (g1 := lookup(coefficient(g, i)$(SUP GF))$GF) =>  equal := 1
        f1 <$PositiveInteger g1 => equal := (-1)
      equal = 1 => true
      false

    initialize() ==
      -- 1) in the case of equal def. polynomials initialize is called only
      --  if one of the rep. types is "normal" and the other one is "polynomial"
      --  we have to compute the basis change matrix 'mat', which i-th
      --  column are the coordinates of a^(q^i), the i-th component of
      --  the normal basis ('a' the root of the def. polynomial and q the
      --  size of the ground field)
      defPol1 =$(SUP GF) defPol2 =>
        -- new code using reducedQPowers
        mat := zero(degree1, degree1)$M
        arr := reducedQPowers(defPol1)$FFPOLY(GF)
        for i in 1..degree1 repeat
          setColumn!(mat, i, vectorise(arr.(i-1), degree1)$SUP(GF))$M
          -- old code
          -- here one of the representation types must be "normal"
          --a := basis()$FFP(GF, defPol1).2  -- the root of the def. polynomial
          --setColumn!(mat, 1, coordinates(a)$FFP(GF, defPol1))$M
          --for i in 2..degree1 repeat
          --  a := a ^$FFP(GF, defPol1) size()$GF
          --  setColumn!(mat, i, coordinates(a)$FFP(GF, defPol1))$M
          --for the direction "normal" -> "polynomial" we have to multiply the
          -- coordinate vector of an element of the normal basis field with
          -- the matrix 'mat'. In this case 'mat' is the correct conversion
          -- matrix for the conversion of F1 to F2, its inverse the correct
          -- inversion matrix for the conversion of F2 to F1
        repType1 = "normal" =>  -- repType2 = "polynomial"
          conMat1to2 := copy(mat)
          conMat2to1 := copy(inverse(mat)$M :: M)
          --we finish the function for one case, hence reset initialization flag
          init? := false
          void()$Void
          -- print("'normal' <=> 'polynomial' matrices initialized"::OUT)
        -- in the other case we have to change the matrices
        -- repType2 = "normal" and repType1 = "polynomial"
        conMat2to1 := copy(mat)
        conMat1to2 := copy(inverse(mat)$M :: M)
        -- print("'normal' <=> 'polynomial' matrices initialized"::OUT)
        --we finish the function for one case, hence reset initialization flag
        init? := false
        void()$Void
      -- 2) in the case of different def. polynomials we have to order the
      --    fields to get the same isomorphism, if the package is called with
      --    the fields F1 and F2 swapped.
      dPbig := defPol2
      rTbig := repType2
      dPsmall := defPol1
      rTsmall := repType1
      degbig := degree2
      degsmall := degree1
      if compare(defPol2, defPol1) then
        degsmall := degree2
        degbig := degree1
        dPbig := defPol1
        rTbig := repType1
        dPsmall := defPol2
        rTsmall := repType2
      -- 3) in every case we need a conversion between the polynomial
      --  represented fields. Therefore we compute 'root' as a root of the
      --  'smaller' def. polynomial in the 'bigger' field.
      --  We compute the matrix 'matsb', which i-th column are the coordinates
      --  of the (i-1)-th power of root, i=1..degsmall. Multiplying a
      --  coordinate vector of an element of the 'smaller' field by this
      --  matrix, we got the coordinates of the corresponding element in the
      --  'bigger' field.
      -- compute the root of dPsmall in the 'big' field
      root := rootOfIrreduciblePoly(dPsmall)$FFPOL2(FFP(GF, dPbig), GF)
      -- set up matrix for polynomial conversion
      matsb := zero(degbig, degsmall)$M
      qsetelt!(matsb, 1, 1, 1$GF)$M
      a := root
      for i in 2..degsmall repeat
        setColumn!(matsb, i, coordinates(a)$FFP(GF, dPbig))$M
        a := a *$FFP(GF, dPbig) root
      --  the conversion from 'big' to 'small': we can't invert matsb
      --  directly, because it has degbig rows and degsmall columns and
      --  may be no square matrix. Therefore we construct a square matrix
      --  mat from degsmall linear independent rows of matsb and invert it.
      --  Now we get the conversion matrix 'matbs' for the conversion from
      --  'big' to 'small' by putting the columns of mat at the indices
      --  of the linear independent rows of matsb to columns of matbs.
      ra : I := 1   -- the rank
      mat : M := transpose(row(matsb, 1))$M -- has already rank 1
      rowind : I := 2
      iVec : Vector I := new(degsmall, 1$I)$(Vector I)
      while ra < degsmall repeat
        if rank(vertConcat(mat, transpose(row(matsb, rowind))$M)$M)$M > ra then
          mat := vertConcat(mat, transpose(row(matsb, rowind))$M)$M
          ra := ra+1
          iVec.ra := rowind
        rowind := rowind + 1
      mat := inverse(mat)$M :: M
      matbs := zero(degsmall, degbig)$M
      for i in 1..degsmall repeat
        setColumn!(matbs, iVec.i, column(mat, i)$M)$M
      -- print(matsb::OUT)
      -- print(matbs::OUT)
      -- 4) if the 'bigger' field is "normal" we have to compose the
      --  polynomial conversion with a conversion from polynomial to normal
      --  between the FFP(GF, dPbig) and FFNBP(GF, dPbig) the 'bigger'
      --  field. Therefore we compute a conversion matrix 'mat' as in 1)
      --  Multiplying with the inverse of 'mat' yields the desired
      --  conversion from polynomial to normal. Multiplying this matrix by
      --  the above computed 'matsb' we got the matrix for converting form
      --  'small polynomial' to 'big normal'.
      -- set up matrix 'mat' for polynomial to normal
      if rTbig = "normal" then
        arr := reducedQPowers(dPbig)$FFPOLY(GF)
        mat := zero(degbig, degbig)$M
        for i in 1..degbig repeat
          setColumn!(mat, i, vectorise(arr.(i-1), degbig)$SUP(GF))$M
        -- old code
        --a := basis()$FFP(GF, dPbig).2  -- the root of the def.Polynomial
        --setColumn!(mat, 1, coordinates(a)$FFP(GF, dPbig))$M
        --for i in 2..degbig repeat
        --  a := a ^$FFP(GF, dPbig) size()$GF
        --  setColumn!(mat, i, coordinates(a)$FFP(GF, dPbig))$M
        -- print(inverse(mat)$M::OUT)
        matsb := (inverse(mat)$M :: M) * matsb
        -- print("inv *.."::OUT)
        matbs := matbs * mat
        -- 5) if the 'smaller' field is "normal" we have first to convert
        --    from 'small normal' to 'small polynomial', that is from
        --    FFNBP(GF, dPsmall) to FFP(GF, dPsmall). Therefore we compute a
        --    conversion matrix 'mat' as in 1). Multiplying with  'mat'
        --    yields the desired conversion from normal to polynomial.
        --    Multiplying the above computed 'matsb' with 'mat' we got the
        --    matrix for converting form 'small normal' to 'big normal'.
      -- set up matrix 'mat' for normal to polynomial
      if rTsmall = "normal" then
        arr := reducedQPowers(dPsmall)$FFPOLY(GF)
        mat := zero(degsmall, degsmall)$M
        for i in 1..degsmall repeat
          setColumn!(mat, i, vectorise(arr.(i-1), degsmall)$SUP(GF))$M
      -- old code
      --b: FFP(GF, dPsmall) := basis()$FFP(GF, dPsmall).2
      --setColumn!(mat, 1, coordinates(b)$FFP(GF, dPsmall))$M
      --for i in 2..degsmall repeat
      --  b := b ^$FFP(GF, dPsmall) size()$GF
      --  setColumn!(mat, i, coordinates(b)$FFP(GF, dPsmall))$M
        -- print(mat::OUT)
        matsb := matsb * mat
        matbs := (inverse(mat) :: M) * matbs
      -- now 'matsb' is the correct conversion matrix for 'small' to 'big'
      -- and 'matbs' the correct one for 'big' to 'small'.
      -- depending on the above ordering the conversion matrices are
      -- initialized
      dPbig =$(SUP GF) defPol2 =>
        conMat1to2 := matsb
        conMat2to1 := matbs
        -- print(conMat1to2::OUT)
        -- print(conMat2to1::OUT)
        -- print("conversion matrices initialized"::OUT)
        --we finish the function for one case, hence reset initialization flag
        init? := false
        void()$Void
      conMat1to2 := matbs
      conMat2to1 := matsb
      -- print(conMat1to2::OUT)
      -- print(conMat2to1::OUT)
      -- print("conversion matrices initialized"::OUT)
      --we finish the function for one case, hence reset initialization flag
      init? := false
      void()$Void


    coerce(x : F1) ==
      inGroundField?(x)$F1 => retract(x)$F1 :: F2
      -- if x is already in GF then we can use a simple coercion
      defPol1 =$(SUP GF) defPol2 => convertWRTsameDefPol12(x)
      convertWRTdifferentDefPol12(x)

    convertWRTsameDefPol12(x : F1)  ==
      repType1 = repType2 => x pretend F2
      -- same ground fields, same defining polynomials, same
      -- representation types --> F1 = F2, x is already in F2
      repType1 = "cyclic" =>
        x = 0$F1 => 0$F2
      -- the SI corresponding to the cyclic representation is the exponent of
      -- the primitiveElement, therefore we exponentiate the primitiveElement
      -- of F2 by it.
        primitiveElement()$F2 ^$F2 (x pretend SI)
      repType2 = "cyclic" =>
        x = 0$F1 => 0$F2
      -- to get the exponent, we have to take the discrete logarithm of the
      -- element in the given field.
        (qconvert(discreteLog(x)$F1)@SI) pretend F2
      -- here one of the representation types is "normal"
      if init? then initialize()
      -- here a conversion matrix is necessary, (see initialize())
      represents(conMat1to2 *$(Matrix GF) coordinates(x)$F1)$F2

    convertWRTdifferentDefPol12(x : F1) ==
      if init? then initialize()
      -- if we want to convert into a 'smaller' field, we have to test,
      -- whether the element is in the subfield of the 'bigger' field, which
      -- corresponds to the 'smaller' field
      if degree1 > degree2 then
        if positiveRemainder(degree2, degree(x)$F1) ~= 0 then
          error "coerce: element doesn't belong to smaller field"
      represents(conMat1to2 *$(Matrix GF) coordinates(x)$F1)$F2

-- the three functions below equal the three functions above up to
-- '1' exchanged by '2' in all domain and variable names


    coerce(x : F2) ==
      inGroundField?(x)$F2 => retract(x)$F2 :: F1
      -- if x is already in GF then we can use a simple coercion
      defPol1 =$(SUP GF) defPol2 => convertWRTsameDefPol21(x)
      convertWRTdifferentDefPol21(x)

    convertWRTsameDefPol21(x : F2)  ==
      repType1 = repType2 => x pretend F1
      -- same ground fields, same defining polynomials,
      -- same representation types --> F1 = F2, that is:
      -- x is already in F1
      repType2 = "cyclic" =>
        x = 0$F2 => 0$F1
        primitiveElement()$F1 ^$F1 (x pretend SI)
      repType1 = "cyclic" =>
        x = 0$F2 => 0$F1
        (qconvert(discreteLog(x)$F2)@SI) pretend F1
      -- here one of the representation types is "normal"
      if init? then initialize()
      represents(conMat2to1 *$(Matrix GF) coordinates(x)$F2)$F1

    convertWRTdifferentDefPol21(x : F2) ==
      if init? then initialize()
      if degree2 > degree1 then
        if positiveRemainder(degree1, degree(x)$F2) ~= 0 then
          error "coerce: element doesn't belong to smaller field"
      represents(conMat2to1 *$(Matrix GF) coordinates(x)$F2)$F1

)abbrev package FFF FiniteFieldFunctions
++ Author: J. Grabmeier, A. Scheerhorn
++ Date Created: 21 March 1991
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ References:
++  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, ATR/5 NP2522.
++ Description:
++  FiniteFieldFunctions(GF) is a package with functions
++  concerning finite extension fields of the finite ground field {\em GF},
++  e.g. Zech logarithms.
++ Keywords: finite field, Zech logarithm, Jacobi logarithm, normal basis

FiniteFieldFunctions(GF) : Exports == Implementation where
  GF  : FiniteFieldCategory  -- the ground field

  PI   ==> PositiveInteger
  NNI  ==> NonNegativeInteger
  I    ==> Integer
  SI   ==> SingleInteger
  SUP  ==> SparseUnivariatePolynomial GF
  V    ==> Vector
  M    ==> Matrix
  L    ==> List
  OUT  ==> OutputForm
  SAE  ==> SimpleAlgebraicExtension
  ARR  ==> PrimitiveArray(SI)
  TERM ==> Record(value : GF, index : SI)
  MM   ==> ModMonic(GF, SUP)
  PF   ==> PrimeField

  Exports ==> with

      createZechTable : SUP -> ARR
        ++ createZechTable(f) generates a Zech logarithm table for the cyclic
        ++ group representation of a extension of the ground field by the
        ++ primitive polynomial {\em f(x)}, i.e. \spad{Z(i)},
        ++ defined by {\em x^Z(i) = 1+x^i} is stored at index i.
        ++ This is needed in particular
        ++ to perform addition of field elements in finite fields represented
        ++ in this way. See \spadtype{FFCGP}, \spadtype{FFCGX}.
      createMultiplicationTable : SUP -> V L TERM
        ++ createMultiplicationTable(f) generates a multiplication
        ++ table for the normal basis of the field extension determined
        ++ by {\em f}. This is needed to perform multiplications
        ++ between elements represented as coordinate vectors to this basis.
        ++ See \spadtype{FFNBP}, \spadtype{FFNBX}.
      createMultiplicationMatrix : V L TERM -> M GF
        ++ createMultiplicationMatrix(m) forms the multiplication table
        ++ {\em m} into a matrix over the ground field.
        -- only useful for the user to visualise the multiplication table
        -- in a nice form
      sizeMultiplication : V L TERM -> NNI
        ++ sizeMultiplication(m) returns the number of entries
        ++ of the multiplication table {\em m}.
        -- the time of the multiplication of field elements depends
        -- on this size
      createLowComplexityTable : PI -> Union(Vector List TERM,"failed")
        ++ createLowComplexityTable(n) tries to find
        ++ a low complexity normal basis of degree {\em n} over {\em GF}
        ++ and returns its multiplication matrix
        ++ Fails, if it does not find a low complexity basis
      createLowComplexityNormalBasis : PI -> Union(SUP, V L TERM)
        ++ createLowComplexityNormalBasis(n) tries to find a
        ++ a low complexity normal basis of degree {\em n} over {\em GF}
        ++ and returns its multiplication matrix
        ++ If no low complexity basis is found it calls
        ++ \spadfunFrom{createNormalPoly}{FiniteFieldPolynomialPackage}(n) to produce a normal
        ++ polynomial of degree {\em n} over {\em GF}

  Implementation ==> add


    createLowComplexityNormalBasis(n) ==
      (u := createLowComplexityTable(n)) case "failed" =>
        createNormalPoly(n)$FiniteFieldPolynomialPackage(GF)
      u

-- try to find a low complexity normal basis multiplication table
-- of the field of extension degree n
-- the algorithm is from:
-- Wassermann A., Konstruktion von Normalbasen,
-- Bayreuther Mathematische Schriften 31 (1989), 1-9.

    createLowComplexityTable(n) ==
      q := size()$GF
      -- this algorithm works only for prime fields
      p := characteristic()$GF
      -- search of a suitable parameter k
      k : NNI := 0
      for i in 1..n-1  while (k = 0) repeat
        if prime?(i*n+1) and not(p = (i*n+1)) then
          primitive?(q::PF(i*n+1))$PF(i*n+1) =>
              a : NNI := 1
              k := i
              t1 : PF(k*n+1) := (q::PF(k*n+1))^n
          gcd(n, a := discreteLog(q::PF(n*i+1))$PF(n*i+1))$I = 1 =>
              k := i
              t1 := primitiveElement()$PF(k*n+1)^n
      k = 0 => "failed"
      -- initialize some start values
      multmat : M PF(p) := zero(n, n)
      p1 := (k*n+1)
      pkn := q::PF(p1)
      t := t1 pretend PF(p1)
      if odd?(k) then
          jt : I := (n quo 2)+1
          vt : I := positiveRemainder((k-a) quo 2, k)+1
        else
          jt : I := 1
          vt : I := (k quo 2)+1
      -- compute matrix
      vec : Vector I := zero(qcoerce(p1))
      for x in 1..k repeat
        for l in 1..n repeat
          vec.((t^(x-1) * pkn^(l-1)) pretend Integer+1) := _
                                            positiveRemainder(l, p1)
      lvj : M I := zero(k::NNI, n)
      for v in 1..k repeat
        for j in 1..n repeat
          if (j ~= jt) or (v ~= vt) then
            help : PF(p1) := t^(v-1)*pkn^(j-1)+1@PF(p1)
            setelt!(lvj, v, j, vec.(help pretend I + 1))
      for j in 1..n repeat
        if j ~= jt then
          for v in 1..k repeat
            lvjh := elt(lvj, v, j)
            setelt!(multmat, j, lvjh, elt(multmat, j, lvjh) + 1)
      for i in 1..n repeat
        setelt!(multmat, jt, i, positiveRemainder(-k, p)::PF(p))
      for v in 1..k repeat
        if v ~= vt then
          lvjh := elt(lvj, v, jt)
          setelt!(multmat, jt, lvjh, elt(multmat, jt, lvjh) + 1)
      -- multmat
      m := nrows(multmat)$(M PF(p))
      multtable : V L TERM := new(m, []$(L TERM))$(V L TERM)
      for i in 1..m repeat
        l : L TERM := []
        v : V PF(p) := row(multmat, i)
        for j in (1::I)..(m::I) repeat
          if (v.j ~= 0) then
            -- take -v.j to get trace 1 instead of -1
            term : TERM := [(convert(-v.j)@I)::GF, qconvert(j - 2)@SI]$TERM
            l := cons(term, l)$(L TERM)
        qsetelt!(multtable, i, copy l)$(V L TERM)
      multtable

    sizeMultiplication(m) ==
      s : NNI := 0
      for i in 1..#m repeat
        s := s + #(m.i)
      s

    createMultiplicationTable(f : SUP) ==
      sizeGF : NNI := size()$GF -- the size of the ground field
      m : PI := qcoerce(degree(f)$SUP)
      m = 1 =>
        [[[-coefficient(f, 0)$SUP, (-1)::SI]$TERM]$(L TERM)]::(V L TERM)
      m1 : I := m-1
      -- initialize basis change matrices
      setPoly(f)$MM
      e := reduce(monomial(1, 1)$SUP)$MM ^ sizeGF
      w := 1$MM
      qpow : PrimitiveArray(MM) := new(m, 0)
      qpow(0) := 1$MM
      for i in 1..m1 repeat
          qpow(i) := (w := w*e)
      qexp : PrimitiveArray(MM) := new(m, 0)
      qexp(0) := reduce(monomial(1, 1)$SUP)$MM
      mat : M(GF) := zero(m, m)$M(GF)
      qsetelt!(mat, 2, 1, 1$GF)$M(GF)
      hm := qpow(1)
      qexp(1) := hm
      setColumn!(mat, 2, Vectorise(hm)$MM)$(M GF)
      for i in 2..m1 repeat
        g := 0$MM
        hs := lift(hm)
        while hs ~= 0 repeat
          g := g + leadingCoefficient(hs) * qpow(degree(hs))
          hs := reductum(hs)
        qexp(i) := hm := g
        setColumn!(mat, i+1, Vectorise(g)$MM)$M(GF)
      -- loop invariant: qexp.i = x^(q^i)
      mat1 := inverse(mat)$M(GF)
      mat1 = "failed" =>
        error "createMultiplicationTable: polynomial must be normal"
      mat := mat1 :: M(GF)
      -- initialize multiplication table
      multtable : V(L(TERM)) := new(m, []$L(TERM))$V(L(TERM))
      for i in 1..m repeat
        l : L(TERM) := []
        v : V(GF) := mat *$M(GF) Vectorise(qexp(i-1) *$MM qexp(0))$MM
        for j in (1::SI)..(m::SI) repeat
          if v(j) ~= 0$GF then
            term : TERM := [v(j), j-(2::SI)]$TERM
            l := cons(term, l)$L(TERM)
        qsetelt!(multtable, i, copy l)$V(L(TERM))
      multtable


    createZechTable(f : SUP) ==
      sizeGF : NNI := size()$GF -- the size of the ground field
      m := degree(f)$SUP::PI
      qm1 : SI := qconvert(sizeGF^m - 1)
      zechlog : ARR := new(((sizeGF ^ m + 1) quo 2)::NNI, -1::SI)$ARR
      helparr : ARR := new(sizeGF ^ m::NNI, 0$SI)$ARR
      primElement := reduce(monomial(1, 1)$SUP)$SAE(GF, SUP, f)
      a := primElement
      for i in 1..qm1-1 repeat
        helparr(lookup(a -$SAE(GF, SUP, f) 1$SAE(GF, SUP, f)_
           )$SAE(GF, SUP, f)) := i::SI
        a := a * primElement
      characteristic() = 2 =>
        a := primElement
        for i in 1..(qm1 quo 2) repeat
            zechlog(i) := helparr(lookup(a)$SAE(GF, SUP, f))
            a := a * primElement
        zechlog
      a := 1$SAE(GF, SUP, f)
      for i in 0..((qm1-2) quo 2) repeat
        zechlog(i) := helparr(lookup(a)$SAE(GF, SUP, f))
        a := a * primElement
      zechlog

    createMultiplicationMatrix(m) ==
      n : NNI := #m
      mat : M GF := zero(n, n)$(M GF)
      for i in 1..n repeat
        for t in m.i repeat
          qsetelt!(mat, i, t.index+2, t.value)
      mat

--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.