axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / INBFF.spad

)abbrev package INBFF InnerNormalBasisFieldFunctions
++ Authors: J.Grabmeier, A.Scheerhorn
++ Date Created: 26.03.1991
++ Date Last Updated: 31 March 1991
++ References:
++ Grab92 Finite Fields in Axiom
++ Lidl83 Finite Field, Encyclopedia of Mathematics and Its Applications
++ Stin90 Some observations on parallel Algorithms for fast exponentiation 
++        in GF(2^n)
++ Itoh88 A fast algorithm for computing multiplicative inverses
++        in GF(2^m) using normal bases
++ Description:
++ InnerNormalBasisFieldFunctions(GF) (unexposed):
++ This package has functions used by
++ every normal basis finite field extension domain.

InnerNormalBasisFieldFunctions(GF) : SIG == CODE 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)
  MM   ==> ModMonic(GF,SUP GF)

  SIG ==> 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 that 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 \axiomFunFrom{**}{DivisionRing}

    "*" :(VGF,VGF) -> VGF
      ++ x*y \undocumented{}
      ++ See \axiomFunFrom{*}{SemiGroup}

    "/" :(VGF,VGF) -> VGF
      ++ x/y \undocumented{}
      ++ See \axiomFunFrom{/}{Field}

    norm :(VGF,PI) -> VGF
      ++ norm(x,n) \undocumented{}
      ++ See \axiomFunFrom{norm}{FiniteAlgebraicExtensionField}

    trace :(VGF,PI) -> VGF
      ++ trace(x,n) \undocumented{}
      ++ See \axiomFunFrom{trace}{FiniteAlgebraicExtensionField}

    inv : VGF -> VGF
      ++ inv x \undocumented{}
      ++ See \axiomFunFrom{inv}{DivisionRing}

    lookup : VGF -> PI
      ++ lookup(x) \undocumented{}
      ++ See \axiomFunFrom{lookup}{Finite}

    normal? : VGF -> Boolean
      ++ normal?(x) \undocumented{}
      ++ See \axiomFunFrom{normal?}{FiniteAlgebraicExtensionField}

    basis : PI -> V VGF
      ++ basis(n) \undocumented{}
      ++ See \axiomFunFrom{basis}{FiniteAlgebraicExtensionField}

    normalElement : PI -> VGF
      ++ normalElement(n) \undocumented{}
      ++ See \axiomFunFrom{normalElement}{FiniteAlgebraicExtensionField}

    minimalPolynomial : VGF -> SUP GF
      ++ minimalPolynomial(x) \undocumented{}
      ++ See \axiomFunFrom{minimalPolynomial}{FiniteAlgebraicExtensionField}

  CODE ==> add

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

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

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

    trGen:GF:=1$GF
    -- controls the imbedding 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
      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 =>
          repSq(a,e)
        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 ((k-1)pretend NNI) +$I_
              ((logqe +$SI k -$SI 1$SI) quo$SI k)::I -$I 2
      (10*mult) >= (15 * (e1::I)) =>
        repSq(a,e)
      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,(e quo 2) @ NNI)
      repSq(a*a,(e quo 2) @ NNI)

    -- 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**(k pretend NNI)
      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:=nil()
      ex:I:=e
      while not(ex = 0) repeat
        l:=concat(l,positiveRemainder(ex,qk) pretend NNI)
        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:GF:=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 devide 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

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

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