From 40b96b5a524ba705251172c030ac9662e7c9dd06 Mon Sep 17 00:00:00 2001 From: Max Horn Date: Tue, 4 Aug 2026 21:29:40 +0200 Subject: [PATCH 1/3] Speed up powering a matrix by a large exponent when its minimal polynomial is small Powering a matrix by a large exponent reduces x^n modulo the characteristic polynomial, after a base change making the matrix block companion. If the minimal polynomial has degree e rather than d, both the Log2(n) polynomial multiplications and the final evaluation work with degree e, which is worth a factor of 2 to 60. Computing the minimal polynomial of a dense matrix costs about as much as the whole characteristic polynomial method, so probe first with a single spun random vector, abandoned once its order polynomial grows past the bound where the minimal polynomial stops being worth using. That order polynomial divides the minimal polynomial, so abandoning settles the question, and the probe does not show up in the timings. The break even points were determined over a grid of dimensions from 60 to 2000, representations GF(2), GF(5), GF(251), GF(257) and GF(3^10), minimal polynomial degrees from d/30 up to d, and exponents from 2^20 to 2^1000. They depend mostly on the representation: 6*e <= d for compressed matrices, 2*e <= d for the rest, where matrix multiplication is far more expensive than polynomial arithmetic. The choice is within a factor of 1.07 of the best of the three methods everywhere measured, and at least 1.5 times faster than before on a third of the cells. Matrix objects without access to their rows, such as those in IsGenericMatrixRep, now work here at all: the helper evaluating a polynomial at a matrix tested mutability by looking at the first row. With the preceding commit they reach the generic methods for the characteristic and minimal polynomial again, so powering them uses the polynomial method rather than failing with "row access unsupported". For this, Matrix_OrderPolynomialInner gained an optional degree bound, and the helpers of POW_MAT_INT moved out of its body to be shared. This commit was prepared with assistance from the AI tool Claude Code (benchmarking, analysis and drafting of the implementation and tests). Co-Authored-By: Claude Opus 5 (1M context) --- lib/matrix.gi | 335 +++++++++++++++++++++++-------------- tst/testinstall/matrix.tst | 121 ++++++++++++++ 2 files changed, 335 insertions(+), 121 deletions(-) diff --git a/lib/matrix.gi b/lib/matrix.gi index 25a7226a65..f709a2afa8 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,25 @@ 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. +## #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) +BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs, bound... ) local d, w, p, one, zero, zeroes, piv, pols, x; + 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"); @@ -489,6 +500,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 +4776,164 @@ end); ## return Value(pol, mat); ## end); ## +# helper function for POW_MAT_INT: build up a semi-echelon basis +BindGlobal("POW_MAT_INT_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); + +# helper function for POW_MAT_INT: compute a base change matrix t such that +# mm := t*m*t^-1 is block triangular with companion matrices along the +# diagonal, and return [ t, t^-1, mm ]. +# +# The rows v_1, ..., v_d of t are built 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 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(m) + local b, t, r, a, ends, images, d, 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 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. + # maybe better start with a random vector? + i := 0; + while Length(t) < d do + i := i + 1; + a := StandardBasisVector(d, m, i); + r := POW_MAT_INT_ADDB(b,a); + if r = true then + repeat + Add(t, a); + a := a*m; + r := POW_MAT_INT_ADDB(b,a); + until r <> true; + # a is the image of the last vector of this chain, and depends on 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; + 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: spin a single random vector and return the +# degree of its order polynomial, or 'fail' if that degree exceeds . +# That order polynomial divides the minimal polynomial of , and equals it +# with high probability. Costs O( * d^2) field operations. +BindGlobal("POW_MAT_INT_PROBE", function(f, mat, bound) + local vec, i, op; + # 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); + op := Matrix_OrderPolynomialInner(f, mat, vec, [], bound); + if op = fail then + return fail; + fi; + return Length(op) - 1; +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, ratio, f, 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 +4951,68 @@ 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); + + # If the minimal polynomial has degree e, both the Log2(n) polynomial + # multiplications and the final evaluation work with degree e instead of + # degree d, which is worth a factor of 2 to 60. The break even point depends + # on the representation: matrix multiplication is cheap for compressed + # matrices and expensive for everything else, so those need 6*e <= d and + # 2*e <= d respectively. + # Computing the minimal polynomial of a dense matrix costs about as much as + # the whole characteristic polynomial method, so first probe with a single + # random vector, which is a fraction of one spinning: its order polynomial + # divides the minimal polynomial, so giving up on it settles the question. + if IsGF2MatrixRep(mat) or Is8BitMatrixRep(mat) then + ratio := 6; + else + ratio := 2; + fi; + e := fail; + if d >= ratio then + e := POW_MAT_INT_PROBE(f, mat, QuoInt(d, ratio)); + fi; + if e <> fail then + pol := MinimalPolynomial(f, mat, 1); + e := DegreeOfLaurentPolynomial(pol); + # the probe only gives a lower bound for e, so check again; if the minimal + # polynomial is large after all, fall through to the method below + if ratio*e <= d then + ind := IndeterminateOfUnivariateRationalFunction(pol); + pol := PowerMod(ind, n, pol); + # Evaluating at mat costs e dense multiplications, the base change + # below turns those into sparse ones but has to be computed. So use mat + # directly while e stays below the cost of that base change, measured in + # dense multiplications -- which over GF(2) is much further. + if IsGF2MatrixRep(mat) then + limit := 125; + else + limit := 12; fi; - od; - pos := PositionNonZero(vv); - if pos <= Length(vv) then - if not IsOne(vv[pos]) then - vv := vv/vv[pos]; + if e <= limit then + return POW_MAT_INT_VALUE(pol, mat); fi; - Add(rows, vv); - Add(pivots, pos); - seb.heads[pos] := len + 1; - return true; - else - return false; + t := POW_MAT_INT_TRAFO(mat); + ti := t[2]; + mm := t[3]; + t := t[1]; + mm := POW_MAT_INT_VALUE(pol, mm); + return ti * mm * t; 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 ]. - # - # 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); + fi; + + t := POW_MAT_INT_TRAFO(mat); ti := t[2]; mm := t[3]; t := t[1]; pol := CharacteristicPolynomial(mm); 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..ec6d3ed132 100644 --- a/tst/testinstall/matrix.tst +++ b/tst/testinstall/matrix.tst @@ -217,3 +217,124 @@ 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 + +# 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 From 8ab7a7eb70a67e6d3f78c80358674207f2c2e2bd Mon Sep 17 00:00:00 2001 From: Max Horn Date: Mon, 17 Aug 2026 09:42:31 +0200 Subject: [PATCH 2/3] Verify a spun order polynomial instead of computing the minimal polynomial POW_MAT_INT probed for a small minimal polynomial and, when the probe succeeded, called MinimalPolynomial. That call dominated everything else: on a 1000 x 1000 permutation matrix over GF(2) with minimal polynomial of degree 4 it took 2.9 s, while reducing modulo the characteristic polynomial -- the method it was supposed to beat -- took 0.29 s. For d = 504 the minimal polynomial method was between 1.15 and 5.4 times slower over GF(2), and up to 1.4 times slower over GF(251). Any annihilating polynomial serves as the modulus, not just the minimal one, so use the probe's own order polynomial once Horner confirms that it kills the matrix, at a cost of e-1 matrix multiplications. Do that at the block companion conjugate rather than at the matrix itself unless e is tiny; the base change is needed by the fallback anyway, so a failed verification wastes nothing. For d = 504 and e from 2 to d/6 this is now faster than the characteristic polynomial method on every measured cell: 1.4 to 115 times over GF(2), 1.8 to 9 times over GF(251). On random matrices, which are cyclic with probability about 1-1/q^3 and hence never take this path, the probe costs 0.5% to 2% of the total. The break even points 6*e <= d and 2*e <= d are inherited from the previous cost model and are now conservative: over GF(2) and GF(251) the order polynomial method still wins at e = d/2. Prepared with assistance from the AI tool Claude Code (benchmarking, analysis and drafting of the implementation and tests). Co-Authored-By: Claude Opus 5 --- lib/matrix.gi | 122 +++++++++++++++++++++++-------------- tst/testinstall/matrix.tst | 18 ++++++ 2 files changed, 95 insertions(+), 45 deletions(-) diff --git a/lib/matrix.gi b/lib/matrix.gi index f709a2afa8..21e3317752 100644 --- a/lib/matrix.gi +++ b/lib/matrix.gi @@ -4897,11 +4897,11 @@ BindGlobal("POW_MAT_INT_VALUE", function(pol, mat) end); # helper function for POW_MAT_INT: spin a single random vector and return the -# degree of its order polynomial, or 'fail' if that degree exceeds . -# That order polynomial divides the minimal polynomial of , and equals it -# with high probability. Costs O( * d^2) field operations. +# coefficients of its order polynomial, or 'fail' if its degree exceeds +# . That order polynomial divides the minimal polynomial of , and +# equals it with high probability. Costs O( * d^2) field operations. BindGlobal("POW_MAT_INT_PROBE", function(f, mat, bound) - local vec, i, op; + 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); @@ -4909,18 +4909,43 @@ BindGlobal("POW_MAT_INT_PROBE", function(f, mat, bound) vec[i] := Random(f); od; MakeImmutable(vec); - op := Matrix_OrderPolynomialInner(f, mat, vec, [], bound); - if op = fail then - return fail; + return Matrix_OrderPolynomialInner(f, mat, vec, [], bound); +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; - return Length(op) - 1; + 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, ratio, f, e, pol, ind, t, ti, mm; + local d, k, limit, ratio, f, op, e, pol, ind, t, ti, mm; d := NrRows(mat); # Decide between repeated squaring (POW_OBJ_INT, about Log2(n) matrix # multiplications) and the methods below, which have a considerable fixed @@ -4957,55 +4982,62 @@ BindGlobal("POW_MAT_INT", function(mat, n) return POW_OBJ_INT(mat, n); fi; - # If the minimal polynomial has degree e, both the Log2(n) polynomial - # multiplications and the final evaluation work with degree e instead of - # degree d, which is worth a factor of 2 to 60. The break even point depends - # on the representation: matrix multiplication is cheap for compressed - # matrices and expensive for everything else, so those need 6*e <= d and - # 2*e <= d respectively. - # Computing the minimal polynomial of a dense matrix costs about as much as - # the whole characteristic polynomial method, so first probe with a single - # random vector, which is a fraction of one spinning: its order polynomial - # divides the minimal polynomial, so giving up on it settles the question. + # 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. The break even point depends on the representation: + # matrix multiplication is cheap for compressed matrices and expensive for + # everything else, so those need 6*e <= d and 2*e <= d respectively. + # Spin a single random vector: its order polynomial divides the minimal + # polynomial, and equals it unless the vector was unlucky, which Horner + # settles for e-1 matrix multiplications. Computing the minimal polynomial + # outright instead costs more than this whole method. if IsGF2MatrixRep(mat) or Is8BitMatrixRep(mat) then ratio := 6; else ratio := 2; fi; - e := fail; + op := fail; if d >= ratio then - e := POW_MAT_INT_PROBE(f, mat, QuoInt(d, ratio)); + op := POW_MAT_INT_PROBE(f, mat, QuoInt(d, ratio)); fi; - if e <> fail then - pol := MinimalPolynomial(f, mat, 1); - e := DegreeOfLaurentPolynomial(pol); - # the probe only gives a lower bound for e, so check again; if the minimal - # polynomial is large after all, fall through to the method below - if ratio*e <= d then - ind := IndeterminateOfUnivariateRationalFunction(pol); - pol := PowerMod(ind, n, pol); - # Evaluating at mat costs e dense multiplications, the base change - # below turns those into sparse ones but has to be computed. So use mat - # directly while e stays below the cost of that base change, measured in - # dense multiplications -- which over GF(2) is much further. - if IsGF2MatrixRep(mat) then - limit := 125; - else - limit := 12; - fi; - if e <= limit then - return POW_MAT_INT_VALUE(pol, mat); - fi; + t := fail; + if op <> fail and 1 < Length(op) then + e := Length(op) - 1; + # 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 computed first; and it is needed by the fallback anyway, so nothing + # is lost if the verification fails. Use mat directly only while 2*e + # stays below the cost of the base change, measured in dense + # multiplications -- which over GF(2) is much further. + if IsGF2MatrixRep(mat) then + limit := 60; + else + limit := 6; + fi; + if limit < e then t := POW_MAT_INT_TRAFO(mat); - ti := t[2]; + # conjugate matrices have the same annihilating polynomials mm := t[3]; - t := t[1]; + else + mm := mat; + fi; + if POW_MAT_INT_ANNIHILATES(op, mm) then + pol := UnivariatePolynomialByCoefficients(ElementsFamily(FamilyObj(f)), + op, 1); + ind := IndeterminateOfUnivariateRationalFunction(pol); + pol := PowerMod(ind, n, pol); mm := POW_MAT_INT_VALUE(pol, mm); - return ti * mm * t; + if t = fail then + return mm; + fi; + return t[2] * mm * t[1]; fi; fi; - t := POW_MAT_INT_TRAFO(mat); + if t = fail then + t := POW_MAT_INT_TRAFO(mat); + fi; ti := t[2]; mm := t[3]; t := t[1]; diff --git a/tst/testinstall/matrix.tst b/tst/testinstall/matrix.tst index ec6d3ed132..59d32325a9 100644 --- a/tst/testinstall/matrix.tst +++ b/tst/testinstall/matrix.tst @@ -291,6 +291,24 @@ gap> [ NrRows(m), DegreeOfLaurentPolynomial(MinimalPolynomial(m)) ]; 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)); From 6eddee38cde21e7fbf74c14902958f88c4f94344 Mon Sep 17 00:00:00 2001 From: Max Horn Date: Mon, 17 Aug 2026 11:58:21 +0200 Subject: [PATCH 3/3] Start the base change at a random vector and keep what it spins POW_MAT_INT spun a random vector to guess the minimal polynomial, threw that spinning away, and then spun standard basis vectors from scratch to build the base change, after which CharacteristicPolynomial spun a third time. Frank Luebeck pointed out that for a cyclic matrix the first spinning already is the base change. So spin the random vector as the first Krylov chain of the base change itself. For a cyclic matrix -- which a random matrix is with probability about 1-1/q^3 -- that one chain spans everything, and its order polynomial is the characteristic polynomial; in general the product of the chains' order polynomials is. Either way the separate CharacteristicPolynomial call is gone, and the guess costs nothing even when it fails. Since nothing is staked on the guess any more, the break even point need no longer be fitted: 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. That replaces the measured 6*e <= d for compressed matrices, widening their range from d/6 to d/2. Timings for n = 2^500+1, milliseconds, best of 5, against the same matrices; e = d is the random matrix case: e=d/12 e=d/6 e=d/4 e=d/2 e=d GF(2) d=504 133 107 106 108 109 before 36 38 49 78 82 after GF(251) d=504 848 847 852 867 888 before 409 460 508 684 855 after GF(2) d=1008 596 627 569 579 589 before 140 211 270 483 491 after GF(251) d=1008 5300 5374 5441 5564 5695 before 3028 3331 3661 4756 5556 after Matrix_OrderPolynomialInner gained an optional argument handing back the Krylov chain it computes anyway, which lets POW_MAT_INT_ADDB -- a second implementation of the same spinning -- be dropped. Prepared with assistance from the AI tool Claude Code (benchmarking, analysis and drafting of the implementation). Co-Authored-By: Claude Opus 5 --- lib/matrix.gi | 210 ++++++++++++++++++++++---------------------------- 1 file changed, 94 insertions(+), 116 deletions(-) diff --git a/lib/matrix.gi b/lib/matrix.gi index 21e3317752..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 @@ -447,6 +447,11 @@ end ); ## 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 ## @@ -454,7 +459,11 @@ end ); #N performance benefit ## BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs, bound... ) - local d, w, p, one, zero, zeroes, piv, pols, x; + 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 @@ -475,6 +484,9 @@ BindGlobal( "Matrix_OrderPolynomialInner", function( fld, mat, vec, vecs, bound. # 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); @@ -4776,80 +4788,64 @@ end); ## return Value(pol, mat); ## end); ## -# helper function for POW_MAT_INT: build up a semi-echelon basis -BindGlobal("POW_MAT_INT_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]; +# 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; - Add(rows, vv); - Add(pivots, pos); - seb.heads[pos] := len + 1; - return true; - else - return false; fi; + return op; end); -# helper function for POW_MAT_INT: compute a base change matrix t such that +# 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, and return [ t, t^-1, mm ]. +# diagonal. # -# The rows v_1, ..., v_d of t are built 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 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(m) - local b, t, r, a, ends, images, d, i, ti, mm, j; +# 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); - 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 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. - # maybe better start with a random vector? i := 0; - while Length(t) < d do + while Length(st.t) < d do i := i + 1; - a := StandardBasisVector(d, m, i); - r := POW_MAT_INT_ADDB(b,a); - if r = true then - repeat - Add(t, a); - a := a*m; - r := POW_MAT_INT_ADDB(b,a); - until r <> true; - # a is the image of the last vector of this chain, and depends on the - # vectors collected so far - Add(ends, Length(t)); - Add(images, a); - fi; + POW_MAT_INT_SPIN(f, m, st, StandardBasisVector(d, m, i)); od; - t := Matrix(t, m); + 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(ends)] do - mm[ends[j]] := images[j] * ti; + 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 ]; @@ -4896,11 +4892,8 @@ BindGlobal("POW_MAT_INT_VALUE", function(pol, mat) return val; end); -# helper function for POW_MAT_INT: spin a single random vector and return the -# coefficients of its order polynomial, or 'fail' if its degree exceeds -# . That order polynomial divides the minimal polynomial of , and -# equals it with high probability. Costs O( * d^2) field operations. -BindGlobal("POW_MAT_INT_PROBE", function(f, mat, bound) +# 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 . @@ -4909,7 +4902,7 @@ BindGlobal("POW_MAT_INT_PROBE", function(f, mat, bound) vec[i] := Random(f); od; MakeImmutable(vec); - return Matrix_OrderPolynomialInner(f, mat, vec, [], bound); + return vec; end); # helper function for POW_MAT_INT: does the polynomial with coefficient list @@ -4945,7 +4938,7 @@ end); # (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, ratio, f, op, e, pol, ind, t, ti, mm; + 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 methods below, which have a considerable fixed @@ -4985,63 +4978,48 @@ BindGlobal("POW_MAT_INT", function(mat, n) # 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. The break even point depends on the representation: - # matrix multiplication is cheap for compressed matrices and expensive for - # everything else, so those need 6*e <= d and 2*e <= d respectively. - # Spin a single random vector: its order polynomial divides the minimal - # polynomial, and equals it unless the vector was unlucky, which Horner - # settles for e-1 matrix multiplications. Computing the minimal polynomial - # outright instead costs more than this whole method. - if IsGF2MatrixRep(mat) or Is8BitMatrixRep(mat) then - ratio := 6; + # 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. + # + # 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 - ratio := 2; + limit := 6; fi; - op := fail; - if d >= ratio then - op := POW_MAT_INT_PROBE(f, mat, QuoInt(d, ratio)); - fi; - t := fail; - if op <> fail and 1 < Length(op) then - e := Length(op) - 1; - # 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 computed first; and it is needed by the fallback anyway, so nothing - # is lost if the verification fails. Use mat directly only while 2*e - # stays below the cost of the base change, measured in dense - # multiplications -- which over GF(2) is much further. - if IsGF2MatrixRep(mat) then - limit := 60; - else - limit := 6; - fi; - if limit < e then - t := POW_MAT_INT_TRAFO(mat); - # conjugate matrices have the same annihilating polynomials - mm := t[3]; - else - mm := mat; - fi; - if POW_MAT_INT_ANNIHILATES(op, mm) then - pol := UnivariatePolynomialByCoefficients(ElementsFamily(FamilyObj(f)), - op, 1); - ind := IndeterminateOfUnivariateRationalFunction(pol); - pol := PowerMod(ind, n, pol); - mm := POW_MAT_INT_VALUE(pol, mm); - if t = fail then - return mm; - fi; - return t[2] * mm * t[1]; - 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; - if t = fail then - t := POW_MAT_INT_TRAFO(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 := POW_MAT_INT_VALUE(pol, mm);