path: root/src/algebra/intclos.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra intclos.spad}
\author{Victor Miller, Barry Trager, Clifton Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{package TRIMAT TriangularMatrixOperations}
<<package TRIMAT TriangularMatrixOperations>>=
)abbrev package TRIMAT TriangularMatrixOperations
++ Fraction free inverses of triangular matrices
++ Author: Victor Miller
++ Date Created:
++ Date Last Updated: 24 Jul 1990
++ Keywords:
++ Examples:
++ References:
++ Description:
++ This package provides functions that compute "fraction-free"
++ inverses of upper and lower triangular matrices over a integral
++ domain.  By "fraction-free inverses" we mean the following:
++ given a matrix B with entries in R and an element d of R such that
++ d * inv(B) also has entries in R, we return d * inv(B).  Thus,
++ it is not necessary to pass to the quotient field in any of our
++ computations.


TriangularMatrixOperations(R,Row,Col,M): Exports == Implementation where
  R   : IntegralDomain
  Row : FiniteLinearAggregate R
  Col : FiniteLinearAggregate R
  M   : MatrixCategory(R,Row,Col)

  Exports ==> with

    UpTriBddDenomInv: (M,R) -> M
      ++ UpTriBddDenomInv(B,d) returns M, where
      ++ B is a non-singular upper triangular matrix and d is an
      ++ element of R such that \spad{M = d * inv(B)} has entries in R.
    LowTriBddDenomInv:(M,R) -> M
      ++ LowTriBddDenomInv(B,d) returns M, where
      ++ B is a non-singular lower triangular matrix and d is an
      ++ element of R such that \spad{M = d * inv(B)} has entries in R.

  Implementation ==> add

    UpTriBddDenomInv(A,denom) ==
      AI := zero(nrows A, nrows A)$M
      offset := minColIndex AI - minRowIndex AI
      for i in minRowIndex AI .. maxRowIndex AI
        for j in minColIndex AI .. maxColIndex AI repeat
          qsetelt_!(AI,i,j,(denom exquo qelt(A,i,j))::R)
      for i in minRowIndex AI .. maxRowIndex AI repeat
        for j in offset + i + 1 .. maxColIndex AI repeat
          qsetelt_!(AI,i,j, - (((+/[qelt(AI,i,k) * qelt(A,k-offset,j)
                                   for k in i+offset..(j-1)])
                                     exquo qelt(A, j-offset, j))::R))
      AI

    LowTriBddDenomInv(A, denom) ==
      AI := zero(nrows A, nrows A)$M
      offset := minColIndex AI - minRowIndex AI
      for i in minRowIndex AI .. maxRowIndex AI
        for j in minColIndex AI .. maxColIndex AI repeat
          qsetelt_!(AI,i,j,(denom exquo qelt(A,i,j))::R)
      for i in minColIndex AI .. maxColIndex AI repeat
        for j in i - offset + 1 .. maxRowIndex AI repeat
          qsetelt_!(AI,j,i, - (((+/[qelt(A,j,k+offset) * qelt(AI,k,i)
                                    for k in i-offset..(j-1)])
                                      exquo qelt(A, j, j+offset))::R))
      AI

@
\section{package IBATOOL IntegralBasisTools}
<<package IBATOOL IntegralBasisTools>>=
)abbrev package IBATOOL IntegralBasisTools
++ Functions common to both integral basis packages
++ Author: Victor Miller, Barry Trager, Clifton Williamson
++ Date Created: 11 April 1990
++ Date Last Updated: 20 September 1994
++ Keywords: integral basis, function field, number field
++ Examples:
++ References:
++ Description:
++ This package contains functions used in the packages
++ FunctionFieldIntegralBasis and NumberFieldIntegralBasis.

