Skip to content
Draft
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
128 changes: 0 additions & 128 deletions lib/matrix.gi
Original file line number Diff line number Diff line change
Expand Up @@ -4587,134 +4587,6 @@ local M, n, p, vars, slackVars, i, id, bestMove,
return [point,value];
end);


# can do better for matrices and large exponents, preliminary improvement,
# will be further improved (FL)
## InstallMethod( \^,
## "for matrices, use char. poly. for large exponents",
## [ IsMatrix, IsPosInt ],
## 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 := CharacteristicPolynomial(mat);
## 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);
##
# 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, addb, trafo, value, t, ti, mm, pol, ind;
d := NrRows(mat);
# finding a better break even point probably also depends on q
if n < 2^QuoInt(3*d,4) or not IsField(DefaultFieldOfMatrix(mat)) 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 returns a base change matrix such that t*m*t^-1 is block triangular
# with companion matrices along the diagonal
# (could/should? be improved to return t^-1, t*m*t^-1 and the
# characteristic polynomial of m at the same time)
trafo := function(m)
local id, b, t, r, a;
id := m^0;
b := rec(vectors := [], pivots := [], heads := []);
t := [];
# maybe better start with a random vector?
for a in id do
r := addb(b,a);
if r = true then
repeat
Add(t, a);
a := a*m;
r := addb(b,a);
until r <> true;
fi;
od;
t := Matrix(t, m);
return t;
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);
ti := t^-1;
mm := t * mat * ti;
pol := CharacteristicPolynomial(mm);
ind := IndeterminateOfUnivariateRationalFunction(pol);
pol := PowerMod(ind, n, pol);
mm := value(pol, mm);
return ti * mm * t;
end);

InstallMethod( \^,
"for matrices, use char. poly. for large exponents",
[ IsMatrixOrMatrixObj, IsPosInt ], POW_MAT_INT );

InstallGlobalFunction(RationalCanonicalFormTransform,function(mat)
local cr,R,x,com,nf,matt,p,i,j,di,d,v;
matt:=TransposedMat(mat);
Expand Down
Loading