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

)abbrev package FFPOLY2 FiniteFieldPolynomialPackage2
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Date Last Updated:
++ 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, Encycoldia of Mathematics and
++  Its Applications, Vol. 20, Cambridge Univ. Press, 1983, ISBN 0 521 30240 4
++  AXIOM Technical Report Series, to appear.
++ 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 -> F
	++ coerce(x) \undocumented{}
      lookup: F -> PositiveInteger
	++ lookup(x) \undocumented{}
      basis: PositiveInteger -> Vector F
	++ basis(n) \undocumented{}
      Frobenius: F -> F
	++ 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
  MM  ==> ModMonic(GF,SUP GF)
  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 positive? degree(trModp)$(SUP F) 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(j pretend PI)$GF::F::(SUP F)))$(SUP F)
            -- make the gcd polynomial monic
            if not one? leadingCoefficient(h)$(SUP 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