IntegralBasisTools(R,UP,F): Exports == Implementation where
  R  : EuclideanDomain with 
	squareFree: $ -> Factored $
		++ squareFree(x) returns a square-free factorisation of x 
  UP : UnivariatePolynomialCategory R
  F  : FramedAlgebra(R,UP)
  Mat ==> Matrix R
  NNI ==> NonNegativeInteger
  Ans ==> Record(basis: Mat, basisDen: R, basisInv:Mat)

  Exports ==> with

    diagonalProduct: Mat -> R
      ++ diagonalProduct(m) returns the product of the elements on the
      ++ diagonal of the matrix m
    matrixGcd: (Mat,R,NNI) -> R
      ++ matrixGcd(mat,sing,n) is \spad{gcd(sing,g)} where \spad{g} is the
      ++ gcd of the entries of the \spad{n}-by-\spad{n} upper-triangular
      ++ matrix \spad{mat}.
    divideIfCan_!: (Matrix R,Matrix R,R,Integer) -> R
      ++ divideIfCan!(matrix,matrixOut,prime,n) attempts to divide the
      ++ entries of \spad{matrix} by \spad{prime} and store the result in
      ++ \spad{matrixOut}.  If it is successful, 1 is returned and if not,
      ++ \spad{prime} is returned.  Here both \spad{matrix} and
      ++ \spad{matrixOut} are \spad{n}-by-\spad{n} upper triangular matrices.
    leastPower: (NNI,NNI) -> NNI
      ++ leastPower(p,n) returns e, where e is the smallest integer
      ++ such that \spad{p **e >= n}
    idealiser: (Mat,Mat) -> Mat
      ++ idealiser(m1,m2) computes the order of an ideal defined by m1 and m2
    idealiser: (Mat,Mat,R) -> Mat
      ++ idealiser(m1,m2,d) computes the order of an ideal defined by m1 and m2
      ++ where d is the known part of the denominator
    idealiserMatrix: (Mat, Mat) -> Mat
      ++ idealiserMatrix(m1, m2) returns the matrix representing the linear
      ++ conditions on the Ring associatied with an ideal defined by m1 and m2.
    moduleSum: (Ans,Ans) -> Ans
      ++ moduleSum(m1,m2) returns the sum of two modules in the framed
      ++ algebra \spad{F}.  Each module \spad{mi} is represented as follows:
      ++ F is a framed algebra with R-module basis \spad{w1,w2,...,wn} and
      ++ \spad{mi} is a record \spad{[basis,basisDen,basisInv]}.  If
      ++ \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ a basis \spad{v1,...,vn} for \spad{mi} is given by
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of 'basis' contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.

  Implementation ==> add
    import ModularHermitianRowReduction(R)
    import TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)

    diagonalProduct m ==
      ans : R := 1
      for i in minRowIndex m .. maxRowIndex m
        for j in minColIndex m .. maxColIndex m repeat
          ans := ans * qelt(m, i, j)
      ans

    matrixGcd(mat,sing,n) ==
      -- note: 'matrix' is upper triangular;
      -- no need to do anything below the diagonal
      d := sing
      for i in 1..n repeat
        for j in i..n repeat
          if not zero?(mij := qelt(mat,i,j)) then d := gcd(d,mij)
          one? d => return d
      d

    divideIfCan_!(matrix,matrixOut,prime,n) ==
    -- note: both 'matrix' and 'matrixOut' will be upper triangular;
    -- no need to do anything below the diagonal
      for i in 1..n repeat
        for j in i..n repeat
          (a := (qelt(matrix,i,j) exquo prime)) case "failed" => return prime
          qsetelt_!(matrixOut,i,j,a :: R)
      1

    leastPower(p,n) ==
      -- efficiency is not an issue here
      e : NNI := 1; q := p
      while q < n repeat (e := e + 1; q := q * p)
      e

    idealiserMatrix(ideal,idealinv) ==
      -- computes the Order of the ideal
      n  := rank()$F
      bigm := zero(n * n,n)$Mat
      mr := minRowIndex bigm; mc := minColIndex bigm
      v := basis()$F
      for i in 0..n-1 repeat
        r := regularRepresentation qelt(v,i + minIndex v)
        m := ideal * r * idealinv
        for j in 0..n-1 repeat
          for k in 0..n-1 repeat
            bigm(j * n + k + mr,i + mc) := qelt(m,j + mr,k + mc)
      bigm

    idealiser(ideal,idealinv) ==
      bigm := idealiserMatrix(ideal, idealinv)
      transpose squareTop rowEch bigm

    idealiser(ideal,idealinv,denom) ==
      bigm := (idealiserMatrix(ideal, idealinv) exquo denom)::Mat
      transpose squareTop rowEchelon(bigm,denom)

    moduleSum(mod1,mod2) ==
      rb1 := mod1.basis; rbden1 := mod1.basisDen; rbinv1 := mod1.basisInv
      rb2 := mod2.basis; rbden2 := mod2.basisDen; rbinv2 := mod2.basisInv
      -- compatibility check: doesn't take much computation time
      (not square? rb1) or (not square? rbinv1) or (not square? rb2) _
        or (not square? rbinv2) =>
        error "moduleSum: matrices must be square"
      ((n := nrows rb1) ~= (nrows rbinv1)) or (n ~= (nrows rb2)) _
        or (n ~= (nrows rbinv2)) =>
        error "moduleSum: matrices of imcompatible dimensions"
      (zero? rbden1) or (zero? rbden2) =>
        error "moduleSum: denominator must be non-zero"
      den := lcm(rbden1,rbden2); c1 := den quo rbden1; c2 := den quo rbden2
      rb := squareTop rowEchelon(vertConcat(c1 * rb1,c2 * rb2),den)
      rbinv := UpTriBddDenomInv(rb,den)
      [rb,den,rbinv]

