1############################################################################# 2## 3## This file is part of GAP, a system for computational discrete algebra. 4## This file's authors include Michael Smith, Alexander Hulpke. 5## 6## Copyright of GAP belongs to its developers, whose names are too numerous 7## to list here. Please refer to the COPYRIGHT file for details. 8## 9## SPDX-License-Identifier: GPL-2.0-or-later 10## 11## This file contains meataxe type routines to compute module homomorphisms 12## for modules that are not necessarily irreducible. They are mainly a 13## conversion of the module routines in the GAP3 package `autag' by the 14## first author. 15## 16 17InstallGlobalFunction(TestModulesFitTogether,function(m1,m2); 18 if m1.field<>m2.field then 19 Error("different fields"); 20 fi; 21 if Length(m1.generators)<>Length(m2.generators) then 22 Error("generators are different lengths"); 23 fi; 24end); 25 26# basis for homomorphism space, efficient only for small dimensions 27BindGlobal("SmalldimHomomorphismsModules",function(m1,m2) 28local f, d1, d2, e, z, g1, g2, r, b, n, a, gp, i, j, k; 29 f:=m1.field; 30 d1:=m1.dimension; 31 d2:=m2.dimension; 32 e:=[]; 33 z:=ListWithIdenticalEntries(d1*d2,Zero(f)); 34 z:=ImmutableVector(f,z); 35 for gp in [1..Length(m1.generators)] do 36 g1:=m1.generators[gp]; 37 g2:=m2.generators[gp]; 38 for i in [1..d1] do 39 for j in [1..d2] do 40 # calculate equation for i-th row, j-th column 41 r:=ShallowCopy(z); 42 # the entry in g*hom is the i-th row of g with the variables in the 43 # j-th column 44 for k in [1..d1] do 45 b:=(k-1)*d2+j; 46 r[b]:=r[b]+g1[i][k]; 47 od; 48 # the entry in hom*g is the variables in the i-th row of hom with the 49 # j-th column of g 50 for k in [1..d2] do 51 b:=(i-1)*d2+k; 52 r[b]:=r[b]-g2[k][j]; 53 od; 54 Add(e,r); 55 od; 56 od; 57 od; 58 n:=NullspaceMat(TransposedMat(e)); 59 b:=[]; 60 for i in n do 61 # convert back to d1 x d2 matrix 62 a:=[]; 63 for j in [1..d1] do 64 Add(a,i{[(j-1)*d2+1..(j-1)*d2+d2]}); 65 od; 66 a:=ImmutableMatrix(f,a); 67 Add(b,a); 68 od; 69 return b; 70end); 71 72 73# the following code is essentially due to Michael Smith 74 75# These routines are designed to accumulate a system of linear equations 76# 77# M_1 X = V_1, M_2 X = V_2 ... M_t X = V_t 78# 79# Where each M_i is an m_i*n matrix, X is the unknown length n vector, and 80# each V is an length m_i vector. The equations can be added as each batch 81# is calculated. Here is some pseudo-code to demonstrate: 82# 83# eqns := newEqns (n, field); 84# i := 1; 85# repeat 86# <calculate M_i and V_i> 87# addEqns(M_i, V_i) 88# increment i; 89# until i > t or eqns.failed; 90# if not eqns.failed then 91# S := solveEqns(eqns); 92# fi; 93# 94# As demonstrated by the example, an early notification of failure is 95# available by checking ".failed". All new equations are sifted with respect 96# to the current set, and only added if they are independent of the current 97# set. If a new equation reduces to the zero row and a nonzero vector 98# entry, then there is no solution and this is immediately returned by 99# setting eqns.failed to true. The function solveEqns has an already 100# triangulised system of equations, so it simply reduces above the pivots 101# and returns the solution vector. 102 103 104BindGlobal("SMTX_AddEqns",function ( eqns, newmat, newvec) 105local n, zero, weights, mat, vec, NextPositionProperty, ReduceRow, t, 106 newweight, newrow, newrhs, i, l, k; 107 108# Add a bunch of equations to the system of equations in <eqns>. Each 109# row of <newmat> is the left-hand side of a new equation, and the 110# corresponding row of <newvec> the right-hand side. Each equation in 111# filtered against the current echelonised system stored in <eqns> and 112# then added if it is independent of the system. As soon as a 113# left-hand side reduces to 0 with a non-zero right-hand side, the flag 114# <eqns.failed> is set. 115 116 Info(InfoMtxHom,6,"addEqns: entering" ); 117 118 n := eqns.dim; 119 zero := Zero(eqns.field); 120 weights := eqns.weights; 121 mat := eqns.mat; 122 vec := eqns.vec; 123 124 NextPositionProperty := function (list, func, start ) 125 local i; 126 for i in [ start .. Length( list ) ] do 127 if func( list[ i ] ) then 128 return i; 129 fi; 130 od; 131 return fail; 132 end; 133 134 # reduce the (lhs,rhs) against the semi-echelonised current matrix, 135 # and return either: (1) the reduced rhs if the lhs reduces to zero, 136 # or (2) a list containing the new echelon weight, the new row and 137 # the new rhs for the system, and the row number that this 138 # equation should placed. 139 ReduceRow := function (lhs, rhs) 140 local lead, i, z; 141 lead := PositionProperty(lhs, i->i<>zero); 142 if lead = fail then 143 return rhs; 144 fi; 145 for i in [1..Length(weights)] do 146 if weights[i] = lead then 147 z := lhs[lead]; 148 lhs := lhs - z * mat[i]; rhs := rhs - z * vec[i]; 149 lead := NextPositionProperty(lhs, i->i<>zero, lead); 150 if lead = fail then 151 return rhs; 152 fi; 153 elif weights[i] > lead then 154 return [lead, lhs, rhs, i]; 155 fi; 156 od; 157 return [lead, lhs, rhs, Length(weights)+1]; 158 end; 159 160 for k in [1..Length(newmat)] do 161 t := ReduceRow(newmat[k], newvec[k]); 162 163 if IsList(t) then 164 # new equation 165 newweight := t[1]; 166 newrow := t[2]; 167 newrhs := t[3]; 168 i := t[4]; # position for new row 169 170 # normalise so that leading entry is 1 171 newrhs := newrhs / newrow[newweight]; 172 newrow := newrow / newrow[newweight]; # NB: in this order 173 174 if i = Length(mat)+1 then 175 # add new equation to end of list 176 Add(mat, newrow); 177 Add(vec, newrhs); 178 Add(weights, newweight); 179 else 180 l := Length(mat); 181 # move down other rows to make space for this new one... 182 mat{[i+1..l+1]} := mat{[i..l]}; 183 vec{[i+1..l+1]} := vec{[i..l]}; 184 # and then slot it in 185 mat[i] := newrow; 186 vec[i] := newrhs; 187 weights{[i+1..l+1]} := weights{[i..l]}; 188 weights[i] := newweight; 189 fi; 190 191 else 192 # no new equation, check whether inconsistent due to 193 # nonzero rhs reduction 194 195 if not IsZero(t) then 196 Info(InfoMtxHom,6,"addEqns: FAIL!" ); 197 eqns.failed := true; 198 return eqns; # return immediately 199 fi; 200 fi; 201 od; 202end); 203 204BindGlobal("SMTX_NewEqns",function (arg) 205local X, n, F, V, eqns; 206 207 if Length(arg) <2 then 208 Error("NewEqns(dim, field) or NewEqns(X, V)"); 209 fi; 210 211 if IsInt(arg[1]) then 212 X := false; 213 n := arg[1]; 214 F := arg[2]; 215 else 216 X := arg[1]; 217 V := arg[2]; 218 n := Length(X[1]); 219 F := Field(X[1][1]); # Note: prime field only 220 fi; 221 222 eqns := rec(); 223 eqns.dim := n; # number of variables 224 eqns.field := F; # field over which the equation hold 225 eqns.mat := []; # left-hand sides of system 226 eqns.weights := []; # echelon weights for lhs matrix 227 eqns.vec := []; # right-hand sides of system 228 eqns.failed := false; # flag to indicate inconsistent system 229 eqns.index := []; # index for row ordering 230 231 if IsMatrix(X) then 232 SMTX_AddEqns(eqns, X, V); 233 fi; 234 235 return eqns; 236end); 237 238BindGlobal("SMTX_KillAbovePivotsEqns",function (eqns) 239# Eliminate entries above pivots. Note that the pivot entries are 240# all 1 courtesy of SMTX_AddEqns. 241 242local m, n, zero, i, c, j, factor; 243 244 Info(InfoMtxHom,6,"killAbovePivotsEqns: entering" ); 245 m := Length(eqns.mat); 246 n := eqns.dim; 247 if m > 0 then 248 zero := Zero(eqns.field); 249 for i in [1..m] do 250 c := eqns.weights[i]; 251 for j in [1..i-1] do 252 if eqns.mat[j][c] <> zero then 253 Info(InfoMtxHom,6,"solveEqns: kill mat[",j,",",c,"]"); 254 factor := eqns.mat[j][c]; 255 eqns.mat[j] := eqns.mat[j] - factor*eqns.mat[i]; 256 eqns.vec[j] := eqns.vec[j] - factor*eqns.vec[i]; 257 fi; 258 od; 259 od; 260 fi; 261 Info(InfoMtxHom,6,"killAbovePivotsEqns: leaving" ); 262end); 263 264BindGlobal("SMTX_NullspaceEqns",function(e) 265 266# Take the matrix stored in equation record <e> and compute a basis 267# for its nullspace, ie x such that mat * x = 0. Note that the 268# vector is on the other side of the matrix from GAP's NullspaceMat. 269# This means we get to skip the Transposing that occurs at the top 270# of that function (a bonus!). 271# 272# This function is a modified version NullspaceMat in matrix.g 273 274local mat, m, n, zero, one, empty, i, k, nullspace, row,mi; 275 276 SMTX_KillAbovePivotsEqns(e); 277 mat := e.mat; 278 279 m := Length(mat); 280 n := e.dim; 281 zero := Zero(e.field); 282 one := One(e.field); 283 284 # insert empty rows to bring the leading term of each row on the diagonal 285 empty := ListWithIdenticalEntries(n,zero); 286 empty := ImmutableVector(e.field,empty); 287 i := 1; 288 while i <= Length(mat) do 289 if i < n and mat[i][i] = zero then 290 mi:=Minimum(Length(mat),n-1); 291 for k in [mi,mi-1..i] do 292 mat[k+1] := mat[k]; 293 od; 294 mat[i] := empty; 295 fi; 296 i := i+1; 297 od; 298 for i in [ Length(mat)+1 .. n ] do 299 mat[i] := empty; 300 od; 301 302 # The following comment from NullspaceMat: 303 # 'mat' now looks like [ [1,2,0,2], [0,0,0,0], [0,0,1,3], [0,0,0,0] ], 304 # and the solutions can be read in those columns with a 0 on the diagonal 305 # by replacing this 0 by a -1, in this example [2,-1,0,0], [2,0,3,-1]. 306 nullspace := []; 307 for k in [1..n] do 308 if mat[k][k] = zero then 309 row := []; 310 for i in [1..k-1] do row[i] := -mat[i][k]; od; 311 row[k] := one; 312 for i in [k+1..n] do row[i] := zero; od; 313 Add( nullspace, row ); 314 fi; 315 od; 316 317 return nullspace; 318end); 319 320 321BindGlobal("EchResidueCoeffs",function (base, ech, v,mode) 322local n, coeffs, x, zero, z, i; 323# 324# Take a semi-ech basis <base>, with ech weights <ech>, and a vector 325# <v> in the subspace spanned by <base>. Returns: 326# if mode>2: 327# a record containing the 328# residue after removing projection of <v> onto subspace spanned by 329# <base>, as well as the coefficients of the linear combination of 330# <base> elements used to obtain the projection. Also return the 331# projection. 332# if mode =1 returns only the coefficients 333# if mode=2 returns only the residue 334 335# Note that the pivots of <base> must be set to 1. 336 337 n:=Length(base); 338 339 if n = 0 then 340 coeffs:=[]; 341 x:=v; 342 else 343 x:=v; 344 zero:=x[1]*0; 345 coeffs:=ListWithIdenticalEntries(n, zero); 346 for i in [1..n] do 347 z:=x[ech[i]]; 348 if z <> zero then 349 x:=x - z * base[i]; 350 coeffs[i]:=z; 351 fi; 352 od; 353 fi; 354 355 if mode=1 then 356 return coeffs; 357 elif mode=2 then 358 return x; 359 else 360 return rec(coeffs:=coeffs, 361 residue:=x, 362 projection:=v - x 363 ); 364 fi; 365end); 366 367BindGlobal("SpinSpaceVector",function (V, U, ech, v,zero) 368local gens, pos, settled, oldlen, i, j; 369 370# Take <U> a semi-ech basis for a submodule of <V>, with ech-weights 371# <ech>, and a vector <v> in <V>. Return a semi-ech basis for the 372# submodule generated by <U> and <v>. 373 374 U:=ShallowCopy(U); 375 ech:=ShallowCopy(ech); 376 gens:=V.generators; 377 378 v:=EchResidueCoeffs(U, ech, v,2); 379 pos:=PositionProperty(v, i->i<>zero); 380 if pos = fail then 381 return U; 382 fi; 383 Add(U, v/v[pos]); Add(ech, pos); 384 385 settled:=Maximum(Length(U),1); # <U> is a submodule 386 repeat 387 oldlen:=Length(U); 388 for i in [settled+1..Length(U)] do 389 for j in [1..Length(gens)] do 390 v:=EchResidueCoeffs(U, ech, (U[i] * gens[j]),2); 391 pos:=PositionProperty(v, i->i<>zero); 392 if pos <> fail then 393 Add(U, v/v[pos]); Add(ech, pos); 394 fi; 395 od; 396 od; 397 settled:=oldlen; 398 until oldlen = Length(U); 399 return U; 400end); 401 402BindGlobal("SpinHomFindVector",function (r) 403local V, nv, W, nw, U, echu, F, matsV, matsW, k, g1, g2, max_stack_len, _t, 404 newstack, v0, extradim, N, count, look_lim, done, grpalg, i, M, pos, A, 405 j,zero; 406 407# <r> contains information about modules <V> and <W>, and a submodule 408# <U> of <V> with semi-ech information <echu>. The routine selects 409# an element of <V> lying outside of <U> that will be used to spin 410# up to a new submodule U'. 411# 412# It returns a list [<v0>, <M>] where <v0> is the element of <V> 413# and <M> is a basis for a submodule of <W> which <v0> must map into 414# under any hom. 415 416 V:=r.V; nv:=V.dimension; 417 W:=r.W; nw:=W.dimension; 418 U:=r.U; echu:=r.echu; 419 F:=V.field; 420 zero:=Zero(F); 421 422 if not IsBound(r.mats) then 423 matsV:=V.generators; 424 matsW:=W.generators; 425 k:=Length(matsV); 426 r.mats:=List([1..k], i -> [matsV[i], matsW[i]]); 427 428 # do preprocessing to make random matrices list in parallel 429 430 for i in [1..10] do 431 g1:=Random(1, k); 432 g2:=g1; 433 while g2 = g1 and Length(r.mats)>1 do 434 g2:=Random(1, k); 435 od; 436 Add(r.mats,[r.mats[g1][1]*r.mats[g2][1], 437 r.mats[g1][2]*r.mats[g2][2]]); 438 k:=k + 1; 439 od; 440 441 r.zero:=[ ImmutableMatrix(F,NullMat(nv,nv,F)), 442 ImmutableMatrix(F,NullMat(nw,nw,F)) ]; 443 444 # we build a stack of good grpalg elements to use for choosing 445 # elements <v0> --- an element <A> in <stack> is of the form: 446 # A[1] = v0 447 # A[2] = grpalg element whose nullspace contains v0 448 # A[3] = Dim(<U,v0>^G)-Dim(U) i.e. increase in dim by adding 449 # <v0> to <U> 450 r.stack:=[]; 451 452 else 453 k:=Length(r.mats); 454 fi; 455 456 max_stack_len:=10; 457 458 # adjust the elements of the stack to account for the larger 459 # submodule <U> we now have 460 _t:=Runtime(); 461 newstack:=[]; 462 for A in r.stack do 463 v0:=A[1]; 464 extradim:=Length(SpinSpaceVector(V, U, echu, v0,zero)) 465 - Length(U); 466 if extradim > 0 then 467 Add(newstack, [v0, A[2], extradim]); 468 fi; 469 od; 470 r.stack:=newstack; 471 Info(InfoMtxHom,2,"stack reduced to length ", Length(r.stack), " (", 472 Runtime()-_t, ")"); 473 474 # <N> contains the nullspace in <V> of a group algebra element --- 475 # initialise it to the empty list for the following repeat loop 476 N:=[]; 477 478 count:=0; 479 look_lim:=5; # give up after this many random grpalg elements 480 481 _t:=Runtime(); 482 483 if Length(r.stack) > 0 then 484 # if we have something left, don't bother generating any new 485 # grpalg elements (?) 486 count:=look_lim + 1; 487 fi; 488 489 done:=false; 490 while count < look_lim and Length(r.stack) < max_stack_len and not done do 491 492 # we look for a while and take the best element found 493 494 # We are looking for an element <v0> of a nullspace that lies 495 # outside of <U> 496 497 repeat 498 # Take a work record <r> containing the information about the two 499 # modules <V> and <W>, and return a random group algebra element 500 # record containing its action on each of the modules. 501 502 # first take two elements of the list and multiply them 503 # together 504 g1:=Random(1, k); 505 repeat 506 g2:=Random(1, k); 507 until g2 <> g1 or Length(r.mats)=1; 508 Add(r.mats,[r.mats[g1][1]*r.mats[g2][1], 509 r.mats[g1][2]*r.mats[g2][2]]); 510 k:=k + 1; 511 512 # Now take a random linear sum of the existing generators as new 513 # generator. Record the sum in coefflist 514 515 grpalg:=ShallowCopy(r.zero); 516 517 for g1 in [1..k] do 518 g2:=Random(F); 519 if not IsZero(g2) then 520 grpalg[1]:=grpalg[1] + g2*r.mats[g1][1]; 521 grpalg[2]:=grpalg[2] + g2*r.mats[g1][2]; 522 fi; 523 od; 524 N:=TriangulizedNullspaceMat(grpalg[1]); 525 count:=count + 1; 526 until Length(N) > 0 or count >= look_lim; 527 528 if Length(N) > 0 then 529 530 # now find best element of <N> for adding to <stack> 531 extradim:=List(N, y -> 532 Length(SpinSpaceVector(V, U, echu, y,zero)) 533 - Length(U)); 534 i:=1; 535 for j in [2..Length(extradim)] do 536 if extradim[j] > extradim[i] then 537 i:=j; 538 fi; 539 od; 540 if extradim[i] > 0 then 541 # exit early if we have found an element that gets use all 542 # of <V> after spinning 543 done:=extradim[i] = nv - Length(U); 544 if done then 545 r.stack:=[[N[i], grpalg, extradim[i]]]; 546 else 547 Add(r.stack, [N[i], grpalg, extradim[i]]); 548 fi; 549 fi; 550 fi; 551 552 od; 553 Info(InfoMtxHom,2,"stack loop done, stack now length ", Length(r.stack), " (", 554 Runtime()-_t, ")"); 555 556 if Length(r.stack) > 0 then 557 # 558 # find best element in r.stack and use it 559 i:=1; 560 for j in [2..Length(r.stack)] do 561 if r.stack[j][3] > r.stack[i][3] then 562 i:=j; 563 fi; 564 od; 565 v0:=r.stack[i][1]; 566 M:=TriangulizedNullspaceMat(r.stack[i][2][2]); 567 568 else 569 570 # we haven't found a good grpalg element, so just choose 571 # something outside of <U> and use it 572 573 Info(InfoMtxHom,1,"too many random grpalg elements..."); 574 M:=IdentityMat(nw,F); 575 pos:=Difference([1..nv], echu)[1]; 576 v0:=ListWithIdenticalEntries(nv,zero); 577 v0[pos]:=One(F); 578 v0:=ImmutableVector(F,v0); 579 fi; 580 return [v0, M]; 581 582end); 583 584# compute a semi-echelonised basis for a matrix algebra 585# If a linearly dependent set of elements is supplied, this 586# routine will trim it down to a basis. 587BindGlobal("SMTX_EcheloniseMats",function (gens, F) 588local n, m, zero, ech, k, i, j, found, l; 589 590 if Length(gens) = 0 then 591 return [ [], [] ]; 592 fi; 593 # copy the list to avoid destroying the original list 594 gens:=List(gens,i->List(i,ShallowCopy)); 595 596 n:=Length(gens[1]); 597 m:=Length(gens[1][1]); 598 zero:=Zero(F); 599 600 ech:=[]; 601 k:=1; 602 603 while k <= Length(gens) do 604 i:=1; j:=1; 605 found:=false; 606 while not found and i <= n do 607 if (gens[k][i][j] <> zero) then 608 found:=true; 609 else 610 j:=j + 1; 611 if (j > m) then 612 j:=1; i:=i + 1; 613 fi; 614 fi; 615 od; 616 617 if found then 618 619 # Now basis element k will have echelonisation index [i,j] 620 Add(ech, [i,j]); 621 622 # First normalise the [i,j] position to 1 623 gens[k]:=gens[k] / gens[k][i][j]; 624 625 # Now zero position [i,j] in all further generators 626 for l in [k+1..Length(gens)] do 627 if (gens[l][i][j] <> zero) then 628 gens[l]:=gens[l] - gens[k] * gens[l][i][j]; 629 fi; 630 od; 631 k:=k + 1; 632 else 633 # no non-zero element found, delete from list 634 gens{[k..Length(gens)-1]}:=gens{[k+1..Length(gens)]}; 635 Unbind(gens[Length(gens)]); 636 # WAS: gens:=gens{ Cat([1..k-1], [k+1..Length(gens)])}; 637 fi; 638 od; 639 return [List(gens,i->ImmutableMatrix(F,i)), ech]; 640end); 641 642# The SpinHom routine in this file was written during August 1996. The 643# basic idea comes from a discussion I had with Charles Leedham-Green early 644# in 1995. He gave me a rough sketch of the algorithm that he and John 645# Cannon developed for Magma. Some details were missing, and this is my 646# attempt at filling in some of them. 647# 648# Many improvements were made on my earlier version, in large part due to a 649# discussion I had with Alice Niemeyer in early 1996. She relayed to me 650# some comments of Klaus Lux on my earlier version. This is a combination 651# of the suggestions of Klaus and Alice and my own ideas. 652# 653# Note: This provides an enormous speed-up on the default GAP routine, 654# and on my own naive intertwining routine, especially when the module is 655# large enough and/or it is irreducible. However, this routine is nowhere 656# near as good as the Magma algorithm, and I do not know how to improve it. 657# 658# The code is heavily commented, and I appreciate suggestions on how to 659# improve it (particularly bits of code). 660BindGlobal("SpinHom",function (V, W) 661local nv, nw, F, zero, zeroW, gV, gW, k, U, echu, r, homs, s, work, ans, v0, 662 M, x, pos, z, echm, t, v, echv, a, u, e, start, oldlen, ag, m, uu, ret, 663 c, s1, X, mat, uuc, uic, newhoms, hom, Uhom, imv0, imv0c, image, i, j, l; 664 665# Compute Hom(V,W) for G-modules <V> and <W>. The algorithm starts with 666# the trivial submodule <U> of <V> for which Hom(U,V) is trivial. It 667# then computes Hom(U',W) for U' a submodule generated by <U> and a 668# single element <v0> in <V>. This U' becomes the next <U> as the process 669# is iterated, ending when <U'> = <V>. The element <v0> is chosen in a 670# nullspace of a group algebra element in order to restrict it possible 671# images in <W>. 672 673 nv:=V.dimension; 674 nw:=W.dimension; 675 676 F:=V.field; 677 if F<>W.field then 678 Error("different fields"); 679 fi; 680 zero:=Zero(F); 681 682 zeroW:=ListWithIdenticalEntries(nw,zero); 683 zeroW:=ImmutableVector(F,zeroW); 684 685 # group generating sets acting on each module 686 gV:=V.generators; 687 gW:=W.generators; 688 689 # <k> is the number of generators of the acting group 690 k:=Length(gV); 691 if k<>Length(gW) then 692 Error("generator lengths"); 693 fi; 694 695 # <U> is the semi-ech basis for the currently known submodule, of 696 # dimension <r> 697 U:=[]; 698 echu:=[]; 699 r:=0; 700 701 # <homs> contains a basis for Hom(U,W), of dimension <s> 702 homs:=[]; 703 s:=0; 704 705 # define a record which stores information about the modules <V>, <W> 706 # and <U> for passing into a routine that selects a new vector <v0> 707 # for spinning up to a larger submodule U'. 708 work:=rec(V:=V, W:=W, U:=U, echu:=echu); 709 710 repeat 711 712 # we loop until <U> is the whole of <V> 713 714 ans:=SpinHomFindVector(work); 715 v0:=ans[1]; 716 M:=ans[2]; 717 718 # find residue of <v0> modulo current submodule <U> 719 x:=EchResidueCoeffs(U, echu, v0,2); 720 721 # normalise <x> (ie get a 1 in leading position) 722 pos:=PositionProperty(x, i->i<>zero); 723 z:=x[pos]; 724 x:=x / z; 725 v0:=v0 / z; 726 727 # we know that <v0> has to map into the subspace <M> of <W>. 728 echm:=List(M, y -> PositionProperty(y, i->i<>zero)); 729 t:=Length(M); 730 731 # now we start building extension of semi-echelonised basis for 732 # the submodule U' generated by <U> and <v0> 733 # 734 # new elements of semi-ech basis will be stored in <v>, with 735 # echelon weights stored in <echv> 736 737 v:=[ x ]; 738 echv:=[ pos ]; 739 740 # we need to keep track of how each new element of the semi-ech 741 # basis was obtained from <v0> --- new basis element <v[i]> will 742 # satisfy: 743 # 744 # v[i] = v0*a[i] + u[i] 745 # 746 # where <a[i]> is an element of the group algebra FG, and <u[i]> is 747 # the element of <U> that was subtracted during semi-ech reduction 748 749 a:=[ M ]; 750 u:=[ x - v0 ]; 751 752 # we will accumulate the homogeneous linear system in <e> 753 # 754 # the first <s> variables are the coefficients of basis elements of 755 # Hom(U,W), which describes how a hom of U' acts on submodule <U> 756 # 757 # the other <t> variables are the coefficients of basis elements of 758 # <M>, which describes the image of <v0> under a hom 759 # 760 e:=SMTX_NewEqns(s + t, F); 761 762 # we will close the submodule by spinning <v0> --- the variable 763 # <start> will trim off the elements of <v> that we have already 764 # used 765 start:=1; 766 767 repeat 768 769 # take an element <v[i]> of <v> and a group generator <g[j]> 770 # and check whether <v[i]^g[j]> is a new basis element. 771 # 772 # if it is, add it to the basis, with its definition. 773 # 774 # if it isn't, we get an equation which an element of Hom(U',W) 775 # must satisfy 776 777 oldlen:=Length(v); 778 779 for i in [start..oldlen] do ### loop on vectors in <v> 780 for j in [1..k] do ### loop on generators of G 781 782 if Length(a[i])=0 then 783 #T: special treatment 0-dimensional 784 ag:=[]; 785 else 786 ag:=a[i] * gW[j]; 787 fi; 788 789 # create new element <x>, with its definition as the 790 # difference between <v0^m> and <uu> in <U>. 791 x:=v[i] * gV[j]; 792 m:=ag; 793 uu:=u[i] * gV[j]; 794 795 ret:=EchResidueCoeffs(U, echu, x,3); 796 x:=ret.residue; 797 uu:=uu - ret.projection; 798 799 # reduce modulo the new semi-ech basis elements in <v>, 800 # storing the coefficients in <c> 801 # 802 c:=ListWithIdenticalEntries(Length(v),zero); 803 for l in [1..Length(v)] do 804 z:=x[echv[l]]; 805 if z <> zero then 806 x:=x - z * v[l]; 807 if Length(m) > 0 then 808 m:=m - z * a[l]; 809 fi; 810 c[l]:=c[l] + z; 811 uu:=uu - z * u[l]; 812 fi; 813 od; 814 c:=ImmutableVector(F,c); 815 816 # Note: at this point, <x> has been reduced modulo the 817 # semi-ech basis <U> union <v>, and that 818 # 819 # x = v0 * a[i] + uu 820 821 pos:=PositionProperty(x, i->i<>zero); 822 if pos <> fail then 823 824 # new semi-ech basis element <x> 825 826 z:=x[pos]; 827 Add(v, x/z); 828 Add(echv, pos); 829 Add(a, m/z); 830 Add(u, uu/z); 831 832 else 833 834 # we get some equations ! 835 836 s1:=Sum([1..Length(v)], y -> c[y] * v[y]); 837 uu:=v[i] * gV[j] - s1; 838 839 X:=NullMat(t, nw, F); 840 for l in [1..Length(v)] do 841 if c[l] <> zero then 842 if Length(X) > 0 then 843 X:=X + c[l] * a[l]; 844 fi; 845 uu:=uu + c[l] * u[l]; 846 fi; 847 od; 848 849 if Length(X) > 0 then 850 X:=X - ag; 851 fi; 852 853 mat:=[]; 854 uuc:=EchResidueCoeffs(U, echu, uu,1); 855 uic:=EchResidueCoeffs(U, echu, u[i],1); 856 for l in [1..s] do 857 Add(mat, uuc * homs[l] - uic * homs[l] * gW[j]); 858 od; 859 Append(mat, X); 860 SMTX_AddEqns(e, TransposedMat(mat), zeroW); 861 fi; 862 od; 863 od; 864 865 start:=oldlen+1; 866 867 # exit when no new elements were added --- i.e. the subspace 868 # is closed under action of G and is therefore a submodule 869 870 until oldlen = Length(v); 871 872 # we have the system of equations, so find its solution space 873 874 ans:=SMTX_NullspaceEqns(e); 875 876 # Now build the homomorphisms 877 878 newhoms:=[]; 879 for i in [1..Length(ans)] do 880 881 # Each row of ans is of the form: 882 # 883 # [ b_1, b_2, ..., b_s, c_1, c_2, ..., c_t ] 884 # 885 # where the action of this hom on <U> is as \Sum{b_l homs[l]} 886 # and the hom sends <v0> to Sum{c_l M[l]} 887 888 hom:=[]; 889 if r > 0 then 890 Uhom:=NullMat(r, nw, F); 891 for l in [1..s] do 892 if ans[i][l] <> zero then 893 Uhom:=Uhom + ans[i][l] * homs[l]; 894 fi; 895 od; 896 for l in [1..r] do 897 Add(hom, Uhom[l]); 898 od; 899 fi; 900 901 imv0:=zeroW * zero; 902 for l in [1..t] do 903 if ans[i][s+l] <> zero then 904 imv0:=imv0 + ans[i][s+l] * M[l]; 905 fi; 906 od; 907 imv0c:=EchResidueCoeffs(M, echm, imv0,1); 908 for l in [1..Length(v)] do 909 image:=imv0c * a[l]; 910 if r > 0 then 911 image:=image + EchResidueCoeffs(U, echu, u[l],1) * Uhom; 912 fi; 913 Add(hom, image); 914 od; 915 hom:=ImmutableMatrix(F,hom); 916 Assert(1,hom<>0*hom); 917 Add(newhoms, hom); 918 od; 919 920 # now update <U> to be the now larger submodule 921 922 Append(U,v); 923 Append(echu, echv); 924 homs:=newhoms; 925 r:=Length(U); 926 s:=Length(homs); 927 928 Info(InfoMtxHom,1,"U is now dimension ", r, " and dim(Hom(U,W)) = ", s); 929 930 until r = nv; # i.e. <U> = <V> 931 932 if Length(homs)=0 then 933 return homs; 934 fi; 935 936 # We must change basis on <V> from <U> to the usual one before returning 937 938 U:=ImmutableMatrix(F,U); 939 return U^-1 * homs; 940 941end); 942 943 944# module isomorphism and decomposition routines 945# 946# These are functions for computing with modules, including: 947# 948# (1) computing a direct sum decomposition of a module into 949# indecomposable summands. 950# 951# (2) deciding module isomorphism using the decomposition. 952# 953# The algorithm for deciding indecomposability is based on the algorithm 954# described by G. Schneider in the Journal of Symbolic Computation, 955# Volume 9, Numbers 5 & 6, 1990 956 957 958# Take a Fitting element and use it to split M into a direct sum 959# of submodules. Return the submodules. 960# r is the rank of a (which migth be known before 961BindGlobal("FittingSplitModule",function (a,r,F) 962local n, ro; 963 964 # do we have a fitting matrix? 965 # a matrix is a fitting matrix if it is singular but not nilpotent. 966 # case 967 n:=Length(a); 968 if r=n or r=0 then 969 # not singular or zero. 970 return fail; 971 fi; 972 # now square repeatedly until the rank stays the same and >0 973 repeat 974 ro:=r; 975 a:=a^2; 976 r:=RankMat(a); 977 until ro=r or r=0; 978 if r=0 then 979 return fail; 980 fi; 981 # otherwise a is a power of a fitting matrix, the space will split in 982 # Kern(a) \oplus Image(a) 983 984 Info(InfoMeatAxe,2,"Decomposition ",r,":",n-r," found"); 985 986 return [ImmutableMatrix(F,BaseMat(a)),NullspaceMat(a)]; 987end); 988 989# Take a module and break it into two pieces if possible. 990# The function searches for a decomposition of the module M while 991# attempting to prove indecomposability at the same time. Of course, 992# only one of these will succeed. 993BindGlobal("ProperModuleDecomp",function (M) 994local proveIndecomposability, addnilpotent, n, F, zero, basis, enddim, 995 echelon, nildim, p, maxorder, maxa, nilbase, nilech, cnt, remain, 996 used, coeffs, a, rk, order, fit, pos, newa, lastdim, i; 997 998 # Check whether we have found the indecomposability proof. That is, 999 # see whether our regular element generates a subalgebra which 1000 # complements the current nilpotent ideal (the approximation to 1001 # radical) 1002 proveIndecomposability:=function () 1003 local maxaord; 1004 # NB: <maxa> is not local 1005 1006 if enddim - nildim = LogInt(maxorder + 1,p) then 1007 # Yes, found the residue field root and proved indecomposability! 1008 maxaord:=Order(maxa); 1009 while maxaord > maxorder do 1010 maxa:=maxa^p; 1011 maxaord:=maxaord / p; 1012 od; 1013 SMTX.SetEndAlgResidue(M, [maxa, maxaord]); 1014 Info(InfoMtxHom,3,"proved ",Length(nilbase)); 1015 SMTX.SetBasisEndomorphismsRadical(M, nilbase); 1016 return true; 1017 fi; 1018 return false; 1019 end; 1020 1021 # take a new nilpotent element and sift against current nilpotent 1022 # ideal basis. If it does not lie in the space spanned so far, 1023 # add it to nilbasis 1024 addnilpotent:=function (a) 1025 local i, r, c, k, done, l; 1026 # NB: <remain> and <nildim> and <cnt> are not local 1027 1028 for i in [1..nildim] do 1029 r:=echelon[nilech[i]][1]; c:=echelon[nilech[i]][2]; 1030 if a[r][c] <> zero then 1031 a:=a - a[r][c] * nilbase[i] / nilbase[i][r][c]; 1032 fi; 1033 od; 1034 1035 # find which echelon index to remove due to this new element 1036 k:=1; done:=false; 1037 while not done and k <= Length(remain) do 1038 l:=remain[k]; 1039 r:=echelon[l][1]; c:=echelon[l][2]; 1040 if a[r][c] <> zero then 1041 done:=true; 1042 else 1043 k:=k + 1; 1044 fi; 1045 od; 1046 1047 if k > Length(remain) then 1048 # in nilpotent ideal already, return 1049 return false; 1050 fi; 1051 1052 # We now know this nilpotent element is a new one 1053 Add(nilbase, a); 1054 1055 # the k-th basis element was used to make the new element a. So 1056 # remove it from future random element calculations 1057 # 1058 Add(nilech, remain[k]); 1059 remain:=Difference(remain, [remain[k]]); 1060 nildim:=nildim + 1; 1061 cnt:=1; 1062 return true; 1063 end; 1064 1065 if not M.IsOverFiniteField then 1066 return Error ("Argument of ProperModuleDecomp is not over a finite field."); 1067 fi; 1068 n:=M.dimension; 1069 F:=M.field; 1070 1071 zero:=Zero(F); 1072 Info(InfoMtxHom,2,"ProperModuleDecomp for module of dimension ", n); 1073 1074 if n = 1 then 1075 # A 1-dimensional module is always indecomposable 1076 Info(InfoMtxHom,3,"1dimensional"); 1077 SMTX.SetEndAlgResidue(M, [[[ PrimitiveElement(F) ]], Size(F) - 1]); 1078 SMTX.SetBasisEndomorphismsRadical(M, []); 1079 return fail; 1080 fi; 1081 1082 basis:=SMTX.BasisModuleEndomorphisms(M); 1083 if Length(basis) = 1 then 1084 # if endomorphism algebra has dimension 1 then indecomposable 1085 #SMTX.SetEndAlgResidueFlag(M, F.root * GModOps.EndAlgBasisFlag(M)[1], F.size - 1); 1086 SMTX.SetEndAlgResidue(M, [PrimitiveElement(F)*basis[1], Size(F) - 1]); 1087 Info(InfoMtxHom,3,"basislength 1"); 1088 SMTX.SetBasisEndomorphismsRadical(M, []); 1089 return fail; 1090 fi; 1091 1092 enddim:=Length(basis); # dim of endo algebra 1093 echelon:=SMTX_EcheloniseMats(basis,F)[2]; # echelon indices for endalg basis 1094 nildim:=0; # dim of current approx to radical 1095 p:=Size(F); 1096 maxorder:=1; # order of largest order regular elmt 1097 # found so far 1098 maxa:=IdentityMat(n,F); # the regular elmt with order maxorder 1099 nilbase:=[]; # basis for approx to radical 1100 nilech:=[]; 1101 1102 cnt:=1; 1103 1104 # We will "quotient" out the nilpotent subspace as we go. The elements 1105 # of remain tell us which (echelonised) basis elements of the 1106 # endomorphism algebra we will take use in our random linear 1107 # combination. 1108 # 1109 remain:=[1..enddim]; 1110 used:=[]; 1111 1112 repeat 1113 # we will loop until too many passes without an improvement in knowledge 1114 repeat 1115 # randomly sample endomorphism algebra 1116 repeat 1117 coeffs:=List([1..enddim], x -> Random(F)); 1118 until ForAny(remain,x->not IsZero(coeffs[x])); 1119 1120 a:=LinearCombination(basis,coeffs); 1121 1122 rk:=RankMat(a); 1123 if rk=n then 1124 # a regular element, check to see whether its order is 1125 # larger than previously known, and if so whether it 1126 # generates the residue field modulo current nilpotent ideal 1127 order:=Order(a); 1128 1129 while (order mod p = 0) do 1130 order:=order / p; 1131 od; 1132 if order > maxorder then 1133 maxorder:=order; 1134 maxa:=a; 1135 if proveIndecomposability() then 1136 return fail; 1137 fi; 1138 cnt:=1; 1139 else 1140 cnt:=cnt + 1; 1141 fi; 1142 else 1143 fit:=FittingSplitModule(a,rk,F); 1144 if fit<>fail then 1145 return fit; 1146 elif addnilpotent(a) then 1147 # new nilpotent element, added to nilbasis. Now close nilbasis to 1148 # basis for an ideal. 1149 1150 # keep a pointer to the first new element added to nilbase 1151 pos:=nildim; # a was just added 1152 1153 # first add powers of a 1154 newa:=a^2; 1155 repeat 1156 lastdim:=nildim; 1157 addnilpotent(newa); 1158 newa:=newa * a; 1159 until lastdim = nildim or IsZero(newa); 1160 1161 # now close nilbase to make ideal basis 1162 repeat 1163 for i in [1..enddim] do 1164 a:=nilbase[pos] * basis[i]; 1165 fit:=FittingSplitModule(a,RankMat(a),F); 1166 if fit <> fail then 1167 return fit; 1168 fi; 1169 addnilpotent(a); 1170 od; 1171 pos:=pos + 1; 1172 until pos = nildim + 1; 1173 fi; 1174 fi; 1175 1176 if proveIndecomposability() then 1177 return fail; 1178 fi; 1179 until (cnt >= 20000); 1180 Error("Unable to ascertain module decomposition within time limits.\n", 1181 "Call `return;' to try again."); 1182 cnt:=0; 1183 until false; 1184end); 1185 1186 1187BindGlobal("SMTX_Indecomposition",function(m) 1188local n, F, stack, i, d, d2, md, b, endo, sel, e1, e2; 1189 if not IsBound(m.indecomposition) then 1190 n:=m.dimension; 1191 F:=m.field; 1192 stack:=[[IdentityMat(n,F),m]]; 1193 i:=1; 1194 while i<=Length(stack) do 1195 d:=ProperModuleDecomp(stack[i][2]); 1196 if d<>fail then 1197 if Length(stack[i][1])<n then 1198 d2:=List(d,j->j*stack[i][1]); 1199 else 1200 d2:=d; 1201 fi; 1202 md:=List(d2,i->SMTX.InducedActionSubmodule(m,i)); 1203 Assert(1,ForAll(md,i->i<>fail)); 1204 # Translate endomorphism rings 1205 b:=Concatenation(d[1],d[2]); # local new basis 1206 # basechange 1207 endo:=List(stack[i][2].basisModuleEndomorphisms, 1208 i->b*i/b); 1209 sel:=[1..Length(d[1])]; 1210 e1:=List(endo,i->i{sel}{sel}); 1211 e1:=SMTX_EcheloniseMats(e1,F)[1]; 1212 Assert(1,ForAll(md[1].generators,i->ForAll(e1,j->i*j=j*i))); 1213 md[1].basisModuleEndomorphisms:=e1; 1214 sel:=[Length(d[1])+1..stack[i][2].dimension]; 1215 e2:=List(endo,i->i{sel}{sel}); 1216 e2:=SMTX_EcheloniseMats(e2,F)[1]; 1217 Assert(1,ForAll(md[2].generators,i->ForAll(e2,j->i*j=j*i))); 1218 md[2].basisModuleEndomorphisms:=e2; 1219 stack[i]:=[d2[1],md[1]]; 1220 Add(stack,[d2[2],md[2]]); 1221 else 1222 SMTX.SetIsIndecomposable(stack[i][2],true); 1223 i:=i+1; 1224 fi; 1225 od; 1226 m.indecomposition:=stack; 1227 fi; 1228 return m.indecomposition; 1229end); 1230 1231SMTX.Indecomposition:=SMTX_Indecomposition; 1232 1233 1234# Check isomorphism of indecomposable modules. 1235# 1236# If they are isomorphic then the homomorphism space between them is a 1237# disguised copy of the endomorphism algebra. This is a local algebra, 1238# and hence all singular elements are nilpotent. Certainly it cannot 1239# have a basis consisting entirely of nilpotent elements (a theorem of 1240# Wedderburn), so at least one basis element for Hom(M1,M2) must be an 1241# isomorphism if they are isomorphic. 1242BindGlobal("IsomIndecModules",function (M1, M2) 1243local base, i,n; 1244 1245 if not (SMTX.IsIndecomposable(M1) and SMTX.IsIndecomposable(M2)) then 1246 Error("IsomIndecModules: requires indecomposable modules"); 1247 fi; 1248 1249 n:=M1.dimension; 1250 # module dimensions certainly must match 1251 if n<>M2.dimension or 1252 1253 # their endomorphism algebras must have same dimension 1254 Length(SMTX.BasisModuleEndomorphisms(M1)) <> 1255 Length(SMTX.BasisModuleEndomorphisms(M2)) or 1256 1257 (SMTX.BasisEndomorphismsRadical(M1)<>fail and 1258 SMTX.BasisEndomorphismsRadical(M2)<>fail and 1259 Length(SMTX.BasisEndomorphismsRadical(M1))<> 1260 Length(SMTX.BasisEndomorphismsRadical(M2)) ) then 1261 return fail; 1262 fi; 1263 # the easy options have run out 1264 1265 # Last case, both modules are idecomposable but not necessarily irreducible. 1266 # In this case, compute Hom and look for isom in the basis. 1267 1268 base:=SMTX.BasisModuleHomomorphisms(M1, M2); 1269 1270 for i in base do 1271 if RankMat(i) = n then 1272 return i; 1273 fi; 1274 od; 1275 1276 return fail; 1277 1278end); 1279 1280BindGlobal("SMTX_HomogeneousComponents",function(m) 1281local d, h, found, i, m1, idx, imgs, hom, j; 1282 d:=SMTX.Indecomposition(m); 1283 h:=[]; 1284 found:=[]; 1285 i:=1; 1286 while Length(found)<Length(d) do 1287 if not i in found then 1288 m1:=d[i][2]; 1289 idx:=[i]; 1290 AddSet(found,i); 1291 imgs:=[]; 1292 for j in [i+1..Length(d)] do 1293 if not j in found and m1.dimension=d[j][2].dimension then 1294 hom:=IsomIndecModules(d[j][2],m1); 1295 if hom<>fail then 1296 Add(idx,j); 1297 AddSet(found,j); 1298 Add(imgs,rec(component:=d[j],isomorphism:=hom^-1)); 1299 fi; 1300 fi; 1301 od; 1302 Add(h,rec(component:=d[i],images:=imgs,indices:=idx)); 1303 fi; 1304 i:=i+1; 1305 od; 1306 return h; 1307end); 1308 1309SMTX.HomogeneousComponents:=SMTX_HomogeneousComponents; 1310 1311 1312# Test for isomorphism of modules. Will return one of: 1313# 1314# (1) the isomorphism as an F-matrix between M1 and M2 1315# (2) fail if the two modules are definitely not isomorphic 1316# 1317# Note that the isomorphism X is such that conjugating each generator 1318# acting on M1 by X gives the corresponding action on M2. Therefore 1319# X^-1 is a matrix whose rows correspond to a new basis of M1 that 1320# duplicates the action of M2 on M1. 1321# 1322# If necessary, uses the decomposition into indecomposable summands. A 1323# homogeneous component is a direct sum of multiple copies of a single 1324# indecomposable summand. The homogeneous components must match between 1325# each module, with their multiplicities. 1326BindGlobal("SMTX_IsomorphismModules",function (M1, M2) 1327local F, n, hc1, hc2, nc, b1, b2, map, remain, j, found, hom, i, k; 1328 1329 TestModulesFitTogether(M1,M2); 1330 F:=M1.field; 1331 n:=M1.dimension; 1332 1333 if n <> M2.dimension then 1334 # Modules have different dimensions 1335 return fail; 1336 elif (SMTX.BasisEndomorphismsRadical(M1)<>fail and 1337 SMTX.BasisEndomorphismsRadical(M2)<>fail and 1338 Length(SMTX.BasisEndomorphismsRadical(M1))<> 1339 Length(SMTX.BasisEndomorphismsRadical(M2)) ) then 1340 # different endomorphism algebra dimensions 1341 return fail; 1342 fi; 1343 1344 hc1:=SMTX.HomogeneousComponents(M1); 1345 hc2:=SMTX.HomogeneousComponents(M2); 1346 1347 nc:=Length(hc1); 1348 if nc <> Length(hc2) then 1349 return fail; 1350 fi; 1351 1352 # build bases that must be mapped to each other iteratively 1353 b1:=[]; 1354 b2:=[]; 1355 map:=[]; 1356 1357 remain:=[1..nc]; 1358 for i in [1..nc] do 1359 j:=1;found:=false; 1360 while j<=nc and not found do 1361 if j in remain and Length(hc1[i].indices)=Length(hc2[j].indices) then 1362 # test: i isomorphic j? 1363 hom:=IsomIndecModules(hc1[i].component[2],hc2[j].component[2]); 1364 if hom<>fail then 1365 # the homogeneous components are isomorphic 1366 found:=true; 1367 Append(b1,hc1[i].component[1]); 1368 Append(b2,hc2[j].component[1]); 1369 Add(map,hom); 1370 for k in [1..Length(hc1[i].images)] do 1371 Append(b1,hc1[i].images[k].component[1]); 1372 Append(b2,hc2[j].images[k].component[1]); 1373 Add(map,hc1[i].images[k].isomorphism^-1*hom* 1374 hc2[j].images[k].isomorphism); 1375 od; 1376 fi; 1377 fi; 1378 j:=j+1; 1379 od; 1380 if found=false then 1381 # one homogeneous component has no image -- the modules cannot be 1382 # isomorphic 1383 return fail; 1384 fi; 1385 od; 1386 b1:=ImmutableMatrix(M1.field,b1); 1387 b2:=ImmutableMatrix(M1.field,b2); 1388 return b1^-1*ImmutableMatrix(M1.field,DirectSumMat(map))*b2; 1389end); 1390 1391SMTX.IsomorphismModules:=SMTX_IsomorphismModules; 1392 1393# Note: matalg is a basis for a nilpotent matrix algebra whose elements 1394# are all in lower diagonal form (zeros on the main diagonal). 1395# 1396# Echelonisation indices are chosen as the earliest non-zero entries 1397# running down diagonals below the main diagonal: 1398# [2,1], [3,2], [4,3], ..., [3,1], [4,2], ..., [n-1,1], [n, 2], [n,1] 1399BindGlobal("SMTX_EcheloniseNilpotentMatAlg",function (matalg, F) 1400local zero, n, flags, base, ech, k, diff, i, j, found, l; 1401 1402 zero:=Zero(F); 1403 n := Length(matalg[1][1]); 1404 flags := NullMat(n,n); 1405 1406 base := matalg; 1407 ech := []; 1408 k := 1; 1409 1410 while k <= Length(base) do 1411 diff := 1; 1412 i := 2; j := i - diff; 1413 found := false; 1414 while not found and diff < n do 1415 if (base[k][i][j] <> zero) and 1416 (flags[i][j] = 0) then 1417 found := true; 1418 else 1419 i := i + 1; 1420 j := i - diff; 1421 if (i > n) then 1422 diff := diff + 1; 1423 i := diff + 1; 1424 j := i - diff; 1425 fi; 1426 fi; 1427 od; 1428 1429 if found then 1430 1431 # Now basis element k will have echelonisation index [i,j] 1432 Add(ech, [i,j]); 1433 1434 # First normalise the [i,j] position to 1 1435 base[k] := base[k] / base[k][i][j]; 1436 1437 # Now zero position [i,j] in all other basis elements 1438 for l in [1..Length(base)] do 1439 if (l <> k) and (base[l][i][j] <> zero) then 1440 base[l] := base[l] - base[k] * base[l][i][j]; 1441 fi; 1442 od; 1443 k := k + 1; 1444 1445 else 1446 # no non-zero element found, delete from list 1447 base := base{ Concatenation([1..k-1], [k+1..Length(base)])}; 1448 fi; 1449 od; 1450 return [base, ech]; 1451end); 1452 1453# compute a change of basis that exhibits the matrix algebra 1454# defined by the basis 'matalg' in triangular form. 1455BindGlobal("SMTX_NilpotentBasis",function (matalg) 1456local decompose, field, Y, mats, newbase; 1457 1458 decompose := function ( m, b ) 1459 local n, subs, vs, vsi,rep, newm,j,ran; 1460 1461 if Length(m) = 0 then 1462 # all action is now zero, so append current full basis and 1463 # finish up 1464 Append(Y, b); 1465 else 1466 1467 n := Length(m[1][1]); 1468 1469 # find the intersection of the nullspaces 1470 subs:=NullspaceMat(m[1]); 1471 for j in [2..Length(m)] do 1472 subs:=SumIntersectionMat(subs,NullspaceMat(m[j]))[2]; 1473 od; 1474 1475 1476 # Use matrix group routine to compute action of nilpotent 1477 # matrices on the quotient vectorspace 1478 vs := BaseSteinitzVectors(IdentityMat(n,field),subs); 1479 vs:=Concatenation(vs.subspace,vs.factorspace); 1480 vs:=ImmutableMatrix(field,vs); 1481 vsi:=vs^-1; 1482 ran:=[Length(subs)+1..n]; 1483 rep:=List(m,i->vs*i*vsi); 1484 rep:=List(rep,i->i{ran}{ran}); 1485 1486 # Take a copy of the non-zero matrices acting on the quotient space 1487 # 1488 newm := Filtered(rep,x->not IsZero(x)); 1489 1490 Append(Y, subs * b); 1491 decompose( newm, vs{ran} * b ); 1492 1493 fi; 1494 end; 1495 1496 # return empty list if empty matrix list 1497 if Length(matalg) = 0 then return []; fi; 1498 1499 field := DefaultField(matalg[1][1]); 1500 1501 Y := []; 1502 1503 decompose( matalg, IdentityMat(Length(matalg[1][1]), field)); 1504 # 1505 # Y is the change of basis matrix 1506 1507 if Length(matalg) > 0 then 1508 mats := Y * matalg / Y; 1509 fi; 1510 # 1511 # mats is now a list of matrices in lower triangular form 1512 1513 # echelonise them along lower diagonals 1514 # 1515 newbase := SMTX_EcheloniseNilpotentMatAlg(mats, field)[1]; 1516 1517 return [newbase, Y]; 1518 1519end); 1520 1521 1522# module automorphism group 1523BindGlobal("SMTX_ModuleAutomorphisms",function(m) 1524 local f, h, hb, hbi, bas, auts, autorder, dim, nb, nbi, r, q, w, Fqr, gl, a, subm, nilbase, homs, sub, endos, au, aubas, it, coeffs, cnt, ol, ind, i, j, g, k; 1525 f:=m.field; 1526 h:=MTX.HomogeneousComponents(m); 1527 # construct basis for each homogeneous component 1528 hb:=[]; 1529 for i in h do 1530 # basis of component 1531 bas:=ShallowCopy(i.component[1]); 1532 for j in i.images do 1533 #Append(bas,LeftQuotient(j.isomorphism,j.component[1])); 1534 Append(bas,j.isomorphism*j.component[1]); 1535 od; 1536 #bas:=MTX.NormedBasisAndBaseChange(bas)[1]; 1537 Add(hb,bas); 1538 od; 1539 1540 # each homogeneous component separately 1541 auts:=[]; 1542 autorder:=1; 1543 for i in [1..Length(h)] do 1544 # basis of component 1545 bas:=hb[i]; 1546 dim:=h[i].component[2].dimension; 1547 nb:=Concatenation(bas,Concatenation(hb{Difference([1..Length(h)],[i])})); 1548 nb:=ImmutableMatrix(f,nb); 1549 nbi:=nb^-1; 1550 1551 # start by building those automorphisms that fix the homogeneous 1552 # components - ie, do not involve maps from M_i to M_j unless 1553 # M_i is the same isomorphism type as M_j 1554 r:=Length(h[i].indices); 1555 # first the subgroup GL(multiplicity, residue field) 1556 q:=SMTX.EndAlgResidue(h[i].component[2]); 1557 w:=q[1]; 1558 q:=q[2]+1; 1559 Fqr:=PrimitiveElement(GF(q)); 1560 gl:=GL(r,q); 1561 autorder:=autorder*Size(gl); 1562 Info(InfoMtxHom,3,"increase by gl",Size(gl)," ",autorder); 1563 for g in GeneratorsOfGroup(gl) do 1564 a:=IdentityMat(m.dimension,f); 1565 for j in [1..r] do 1566 for k in [1..r] do 1567 if IsZero(g[j][k]) then 1568 subm:=w*0; 1569 else 1570 subm:=w^LogFFE(g[j][k],Fqr); 1571 fi; 1572 a{[(j-1)*dim+1..j*dim]}{[(k-1)*dim+1..k*dim]}:=subm; 1573 od; 1574 od; 1575 a:=nbi*a*nb; 1576 Assert(1,ForAll(m.generators,i->i*a=a*i)); 1577 Add(auts,a); 1578 od; 1579 1580 # now the subgroup { I + Y | Y in S } where S generates the radical 1581 # of the endomorphism algebra as a circle group 1582 nilbase:=SMTX.BasisEndomorphismsRadical(h[i].component[2]); 1583 if Length(nilbase)>0 then 1584 nilbase:=SMTX_NilpotentBasis(nilbase); 1585 nilbase:=nilbase[2]^-1*nilbase[1]*nilbase[2]; 1586 fi; 1587 a:=(Size(f)^Length(nilbase))^(r^2); 1588 autorder := autorder * a; 1589 Info(InfoMtxHom,3,"increase by radical",a," ",autorder); 1590 1591 for j in nilbase do; 1592 a:=IdentityMat(m.dimension,f); 1593 subm:=IdentityMat(dim,f)+j; 1594 a{[1..dim]}{[1..dim]}:=subm; 1595 a:=nbi*a*nb; 1596 Assert(1,ForAll(m.generators,i->i*a=a*i)); 1597 Add(auts,a); 1598 od; 1599 1600 # Now the automorphisms that act trivially when restricted to 1601 # each homogeneous component, but which include action between 1602 # homogeneous components via elements of Hom(M_i, M_j) 1603 for j in [1..Length(h)] do 1604 if i <> j then 1605 homs:=SMTX.BasisModuleHomomorphisms(h[i].component[2], 1606 h[j].component[2]); 1607 if Length(homs) > 0 then 1608 hbi:=0; 1609 for k in [1..j-1] do 1610 hbi:=hbi+Length(hb[k]); 1611 od; 1612 if i>j then 1613 hbi:=hbi+Length(hb[i]); 1614 fi; 1615 hbi:=hbi+[1..h[j].component[2].dimension]; 1616 1617 a:=(Size(f)^Length(homs))^(r*Length(h[j].indices)); 1618 autorder:=autorder*a; 1619 Info(InfoMtxHom,3,"increase by mixing ",j,":",a," ",autorder); 1620 for k in homs do 1621 a:=IdentityMat(m.dimension,f); 1622 a{[1..dim]}{hbi}:=k; 1623 a:=nbi*a*nb; 1624 Assert(1,ForAll(m.generators,i->i*a=a*i)); 1625 Add(auts,a); 1626 od; 1627 fi; 1628 fi; 1629 od; 1630 od; 1631 1632 if Length(auts)=0 then 1633 return Group(auts,IdentityMat(m.dimension,f)); 1634 else 1635 a:=Group(auts); 1636 Assert(1,Size(a)=autorder); 1637 SetSize(a,autorder); 1638 return a; 1639 fi; 1640 1641end); 1642 1643SMTX.ModuleAutomorphisms:=SMTX_ModuleAutomorphisms; 1644 1645SMTX.SetIsIndecomposable:=function(m,b) 1646 m.isIndecomposable:=b; 1647end; 1648 1649SMTX.HasIsIndecomposable:=function(m) 1650 return IsBound(m.isIndecomposable); 1651end; 1652 1653SMTX.IsIndecomposable:=function(m) 1654 if not SMTX.HasIsIndecomposable(m) then 1655 m.isIndecomposable:=Length(SMTX.Indecomposition(m))=1; 1656 fi; 1657 return m.isIndecomposable; 1658end; 1659 1660 1661SMTX_BasisModuleHomomorphisms:=function(m1,m2) 1662local b; 1663 TestModulesFitTogether(m1,m2); 1664 if m1.dimension>5 then 1665 b:= SpinHom(m1,m2); 1666 Assert(1,Length(b)=Length(SmalldimHomomorphismsModules(m1,m2))); 1667 else 1668 b:= SmalldimHomomorphismsModules(m1,m2); 1669 fi; 1670 Assert(1,ForAll([1..Length(m1.generators)], 1671 i->ForAll(b,j->m1.generators[i]*j=j*m2.generators[i]))); 1672 return b; 1673end; 1674 1675SMTX.BasisModuleHomomorphisms:=SMTX_BasisModuleHomomorphisms; 1676 1677SMTX_BasisModuleEndomorphisms:=function(m) 1678 if not IsBound(m.basisModuleEndomorphisms) then 1679 m.basisModuleEndomorphisms:=Immutable(SMTX.BasisModuleHomomorphisms(m,m)); 1680 fi; 1681 return m.basisModuleEndomorphisms; 1682end; 1683 1684SMTX.BasisModuleEndomorphisms:=SMTX_BasisModuleEndomorphisms; 1685 1686SMTX.SetBasisEndomorphismsRadical:=SMTX.Setter("basisEndoRad"); 1687SMTX.BasisEndomorphismsRadical:=SMTX.Getter("basisEndoRad"); 1688 1689SMTX.SetEndAlgResidue:=SMTX.Setter("endAlgResidue"); 1690SMTX.EndAlgResidue:=SMTX.Getter("endAlgResidue"); 1691 1692if IsHPCGAP then 1693 MakeReadOnlyObj(SMTX); 1694fi; 1695