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

)abbrev domain SAE SimpleAlgebraicExtension
++ Author: Barry Trager, Manuel Bronstein, Clifton Williamson
++ Date Created: 1986
++ Date Last Updated: 9 May 1994
++ References:
++ Grab92 Finite Fields in Axiom
++ Description:
++ Algebraic extension of a ring by a single polynomial.
++ Domain which represents simple algebraic extensions of arbitrary
++ rings. The first argument to the domain, R, is the underlying ring,
++ the second argument is a domain of univariate polynomials over K,
++ while the last argument specifies the defining minimal polynomial.
++ The elements of the domain are canonically represented as polynomials
++ of degree less than that of the minimal polynomial with coefficients
++ in R. The second argument is both the type of the third argument and
++ the underlying representation used by \spadtype{SAE} itself.

SimpleAlgebraicExtension(R,UP,M) : SIG == CODE where
  R : CommutativeRing
  UP : UnivariatePolynomialCategory R
  M : UP

  SIG ==> MonogenicAlgebra(R, UP)

  CODE ==> add

    --sqFr(pb): FactorS(Poly) from UnivPolySquareFree(Poly)

    --degree(M) > 0 and M must be monic if R is not a field.
    if (r := recip leadingCoefficient M) case "failed" then
                                    error "Modulus cannot be made monic"
    Rep := UP
    x,y :$
    c: R

    mkDisc   : Boolean -> Void

    mkDiscMat: Boolean -> Void

    M   := r::R * M

    d   := degree M

    d1  := subtractIfCan(d,1)::NonNegativeInteger

    discmat:Matrix(R) := zero(d, d)

    nodiscmat?:Reference(Boolean) := ref true

    disc:Reference(R) := ref 0

    nodisc?:Reference(Boolean) := ref true

    bsis := [monomial(1, i)$Rep for i in 0..d1]$Vector(Rep)

    if R has Finite then

         size == size$R ** d

         random == represents([random()$R for i in 0..d1])

    0 == 0$Rep

    1 == 1$Rep

    c * x == c *$Rep x

    n:Integer * x == n *$Rep x

    coerce(n:Integer):$   == coerce(n)$Rep

    coerce(c) == monomial(c,0)$Rep

    coerce(x):OutputForm == coerce(x)$Rep

    lift(x) == x pretend Rep

    reduce(p:UP):$ == (monicDivide(p,M)$Rep).remainder

    x = y == x =$Rep y

    x + y == x +$Rep y

    - x == -$Rep x

    x * y == reduce((x *$Rep y) pretend UP)

    coordinates(x) == [coefficient(lift(x),i) for i in 0..d1]

    represents(vect) == +/[monomial(vect.(i+1),i) for i in 0..d1]

    definingPolynomial()  == M

    characteristic()      == characteristic()$R

    rank()                == d::PositiveInteger

    basis()               == copy(bsis@Vector(Rep) pretend Vector($))

    if R has Field then

      minimalPolynomial x == squareFreePart characteristicPolynomial x

    if R has Field then

      coordinates(x:$,bas: Vector $) ==
        (m := inverse transpose coordinates bas) case "failed" =>
          error "coordinates: second argument must be a basis"
        (m :: Matrix R) * coordinates(x)

    else if R has IntegralDomain then

      coordinates(x:$,bas: Vector $) ==
        -- we work over the quotient field of R to invert a matrix
        qf := Fraction R
        imatqf := InnerMatrixQuotientFieldFunctions(R,Vector R,Vector R,_
                   Matrix R,qf,Vector qf,Vector qf,Matrix qf)
        mat := transpose coordinates bas
        (m := inverse(mat)$imatqf) case "failed" =>
          error "coordinates: second argument must be a basis"
        coordsQF: Vector qf := 
          map(y +-> y::qf,coordinates x)$VectorFunctions2(R,qf)
        -- here are the coordinates as elements of the quotient field:
        vecQF := (m :: Matrix qf) * coordsQF
        vec : Vector R := new(d,0)
        for i in 1..d repeat
          xi := qelt(vecQF,i)
          denom(xi) = 1 => qsetelt_!(vec,i,numer xi)
          error "coordinates: coordinates are not integral over ground ring"
        vec

    reducedSystem(m:Matrix $):Matrix(R) ==
      reducedSystem(map(lift, m)$MatrixCategoryFunctions2($, Vector $,
               Vector $, Matrix $, UP, Vector UP, Vector UP, Matrix UP))

    reducedSystem(m:Matrix $, v:Vector $):Record(mat:Matrix R,vec:Vector R) ==
      reducedSystem(map(lift, m)$MatrixCategoryFunctions2($, Vector $,
               Vector $, Matrix $, UP, Vector UP, Vector UP, Matrix UP),
                                    map(lift, v)$VectorFunctions2($, UP))

    discriminant() ==
      if nodisc?() then mkDisc false
      disc()

    mkDisc b ==
      nodisc?() := b
      disc() := discriminant M
      void

    traceMatrix() ==
      if nodiscmat?() then mkDiscMat false
      discmat

    mkDiscMat b ==
      nodiscmat?() := b
      mr := minRowIndex discmat; mc := minColIndex discmat
      for i in 0..d1 repeat
        for j in 0..d1 repeat
          qsetelt_!(discmat,mr + i,mc + j,trace reduce monomial(1,i + j))
      void

    trace x ==          --this could be coded perhaps more efficiently
      xn := x;  ans := coefficient(lift xn, 0)
      for n in 1..d1 repeat
        (xn := generator() * xn;  ans := coefficient(lift xn, n) + ans)
      ans

    if R has Finite then

       index k ==
         i:Integer := k rem size()
         p:Integer := size()$R
         ans:$ := 0
         for j in 0.. while i > 0 repeat
           h := i rem p
           -- index(p) = 0$R
           if h ^= 0 then
             -- here was a bug: "index" instead of
             -- "coerce", otherwise it wouldn't work for
             -- Rings R where "coerce: I-> R" is not surjective
             a := index(h :: PositiveInteger)$R
             ans := ans + reduce monomial(a, j)
           i := i quo p
         ans

       lookup(z : $) : PositiveInteger ==
         -- z = index lookup z, n = lookup index n
         -- the answer is merely the Horner evaluation of the
         -- representation with the size of R (as integers).
         zero?(z) => size()$$ pretend PositiveInteger
         p  :            Integer := size()$R
         co :            Integer := lookup(leadingCoefficient z)$R
         n  : NonNegativeInteger := degree(z)
         while not zero?(z := reductum z) repeat
          co := co * p ** ((n - (n := degree z)) pretend
            NonNegativeInteger) + lookup(leadingCoefficient z)$R
         n = 0 => co pretend PositiveInteger
         (co * p ** n) pretend PositiveInteger