@
\section{package FFINTBAS FunctionFieldIntegralBasis}
<<package FFINTBAS FunctionFieldIntegralBasis>>=
)abbrev package FFINTBAS FunctionFieldIntegralBasis
++ Integral bases for function fields of dimension one
++ Author: Victor Miller
++ Date Created: 9 April 1990
++ Date Last Updated: 20 September 1994
++ Keywords:
++ Examples:
++ References:
++ Description:
++ In this package R is a Euclidean domain and F is a framed algebra
++ over R.  The package provides functions to compute the integral
++ closure of R in the quotient field of F.  It is assumed that
++ \spad{char(R/P) = char(R)} for any prime P of R.  A typical instance of
++ this is when \spad{R = K[x]} and F is a function field over R.


FunctionFieldIntegralBasis(R,UP,F): Exports == Implementation where
  R  : EuclideanDomain with 
	squareFree: $ -> Factored $
		++ squareFree(x) returns a square-free factorisation of x 
  UP : UnivariatePolynomialCategory R
  F  : FramedAlgebra(R,UP)

  I   ==> Integer
  Mat ==> Matrix R
  NNI ==> NonNegativeInteger

  Exports ==> with
    integralBasis : () -> Record(basis: Mat, basisDen: R, basisInv:Mat)
      ++ \spad{integralBasis()} returns a record
      ++ \spad{[basis,basisDen,basisInv]} containing information regarding
      ++ the integral closure of R in the quotient field of F, where
      ++ F is a framed algebra with R-module basis \spad{w1,w2,...,wn}.
      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
    localIntegralBasis : R -> Record(basis: Mat, basisDen: R, basisInv:Mat)
      ++ \spad{integralBasis(p)} returns a record
      ++ \spad{[basis,basisDen,basisInv]} containing information regarding
      ++ the local integral closure of R at the prime \spad{p} in the quotient
      ++ field of F, where F is a framed algebra with R-module basis
      ++ \spad{w1,w2,...,wn}.
      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the local integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.

  Implementation ==> add
    import IntegralBasisTools(R, UP, F)
    import ModularHermitianRowReduction(R)
    import TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)

    squaredFactors: R -> R
    squaredFactors px ==
           */[(if ffe.exponent > 1 then ffe.factor else 1$R)
                for ffe in factors squareFree px]

    iIntegralBasis: (Mat,R,R) -> Record(basis: Mat, basisDen: R, basisInv:Mat)
    iIntegralBasis(tfm,disc,sing) ==
      -- tfm = trace matrix of current order
      n := rank()$F; tfm0 := copy tfm; disc0 := disc
      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
      -- rb    = basis matrix of current order
      -- rbinv = inverse basis matrix of current order
      -- these are wrt the original basis for F
      rbden : R := 1; index : R := 1; oldIndex : R := 1
      -- rbden = denominator for current basis matrix
      -- index = index of original order in current order
      not sizeLess?(1, sing) => [rb, rbden, rbinv]
      repeat
        -- compute the p-radical
        idinv := transpose squareTop rowEchelon(tfm, sing)
        -- [u1,..,un] are the coordinates of an element of the p-radical
        -- iff [u1,..,un] * idinv is in sing * R^n
        id := rowEchelon LowTriBddDenomInv(idinv, sing)
        -- id = basis matrix of the p-radical
        idinv := UpTriBddDenomInv(id, sing)
        -- id * idinv = sing * identity
        -- no need to check for inseparability in this case
        rbinv := idealiser(id * rb, rbinv * idinv, sing * rbden)
        index := diagonalProduct rbinv
        rb := rowEchelon LowTriBddDenomInv(rbinv, rbden * sing)
        g := matrixGcd(rb,sing,n)
        if sizeLess?(1,g) then rb := (rb exquo g) :: Mat
        rbden := rbden * (sing quo g)
        rbinv := UpTriBddDenomInv(rb, rbden)
        disc := disc0 quo (index * index)
        indexChange := index quo oldIndex; oldIndex := index
        sing := gcd(indexChange, squaredFactors disc)
        not sizeLess?(1, sing) => return [rb, rbden, rbinv]
        tfm := ((rb * tfm0 * transpose rb) exquo (rbden * rbden)) :: Mat

    integralBasis() ==
      n := rank()$F; p := characteristic$F
      (not zero? p) and (n >= p) =>
        error "integralBasis: possible wild ramification"
      tfm := traceMatrix()$F; disc := determinant tfm
      sing := squaredFactors disc    -- singularities of relative Spec
      iIntegralBasis(tfm,disc,sing)

    localIntegralBasis prime ==
      n := rank()$F; p := characteristic$F
      (not zero? p) and (n >= p) =>
        error "integralBasis: possible wild ramification"
      tfm := traceMatrix()$F; disc := determinant tfm
      (disc exquo (prime * prime)) case "failed" =>
        [scalarMatrix(n,1),1,scalarMatrix(n,1)]
      iIntegralBasis(tfm,disc,prime)

