Initial population.

1 files changed, 816 insertions, 0 deletions

diff --git a/src/algebra/intclos.spad.pamphlet b/src/algebra/intclos.spad.pamphlet
new file mode 100644
index 00000000..4470fc41
--- /dev/null
+++ b/src/algebra/intclos.spad.pamphlet
@@ -0,0 +1,816 @@
+\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 = 1) => 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]
+      (sing = 1) 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]
+        (indexChange = 1) => 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 =>
+        (indexChange = 1) 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}