diff --git a/lib/matrix.gi b/lib/matrix.gi index 25a7226a65..8a8292d7b2 100644 --- a/lib/matrix.gi +++ b/lib/matrix.gi @@ -428,7 +428,7 @@ end ); ############################################################################# ## -#F Matrix_OrderPolynomialInner( , , , ) +#F Matrix_OrderPolynomialInner( , , , [, [, ]] ) ## ## Returns the coefficients of the order polynomial of at ## modulo . No conversions are attempted on or @@ -441,14 +441,34 @@ end ); ## The result, and any vectors added to are compressed ## and immutable ## +## If the optional argument is given, the computation is abandoned +## as soon as it is clear that the order polynomial has degree larger than +## , and 'fail' is returned. In that case has still +## been extended by the images computed so far, so a caller that wants to +## reuse it must pass in a copy. +## +## If the optional argument is given, the images of under +## powers of are appended to it, up to and including the first one +## lying in the span of its predecessors and . So for an order +## polynomial of degree e it receives e+1 vectors. +## #N In characteristic zero, or for structured sparse matrices, the naive #N Gaussian elimination here may not be optimal ## #N Shift to using ClearRow once we have kernel methods that give a #N performance benefit ## -BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs) - local d, w, p, one, zero, zeroes, piv, pols, x; +BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs, bound... ) + local d, w, p, one, zero, zeroes, piv, pols, x, krylov; + krylov := fail; + if 1 < Length(bound) then + krylov := bound[2]; + fi; + if Length(bound) = 0 then + bound := infinity; + else + bound := bound[1]; + fi; Info(InfoMatrix,2,"Order Polynomial Inner on ",NrRows(mat), " x ",NrCols(mat)," matrix over ",fld," with ", Number(vecs)," basis vectors already given"); @@ -464,6 +484,9 @@ BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs) # when we succeed, we know the order polynomial repeat + if krylov <> fail then + Add(krylov, vec); + fi; w := ShallowCopy(vec); p := ShallowCopy(zeroes); Add(p,one); @@ -489,6 +512,11 @@ BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs) # if piv <=d then + # the order polynomial has degree at least Length(zeroes)+1 + if Length(zeroes) >= bound then + Info(InfoMatrix,2,"Order Polynomial exceeds bound ",bound); + return fail; + fi; x := Inverse(w[piv]); MultVector(p, x); MakeImmutable(p); @@ -4760,24 +4788,170 @@ end); ## return Value(pol, mat); ## end); ## +# helper function for POW_MAT_INT: spin under , extending the state +# by the resulting Krylov chain, and return the coefficients of its order +# polynomial modulo what already spans. +# +# collects the chains in st.t, the positions where they end in st.ends, +# the image of the last vector of each chain in st.images, and the product of +# the order polynomials -- the characteristic polynomial of , once the +# chains span everything -- in st.cp. +BindGlobal("POW_MAT_INT_SPIN", function(f, m, st, vec) + local krylov, op, e; + krylov := []; + op := Matrix_OrderPolynomialInner(f, m, vec, st.vecs, infinity, krylov); + e := Length(op) - 1; + if 0 < e then + Append(st.t, krylov{[1..e]}); + Add(st.ends, Length(st.t)); + # the last image is a linear combination of the vectors collected so far + Add(st.images, Last(krylov)); + if st.cp = fail then + st.cp := op; + else + st.cp := ProductCoeffs(st.cp, op); + fi; + fi; + return op; +end); + +# helper function for POW_MAT_INT: complete the spinning begun in and +# return [ t, t^-1, mm ], where t is a base change matrix such that +# mm := t*m*t^-1 is block triangular with companion matrices along the +# diagonal. +# +# The rows v_1, ..., v_d of t are the spun Krylov chains, so that within one +# chain we have v_{i+1} = v_i*m. Now the i-th row of mm is the coordinate +# vector of v_i*m with respect to v_1, ..., v_d, hence it is the standard +# basis vector e_{i+1} for all i inside a chain, and only at the end of a +# chain is there anything to compute. So instead of multiplying out t*m*t^-1, +# which costs two matrix multiplications, we assemble mm from the images the +# spinning has produced anyway, using one vector-matrix product per chain. +BindGlobal("POW_MAT_INT_TRAFO", function(f, m, st) + local d, i, t, ti, mm, j; + d := NrRows(m); + # Spin up standard basis vectors, created one at a time as they are needed, + # until they span the whole space. Stopping there matters: once the basis is + # complete, every further vector would still be reduced against all of it, + # which for a cyclic matrix is as much work again as the spinning itself. + i := 0; + while Length(st.t) < d do + i := i + 1; + POW_MAT_INT_SPIN(f, m, st, StandardBasisVector(d, m, i)); + od; + Assert(2, Length(st.cp) = d+1); + t := Matrix(st.t, m); + ti := t^-1; + # all rows but those ending a chain are standard basis vectors + mm := List([2..d], k -> StandardBasisVector(d, m, k)); + for j in [1..Length(st.ends)] do + mm[st.ends[j]] := st.images[j] * ti; + od; + mm := Matrix(mm, m); + return [ t, ti, mm ]; +end); + +# helper function for POW_MAT_INT: evaluate at . +# Compared to the standard method, we avoid some zero or identity matrices +# and we multiply with mat from left to take advantage of sparseness of mat +BindGlobal("POW_MAT_INT_VALUE", function(pol, mat) + local f, c, i, val, j; + f := CoefficientsOfLaurentPolynomial(pol); + c := f[1]; + i := Length(c); + if i = 0 then + return 0*mat; + fi; + if i = 1 then + val := POW_OBJ_INT(mat, f[2]); + return c[1] * val; + fi; + val := c[i] * mat; + if IsMatrixObj(val) and not IsRowListMatrix(val) then + # no access to the rows, so mutability is that of the matrix itself + if not IsMutable(val) then + val := MutableCopyMatrix(val); + fi; + elif not IsMutable(val[1]) then + val := MutableCopyMatrix(val); + fi; + i := i-1; + for j in [1..NrRows(mat)] do + val[j,j] := val[j,j]+c[i]; + od; + while 1 < i do + val := mat * val; + i := i - 1; + for j in [1..NrRows(mat)] do + val[j,j] := val[j,j]+c[i]; + od; + od; + if 0 <> f[2] then + val := val * POW_OBJ_INT(mat, f[2]); + fi; + return val; +end); + +# helper function for POW_MAT_INT: a random vector in the row space of . +BindGlobal("POW_MAT_INT_RANDOMVEC", function(f, mat) + local vec, i; + # a fixed vector has systematic bad cases: the all ones vector is fixed by + # every permutation matrix. ZeroVector to match the representation of . + vec := ZeroVector(NrCols(mat), mat); + for i in [1..Length(vec)] do + vec[i] := Random(f); + od; + MakeImmutable(vec); + return vec; +end); + +# helper function for POW_MAT_INT: does the polynomial with coefficient list +# annihilate ? Horner, costing Length()-2 matrix +# multiplications. +BindGlobal("POW_MAT_INT_ANNIHILATES", function(coeffs, mat) + local i, val, j; + i := Length(coeffs); + val := coeffs[i] * mat; + if IsMatrixObj(val) and not IsRowListMatrix(val) then + # no access to the rows, so mutability is that of the matrix itself + if not IsMutable(val) then + val := MutableCopyMatrix(val); + fi; + elif not IsMutable(val[1]) then + val := MutableCopyMatrix(val); + fi; + i := i-1; + for j in [1..NrRows(mat)] do + val[j,j] := val[j,j]+coeffs[i]; + od; + while 1 < i do + val := mat * val; + i := i - 1; + for j in [1..NrRows(mat)] do + val[j,j] := val[j,j]+coeffs[i]; + od; + od; + return IsZero(val); +end); + # next iteration, conjugate matrix such that it is often very sparse # (a companion matrix), could still be improved, maybe with kernel functions # for compact matrices (FL) BindGlobal("POW_MAT_INT", function(mat, n) - local d, k, limit, f, addb, trafo, value, t, ti, mm, pol, ind; + local d, k, limit, f, fam, st, op, e, pol, ind, t, ti, mm; d := NrRows(mat); # Decide between repeated squaring (POW_OBJ_INT, about Log2(n) matrix - # multiplications) and the method below, which has a considerable fixed - # overhead (base change, characteristic polynomial, about d matrix - # multiplications) but afterwards only needs about Log2(n) polynomial - # multiplications modulo the characteristic polynomial, which are much - # cheaper than matrix multiplications when d is large. + # multiplications) and the methods below, which have a considerable fixed + # overhead (base change, characteristic or minimal polynomial, about d + # matrix multiplications) but afterwards only need about Log2(n) polynomial + # multiplications modulo that polynomial, which are much cheaper than + # matrix multiplications when d is large. # The break even points below were determined experimentally, on the basis # that both costs grow linearly in Log2(n) for a fixed matrix; see the # discussion in https://github.com/gap-system/gap/pull/6293 for details. # They depend on the representation: for compressed matrices over small # finite fields multiplication is very fast compared to the (partially - # interpreted) overhead of the method below, and for compressed matrices + # interpreted) overhead of the methods below, and for compressed matrices # over GF(2) it is so fast that the overhead only pays off for huge # exponents. k := LogInt(n, 2); @@ -4795,131 +4969,60 @@ BindGlobal("POW_MAT_INT", function(mat, n) if k < limit then return POW_OBJ_INT(mat, n); fi; - # the method below requires the entries to lie in a field + # the methods below require the entries to lie in a field f := DefaultFieldOfMatrix(mat); if f = fail or not IsField(f) then return POW_OBJ_INT(mat, n); fi; - # helper function to build up a semi-echelon basis - addb := function(seb, v) - local rows, pivots, len, vv, c, pos, i; - rows := seb.vectors; - pivots := seb.pivots; - len := Length(rows); - vv := ShallowCopy(v); - for i in [1..len] do - c := vv[pivots[i]]; - if not IsZero(c) then - AddRowVector(vv, rows[i], -c); - fi; - od; - pos := PositionNonZero(vv); - if pos <= Length(vv) then - if not IsOne(vv[pos]) then - vv := vv/vv[pos]; - fi; - Add(rows, vv); - Add(pivots, pos); - seb.heads[pos] := len + 1; - return true; - else - return false; - fi; - end; - # This computes a base change matrix t such that mm := t*m*t^-1 is block - # triangular with companion matrices along the diagonal, and returns the - # triple [ t, t^-1, mm ]. + + # Any polynomial annihilating mat can serve as the modulus, not just the + # characteristic one. If mat has one of degree e, both the Log2(n) + # polynomial multiplications and the final evaluation work with degree e + # instead of degree d. Verifying a candidate costs e-1 of the matrix + # multiplications that the evaluation then saves d-e of, so it is worth + # trying as soon as 2*e <= d. # - # The rows v_1, ..., v_d of t are obtained by spinning up standard basis - # vectors, so that within one such Krylov chain we have v_{i+1} = v_i*m. - # Now the i-th row of mm is the coordinate vector of v_i*m with respect to - # v_1, ..., v_d, hence it is the standard basis vector e_{i+1} for every i - # inside a chain, and only at the end of a chain is there anything to - # compute -- and the image needed there is exactly the vector on which the - # spinning stopped. So rather than multiplying out t*m*t^-1, which costs - # two matrix multiplications, we assemble mm from what the spinning has - # produced anyway, using one vector-matrix product per chain. - trafo := function(m) - local d, b, t, r, a, ends, images, i, ti, mm, j; - d := NrRows(m); - b := rec(vectors := [], pivots := [], heads := []); - t := []; - ends := []; - images := []; - # Spin up standard basis vectors, created one at a time as they are - # needed, until they span the whole space. Stopping as soon as that - # happens matters: any further vector would still be reduced against the - # complete basis, which for a cyclic matrix amounts to as much work again - # as the spinning itself. - # maybe better start with a random vector? - i := 0; - while Length(t) < d do - i := i + 1; - a := StandardBasisVector(d, m, i); - r := addb(b,a); - if r = true then - repeat - Add(t, a); - a := a*m; - r := addb(b,a); - until r <> true; - # a is the image of the last vector of this chain, and is a linear - # combination of the vectors collected so far - Add(ends, Length(t)); - Add(images, a); - fi; - od; - t := Matrix(t, m); - ti := t^-1; - # all rows but those ending a chain are standard basis vectors - mm := List([2..d], k -> StandardBasisVector(d, m, k)); - for j in [1..Length(ends)] do - mm[ends[j]] := images[j] * ti; - od; - return [ t, ti, Matrix(mm, m) ]; - end; - # compared to standard method, we avoid some zero or identity matrices - # and we multiply with mat from left to take advantage of sparseness of mat - value := function(pol, mat) - local f, c, i, val, j; - f := CoefficientsOfLaurentPolynomial(pol); - c := f[1]; - i := Length(c); - if i = 0 then - return 0*mat; - fi; - if i = 1 then - val := POW_OBJ_INT(mat, f[2]); - return c[1] * val; - fi; - val := c[i] * mat; - if not IsMutable(val[1]) then - val := MutableCopyMatrix(val); - fi; - i := i-1; - for j in [1..NrRows(mat)] do - val[j,j] := val[j,j]+c[i]; - od; - while 1 < i do - val := mat * val; - i := i - 1; - for j in [1..NrRows(mat)] do - val[j,j] := val[j,j]+c[i]; - od; - od; - if 0 <> f[2] then - val := val * POW_OBJ_INT(mat, f[2]); - fi; - return val; - end; - t := trafo(mat); + # Start the base change at a random vector rather than at e_1. Its order + # polynomial divides the minimal polynomial, and equals it unless the vector + # was unlucky; computing the minimal polynomial outright instead costs more + # than this whole method. Nothing is staked on the guess: for a cyclic + # matrix -- which a random matrix is with probability about 1-1/q^3 -- that + # one chain already is the whole base change, and its order polynomial the + # characteristic polynomial, so both methods share the same first step. + st := rec(vecs := [], t := [], ends := [], images := [], cp := fail); + op := POW_MAT_INT_SPIN(f, mat, st, POW_MAT_INT_RANDOMVEC(f, mat)); + e := Length(op) - 1; + fam := ElementsFamily(FamilyObj(f)); + # Verifying and evaluating at mat costs 2*e dense multiplications. The base + # change below makes both sparse and thus much cheaper, but has to be + # completed first. So use mat directly only while 2*e stays below the cost + # of that completion, measured in dense multiplications -- which over GF(2) + # is much further. + if IsGF2MatrixRep(mat) then + limit := 60; + else + limit := 6; + fi; + if 0 < e and e <= limit and POW_MAT_INT_ANNIHILATES(op, mat) then + pol := UnivariatePolynomialByCoefficients(fam, op, 1); + ind := IndeterminateOfUnivariateRationalFunction(pol); + return POW_MAT_INT_VALUE(PowerMod(ind, n, pol), mat); + fi; + + t := POW_MAT_INT_TRAFO(f, mat, st); ti := t[2]; mm := t[3]; t := t[1]; - pol := CharacteristicPolynomial(mm); + # conjugate matrices have the same annihilating polynomials, and verifying + # at the sparse mm is much cheaper than at mat + if 0 < e and 2*e <= d and POW_MAT_INT_ANNIHILATES(op, mm) then + pol := UnivariatePolynomialByCoefficients(fam, op, 1); + else + pol := UnivariatePolynomialByCoefficients(fam, st.cp, 1); + fi; ind := IndeterminateOfUnivariateRationalFunction(pol); pol := PowerMod(ind, n, pol); - mm := value(pol, mm); + mm := POW_MAT_INT_VALUE(pol, mm); return ti * mm * t; end); diff --git a/tst/testinstall/matrix.tst b/tst/testinstall/matrix.tst index 9876343a2d..59d32325a9 100644 --- a/tst/testinstall/matrix.tst +++ b/tst/testinstall/matrix.tst @@ -217,3 +217,142 @@ gap> NumberColumns(m); 2 gap> DimensionsMat(m); [ 2, 2 ] + +# +# powering matrices by large exponents: POW_MAT_INT dispatches between +# repeated squaring, reduction modulo the characteristic polynomial and +# reduction modulo the minimal polynomial +# +gap> conj := function(m, F) +> local d, t, i; +> d := NrRows(m); +> t := List(IdentityMat(d, F), ShallowCopy); +> for i in [1..d-1] do t[i][i+1] := One(F); od; +> t := ImmutableMatrix(F, t); +> return ImmutableMatrix(F, t * m * t^-1); +> end;; +gap> blockdiag := function(b, nb, F) +> local e, d, m, i; +> e := NrRows(b); d := e*nb; m := NullMat(d, d, F); +> for i in [1..nb] do m{[(i-1)*e+1..i*e]}{[(i-1)*e+1..i*e]} := b; od; +> return conj(ImmutableMatrix(F, m), F); +> end;; +gap> comp := function(coeffs, F) +> local d, m, i; +> d := Length(coeffs) - 1; m := NullMat(d, d, F); +> for i in [1..d-1] do m[i][i+1] := One(F); od; +> for i in [1..d] do m[d][i] := -coeffs[i]/Last(coeffs); od; +> return ImmutableMatrix(F, m); +> end;; +gap> check := m -> ForAll([64, 100, 300], +> k -> POW_MAT_INT(m, 2^k+1) = POW_OBJ_INT(m, 2^k+1));; + +# the bounded variant of the spinning used to detect a small minimal polynomial +gap> F := GF(5);; b := comp([1,1,1,1,1,1]*One(F), F);; +gap> DegreeOfLaurentPolynomial(MinimalPolynomial(F, b, 1)); +5 +gap> v := ImmutableVector(F, One(F)*[1,0,0,0,0]);; +gap> Length(Matrix_OrderPolynomialInner(F, b, v, [])) - 1; +5 +gap> Matrix_OrderPolynomialInner(F, b, v, [], 4); +fail +gap> Length(Matrix_OrderPolynomialInner(F, b, v, [], 5)) - 1; +5 + +# generic matrices: minimal polynomial equals characteristic polynomial, +# so the probe fails and the characteristic polynomial is used +gap> m := conj(comp(One(GF(2))*[1,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1], GF(2)), GF(2));; +gap> DegreeOfLaurentPolynomial(MinimalPolynomial(m)) = NrRows(m); +true +gap> check(m); +true +gap> m := conj(comp(One(GF(5))*[2,1,0,0,3,0,0,1,0,0,0,4,0,0,0,0,1], GF(5)), GF(5));; +gap> DegreeOfLaurentPolynomial(MinimalPolynomial(m)) = NrRows(m); +true +gap> check(m); +true + +# small minimal polynomial, evaluated at the matrix itself +gap> m := blockdiag(comp(One(GF(5))*[3,1,4,1], GF(5)), 20, GF(5));; +gap> [ NrRows(m), DegreeOfLaurentPolynomial(MinimalPolynomial(m)) ]; +[ 60, 3 ] +gap> check(m); +true +gap> m := blockdiag(comp(One(GF(251))*[7,1,4,1,1], GF(251)), 15, GF(251));; +gap> [ NrRows(m), DegreeOfLaurentPolynomial(MinimalPolynomial(m)) ]; +[ 60, 4 ] +gap> check(m); +true + +# small minimal polynomial, evaluated at the block companion form +gap> m := blockdiag(comp(One(GF(5))*[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1], GF(5)), 8, GF(5));; +gap> [ NrRows(m), DegreeOfLaurentPolynomial(MinimalPolynomial(m)) ]; +[ 120, 15 ] +gap> check(m); +true + +# a polynomial annihilates a matrix iff it is a multiple of its minimal +# polynomial; a proper divisor of the latter does not +gap> b := comp(One(GF(5))*[1,1,1,1,1,1], GF(5));; +gap> POW_MAT_INT_ANNIHILATES(One(GF(5))*[1,1,1,1,1,1], b); +true +gap> POW_MAT_INT_ANNIHILATES(One(GF(5))*[1,1,1,1,1,1,1], b); +false +gap> POW_MAT_INT_ANNIHILATES(One(GF(5))*[4,1], b); +false + +# a tiny minimal polynomial next to a large dimension: here computing the +# minimal polynomial outright would cost far more than the whole method +gap> m := ImmutableMatrix(GF(2), PermutationMat((1,2,3,4), 120, GF(2)));; +gap> [ NrRows(m), DegreeOfLaurentPolynomial(MinimalPolynomial(m)) ]; +[ 120, 4 ] +gap> check(m); +true + +# scalar and unipotent matrices +gap> m := ImmutableMatrix(GF(9), Z(9)*IdentityMat(20, GF(9)));; +gap> DegreeOfLaurentPolynomial(MinimalPolynomial(m)); +1 +gap> check(m); +true +gap> m := List(IdentityMat(40, GF(2)), ShallowCopy);; +gap> for i in [1..20] do m[i][i+20] := One(GF(2)); od; +gap> m := ImmutableMatrix(GF(2), m);; +gap> DegreeOfLaurentPolynomial(MinimalPolynomial(m)); +2 +gap> check(m); +true + +# entries not in a field, or in a ring of characteristic 0: repeated squaring. +# The exponents are large enough that the polynomial methods would be used if +# the base domain were a field. +gap> POW_MAT_INT([[1,1],[0,1]], 2^100+1) = [[1, 2^100+1],[0,1]]; +true +gap> m := One(Integers mod 6) * [[1,1],[0,1]];; +gap> POW_MAT_INT(m, 2^4100+1) = POW_OBJ_INT(m, 2^4100+1); +true +gap> m := List([1..4], i -> List([1..4], j -> Random(Integers mod 8)));; +gap> IsField(DefaultFieldOfMatrix(m)); +false +gap> POW_MAT_INT(m, 2^2050+1) = POW_OBJ_INT(m, 2^2050+1); +true + +# matrix objects which do not offer access to their rows are handled as well, +# both by the characteristic and by the minimal polynomial method +gap> m := Matrix(IsGenericMatrixRep, GF(257), +> List([1..16], i -> List([1..16], j -> Random(GF(257)))));; +gap> IsRowListMatrix(m); +false +gap> POW_MAT_INT(m, 2^514+1) = POW_OBJ_INT(m, 2^514+1); +true +gap> m := Matrix(IsGenericMatrixRep, GF(257), +> Unpack(blockdiag(comp(One(GF(257))*[3,1,4,1], GF(257)), 8, GF(257))));; +gap> DegreeOfLaurentPolynomial(MinimalPolynomial(GF(257), m, 1)); +3 +gap> POW_MAT_INT(m, 2^514+1) = POW_OBJ_INT(m, 2^514+1); +true + +# matrix objects +gap> m := Matrix(IsPlistMatrixRep, GF(5), Unpack(blockdiag(comp(One(GF(5))*[3,1,4,1], GF(5)), 20, GF(5))));; +gap> POW_MAT_INT(m, 2^500+1) = POW_OBJ_INT(m, 2^500+1); +true