@
\section{package WFFINTBS WildFunctionFieldIntegralBasis}
<<package WFFINTBS WildFunctionFieldIntegralBasis>>=
)abbrev package WFFINTBS WildFunctionFieldIntegralBasis
++ Authors: Victor Miller, Clifton Williamson
++ Date Created: 24 July 1991
++ Date Last Updated: 20 September 1994
++ Basic Operations: integralBasis, localIntegralBasis
++ Related Domains: IntegralBasisTools(R,UP,F),
++   TriangularMatrixOperations(R,Vector R,Vector R,Matrix R)
++ Also See: FunctionFieldIntegralBasis, NumberFieldIntegralBasis
++ AMS Classifications:
++ Keywords: function field, integral basis
++ Examples:
++ References:
++ Description:
++ In this package K is a finite field, R is a ring of univariate
++ polynomials over K, and F is a framed algebra over R.  The package
++ provides a function to compute the integral closure of R in the quotient
++ field of F as well as a function to compute a "local integral basis"
++ at a specific prime.

WildFunctionFieldIntegralBasis(K,R,UP,F): Exports == Implementation where
  K  : FiniteFieldCategory
  --K  : Join(Field,Finite)
  R  : UnivariatePolynomialCategory K
  UP : UnivariatePolynomialCategory R
  F  : FramedAlgebra(R,UP)

  I       ==> Integer
  Mat     ==> Matrix R
  NNI     ==> NonNegativeInteger
  SAE     ==> SimpleAlgebraicExtension
  RResult ==> Record(basis: Mat, basisDen: R, basisInv:Mat)
  IResult ==> Record(basis: Mat, basisDen: R, basisInv:Mat,discr: R)
  MATSTOR ==> StorageEfficientMatrixOperations

  Exports ==> with
    integralBasis : () -> RResult
      ++ \spad{integralBasis()} returns a record
      ++ \spad{[basis,basisDen,basisInv]} containing information regarding
      ++ the integral closure of R in the quotient field of F, where
      ++ F is a framed algebra with R-module basis \spad{w1,w2,...,wn}.
      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
    localIntegralBasis : R -> RResult
      ++ \spad{integralBasis(p)} returns a record
      ++ \spad{[basis,basisDen,basisInv]} containing information regarding
      ++ the local integral closure of R at the prime \spad{p} in the quotient
      ++ field of F, where F is a framed algebra with R-module basis
      ++ \spad{w1,w2,...,wn}.
      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the local integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.

  Implementation ==> add
    import IntegralBasisTools(R, UP, F)
    import ModularHermitianRowReduction(R)
    import TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)
    import DistinctDegreeFactorize(K,R)

    listSquaredFactors: R -> List R
    listSquaredFactors px ==
      -- returns a list of the factors of px which occur with
      -- exponent > 1
      ans : List R := empty()
      factored := factor(px)$DistinctDegreeFactorize(K,R)
      for f in factors(factored) repeat
        if f.exponent > 1 then ans := concat(f.factor,ans)
      ans

    iLocalIntegralBasis: (Vector F,Vector F,Matrix R,Matrix R,R,R) -> IResult
    iLocalIntegralBasis(bas,pows,tfm,matrixOut,disc,prime) ==
      n := rank()$F; standardBasis := basis()$F
      -- 'standardBasis' is the basis for F as a FramedAlgebra;
      -- usually this is [1,y,y**2,...,y**(n-1)]
      p2 := prime * prime; sae := SAE(K,R,prime)
      p := characteristic$F; q := size()$sae
      lp := leastPower(q,n)
      rb := scalarMatrix(n,1); rbinv := scalarMatrix(n,1)
      -- rb    = basis matrix of current order
      -- rbinv = inverse basis matrix of current order
      -- these are wrt the orginal basis for F
      rbden : R := 1; index : R := 1; oldIndex : R := 1
      -- rbden = denominator for current basis matrix
      -- index = index of original order in current order
      repeat
        -- pows = [(w1 * rbden) ** q,...,(wn * rbden) ** q], where
        -- bas = [w1,...,wn] is 'rbden' times the basis for the order B = 'rb'
        for i in 1..n repeat
          bi : F := 0
          for j in 1..n repeat
            bi := bi + qelt(rb,i,j) * qelt(standardBasis,j)
          qsetelt_!(bas,i,bi)
          qsetelt_!(pows,i,bi ** p)
        coor0 := transpose coordinates(pows,bas)
        denPow := rbden ** ((p - 1) :: NNI)
        (coMat0 := coor0 exquo denPow) case "failed" =>
          error "can't happen"
        -- the jth column of coMat contains the coordinates of (wj/rbden)**q
        -- with respect to the basis [w1/rbden,...,wn/rbden]
        coMat := coMat0 :: Matrix R
        -- the ith column of 'pPows' contains the coordinates of the pth power
        -- of the ith basis element for B/prime.B over 'sae' = R/prime.R
        pPows := map(reduce,coMat)$MatrixCategoryFunctions2(R,Vector R,
                    Vector R,Matrix R,sae,Vector sae,Vector sae,Matrix sae)
        -- 'frob' will eventually be the Frobenius matrix for B/prime.B over
        -- 'sae' = R/prime.R; at each stage of the loop the ith column will
        -- contain the coordinates of p^k-th powers of the ith basis element
        frob := copy pPows; tmpMat : Matrix sae := new(n,n,0)
        for r in 2..leastPower(p,q) repeat
          for i in 1..n repeat for j in 1..n repeat
            qsetelt_!(tmpMat,i,j,qelt(frob,i,j) ** p)
          times_!(frob,pPows,tmpMat)$MATSTOR(sae)
        frobPow := frob ** lp
        -- compute the p-radical
        ns := nullSpace frobPow
        for i in 1..n repeat for j in 1..n repeat qsetelt_!(tfm,i,j,0)
        for vec in ns for i in 1.. repeat
          for j in 1..n repeat
            qsetelt_!(tfm,i,j,lift qelt(vec,j))
        id := squareTop rowEchelon(tfm,prime)
        -- id = basis matrix of the p-radical
        idinv := UpTriBddDenomInv(id, prime)
        -- id * idinv = prime * identity
        -- no need to check for inseparability in this case
        rbinv := idealiser(id * rb, rbinv * idinv, prime * rbden)
        index := diagonalProduct rbinv
        rb := rowEchelon LowTriBddDenomInv(rbinv,rbden * prime)
        if divideIfCan_!(rb,matrixOut,prime,n) = 1
          then rb := matrixOut
          else rbden := rbden * prime
        rbinv := UpTriBddDenomInv(rb,rbden)
        indexChange := index quo oldIndex
        oldIndex := index
        disc := disc quo (indexChange * indexChange)
        (not sizeLess?(1,indexChange)) or ((disc exquo p2) case "failed") =>
          return [rb, rbden, rbinv, disc]

    integralBasis() ==
      traceMat := traceMatrix()$F; n := rank()$F
      disc := determinant traceMat        -- discriminant of current order
      zero? disc => error "integralBasis: polynomial must be separable"
      singList := listSquaredFactors disc -- singularities of relative Spec
      runningRb := scalarMatrix(n,1); runningRbinv := scalarMatrix(n,1)
      -- runningRb    = basis matrix of current order
      -- runningRbinv = inverse basis matrix of current order
      -- these are wrt the original basis for F
      runningRbden : R := 1
      -- runningRbden = denominator for current basis matrix
      empty? singList => [runningRb, runningRbden, runningRbinv]
      bas : Vector F := new(n,0); pows : Vector F := new(n,0)
      -- storage for basis elements and their powers
      tfm : Matrix R := new(n,n,0)
      -- 'tfm' will contain the coordinates of a lifting of the kernel
      -- of a power of Frobenius
      matrixOut : Matrix R := new(n,n,0)
      for prime in singList repeat
        lb := iLocalIntegralBasis(bas,pows,tfm,matrixOut,disc,prime)
        rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
        disc := lb.discr
        -- update 'running integral basis' if newly computed
        -- local integral basis is non-trivial
        if sizeLess?(1,rbden) then
          mat := vertConcat(rbden * runningRb,runningRbden * rb)
          runningRbden := runningRbden * rbden
          runningRb := squareTop rowEchelon(mat,runningRbden)
          runningRbinv := UpTriBddDenomInv(runningRb,runningRbden)
      [runningRb, runningRbden, runningRbinv]

    localIntegralBasis prime ==
      traceMat := traceMatrix()$F; n := rank()$F
      disc := determinant traceMat        -- discriminant of current order
      zero? disc => error "localIntegralBasis: polynomial must be separable"
      (disc exquo (prime * prime)) case "failed" =>
        [scalarMatrix(n,1), 1, scalarMatrix(n,1)]
      bas : Vector F := new(n,0); pows : Vector F := new(n,0)
      -- storage for basis elements and their powers
      tfm : Matrix R := new(n,n,0)
      -- 'tfm' will contain the coordinates of a lifting of the kernel
      -- of a power of Frobenius
      matrixOut : Matrix R := new(n,n,0)
      lb := iLocalIntegralBasis(bas,pows,tfm,matrixOut,disc,prime)
      [lb.basis, lb.basisDen, lb.basisInv]

