--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 FFF FiniteFieldFunctions
++ Author: J. Grabmeier, A. Scheerhorn
++ Date Created: 21 March 1991
++ Date Last Updated: 31 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
        ++ \axiomFunFrom{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::(V L TERM)

-- 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
      a:NNI
      t1: PF(k*n+1)       -- all that matters is the syntax of the type
      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 := 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(p1 pretend NNI)
      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,nil()$(L TERM))$(V L TERM)
      for i in 1..m repeat
        l:L TERM:=nil()$(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,(j-2) pretend 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:=degree(f)$SUP pretend PI
      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)
      -- qpow.i = x**(i*q)
      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)
      h:=qpow.1
      qexp.1:=h
      setColumn!(mat,2,Vectorise(h)$MM)$(M GF)
      for i in 2..m1 repeat
        g:=0$MM
        while h ~= 0 repeat
          g:=g + leadingCoefficient(h) * qpow.degree(h)$MM
          h:=reductum(h)$MM
        qexp.i:=g
        setColumn!(mat,i+1,Vectorise(h:=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,nil()$(L TERM))$(V L TERM)
      for i in 1..m repeat
        l:L TERM:=nil()$(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:=(sizeGF ** m -1) pretend SI
      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$GF = 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

