Improve heuristic for choosing between matrix powering algorithms

[fingolfin]

(Replaces the original "Remove POW_MAT_INT" proposal, following the discussion below.)

Powering a matrix by a positive integer is handled by POW_MAT_INT, which chooses between
plain repeated squaring (POW_OBJ_INT, about Log2(n) matrix multiplications) and a method
based on computing x^n modulo the characteristic polynomial (considerable fixed overhead,
plus about Log2(n) much cheaper polynomial multiplications). The old rule chose the latter
whenever Log2(n) >= 3d/4, which is far too eager for small matrices (up to 30x slower than
repeated squaring, see the timings in the discussion) and at the same time too conservative
for large compressed matrices.

This PR replaces it by representation dependent break even points determined experimentally:
both algorithms were measured over a grid of dimensions d in 2..200, fields
GF(2), GF(5), GF(251), GF(257), GF(3^10) plus integer matrices, and exponent bit lengths
k = Log2(n) in 8..1600. Since for a fixed matrix the costs of both algorithms grow linearly
in k, the break even point of each cell can be computed reliably from linear fits. Findings:

• compressed 8-bit matrices (Is8BitMatrixRep): break even falls from k ~ 330 at d = 5
  to k ~ 10 at d = 200; modelled as k >= 8 + 768/d
• compressed GF(2) matrices (IsGF2MatrixRep): multiplication is so fast that the break even
  point stays roughly flat at k ~ 47..90 for all measured d; modelled as k >= 64
• everything else (plain lists, matrix objects, characteristic 0): polynomial arithmetic is
  nearly as expensive as small matrix products; the characteristic polynomial method never
  wins below d ~ 20 and the break even point is k ~ 60 for d >= 50; modelled as
  k >= Maximum(64, 8192/d)

On all measured cells the new rule is within a factor of ~1.9 of the optimal choice
(and optimal on most), compared to a factor of up to 30 for the old rule. In particular:

• the cases from the original PR description now run at POW_OBJ_INT speed, e.g.
  powering in GL(2,5) by 3*10^7 (10000 times) takes 52ms instead of 507ms
• @frankluebeck's use case keeps its speedup: powering an element of GL(200,251) by
  251^200 still uses the characteristic polynomial method and remains ~40x faster than
  repeated squaring
• the near-break-even case GL(20, 3^10) with e = 10^7 picks repeated squaring
  (~3% slower than optimal in my measurements)

The decision is based on the actual representation of the matrix, which is cheaper to test
than DefaultFieldOfMatrix (that call alone previously dominated the runtime when powering
small matrices by small exponents) and matches what actually determines the arithmetic cost.
DefaultFieldOfMatrix is only computed once the exponent is known to be large enough, to
verify that the entries lie in a field as required by the characteristic polynomial method;
matrices over other rings fall back to repeated squaring as before.

Not addressed here, possible follow-ups: using the minimal instead of the characteristic
polynomial (POW_BY_MINPOL as sketched by @frankluebeck above, with PowerModCoeffs), and a
kernel implementation for compressed matrices.

AI disclosure: this PR was prepared with assistance from the AI tool Claude Code
(benchmarking, analysis, drafting of the heuristic and of this text); the commit carries a
corresponding Co-authored-by: line.

[frankluebeck]

I do not approve this pull request. This code should not be removed. But the heuristics for choosing the best method in each case can certainly be improved. Certainly in the tiny cases you mention above.

POW_MAT_INT was written to speed up typical applications in group recognition (recog package).

In that context typical calls are for matrices in GL(n,q) and exponents in the order of q^n, e.g.,

gap> x := PseudoRandom(GL(200,251));;
gap> e := 251^200;;
gap> for i in [1..10] do x:=POW_OBJ_INT(x,e); od; time;
84571
gap> for i in [1..10] do x:=POW_MAT_INT(x,e); od; time;
1142

For compressed matrices over small finite fields the matrix multiplication is quite fast, and the overhead in POW_MAT_INT will only beat POW_OBJ_INT in somehow large cases.

For non-compressed matrices the break even point is much lower, e.g.,

gap> x := PseudoRandom(GL(20,3^10));;
gap> e := 10000000;;
gap> for i in [1..100] do x:=POW_OBJ_INT(x,e); od; time;
298
gap> for i in [1..100] do x:=POW_MAT_INT(x,e); od; time;
288
gap> for i in [1..100] do x:=POW_BY_MINPOL(x,e); od; time;   # see below
279

In an attempt to improve the heuristics it may also be useful to reenable the method which is currently commented out (many matrices have efficient methods for characteristic/minimal polynomial and for large exponents one can reduce to n multiplications):

POW_BY_MINPOL := function(mat, n)
  local pol, indet;
  # generic method for small n, break even point probably a bit lower,
  # needs rethinking and some experiments.
  if n < 2^Length(mat) then
    return POW_OBJ_INT(mat, n);
  fi;
  pol := MinimalPolynomial(mat); 
  # could be further improved using PowerModCoeffs:
  indet := IndeterminateOfUnivariateRationalFunction(pol);
  # cost of this needs to be investigated
  pol := PowerMod(indet, n, pol);
  # now we are sure that we need at most Length(mat) matrix multiplications
  return Value(pol, mat);
end;

As far as I remember, after introducing POW_MAT_INT there was the idea to implement a kernel version of this for compressed matrices. Therefore we did not work hard to find a good heuristic. But so far nobody implemented this ...

[fingolfin]

@frankluebeck thanks for the additional data points. That might help coming up with a better heuristic for when to use which algorithm,

[codecov]

Codecov Report

❌ Patch coverage is 83.33333% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 79.04%. Comparing base (d2134de) to head (eae5041).

Files with missing lines	Patch %	Lines
lib/matrix.gi	83.33%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #6293      +/-   ##
==========================================
- Coverage   79.08%   79.04%   -0.04%     
==========================================
  Files         685      685              
  Lines      293699   293705       +6     
  Branches     8664     8641      -23     
==========================================
- Hits       232264   232158     -106     
- Misses      59635    59745     +110     
- Partials     1800     1802       +2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
• 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

[fingolfin]

I've reworked this PR to keep both implementations but hopefully use a much better heuristic for the cross over points.

[fingolfin]

(And when I write "I" then I really mean me together with Claude Code)

[frankluebeck]

Thanks. This improment works as described.