@
\section{package NFINTBAS NumberFieldIntegralBasis}
<<package NFINTBAS NumberFieldIntegralBasis>>=
)abbrev package NFINTBAS NumberFieldIntegralBasis
++ Author: Victor Miller, Clifton Williamson
++ Date Created: 9 April 1990
++ Date Last Updated: 20 September 1994
++ Basic Operations: discriminant, integralBasis
++ Related Domains: IntegralBasisTools, TriangularMatrixOperations
++ Also See: FunctionFieldIntegralBasis, WildFunctionFieldIntegralBasis
++ AMS Classifications:
++ Keywords: number field, integral basis, discriminant
++ Examples:
++ References:
++ Description:
++   In this package F is a framed algebra over the integers (typically
++   \spad{F = Z[a]} for some algebraic integer a).  The package provides
++   functions to compute the integral closure of Z in the quotient
++   quotient field of F.
NumberFieldIntegralBasis(UP,F): Exports == Implementation where
  UP : UnivariatePolynomialCategory Integer
  F  : FramedAlgebra(Integer,UP)

  FR  ==> Factored Integer
  I   ==> Integer
  Mat ==> Matrix I
  NNI ==> NonNegativeInteger
  Ans ==> Record(basis: Mat, basisDen: I, basisInv:Mat,discr: I)

  Exports ==> with
    discriminant: () -> Integer
      ++ \spad{discriminant()} returns the discriminant of the integral
      ++ closure of Z in the quotient field of the framed algebra F.
    integralBasis : () -> Record(basis: Mat, basisDen: I, basisInv:Mat)
      ++ \spad{integralBasis()} returns a record
      ++ \spad{[basis,basisDen,basisInv]}
      ++ containing information regarding the integral closure of Z in the
      ++ quotient field of F, where F is a framed algebra with Z-module
      ++ basis \spad{w1,w2,...,wn}.
      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
    localIntegralBasis : I -> Record(basis: Mat, basisDen: I, basisInv:Mat)
      ++ \spad{integralBasis(p)} returns a record
      ++ \spad{[basis,basisDen,basisInv]} containing information regarding
      ++ the local integral closure of Z at the prime \spad{p} in the quotient
      ++ field of F, where F is a framed algebra with Z-module basis
      ++ \spad{w1,w2,...,wn}.
      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
      ++ the \spad{i}th element of the integral basis is
      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
      ++ respect to the basis \spad{v1,...,vn}: if \spad{basisInv} is the
      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
      ++ \spad{wi = sum(bij * vj, j = 1..n)}.

  Implementation ==> add
    import IntegralBasisTools(I, UP, F)
    import ModularHermitianRowReduction(I)
    import TriangularMatrixOperations(I, Vector I, Vector I, Matrix I)

    frobMatrix              : (Mat,Mat,I,NNI) -> Mat
    wildPrimes              : (FR,I) -> List I
    tameProduct             : (FR,I) -> I
    iTameLocalIntegralBasis : (Mat,I,I) -> Ans
    iWildLocalIntegralBasis : (Mat,I,I) -> Ans

    frobMatrix(rb,rbinv,rbden,p) ==
      n := rank()$F; b := basis()$F
      v : Vector F := new(n,0)
      for i in minIndex(v)..maxIndex(v)
        for ii in minRowIndex(rb)..maxRowIndex(rb) repeat
          a : F := 0
          for j in minIndex(b)..maxIndex(b)
            for jj in minColIndex(rb)..maxColIndex(rb) repeat
              a := a + qelt(rb,ii,jj) * qelt(b,j)
          qsetelt_!(v,i,a**p)
      mat := transpose coordinates v
      ((transpose(rbinv) * mat) exquo (rbden ** p)) :: Mat

    wildPrimes(factoredDisc,n) ==
      -- returns a list of the primes <=n which divide factoredDisc to a
      -- power greater than 1
      ans : List I := empty()
      for f in factors(factoredDisc) repeat
        if f.exponent > 1 and f.factor <= n then ans := concat(f.factor,ans)
      ans

    tameProduct(factoredDisc,n) ==
      -- returns the product of the primes > n which divide factoredDisc
      -- to a power greater than 1
      ans : I := 1
      for f in factors(factoredDisc) repeat
        if f.exponent > 1 and f.factor > n then ans := f.factor * ans
      ans

    integralBasis() ==
      traceMat := traceMatrix()$F; n := rank()$F
      disc := determinant traceMat  -- discriminant of current order
      disc0 := disc                 -- this is disc(F)
      factoredDisc := factor(disc0)$IntegerFactorizationPackage(Integer)
      wilds := wildPrimes(factoredDisc,n)
      sing := tameProduct(factoredDisc,n)
      runningRb := scalarMatrix(n, 1); runningRbinv := scalarMatrix(n, 1)
      -- runningRb    = basis matrix of current order
      -- runningRbinv = inverse basis matrix of current order
      -- these are wrt the original basis for F
      runningRbden : I := 1
      -- runningRbden = denominator for current basis matrix
      one? sing and empty? wilds => [runningRb, runningRbden, runningRbinv]
      -- id = basis matrix of the ideal (p-radical) wrt current basis
      matrixOut : Mat := scalarMatrix(n,0)
      for p in wilds repeat
        lb := iWildLocalIntegralBasis(matrixOut,disc,p)
        rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
        disc := lb.discr
        -- update 'running integral basis' if newly computed
        -- local integral basis is non-trivial
        if sizeLess?(1,rbden) then
          mat := vertConcat(rbden * runningRb,runningRbden * rb)
          runningRbden := runningRbden * rbden
          runningRb := squareTop rowEchelon(mat,runningRbden)
          runningRbinv := UpTriBddDenomInv(runningRb,runningRbden)
      lb := iTameLocalIntegralBasis(traceMat,disc,sing)
      rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
      disc := lb.discr
      -- update 'running integral basis' if newly computed
      -- local integral basis is non-trivial
      if sizeLess?(1,rbden) then
        mat := vertConcat(rbden * runningRb,runningRbden * rb)
        runningRbden := runningRbden * rbden
        runningRb := squareTop rowEchelon(mat,runningRbden)
        runningRbinv := UpTriBddDenomInv(runningRb,runningRbden)
      [runningRb,runningRbden,runningRbinv]

    localIntegralBasis p ==
      traceMat := traceMatrix()$F; n := rank()$F
      disc := determinant traceMat  -- discriminant of current order
      (disc exquo (p*p)) case "failed" =>
        [scalarMatrix(n, 1), 1, scalarMatrix(n, 1)]
      lb :=
        p > rank()$F =>
          iTameLocalIntegralBasis(traceMat,disc,p)
        iWildLocalIntegralBasis(scalarMatrix(n,0),disc,p)
      [lb.basis,lb.basisDen,lb.basisInv]

    iTameLocalIntegralBasis(traceMat,disc,sing) ==
      n := rank()$F; disc0 := disc
      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
      -- rb    = basis matrix of current order
      -- rbinv = inverse basis matrix of current order
      -- these are wrt the original basis for F
      rbden : I := 1; index : I := 1; oldIndex : I := 1
      -- rbden = denominator for current basis matrix
      -- id = basis matrix of the ideal (p-radical) wrt current basis
      tfm := traceMat
      repeat
        -- compute the p-radical = p-trace-radical
        idinv := transpose squareTop rowEchelon(tfm,sing)
        -- [u1,..,un] are the coordinates of an element of the p-radical
        -- iff [u1,..,un] * idinv is in p * Z^n
        id := rowEchelon LowTriBddDenomInv(idinv, sing)
        -- id = basis matrix of the p-radical
        idinv := UpTriBddDenomInv(id, sing)
        -- id * idinv = sing * identity
        -- no need to check for inseparability in this case
        rbinv := idealiser(id * rb, rbinv * idinv, sing * rbden)
        index := diagonalProduct rbinv
        rb := rowEchelon LowTriBddDenomInv(rbinv, sing * rbden)
        g := matrixGcd(rb,sing,n)
        if sizeLess?(1,g) then rb := (rb exquo g) :: Mat
        rbden := rbden * (sing quo g)
        rbinv := UpTriBddDenomInv(rb, rbden)
        disc := disc0 quo (index * index)
        indexChange := index quo oldIndex; oldIndex := index
        one? indexChange => return [rb, rbden, rbinv, disc]
        tfm := ((rb * traceMat * transpose rb) exquo (rbden * rbden)) :: Mat

    iWildLocalIntegralBasis(matrixOut,disc,p) ==
      n := rank()$F; disc0 := disc
      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
      -- rb    = basis matrix of current order
      -- rbinv = inverse basis matrix of current order
      -- these are wrt the original basis for F
      rbden : I := 1; index : I := 1; oldIndex : I := 1
      -- rbden = denominator for current basis matrix
      -- id = basis matrix of the ideal (p-radical) wrt current basis
      p2 := p * p; lp := leastPower(p::NNI,n)
      repeat
        tfm := frobMatrix(rb,rbinv,rbden,p::NNI) ** lp
        -- compute Rp = p-radical
        idinv := transpose squareTop rowEchelon(tfm, p)
        -- [u1,..,un] are the coordinates of an element of Rp
        -- iff [u1,..,un] * idinv is in p * Z^n
        id := rowEchelon LowTriBddDenomInv(idinv,p)
        -- id = basis matrix of the p-radical
        idinv := UpTriBddDenomInv(id,p)
        -- id * idinv = p * identity
        -- no need to check for inseparability in this case
        rbinv := idealiser(id * rb, rbinv * idinv, p * rbden)
        index := diagonalProduct rbinv
        rb := rowEchelon LowTriBddDenomInv(rbinv, p * rbden)
        if divideIfCan_!(rb,matrixOut,p,n) = 1
          then rb := matrixOut
          else rbden := p * rbden
        rbinv := UpTriBddDenomInv(rb, rbden)
        indexChange := index quo oldIndex; oldIndex := index
        disc := disc quo (indexChange * indexChange)
        one? indexChange or gcd(p2,disc) ~= p2 =>
          return [rb, rbden, rbinv, disc]

    discriminant() ==
      disc := determinant traceMatrix()$F
      intBas := integralBasis()
      rb := intBas.basis; rbden := intBas.basisDen
      index := ((rbden ** rank()$F) exquo (determinant rb)) :: Integer
      (disc exquo (index * index)) :: Integer

@
\section{License}
<<license>>=
--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.
@
<<*>>=
<<license>>

<<package TRIMAT TriangularMatrixOperations>>
<<package IBATOOL IntegralBasisTools>>
<<package FFINTBAS FunctionFieldIntegralBasis>>
<<package WFFINTBS WildFunctionFieldIntegralBasis>>
<<package NFINTBAS NumberFieldIntegralBasis>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}