1)abbrev package TRIMAT TriangularMatrixOperations
2++ Fraction free inverses of triangular matrices
3++ Author: Victor Miller
4++ Keywords:
5++ Examples:
6++ References:
7++ Description:
8++ This package provides functions that compute "fraction-free"
9++ inverses of upper and lower triangular matrices over a integral
10++ domain.  By "fraction-free inverses" we mean the following:
11++ given a matrix B with entries in R and an element d of R such that
12++ d * inv(B) also has entries in R, we return d * inv(B).  Thus,
13++ it is not necessary to pass to the quotient field in any of our
14++ computations.
15
16
17TriangularMatrixOperations(R, Row, Col, M) : Exports == Implementation where
18  R   : IntegralDomain
19  Row : FiniteLinearAggregate R
20  Col : FiniteLinearAggregate R
21  M   : MatrixCategory(R, Row, Col)
22
23  Exports ==> with
24
25    UpTriBddDenomInv : (M, R) -> M
26      ++ UpTriBddDenomInv(B, d) returns M, where
27      ++ B is a non-singular upper triangular matrix and d is an
28      ++ element of R such that \spad{M = d * inv(B)} has entries in R.
29    LowTriBddDenomInv : (M, R) -> M
30      ++ LowTriBddDenomInv(B, d) returns M, where
31      ++ B is a non-singular lower triangular matrix and d is an
32      ++ element of R such that \spad{M = d * inv(B)} has entries in R.
33
34  Implementation ==> add
35
36    UpTriBddDenomInv(A, denom) ==
37      AI := zero(nrows A, nrows A)$M
38      offset := minColIndex AI - minRowIndex AI
39      for i in minRowIndex AI .. maxRowIndex AI
40        for j in minColIndex AI .. maxColIndex AI repeat
41          qsetelt!(AI, i, j, (denom exquo qelt(A, i, j))::R)
42      for i in minRowIndex AI .. maxRowIndex AI repeat
43        for j in offset + i + 1 .. maxColIndex AI repeat
44          qsetelt!(AI, i, j, - (((+/[qelt(AI, i, k) * qelt(A, k-offset, j)
45                                   for k in i+offset..(j-1)])
46                                     exquo qelt(A, j-offset, j))::R))
47      AI
48
49    LowTriBddDenomInv(A, denom) ==
50      AI := zero(nrows A, nrows A)$M
51      offset := minColIndex AI - minRowIndex AI
52      for i in minRowIndex AI .. maxRowIndex AI
53        for j in minColIndex AI .. maxColIndex AI repeat
54          qsetelt!(AI, i, j, (denom exquo qelt(A, i, j))::R)
55      for i in minColIndex AI .. maxColIndex AI repeat
56        for j in i - offset + 1 .. maxRowIndex AI repeat
57          qsetelt!(AI, j, i, - (((+/[qelt(A, j, k+offset) * qelt(AI, k, i)
58                                    for k in i-offset..(j-1)])
59                                      exquo qelt(A, j, j+offset))::R))
60      AI
61
62)abbrev package IBATOOL IntegralBasisTools
63++ Functions common to both integral basis packages
64++ Author: Victor Miller, Barry Trager, Clifton Williamson
65++ Date Created: 11 April 1990
66++ Keywords: integral basis, function field, number field
67++ Examples:
68++ References:
69++ Description:
70++ This package contains functions used in the packages
71++ FunctionFieldIntegralBasis and NumberFieldIntegralBasis.
72
73IntegralBasisTools(R, UP, F) : Exports == Implementation where
74  R  : EuclideanDomain with
75        squareFree : % -> Factored %
76                ++ squareFree(x) returns a square-free factorisation of x
77  UP : UnivariatePolynomialCategory R
78  F  : FramedAlgebra(R, UP)
79  Mat ==> Matrix R
80  NNI ==> NonNegativeInteger
81  Ans ==> Record(basis : Mat, basisDen : R, basisInv : Mat)
82
83  Exports ==> with
84
85    diagonalProduct : Mat -> R
86      ++ diagonalProduct(m) returns the product of the elements on the
87      ++ diagonal of the matrix m
88    matrixGcd : (Mat, R, NNI) -> R
89      ++ matrixGcd(mat, sing, n) is \spad{gcd(sing, g)} where \spad{g} is the
90      ++ gcd of the entries of the \spad{n}-by-\spad{n} upper-triangular
91      ++ matrix \spad{mat}.
92    divideIfCan! : (Matrix R, Matrix R, R, Integer) -> R
93      ++ divideIfCan!(matrix, matrixOut, prime, n) attempts to divide the
94      ++ entries of \spad{matrix} by \spad{prime} and store the result in
95      ++ \spad{matrixOut}.  If it is successful, 1 is returned and if not,
96      ++ \spad{prime} is returned.  Here both \spad{matrix} and
97      ++ \spad{matrixOut} are \spad{n}-by-\spad{n} upper triangular matrices.
98    leastPower : (NNI, NNI) -> NNI
99      ++ leastPower(p, n) returns e, where e is the smallest integer
100      ++ such that \spad{p ^e >= n}
101    idealiser : (Mat, Mat) -> Mat
102      ++ idealiser(m1, m2) computes the order of an ideal defined by m1 and m2
103    idealiser : (Mat, Mat, R) -> Mat
104      ++ idealiser(m1, m2, d) computes the order of an ideal defined by m1 and m2
105      ++ where d is the known part of the denominator
106    idealiserMatrix : (Mat, Mat) -> Mat
107      ++ idealiserMatrix(m1, m2) returns the matrix representing the linear
108      ++ conditions on the Ring associated with an ideal defined by m1 and m2.
109    moduleSum : (Ans, Ans) -> Ans
110      ++ moduleSum(m1, m2) returns the sum of two modules in the framed
111      ++ algebra \spad{F}.  Each module \spad{mi} is represented as follows:
112      ++ F is a framed algebra with R-module basis \spad{w1, w2, ..., wn} and
113      ++ \spad{mi} is a record \spad{[basis, basisDen, basisInv]}.  If
114      ++ \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
115      ++ a basis \spad{v1, ..., vn} for \spad{mi} is given by
116      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
117      ++ \spad{i}th row of 'basis' contains the coordinates of the
118      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
119      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
120      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
121      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
122      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
123
124  Implementation ==> add
125    import from ModularHermitianRowReduction(R)
126    import from TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)
127
128    diagonalProduct m ==
129      ans : R := 1
130      for i in minRowIndex m .. maxRowIndex m
131        for j in minColIndex m .. maxColIndex m repeat
132          ans := ans * qelt(m, i, j)
133      ans
134
135    matrixGcd(mat, sing, n) ==
136      -- note: 'matrix' is upper triangular;
137      -- no need to do anything below the diagonal
138      d := sing
139      for i in 1..n repeat
140        for j in i..n repeat
141          if not zero?(mij := qelt(mat, i, j)) then d := gcd(d, mij)
142          (d = 1) => return d
143      d
144
145    divideIfCan!(matrix, matrixOut, prime, n) ==
146    -- note: both 'matrix' and 'matrixOut' will be upper triangular;
147    -- no need to do anything below the diagonal
148      for i in 1..n repeat
149        for j in i..n repeat
150          (a := (qelt(matrix,i,j) exquo prime)) case "failed" => return prime
151          qsetelt!(matrixOut, i, j, a :: R)
152      1
153
154    leastPower(p, n) ==
155      -- efficiency is not an issue here
156      e : NNI := 1; q := p
157      while q < n repeat (e := e + 1; q := q * p)
158      e
159
160    idealiserMatrix(ideal, idealinv) ==
161      -- computes the Order of the ideal
162      n  := rank()$F
163      bigm := zero(n * n, n)$Mat
164      mr := minRowIndex bigm; mc := minColIndex bigm
165      v := basis()$F
166      for i in 0..n-1 repeat
167        r := transpose regularRepresentation qelt(v, i + minIndex v)
168        m := ideal * r * idealinv
169        for j in 0..n-1 repeat
170          for k in 0..n-1 repeat
171            bigm(j * n + k + mr, i + mc) := qelt(m, j + mr, k + mc)
172      bigm
173
174    idealiser(ideal, idealinv) ==
175      bigm := idealiserMatrix(ideal, idealinv)
176      transpose squareTop rowEch bigm
177
178    idealiser(ideal, idealinv, denom) ==
179      bigm := (idealiserMatrix(ideal, idealinv) exquo denom)::Mat
180      transpose squareTop rowEchelon(bigm, denom)
181
182    moduleSum(mod1, mod2) ==
183      rb1 := mod1.basis; rbden1 := mod1.basisDen; rbinv1 := mod1.basisInv
184      rb2 := mod2.basis; rbden2 := mod2.basisDen; rbinv2 := mod2.basisInv
185      -- compatibility check: doesn't take much computation time
186      (not square? rb1) or (not square? rbinv1) or (not square? rb2) _
187        or (not square? rbinv2) =>
188        error "moduleSum: matrices must be square"
189      ((n := nrows rb1) ~= (nrows rbinv1)) or (n ~= (nrows rb2)) _
190        or (n ~= (nrows rbinv2)) =>
191        error "moduleSum: matrices of incompatible dimensions"
192      (zero? rbden1) or (zero? rbden2) =>
193        error "moduleSum: denominator must be non-zero"
194      den := lcm(rbden1, rbden2); c1 := den quo rbden1; c2 := den quo rbden2
195      rb := squareTop rowEchelon(vertConcat(c1 * rb1, c2 * rb2), den)
196      rbinv := UpTriBddDenomInv(rb, den)
197      [rb, den, rbinv]
198
199)abbrev package FFINTBAS FunctionFieldIntegralBasis
200++ Integral bases for function fields of dimension one
201++ Author: Victor Miller
202++ Date Created: 9 April 1990
203++ Keywords:
204++ Examples:
205++ References:
206++ Description:
207++ In this package R is a Euclidean domain and F is a framed algebra
208++ over R.  The package provides functions to compute the integral
209++ closure of R in the quotient field of F.  It is assumed that
210++ \spad{char(R/P) = char(R)} for any prime P of R.  A typical instance of
211++ this is when \spad{R = K[x]} and F is a function field over R.
212
213
214FunctionFieldIntegralBasis(R, UP, F) : Exports == Implementation where
215  R  : EuclideanDomain with
216        squareFree : % -> Factored %
217                ++ squareFree(x) returns a square-free factorisation of x
218  UP : UnivariatePolynomialCategory R
219  F  : FramedAlgebra(R, UP)
220
221  I   ==> Integer
222  Mat ==> Matrix R
223  NNI ==> NonNegativeInteger
224
225  Exports ==> with
226    integralBasis : () -> Record(basis : Mat, basisDen : R, basisInv : Mat)
227      ++ \spad{integralBasis()} returns a record
228      ++ \spad{[basis, basisDen, basisInv]} containing information regarding
229      ++ the integral closure of R in the quotient field of F, where
230      ++ F is a framed algebra with R-module basis \spad{w1, w2, ..., wn}.
231      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
232      ++ the \spad{i}th element of the integral basis is
233      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
234      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
235      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
236      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
237      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
238      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
239      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
240    localIntegralBasis : R -> Record(basis : Mat, basisDen : R, basisInv : Mat)
241      ++ \spad{integralBasis(p)} returns a record
242      ++ \spad{[basis, basisDen, basisInv]} containing information regarding
243      ++ the local integral closure of R at the prime \spad{p} in the quotient
244      ++ field of F, where F is a framed algebra with R-module basis
245      ++ \spad{w1, w2, ..., wn}.
246      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
247      ++ the \spad{i}th element of the local integral basis is
248      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
249      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
250      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
251      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
252      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
253      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
254      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
255
256  Implementation ==> add
257    import from IntegralBasisTools(R, UP, F)
258    import from ModularHermitianRowReduction(R)
259    import from TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)
260
261    squaredFactors : R -> R
262    squaredFactors px ==
263           */[(if ffe.exponent > 1 then ffe.factor else 1$R)
264                for ffe in factorList squareFree px]
265
266    iIntegralBasis : (Mat, R, R) -> Record(basis : Mat, basisDen : R, basisInv : Mat)
267    iIntegralBasis(tfm, disc, sing) ==
268      -- tfm = trace matrix of current order
269      n := rank()$F; tfm0 := copy tfm; disc0 := disc
270      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
271      -- rb    = basis matrix of current order
272      -- rbinv = inverse basis matrix of current order
273      -- these are wrt the original basis for F
274      rbden : R := 1; index : R := 1; oldIndex : R := 1
275      -- rbden = denominator for current basis matrix
276      -- index = index of original order in current order
277      not sizeLess?(1, sing) => [rb, rbden, rbinv]
278      repeat
279        -- compute the p-radical
280        idinv := transpose squareTop rowEchelon(tfm, sing)
281        -- [u1, .., un] are the coordinates of an element of the p-radical
282        -- iff [u1, .., un] * idinv is in sing * R^n
283        id := rowEchelon LowTriBddDenomInv(idinv, sing)
284        -- id = basis matrix of the p-radical
285        idinv := UpTriBddDenomInv(id, sing)
286        -- id * idinv = sing * identity
287        -- no need to check for inseparability in this case
288        rbinv := idealiser(id * rb, rbinv * idinv, sing * rbden)
289        index := diagonalProduct rbinv
290        rb := rowEchelon LowTriBddDenomInv(rbinv, rbden * sing)
291        g := matrixGcd(rb, sing, n)
292        if sizeLess?(1, g) then rb := (rb exquo g) :: Mat
293        rbden := rbden * (sing quo g)
294        rbinv := UpTriBddDenomInv(rb, rbden)
295        disc := disc0 quo (index * index)
296        indexChange := index quo oldIndex; oldIndex := index
297        sing := gcd(indexChange, squaredFactors disc)
298        not sizeLess?(1, sing) => return [rb, rbden, rbinv]
299        tfm := ((rb * tfm0 * transpose rb) exquo (rbden * rbden)) :: Mat
300
301    integralBasis() ==
302      n := rank()$F; p := characteristic()$F
303      (not zero? p) and (n >= p) =>
304        error "integralBasis: possible wild ramification"
305      tfm := traceMatrix()$F; disc := determinant tfm
306      sing := squaredFactors disc    -- singularities of relative Spec
307      iIntegralBasis(tfm, disc, sing)
308
309    localIntegralBasis prime ==
310      n := rank()$F; p := characteristic()$F
311      (not zero? p) and (n >= p) =>
312        error "integralBasis: possible wild ramification"
313      tfm := traceMatrix()$F; disc := determinant tfm
314      (disc exquo (prime * prime)) case "failed" =>
315        [scalarMatrix(n, 1), 1, scalarMatrix(n, 1)]
316      iIntegralBasis(tfm, disc, prime)
317
318)abbrev package WFFINTBS WildFunctionFieldIntegralBasis
319++ Authors: Victor Miller, Clifton Williamson
320++ Date Created: 24 July 1991
321++ Basic Operations: integralBasis, localIntegralBasis
322++ Related Domains: IntegralBasisTools(R, UP, F),
323++   TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)
324++ Also See: FunctionFieldIntegralBasis, NumberFieldIntegralBasis
325++ AMS Classifications:
326++ Keywords: function field, integral basis
327++ Examples:
328++ References:
329++ Description:
330++ In this package K is a finite field, R is a ring of univariate
331++ polynomials over K, and F is a framed algebra over R.  The package
332++ provides a function to compute the integral closure of R in the quotient
333++ field of F as well as a function to compute a "local integral basis"
334++ at a specific prime.
335
336WildFunctionFieldIntegralBasis(K, R, UP, F) : Exports == Implementation where
337  K  : FiniteFieldCategory
338  --K  : Join(Field, Finite)
339  R  : UnivariatePolynomialCategory K
340  UP : UnivariatePolynomialCategory R
341  F  : FramedAlgebra(R, UP)
342
343  I       ==> Integer
344  Mat     ==> Matrix R
345  NNI     ==> NonNegativeInteger
346  SAE     ==> SimpleAlgebraicExtension
347  RResult ==> Record(basis : Mat, basisDen : R, basisInv : Mat)
348  IResult ==> Record(basis : Mat, basisDen : R, basisInv : Mat, discr : R)
349  MATSTOR ==> StorageEfficientMatrixOperations
350
351  Exports ==> with
352    integralBasis : () -> RResult
353      ++ \spad{integralBasis()} returns a record
354      ++ \spad{[basis, basisDen, basisInv]} containing information regarding
355      ++ the integral closure of R in the quotient field of F, where
356      ++ F is a framed algebra with R-module basis \spad{w1, w2, ..., wn}.
357      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
358      ++ the \spad{i}th element of the integral basis is
359      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
360      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
361      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
362      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
363      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
364      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
365      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
366    localIntegralBasis : R -> RResult
367      ++ \spad{integralBasis(p)} returns a record
368      ++ \spad{[basis, basisDen, basisInv]} containing information regarding
369      ++ the local integral closure of R at the prime \spad{p} in the quotient
370      ++ field of F, where F is a framed algebra with R-module basis
371      ++ \spad{w1, w2, ..., wn}.
372      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
373      ++ the \spad{i}th element of the local integral basis is
374      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
375      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
376      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
377      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
378      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
379      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
380      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
381
382  Implementation ==> add
383    import from IntegralBasisTools(R, UP, F)
384    import from ModularHermitianRowReduction(R)
385    import from TriangularMatrixOperations(R, Vector R, Vector R, Matrix R)
386    import from DistinctDegreeFactorize(K, R)
387
388    listSquaredFactors : R -> List R
389    listSquaredFactors px ==
390      -- returns a list of the factors of px which occur with
391      -- exponent > 1
392      ans : List R := empty()
393      factored := factor(px)$DistinctDegreeFactorize(K, R)
394      for f in factorList(factored) repeat
395        if f.exponent > 1 then ans := concat(f.factor, ans)
396      ans
397
398    iLocalIntegralBasis : (Vector F, Vector F, Matrix R, Matrix R, R, R) -> IResult
399    iLocalIntegralBasis(bas, pows, tfm, matrixOut, disc, prime) ==
400      n := rank()$F; standardBasis := basis()$F
401      -- 'standardBasis' is the basis for F as a FramedAlgebra;
402      -- usually this is [1, y, y^2, ..., y^(n-1)]
403      p2 := prime * prime; sae := SAE(K, R, prime)
404      p := characteristic()$F; q := size()$sae
405      lp := leastPower(q, n)
406      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
407      -- rb    = basis matrix of current order
408      -- rbinv = inverse basis matrix of current order
409      -- these are wrt the original basis for F
410      rbden : R := 1; index : R := 1; oldIndex : R := 1
411      -- rbden = denominator for current basis matrix
412      -- index = index of original order in current order
413      repeat
414        -- pows = [(w1 * rbden) ^ q, ..., (wn * rbden) ^ q], where
415        -- bas = [w1, ..., wn] is 'rbden' times the basis for the order B = 'rb'
416        for i in 1..n repeat
417          bi : F := 0
418          for j in 1..n repeat
419            bi := bi + qelt(rb, i, j) * qelt(standardBasis, j)
420          qsetelt!(bas, i, bi)
421          qsetelt!(pows, i, bi ^ p)
422        coor0 := transpose coordinates(pows, bas)
423        denPow := rbden ^ ((p - 1) :: NNI)
424        (coMat0 := coor0 exquo denPow) case "failed" =>
425          error "can't happen"
426        -- the jth column of coMat contains the coordinates of (wj/rbden)^q
427        -- with respect to the basis [w1/rbden, ..., wn/rbden]
428        coMat := coMat0 :: Matrix R
429        -- the ith column of 'pPows' contains the coordinates of the pth power
430        -- of the ith basis element for B/prime.B over 'sae' = R/prime.R
431        pPows : Matrix sae :=
432                map(reduce, coMat)$MatrixCategoryFunctions2(R, Vector R,
433                    Vector R, Matrix R, sae, Vector sae, Vector sae, Matrix sae)
434        -- 'frob' will eventually be the Frobenius matrix for B/prime.B over
435        -- 'sae' = R/prime.R; at each stage of the loop the ith column will
436        -- contain the coordinates of p^k-th powers of the ith basis element
437        frob := copy pPows; tmpMat : Matrix sae := new(n, n, 0)
438        for r in 2..leastPower(p, q) repeat
439          for i in 1..n repeat for j in 1..n repeat
440            qsetelt!(tmpMat, i, j, qelt(frob, i, j) ^ p)
441          times!(frob, pPows, tmpMat)$MATSTOR(sae)
442        frobPow := frob ^ lp
443        -- compute the p-radical
444        ns := nullSpace frobPow
445        for i in 1..n repeat for j in 1..n repeat qsetelt!(tfm, i, j, 0)
446        for vec in ns for i in 1.. repeat
447          for j in 1..n repeat
448            qsetelt!(tfm, i, j, lift qelt(vec, j))
449        id := squareTop rowEchelon(tfm, prime)
450        -- id = basis matrix of the p-radical
451        idinv := UpTriBddDenomInv(id, prime)
452        -- id * idinv = prime * identity
453        -- no need to check for inseparability in this case
454        rbinv := idealiser(id * rb, rbinv * idinv, prime * rbden)
455        index := diagonalProduct rbinv
456        rb := rowEchelon LowTriBddDenomInv(rbinv, rbden * prime)
457        if divideIfCan!(rb, matrixOut, prime, n) = 1
458          then rb := matrixOut
459          else rbden := rbden * prime
460        rbinv := UpTriBddDenomInv(rb, rbden)
461        indexChange := index quo oldIndex
462        oldIndex := index
463        disc := disc quo (indexChange * indexChange)
464        (not sizeLess?(1,indexChange)) or ((disc exquo p2) case "failed") =>
465          return [rb, rbden, rbinv, disc]
466
467    integralBasis() ==
468      traceMat := traceMatrix()$F; n := rank()$F
469      disc := determinant traceMat        -- discriminant of current order
470      zero? disc => error "integralBasis: polynomial must be separable"
471      singList := listSquaredFactors disc -- singularities of relative Spec
472      runningRb := scalarMatrix(n, 1); runningRbinv := scalarMatrix(n, 1)
473      -- runningRb    = basis matrix of current order
474      -- runningRbinv = inverse basis matrix of current order
475      -- these are wrt the original basis for F
476      runningRbden : R := 1
477      -- runningRbden = denominator for current basis matrix
478      empty? singList => [runningRb, runningRbden, runningRbinv]
479      bas : Vector F := new(n, 0); pows : Vector F := new(n, 0)
480      -- storage for basis elements and their powers
481      tfm : Matrix R := new(n, n, 0)
482      -- 'tfm' will contain the coordinates of a lifting of the kernel
483      -- of a power of Frobenius
484      matrixOut : Matrix R := new(n, n, 0)
485      for prime in singList repeat
486        lb := iLocalIntegralBasis(bas, pows, tfm, matrixOut, disc, prime)
487        rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
488        disc := lb.discr
489        -- update 'running integral basis' if newly computed
490        -- local integral basis is non-trivial
491        if sizeLess?(1, rbden) then
492          mat := vertConcat(rbden * runningRb, runningRbden * rb)
493          runningRbden := runningRbden * rbden
494          runningRb := squareTop rowEchelon(mat, runningRbden)
495          runningRbinv := UpTriBddDenomInv(runningRb, runningRbden)
496      [runningRb, runningRbden, runningRbinv]
497
498    localIntegralBasis prime ==
499      traceMat := traceMatrix()$F; n := rank()$F
500      disc := determinant traceMat        -- discriminant of current order
501      zero? disc => error "localIntegralBasis: polynomial must be separable"
502      (disc exquo (prime * prime)) case "failed" =>
503        [scalarMatrix(n, 1), 1, scalarMatrix(n, 1)]
504      bas : Vector F := new(n, 0); pows : Vector F := new(n, 0)
505      -- storage for basis elements and their powers
506      tfm : Matrix R := new(n, n, 0)
507      -- 'tfm' will contain the coordinates of a lifting of the kernel
508      -- of a power of Frobenius
509      matrixOut : Matrix R := new(n, n, 0)
510      lb := iLocalIntegralBasis(bas, pows, tfm, matrixOut, disc, prime)
511      [lb.basis, lb.basisDen, lb.basisInv]
512
513)abbrev package NFINTBAS NumberFieldIntegralBasis
514++ Author: Victor Miller, Clifton Williamson
515++ Date Created: 9 April 1990
516++ Basic Operations: discriminant, integralBasis
517++ Related Domains: IntegralBasisTools, TriangularMatrixOperations
518++ Also See: FunctionFieldIntegralBasis, WildFunctionFieldIntegralBasis
519++ AMS Classifications:
520++ Keywords: number field, integral basis, discriminant
521++ Examples:
522++ References:
523++ Description:
524++   In this package F is a framed algebra over the integers (typically
525++   \spad{F = Z[a]} for some algebraic integer a).  The package provides
526++   functions to compute the integral closure of Z in the quotient
527++   field of F.
528NumberFieldIntegralBasis(UP, F) : Exports == Implementation where
529  UP : UnivariatePolynomialCategory Integer
530  F  : FramedAlgebra(Integer, UP)
531
532  FR  ==> Factored Integer
533  I   ==> Integer
534  Mat ==> Matrix I
535  NNI ==> NonNegativeInteger
536  Ans ==> Record(basis : Mat, basisDen : I, basisInv : Mat, discr : I)
537
538  Exports ==> with
539    discriminant : () -> Integer
540      ++ \spad{discriminant()} returns the discriminant of the integral
541      ++ closure of Z in the quotient field of the framed algebra F.
542    integralBasis : () -> Record(basis : Mat, basisDen : I, basisInv : Mat)
543      ++ \spad{integralBasis()} returns a record
544      ++ \spad{[basis, basisDen, basisInv]}
545      ++ containing information regarding the integral closure of Z in the
546      ++ quotient field of F, where F is a framed algebra with Z-module
547      ++ basis \spad{w1, w2, ..., wn}.
548      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
549      ++ the \spad{i}th element of the integral basis is
550      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
551      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
552      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
553      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
554      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
555      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
556      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
557    localIntegralBasis : I -> Record(basis : Mat, basisDen : I, basisInv : Mat)
558      ++ \spad{integralBasis(p)} returns a record
559      ++ \spad{[basis, basisDen, basisInv]} containing information regarding
560      ++ the local integral closure of Z at the prime \spad{p} in the quotient
561      ++ field of F, where F is a framed algebra with Z-module basis
562      ++ \spad{w1, w2, ..., wn}.
563      ++ If \spad{basis} is the matrix \spad{(aij, i = 1..n, j = 1..n)}, then
564      ++ the \spad{i}th element of the integral basis is
565      ++ \spad{vi = (1/basisDen) * sum(aij * wj, j = 1..n)}, i.e. the
566      ++ \spad{i}th row of \spad{basis} contains the coordinates of the
567      ++ \spad{i}th basis vector.  Similarly, the \spad{i}th row of the
568      ++ matrix \spad{basisInv} contains the coordinates of \spad{wi} with
569      ++ respect to the basis \spad{v1, ..., vn}: if \spad{basisInv} is the
570      ++ matrix \spad{(bij, i = 1..n, j = 1..n)}, then
571      ++ \spad{wi = sum(bij * vj, j = 1..n)}.
572
573  Implementation ==> add
574    import from IntegralBasisTools(I, UP, F)
575    import from ModularHermitianRowReduction(I)
576    import from TriangularMatrixOperations(I, Vector I, Vector I, Matrix I)
577
578    frobMatrix              : (Mat, Mat, I, NNI) -> Mat
579    wildPrimes              : (FR, I) -> List I
580    tameProduct             : (FR, I) -> I
581    iTameLocalIntegralBasis : (Mat, I, I) -> Ans
582    iWildLocalIntegralBasis : (Mat, I, I) -> Ans
583
584    frobMatrix(rb, rbinv, rbden, p) ==
585      n := rank()$F; b := basis()$F
586      v : Vector F := new(n, 0)
587      for i in minIndex(v)..maxIndex(v)
588        for ii in minRowIndex(rb)..maxRowIndex(rb) repeat
589          a : F := 0
590          for j in minIndex(b)..maxIndex(b)
591            for jj in minColIndex(rb)..maxColIndex(rb) repeat
592              a := a + qelt(rb, ii, jj) * qelt(b, j)
593          qsetelt!(v, i, a^p)
594      mat := transpose coordinates v
595      ((transpose(rbinv) * mat) exquo (rbden ^ p)) :: Mat
596
597    wildPrimes(factoredDisc, n) ==
598      -- returns a list of the primes <=n which divide factoredDisc to a
599      -- power greater than 1
600      ans : List I := empty()
601      for f in factorList(factoredDisc) repeat
602        if f.exponent > 1 and f.factor <= n then ans := concat(f.factor, ans)
603      ans
604
605    tameProduct(factoredDisc, n) ==
606      -- returns the product of the primes > n which divide factoredDisc
607      -- to a power greater than 1
608      ans : I := 1
609      for f in factorList(factoredDisc) repeat
610        if f.exponent > 1 and f.factor > n then ans := f.factor * ans
611      ans
612
613    integralBasis() ==
614      traceMat := traceMatrix()$F; n := rank()$F
615      disc := determinant traceMat  -- discriminant of current order
616      disc0 := disc                 -- this is disc(F)
617      factoredDisc := factor(disc0)$IntegerFactorizationPackage(Integer)
618      wilds := wildPrimes(factoredDisc, n)
619      sing := tameProduct(factoredDisc, n)
620      runningRb := scalarMatrix(n, 1); runningRbinv := scalarMatrix(n, 1)
621      -- runningRb    = basis matrix of current order
622      -- runningRbinv = inverse basis matrix of current order
623      -- these are wrt the original basis for F
624      runningRbden : I := 1
625      -- runningRbden = denominator for current basis matrix
626      (sing = 1) and empty? wilds => [runningRb, runningRbden, runningRbinv]
627      -- id = basis matrix of the ideal (p-radical) wrt current basis
628      matrixOut : Mat := scalarMatrix(n, 0)
629      for p in wilds repeat
630        lb := iWildLocalIntegralBasis(matrixOut, disc, p)
631        rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
632        disc := lb.discr
633        -- update 'running integral basis' if newly computed
634        -- local integral basis is non-trivial
635        if sizeLess?(1, rbden) then
636          mat := vertConcat(rbden * runningRb, runningRbden * rb)
637          runningRbden := runningRbden * rbden
638          runningRb := squareTop rowEchelon(mat, runningRbden)
639          runningRbinv := UpTriBddDenomInv(runningRb, runningRbden)
640      lb := iTameLocalIntegralBasis(traceMat, disc, sing)
641      rb := lb.basis; rbinv := lb.basisInv; rbden := lb.basisDen
642      disc := lb.discr
643      -- update 'running integral basis' if newly computed
644      -- local integral basis is non-trivial
645      if sizeLess?(1, rbden) then
646        mat := vertConcat(rbden * runningRb, runningRbden * rb)
647        runningRbden := runningRbden * rbden
648        runningRb := squareTop rowEchelon(mat, runningRbden)
649        runningRbinv := UpTriBddDenomInv(runningRb, runningRbden)
650      [runningRb, runningRbden, runningRbinv]
651
652    localIntegralBasis p ==
653      traceMat := traceMatrix()$F; n := rank()$F
654      disc := determinant traceMat  -- discriminant of current order
655      (disc exquo (p*p)) case "failed" =>
656        [scalarMatrix(n, 1), 1, scalarMatrix(n, 1)]
657      lb :=
658        p > rank()$F =>
659          iTameLocalIntegralBasis(traceMat, disc, p)
660        iWildLocalIntegralBasis(scalarMatrix(n, 0), disc, p)
661      [lb.basis, lb.basisDen, lb.basisInv]
662
663    iTameLocalIntegralBasis(traceMat, disc, sing) ==
664      n := rank()$F; disc0 := disc
665      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
666      -- rb    = basis matrix of current order
667      -- rbinv = inverse basis matrix of current order
668      -- these are wrt the original basis for F
669      rbden : I := 1; index : I := 1; oldIndex : I := 1
670      -- rbden = denominator for current basis matrix
671      -- id = basis matrix of the ideal (p-radical) wrt current basis
672      tfm := traceMat
673      repeat
674        -- compute the p-radical = p-trace-radical
675        idinv := transpose squareTop rowEchelon(tfm, sing)
676        -- [u1, .., un] are the coordinates of an element of the p-radical
677        -- iff [u1, .., un] * idinv is in p * Z^n
678        id := rowEchelon LowTriBddDenomInv(idinv, sing)
679        -- id = basis matrix of the p-radical
680        idinv := UpTriBddDenomInv(id, sing)
681        -- id * idinv = sing * identity
682        -- no need to check for inseparability in this case
683        rbinv := idealiser(id * rb, rbinv * idinv, sing * rbden)
684        index := diagonalProduct rbinv
685        rb := rowEchelon LowTriBddDenomInv(rbinv, sing * rbden)
686        g := matrixGcd(rb, sing, n)
687        if sizeLess?(1, g) then rb := (rb exquo g) :: Mat
688        rbden := rbden * (sing quo g)
689        rbinv := UpTriBddDenomInv(rb, rbden)
690        disc := disc0 quo (index * index)
691        indexChange := index quo oldIndex; oldIndex := index
692        (indexChange = 1) => return [rb, rbden, rbinv, disc]
693        tfm := ((rb * traceMat * transpose rb) exquo (rbden * rbden)) :: Mat
694
695    iWildLocalIntegralBasis(matrixOut, disc, p) ==
696      n := rank()$F; disc0 := disc
697      rb := scalarMatrix(n, 1); rbinv := scalarMatrix(n, 1)
698      -- rb    = basis matrix of current order
699      -- rbinv = inverse basis matrix of current order
700      -- these are wrt the original basis for F
701      rbden : I := 1; index : I := 1; oldIndex : I := 1
702      -- rbden = denominator for current basis matrix
703      -- id = basis matrix of the ideal (p-radical) wrt current basis
704      p2 := p * p; lp := leastPower(p::NNI, n)
705      repeat
706        tfm := frobMatrix(rb, rbinv, rbden, p::NNI) ^ lp
707        -- compute Rp = p-radical
708        idinv := transpose squareTop rowEchelon(tfm, p)
709        -- [u1, .., un] are the coordinates of an element of Rp
710        -- iff [u1, .., un] * idinv is in p * Z^n
711        id := rowEchelon LowTriBddDenomInv(idinv, p)
712        -- id = basis matrix of the p-radical
713        idinv := UpTriBddDenomInv(id, p)
714        -- id * idinv = p * identity
715        -- no need to check for inseparability in this case
716        rbinv := idealiser(id * rb, rbinv * idinv, p * rbden)
717        index := diagonalProduct rbinv
718        rb := rowEchelon LowTriBddDenomInv(rbinv, p * rbden)
719        if divideIfCan!(rb, matrixOut, p, n) = 1
720          then rb := matrixOut
721          else rbden := p * rbden
722        rbinv := UpTriBddDenomInv(rb, rbden)
723        indexChange := index quo oldIndex; oldIndex := index
724        disc := disc quo (indexChange * indexChange)
725        (indexChange = 1) or gcd(p2, disc) ~= p2 =>
726          return [rb, rbden, rbinv, disc]
727
728    discriminant() ==
729      disc := determinant traceMatrix()$F
730      intBas := integralBasis()
731      rb := intBas.basis; rbden := intBas.basisDen
732      index := ((rbden ^ rank()$F) exquo (determinant rb)) :: Integer
733      (disc exquo (index * index)) :: Integer
734
735--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
736--All rights reserved.
737--
738--Redistribution and use in source and binary forms, with or without
739--modification, are permitted provided that the following conditions are
740--met:
741--
742--    - Redistributions of source code must retain the above copyright
743--      notice, this list of conditions and the following disclaimer.
744--
745--    - Redistributions in binary form must reproduce the above copyright
746--      notice, this list of conditions and the following disclaimer in
747--      the documentation and/or other materials provided with the
748--      distribution.
749--
750--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
751--      names of its contributors may be used to endorse or promote products
752--      derived from this software without specific prior written permission.
753--
754--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
755--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
756--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
757--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
758--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
759--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
760--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
761--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
762--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
763--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
764--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
765