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