1############################################################################# 2## 3#W recograt.gd matgrp package Alexander Hulpke 4## 5## 6#Y Copyright (C) 2014-18, Alexander Hulpke 7## 8## basic setup for matrix fitting free. 9## 10 11SetInfoLevel(InfoRecog,0); # recog will print status messages otherwise 12 13# mod to genss -- rings 14FindHomMethodsMatrix.Nonfield := function(ri, G) 15 if IsBound(ri!.ring) and not IsBound(ri!.field) then 16 Error("hereIAm"); 17 fi; 18 return false; 19end; 20 21AddMethod( FindHomDbMatrix, FindHomMethodsMatrix.Nonfield, 22 5100, "Nonfield", 23 "catch matrix groups defined over nonfield rings" ); 24 25OnSubmoduleCosets:=function(cset,g) 26local v; 27 return [SiftedVector(cset[2],cset[1]*g),cset[2]]; 28end; 29 30MakeSubmoduleCosetAction:=function(basis) 31 return function(vec,g) 32 return SiftedVector(basis,vec*g); 33 end; 34end; 35 36MakeSubmoduleColineAction:=function(basis) 37 return function(vec,g) 38 local c; 39 vec:=SiftedVector(basis,vec*g); 40 c:=PositionNonZero(vec); 41 if c<=Length(vec) then 42 vec := Inverse( vec[c] ) * vec; 43 fi; 44 return vec; 45 end; 46end; 47 48FUNCSPACEHASH:=[]; 49# mod to genss -- use submodules and quotients for base points 50 51MSSFBPC:=function( grp, opt, ii, parentS ) # owf fct to call easily 52local F,dim,orb,orbs,i,fct, 53mo,cs,j,k,dims,bas,basc,basinv,nb,lastdim,cand,fcand,sel,limit,trysel,submodule; 54 55 trysel:=function(recsub,recfac) 56 local lgens; 57 # nor the trivial action 58 if ForAny(mo,x->Order(x{sel}{sel})>1) then 59 Info(InfoFFMat,2,"range ",sel," have ",Length(cand.points)); 60 lgens:=List(mo,x->x{sel}{sel}); 61 fcand:=FindBasePointCandidates(Group(lgens),opt,ii, 62 false:Subrecurse:=recsub,Facrecurse:=recfac); 63 Info( InfoGenSS, 3, "Subfactor module of range ",sel,", ",Length(fcand.ops), 64 " candidates"); 65 for k in [1..Length(fcand.ops)] do 66 if ForAll(lgens,x->fcand.ops[k](fcand.points[k],x)=fcand.points[k]) then 67 Info(InfoFFMat,2,"Ignoring fixed element for base"); 68 69 elif fcand.ops[k]=OnRight or fcand.ops[k]=OnPoints or fcand.ops[k]=OnLines then 70 nb:=fcand.points[k]*bas{sel}; 71 if lastdim=0 then 72 # proper subspace -- just vectors 73 Add(cand.points,nb); 74 Add(cand.ops,fcand.ops[k]); 75 elif true then 76 # # action on cosets 77 submodule:=SemiEchelonBasis(VectorSpace(F,bas{[1..lastdim]},Zero(bas[1]))); 78 nb:=SiftedVector(submodule,nb); 79 if fcand.ops[k]=OnLines then 80 fct:=MakeSubmoduleColineAction(submodule); 81 Add(FUNCSPACEHASH,[fct,submodule]); 82 else 83 fct:=MakeSubmoduleCosetAction(submodule); 84 Add(FUNCSPACEHASH,[fct,submodule]); 85 fi; 86 Add(cand.points,nb); 87 Add(cand.ops,fct); 88 elif Length(sel)=1 then 89 # TODO: 1-dim factor -- need to do cosets 90 Info(InfoWarning,1,"Case not yet implemented"); 91 else 92 # subfactor -- take subspace preimage 93 nb:=OnSubspacesByCanonicalBasis(Concatenation(bas{[1..lastdim]},[nb]), 94 One(grp)); 95 Add(cand.points,nb); 96 Add(cand.ops,OnSubspacesByCanonicalBasis); 97 fi; 98 elif ForAny(FUNCSPACEHASH,x->x[1]=fcand.ops[k]) then 99 fct:=First(FUNCSPACEHASH,x->x[1]=fcand.ops[k]); 100 submodule:=SemiEchelonBasis(VectorSpace(F, 101 Concatenation(bas{[1..lastdim]}, 102 BasisVectors(fct[2])*bas{sel}))); 103 Add(cand.points,fcand.points[k]*bas{sel}); 104 fct:=MakeSubmoduleCosetAction(submodule); 105 Add(FUNCSPACEHASH,[fct,submodule]); 106 Add(cand.ops,fct); 107 Info(InfoFFMat,2,"ACTPOP"); 108 elif fcand.ops[k]=OnSubspacesByCanonicalBasis then 109 nb:=fcand.points[k]*bas{sel}; 110 nb:=OnSubspacesByCanonicalBasis(Concatenation(bas{[1..lastdim]},nb), 111 One(grp)); 112 Add(cand.points,nb); 113 Add(cand.ops,OnSubspacesByCanonicalBasis); 114 else 115 Info(InfoWarning,1,"Action not recognized"); 116 fi; 117 od; 118 return true; 119 else 120 return false; 121 fi; 122 end; 123 124 if IsBound(opt.VeryShortOrbLimit) then 125 limit:=2*opt.VeryShortOrbLimit; 126 else 127 limit:=10^4; 128 fi; 129 130 F := DefaultFieldOfMatrixGroup(grp); 131 132 # don't bother in small cases 133 if Size(F)^Length(One(grp))<=opt.ShortOrbitsOrbLimit then 134 TryNextMethod(); 135 fi; 136 137 cs:=GeneratorsOfGroup(grp); 138 if ForAny(cs,IsObjWithMemory) then 139 cs:=Concatenation(Filtered(cs,x->not IsObjWithMemory(x)), 140 List(Filtered(cs,x->IsObjWithMemory(x)), 141 x->x!.el)); 142 fi; 143 144 cand:=rec(ops:=[],points:=[],used:=0); 145 mo:=GModuleByMats(cs,F); 146 dim:=mo.dimension; 147 if MTX.IsIrreducible(mo) or Size(F)^dim<=limit then 148 TryNextMethod(); 149 fi; 150 151 # build new basis corresponding to comp ser. 152 cs:=MTX.BasesCompositionSeries(mo); 153 Info(InfoFFMat,2,"dims=",List(cs,Length)); 154 dims:=List(cs,Length); 155 bas:=[]; 156 basc:=[]; 157 for j in [2..Length(cs)] do 158 nb:=BaseSteinitzVectors(cs[j],basc); 159 Append(bas,nb.factorspace); 160 basc:=Concatenation(nb.subspace,nb.factorspace); 161 Sort(basc,function(a,b) return PositionNonZero(a)<PositionNonZero(b);end); 162 od; 163 basinv:=bas^-1; 164 165 mo:=List(mo.generators,x->bas*x*basinv); 166 167 # now step up in sizes indicated by short orbit lengths 168 lastdim:=0; 169 j:=2; 170 if ValueOption("Subrecurse")<>false then 171 while j<=Length(dims) and Length(cand.points)<=5 do 172 # don't bother is the space is too small 173 if (j=Length(dims) and lastdim>0) or 174 (j<Length(dims) and Size(F)^(dims[j+1]-lastdim)>limit) then 175 sel:=[lastdim+1..dims[j]]; 176 if trysel(false,true) then 177 # we tried a space 178 lastdim:=dims[j]; 179 fi; 180 181 fi; 182 183 j:=j+1; 184 od; 185 fi; 186 187 if lastdim=0 and ValueOption("Facrecurse")<>false then 188 189 # all but last submodule was trivial -- step down on factors 190 j:=Length(dims)-1; 191 while j>0 and Length(cand.points)<=5 do 192 lastdim:=dims[j]-1; 193 if (j>1 and Size(F)^(dim-dims[j-1])>limit) then 194 sel:=[lastdim+1..dim]; 195 if trysel(false,false) then 196 j:=0; 197 fi; 198 fi; 199 j:=j-1; 200 od; 201 202 fi; 203 204 if Length(cand.points)>0 then return cand; fi; 205 206 # refer to parent 207 if IsBound(opt.PCand) then 208 # try which ones are short 209 sel:=[]; 210 for i in [1..Length(opt.PCand.points)] do 211 orb:=[opt.PCand.points[i]]; 212 mo:=opt.PCand.ops[i]; 213 orbs:=Set(orb); # as short, set is fine. 214 j:=1; 215 while j<=Length(orb) and Length(orb)<=limit do 216 for k in GeneratorsOfGroup(grp) do 217 cs:=mo(orb[j],k); 218 if not cs in orbs then 219 Add(orb,cs); 220 AddSet(orbs,cs); 221 fi; 222 od; 223 j:=j+1; 224 od; 225 if Length(orb)<=limit then 226 Add(sel,i); 227 Add(cand.points,orb[1]); 228 Add(cand.ops,mo); 229 fi; 230 od; 231 Info(InfoFFMat,2,"Selected ",Length(sel)," of ",Length(opt.PCand.points)," group basis vectors"); 232 if Length(cand.points)>0 then return cand; fi; 233 fi; 234 235 TryNextMethod(); 236end; 237 238InstallMethod(FindBasePointCandidates, 239 "for reducible matrix group over a FF, use submodules and quotients", 240 [ IsGroup and IsMatrixGroup and IsFinite, IsRecord, IsInt, IsObject ], 21, 241 MSSFBPC); 242 243GetInformationFromRecog:=function(recog) 244local treerecurse,n,factors,homs,leafgens,niceranges,genum,sz,leafs,g, 245 minranges,mingens,permap,furtherhom,fn,extragens,extranum,extra,allgens, 246 nicegp,map,stc,orbit,gd,cnt; 247 248 # the main worker function -- traverse the tree. 249 treerecurse:=function(r,h) 250 local 251 hom,f,k,sel,u,i,j,x,cs,nsol,xf,bigger,nicehom,ngi,stbc,opt,act,ng,dom,v; 252 if IsLeaf(r) then 253 f:=Length(NiceGens(r)); 254 Info(InfoFFMat,2,"Leaf",Size(r)," ",genum," ",f); 255 k:=[genum+1..genum+f]; 256 if Size(r)=1 then 257 Info(InfoFFMat,2,"ignoring trivial factor"); 258 elif Length(Set(Factors(Size(r))))<3 then 259 Info(InfoFFMat,2,"ignoring solvable factor of order ",Size(r)); 260 261 # store info 262 n:=n+Length(Factors(Size(r))); 263 SetIsSolvableGroup(Grp(r),true); 264 permap[n]:=IdentityMapping(Grp(r)); 265 leafs[n]:=r; 266 homs[n]:=h; 267 factors[n]:=false; 268 sz[n]:=Size(r); 269 leafgens[n]:=f; 270 271 sel:=[1..Length(NiceGens(r))]; 272 mingens[n]:=sel; 273 niceranges[n]:=k; 274 minranges[n]:=k{sel}; 275 276 else 277 nicehom:=false; 278 if IsPermGroup(Grp(r)) then 279 nicehom:=IdentityMapping(Grp(r)); 280 elif IsBound(r!.stabilizerchain) then 281 #use ActionOnOrbit? 282 stbc:=BaseStabilizerChain(r!.stabilizerchain); 283 act:=stbc.ops[1]; 284 if ForAll(stbc.ops,x->x=act) then 285 # all the same action -- union 286 orbit:=[]; 287 for j in [1..Length(stbc.points)] do 288 if not stbc.points[j] in orbit then 289 orbit:=Union(orbit,Orbit(Grp(r),stbc.points[j],stbc.ops[j])); 290 fi; 291 od; 292 if Length(orbit)>1 and Length(orbit)<50000 then 293 nicehom:=ActionHomomorphism(Grp(r),orbit,stbc.ops[1],"surjective"); 294 Info(InfoFFMat,2,"Got degree ",Length(orbit)); 295 fi; 296 fi; 297 fi; 298 299 300 if nicehom=false then 301# Hasfhmethsel, success method: StabilizerChain 302 303 if IsBound(r!.projective) and r!.projective then 304 g:=Grp(r); 305 v:=OnLines(g.1[1],One(g)); 306 dom:=ShallowCopy(Orbit(g,v,OnLines)); 307 nicehom:=ActionHomomorphism(g,dom,OnLines,"surjective"); 308 while Size(Image(nicehom))<Size(r) do 309 cnt:=0; 310 repeat 311 v:=OnLines(Random(DefaultFieldOfMatrixGroup(g)^Length(One(g))), 312 One(g)); 313 cnt:=cnt+1; 314 if cnt>1000 then Error("no vector found");fi; 315 until not v in dom; 316 Append(dom,Orbit(g,v,OnLines)); 317 nicehom:=ActionHomomorphism(g,dom,OnLines,"surjective"); 318 od; 319 else 320 nicehom:=IsomorphismPermGroup(Grp(r)); 321 fi; 322 fi; 323 g:=Image(nicehom); 324 if Size(g)<>Size(r) then 325 Error("some discrepancy happened"); 326 fi; 327 328 #test!! 329 # .isknownsimple, .isknownalmostsimple 330 if IsBound(r!.comment) and r!.comment[1]<>'_' then 331 SetIsSimpleGroup(g,true); 332 gd:=g; 333 elif not IsSolvableGroup(g) then 334 gd:=PerfectResiduum(g); 335 else 336 gd:=fail; 337 fi; 338 339 if gd<>fail then 340 if NrMovedPoints(g)> SufficientlySmallDegreeSimpleGroupOrder(Size(gd)) 341 then 342 nicehom:=nicehom*SmallerDegreePermutationRepresentation(g); 343 gd:=Image(nicehom); 344 Info(InfoFFMat,2,"Improved degree ",NrMovedPoints(g),"->", 345 NrMovedPoints(gd)); 346 if HasIsSimpleGroup(g) and IsSimpleGroup(g) then 347 SetIsSimpleGroup(gd,true); 348 SetIsSolvableGroup(gd,false); 349 fi; 350 g:=gd; 351 fi; 352 fi; 353 354 ## indicate unsolvable 355 #nsol:=IsSimpleGroup(g) or 356 # (Length(Set(Factors(Size(r))))>2 and not IsSolvableGroup(g)); 357 358 359 if IsSolvableGroup(g) then 360 Info(InfoFFMat,2,"Ignoring solvable factor of order ",Size(g)); 361 n:=n+Length(Factors(Size(g))); 362 elif not IsSimpleGroup(g) then 363 Info(InfoFFMat,2,"doing size",Size(g)," ",IsSimpleGroup(g)); 364 365 # not simple -- split 366 cs:=CompositionSeries(g); 367 ng:=List(NiceGens(r),x->ImagesRepresentative(nicehom,x)); 368 for i in [2..Length(cs)] do 369 extra:=[]; 370 n:=n+1; 371 # test generators 372 u:=cs[i]; 373 sel:=[]; 374 for j in [1..f] do 375 x:=ng[j]; 376 if x in cs[i-1] and not x in u then 377 Add(sel,j); 378 u:=ClosureSubgroup(u,x); 379 fi; 380 od; 381 382 nicegp:=Group(ng); 383 map:=EpimorphismFromFreeGroup(nicegp); 384 385 if Size(u)<Size(cs[i-1]) then 386 Info(InfoFFMat,2,"cannot compatibilize generators -- add extras"); 387 while Size(u)<Size(cs[i-1]) do 388 x:=First(GeneratorsOfGroup(cs[i-1]),x->not x in u); 389 u:=ClosureSubgroup(u,x); 390 391 # decompose into word 392 xf:=Factorization(nicegp,x); 393 if xf=fail then Error("factorization error");fi; 394 x:=MappedWord(xf,MappingGeneratorsImages(map)[1],allgens{k}); 395 396 extragens[extranum]:=x; 397 Add(extra,extranum); 398 extranum:=extranum+1; 399 od; 400 fi; 401 hom:=NaturalHomomorphismByNormalSubgroup(cs[i-1],cs[i]); 402 if not IsAbelian(Image(hom)) then 403 permap[n]:=IsomorphismPermGroup(Image(hom)); 404 fi; 405 406 fn:=fn+1; 407 if Size(cs[i])=1 then 408 furtherhom[fn]:=nicehom; 409 else 410 furtherhom[fn]:=nicehom*hom; 411 fi; 412 leafs[n]:=false; 413 homs[n]:=Concatenation(h,[fn]); 414 factors[n]:=Image(hom); 415 sz[n]:=Size(Image(hom)); 416 leafgens[n]:=false; 417 niceranges[n]:=false; 418 minranges[n]:=Concatenation(k{sel},extra); 419 od; 420 421 else 422 # found a simple case 423 n:=n+1; 424 425 permap[n]:=nicehom; 426 # get smaller generating set 427 sel:=[]; 428 u:=TrivialSubgroup(Image(nicehom)); 429 for i in [1..Length(NiceGens(r))] do 430 ngi:=ImagesRepresentative(nicehom,NiceGens(r)[i]); 431 if not ngi in u then 432 u:=ClosureGroup(u,ngi); 433 Add(sel,i); 434 fi; 435 od; 436 437 leafs[n]:=r; 438 homs[n]:=h; 439 factors[n]:=g; 440 sz[n]:=Size(r); 441 leafgens[n]:=f; 442 443 mingens[n]:=sel; 444 niceranges[n]:=k; 445 minranges[n]:=k{sel}; 446 447 fi; 448 449 fi; 450 genum:=genum+f; 451 else 452 #hom:=Homom(r); 453 f:=RIFac(r); 454 treerecurse(f,Concatenation(h,[1])); 455 k:=RIKer(r); 456 if k<>fail then 457 treerecurse(k,Concatenation(h,[0])); 458 fi; 459 fi; 460 461 end; 462 463 extragens:=[]; 464 furtherhom:=[]; 465 fn:=2; 466 n:=0; 467 genum:=0; 468 factors:=[]; 469 homs:=[]; 470 leafgens:=[]; 471 leafs:=[]; 472 niceranges:=[]; 473 minranges:=[]; 474 mingens:=[]; 475 sz:=[]; 476 permap:=[]; 477 allgens:=NiceGens(recog); 478 479 extranum:=Length(allgens)+1; 480 treerecurse(recog,[]); 481 482 return rec( 483 group:=Grp(recog), 484 n:=n, 485 genum:=Length(allgens), 486 recog:=recog, 487 leafs:=leafs, 488 factors:=factors, 489 furtherhom:=furtherhom, 490 sz:=sz, 491 homs:=homs, 492 leafgens:=leafgens, 493 minranges:=minranges, 494 mingens:=mingens, 495 extragens:=extragens, 496 permap:=permap, 497 niceranges:=niceranges); 498end; 499 500# assume that hom i can be applied 501CSIImageHomNr:=function(csi,n,x) 502local r,h,i; 503 r:=csi.recog; 504 h:=csi.homs[n]; 505 for i in h do 506 if i=0 then 507 r:=RIKer(r); 508 elif i=1 then 509 x:=ImageElm(Homom(r),x); 510 r:=RIFac(r); 511 else 512 x:=ImageElm(csi.furtherhom[i],x); 513 fi; 514 od; 515 return x; 516end; 517 518# since nice ranges can contain negatives, we cannot simply sublist 519CSINiceGens:=function(csi,a) 520 if IsList(a) then 521 return List(a,x->CSINiceGens(csi,x)); 522 elif a<=csi.genum then 523 return NiceGens(csi.recog)[a]; 524 else 525 return csi.extragens[a]; 526 fi; 527end; 528 529CSIDepthElm:=function(csi,x) 530local i; 531 i:=1; 532 while i<=csi.n and 533 (not IsBound(csi.leafs[i]) or # early skipped multiple solvable factor 534 ( 535 (IsBool(csi.leafs[i]) or not (IsBound(csi.leafs[i]!.projective) and csi.leafs[i]!.projective) 536 or csi.leafs[i]!.projective=false) 537 and Order(CSIImageHomNr(csi,i,x))=1) 538 or 539 (not IsBool(csi.leafs[i]) and (IsBound(csi.leafs[i]!.projective) and csi.leafs[i]!.projective) and 540 Order(ImagesRepresentative(csi.permap[i],CSIImageHomNr(csi,i,x)))=1)) do 541 i:=i+1; 542 od; 543 return i; 544end; 545 546# to test whether an element act trivially projectively, it is not enough to 547# test on base, but also need to have sum of base 548# elm must be an element that is already image 549CSIProjectiveBases:=function(csi,i,elm) 550local m; 551 if not IsBound(csi.projbases) then 552 csi.projbases:=[]; 553 fi; 554 if not IsBound(csi.projbases[i]) then 555 m:=ShallowCopy(One(elm)); 556 Add(m,Sum(m)); 557 csi.projbases[i]:=m; 558 fi; 559 return csi.projbases[i]; 560end; 561 562CSIDepthElm:=function(csi,x) 563local i,ximg; 564 i:=1; 565 while i<=csi.n do 566 if not IsBound(csi.leafs[i]) or 567 ((IsBool(csi.leafs[i]) or not (IsBound(csi.leafs[i]!.projective) and csi.leafs[i]!.projective) 568 or csi.leafs[i]!.projective=false) 569 and Order(CSIImageHomNr(csi,i,x))=1) then 570 i:=i+1; 571 elif (not IsBool(csi.leafs[i]) and IsBound(csi.leafs[i]!.projective) and 572 csi.leafs[i]!.projective) then 573 ximg:=CSIImageHomNr(csi,i,x); 574 if ForAll(CSIProjectiveBases(csi,i,ximg),z->OnLines(z,ximg)=z) then 575 i:=i+1; 576 else 577 return i; 578 fi; 579 else 580 return i; 581 fi; 582 od; 583 return i; 584end; 585 586CSIAelement:=function(a,localgens,l) 587local limg,rep; 588 limg:=List(l,x->Image(a[2],x)); 589 rep:= RepresentativeAction(a[1],localgens,limg,OnTuples); 590 if rep=fail then 591 Error("no rep"); 592 fi; 593 return rep; 594end; 595 596FindAct:=function(csi) 597local csinice,dci,pool, 598act,n,i,j,k,l,gens,lgens,c,d,genims,gp,hom,auts,isoms,pools,poolnum,a,x, 599isom,diso,conj,dgp,imgdepth,perms,kn,process,ii,goodgens,conjgens,genimgs, 600doesaut,biggens,wrimages,m,w,e,poolimggens,img,localgens,dfgens,wrs,dfimgs,b,perm,lims,map,aels,wrsoc,dfs; 601 602 603 # get a list of generator images and construct the appropriate element of 604 # a[1] inducing these generator images by conjugation 605 aels:=[]; 606 607 # dci is decomposition information needed for restriction to proper 608 # subgroups 609 dci:=rec(csi:=csi); 610 611 612 # isoms[i] An isomorphism from the first group in its pool to group i 613 614 goodgens:=[]; 615 poolimggens:=[]; 616 genimgs:=List([1..Maximum(Union(csi.minranges))],x->[]); 617 n:=csi.n; 618 csinice:=NiceGens(csi.recog); 619 act:=[]; 620 auts:=[]; 621 pools:=[]; 622 perms:=[]; 623 isoms:=[]; 624 n:=csi.n; 625 626 doesaut:=[]; 627 628 process:=[n]; 629 ii:=1; 630 while ii<=n do 631 # do we need to feed in another layer? 632 if Length(process)<ii then 633 Add(process,First([n,n-1..1],x->not x in process)); 634 fi; 635 i:=process[ii]; 636 637 if IsBound(csi.permap[i]) and not IsSolvableGroup(Image(csi.permap[i])) then 638 639 if not IsBound(goodgens[i]) then 640 # no generators set yet -- just take new ones 641 #lgens:=csinice{csi.minranges[i]}; 642 lgens:=CSINiceGens(csi,csi.minranges[i]); 643 goodgens[i]:=lgens; 644 else 645 lgens:=goodgens[i]; 646 fi; 647 648 gp:=Image(csi.permap[i]); 649 act[i]:=[]; 650 gens:=List(lgens,x->ImagesRepresentative(csi.permap[i],CSIImageHomNr(csi,i,x))); 651 652# for x in [1..Length(gens)] do 653# j:=CSIImageHomNr(csi,i,csinice[csi.minranges[i][x]]); 654# k:=Image(csi.permap[i],j); 655# if k<>gens[x] then 656# Error("GENS"); 657# fi; 658# od; 659 660 # nonsolvable factor -- Is it in pools? 661 poolnum:=First([1..Length(pools)],x->i in pools[x]); 662 # pools will always be joined TO the current one, so isom will not 663 # have to change on the way. 664 if poolnum=fail then 665 Add(pools,[i]); 666 poolnum:=Length(pools); 667 auts[poolnum]:=[]; 668 isoms[i]:=false; # representative 669 isom:=IdentityMapping(gp); 670 Add(poolimggens,gens); 671 elif i<>pools[poolnum][1] then 672 isom:=isoms[i]; # iso 673 else 674 isom:=IdentityMapping(gp); # first one -- iso is trivial 675 fi; 676 677 #find out what the previous factors do on it 678 for j in Filtered([1..i-1],x->IsBound(csi.minranges[x])) do 679 Info(InfoFFMat,2,"Try ",i," mapped by ",j); 680 681 for kn in csi.minranges[j] do 682 #k:=csinice[kn]; 683 k:=CSINiceGens(csi,kn); 684 if not IsBound(perms[kn]) then 685 perms[kn]:=[]; 686 genimgs[kn]:=[]; 687 doesaut[kn]:=[]; 688 fi; 689 690 conjgens:=[]; 691 genims:=[]; 692 imgdepth:=fail; 693 # calculate images 694 for l in lgens do 695 c:=l^k; # conjugate 696 Add(conjgens,c); 697 d:=CSIDepthElm(csi,c); 698 if imgdepth=fail then 699 # new image -- isomorphism 700 imgdepth:=d; 701 perms[kn][i]:=d; # record permutations 702 elif imgdepth<>d then 703 # this cannot happen 704 Error("incompatible depths"); 705 fi; 706 Add(genims,ImagesRepresentative(csi.permap[d],CSIImageHomNr(csi,d,c))); 707 od; 708 709 dgp:=Image(csi.permap[d]); 710 if AssertionLevel()>2 then 711 hom:=GroupHomomorphismByImages(gp,dgp,gens,genims); 712 else 713 hom:=GroupHomomorphismByImagesNC(gp,dgp,gens,genims); 714 fi; 715 if hom=fail then Error("should not happen");fi; 716 717 if d=i then 718 # pull generator images back in original group 719 genimgs[kn][i]:=List(genims,x->PreImagesRepresentative(isom,x)); 720 721#if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then 722# Error("imerrD"); 723#fi; 724 725 # component is not moved we just need to conjugate with isom 726 hom:=isom*hom*InverseGeneralMapping(isom); 727 SetIsBijective(hom,true); 728 if not IsInnerAutomorphism(hom) then 729 Add(auts[poolnum],hom); 730 AddSet(doesaut[kn],i); 731 elif not IsOne(hom) then 732 # not automorphism, but still acting 733 AddSet(doesaut[kn],i); 734 fi; 735 elif d in pools[poolnum] then 736 # we know already that the groups are isomorphic -- just store 737 # the new automorphism 738 739 # isomorphism is always to the first element in the pool 740 if d=pools[poolnum][1] then 741 # the group is the normal form -- we do not need to translate 742 diso:=IdentityMapping(dgp); 743 else 744 # isomorphism to canonical form 745 diso:=InverseGeneralMapping(isoms[d]); 746 fi; 747 genimgs[kn][i]:=List(genims,x->ImagesRepresentative(diso,x)); 748 749# if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then 750# Error("imerrC"); 751# fi; 752 hom:=isom*hom*diso; 753 SetIsBijective(hom,true); 754 if not IsInnerAutomorphism(hom) then 755 Add(auts[poolnum],hom); 756 AddSet(doesaut[kn],i); 757 fi; 758 else 759 # isomorphism wasn't known yet -- we need to join pools 760 # hom is the joiner, isom*hom*isoms[d]^-1 the conjugator of the 761 # canonical map 762 763 # the pool we add. We always join the image pool to the current 764 # one. 765 a:=First([1..Length(pools)],x->d in pools[x]); 766 if a=fail then 767 # we did not yet know layer d -- add it 768 Info(InfoFFMat,2,"Found Layer ",d); 769 Add(pools[poolnum],d); 770 isoms[d]:=isom*hom; 771 Assert(3,IsBijective(isoms[d])); 772 # newly included component -- the generators *are* simply the 773 # images 774 genimgs[kn][i]:=List(genims, 775 x->PreImagesRepresentative(isoms[d],x)); 776 777#if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then 778 #Error("imerrB"); 779#fi; 780 goodgens[d]:=conjgens; 781 Add(process,d); 782 else 783 Error("This cannot happen??"); 784 # we knew the layer and join pools 785 Info(InfoFFMat,2,"Join ",pools[a]," to ",pools[poolnum]); 786 conj:=isom*hom; 787 if not d=pools[a][1] then 788 # translate to normal form 789 conj:=conj*InverseGeneralMapping(isoms[d]); 790 fi; 791 792 # join pools 793 for x in pools[a] do 794 Add(pools[poolnum],x); 795 isoms[x]:=conj*isoms[x]; # reroot isomorphism 796 od; 797 798 # translate automorphisms 799 for x in auts[a] do 800 Add(auts[poolnum],conj*x*InverseGeneralMapping(conj)); 801 od; 802 803 # delete old 804 pools[a]:=[]; 805 auts[a]:=[]; 806 fi; 807 fi; 808 od; 809 od; 810 fi; 811 ii:=ii+1; 812 od; 813 814 wrs:=[]; 815 wrsoc:=[]; 816 dfgens:=[]; 817 dfimgs:=[]; 818 819#if ForAny(Flat(genimgs),x->not x in Image(csi.permap[pools[poolnum][1]])) then 820# Error("imerrA"); 821#fi; 822 823 if Length(pools)=0 then 824 # solvable 825 d:=Group(()); 826 a:=GeneratorsOfGroup(csi.group); 827 b:=List(a,x->One(d)); 828 RUN_IN_GGMBI:=true; # hack to skip Nice treatment 829 hom:=GroupHomomorphismByImagesNC(csi.group,d,a,b); 830 RUN_IN_GGMBI:=false; 831 dci:=rec(isTrivial:=true); 832 SetRecogDecompinfoHomomorphism(hom,dci); 833 SetImagesSource(hom,Group(())); 834 return hom; 835 fi; 836 837 dci.pools:=pools; 838 dci.wreathemb:=[]; 839 dci.goodgens:=goodgens; 840 dci.isoms:=isoms; 841 dci.poollocalgens:=[]; 842 843 # now build each group 844 for i in [1..Length(pools)] do 845 pool:=pools[i]; 846 m:=Length(pools[i]); 847 a:=Image(csi.permap[pools[i][1]]); 848 a:=[a,Concatenation(auts[i], 849 List(GeneratorsOfGroup(a),x->InnerAutomorphism(a,x)))]; 850 a:=AutomorphismRepresentingGroup(a[1],a[2]); 851 aels[i]:=a; 852 localgens:=List(poolimggens[i],x->Image(a[2],x)); 853 dci.poollocalgens[i]:=localgens; 854 b:=SymmetricGroup(m); 855 w:=WreathProduct(a[1],b); 856 e:=List([1..m+1],x->Embedding(w,x)); 857 dci.wreathemb[i]:=e; 858 biggens:=[]; 859 wrimages:=[]; 860 861 for j in Filtered([1..n],x->IsBound(csi.minranges[x])) do 862 c:=csi.minranges[j]; 863 Info(InfoFFMat,2,"Images for ",j," ",c); 864 865 # does any of the elements actually do something -- if so, what? 866 867# if true or ForAny(c,x->IsBound(doesaut[x]) 868# and ForAny(pools[i],y->y in doesaut[x])) then 869 Info(InfoFFMat,2,c," acts ",pools[i],j); 870 for kn in c do 871 #Add(biggens,csinice[kn]); 872 Add(biggens,CSINiceGens(csi,kn)); 873 img:=One(w); 874 for l in [1..Length(pools[i])] do 875 if IsBound(genimgs[kn][pools[i][l]]) then 876 d:=Image(e[l],CSIAelement(a,localgens,genimgs[kn][pools[i][l]])); 877 else 878 if pools[i][l]=j then 879 d:=List(goodgens[pools[i][l]], 880 x->ImagesRepresentative(csi.permap[pools[i][l]], 881 CSIImageHomNr(csi,pools[i][l],x^CSINiceGens(csi,kn)))); 882 if isoms[j]<>false then 883 d:=List(d,x->PreImagesRepresentative(isoms[j],x)); 884 fi; 885 d:=Image(e[l],CSIAelement(a,localgens,d)); 886 else 887 # component is not acted on as it lies higher 888 d:=One(a[1]); 889 fi; 890 fi; 891 img:=img*d; 892 od; 893 # is it a nontrivial permutation 894 if IsBound(perms[kn]) and ForAny([1..Length(perms[kn])], 895 x->IsBound(perms[kn][x]) and perms[kn][x]<>x) then 896 # fill up fixed points 897 for d in pool do 898 if not IsBound(perms[kn][d]) then perms[kn][d]:=d;fi; 899 od; 900 901 # permuting part 902 d:=PermList(List(perms[kn]{pool},x->Position(pool,x))); 903 img:=img*Image(e[Length(pool)+1],d); 904 fi; 905 906 Add(wrimages,img); 907 od; 908 909 od; 910 Add(wrs,w); 911 a:=List([1..m],x->PerfectResiduum(Image(Embedding(w,x)))); 912 b:=TrivialSubgroup(w); 913 for l in a do 914 b:=ClosureGroup(b,l); 915 od; 916 SetDirectFactorsFittingFreeSocle(w,a); 917 Add(wrsoc,b); # socle 918 Add(dfgens,biggens); 919 Add(dfimgs,wrimages); 920 od; 921 922 d:=DirectProduct(wrs); 923 dci.dirprod:=d; 924 dci.embeddings:=List([1..Length(pools)],x->Embedding(d,x)); 925 dci.aels:=aels; 926 a:=List([1..Length(pools)],x->List(DirectFactorsFittingFreeSocle(wrs[x]), 927 y->Image(dci.embeddings[x],y))); 928 dfs:=Concatenation(a); 929 SetDirectFactorsFittingFreeSocle(d,dfs); 930 a:=dfgens[1]; 931 b:=List(a,x->One(d)); 932 w:=TrivialSubgroup(d); 933 for i in [1..Length(pools)] do 934 e:=dci.embeddings[i]; 935 w:=ClosureGroup(w,Image(e,GeneratorsOfGroup(wrsoc[i]))); 936 for j in [1..Length(a)] do 937 Assert(1,dfgens[i][j]=a[j]); 938 939 b[j]:=b[j]*Image(e,dfimgs[i][j]); 940 od; 941 od; 942 943 if IsPermGroup(csi.group) and AssertionLevel()>2 then 944 hom:=GroupHomomorphismByImages(csi.group,d,a,b); 945 else 946 RUN_IN_GGMBI:=true; # hack to skip Nice treatment 947 hom:=GroupHomomorphismByImagesNC(csi.group,d,a,b); 948 RUN_IN_GGMBI:=false; 949 fi; 950 if hom=fail then 951 Error("fail!"); 952 fi; 953 b:=Subgroup(d,b); 954 #SetSocle(b,w); 955 dci.socle:=w; 956 dfs:=Filtered(dfs,x->IsSubset(b,x)); 957 SetDirectFactorsFittingFreeSocle(w,dfs); 958 SetDirectFactorsFittingFreeSocle(b,dfs); 959 SetImagesSource(hom,b); 960 SetRecogDecompinfoHomomorphism(hom,dci); 961 return hom; 962 963end; 964 965 966CSIElementAct:=function(dci,elm) 967local csi,pools,result,i,e,a,img,b,dp,d,kn,perm,l; 968 csi:=dci.csi; 969 pools:=dci.pools; 970 result:=One(dci.dirprod); 971 for i in [1..Length(pools)] do 972 e:=dci.wreathemb[i]; 973 a:=dci.aels[i]; 974 img:=One(Image(e[1])); 975 perm:=[]; 976 for l in [1..Length(pools[i])] do 977 b:=List(dci.goodgens[pools[i][l]],x->x^elm); 978 dp:=CSIDepthElm(csi,b[1]); 979 Assert(2,ForAll(b,x->CSIDepthElm(csi,x)=dp)); 980 perm[l]:=Position(pools[i],dp); 981 982 d:=List(dci.goodgens[pools[i][l]], 983 x->ImagesRepresentative(csi.permap[dp], 984 CSIImageHomNr(csi,dp,x^elm))); 985 if dci.isoms[dp]<>false then 986 d:=List(d,x->PreImagesRepresentative(dci.isoms[dp],x)); 987 fi; 988 d:=Image(e[l],CSIAelement(a,dci.poollocalgens[i],d)); 989 img:=img*d; 990 od; 991 992 # permuting part 993 d:=PermList(perm); 994 img:=img*Image(e[Length(pools[i])+1],d); 995 result:=result*Image(dci.embeddings[i],img); 996 997 od; 998 return result; 999end; 1000 1001############################################################################# 1002## 1003#M ImagesRepresentative(<hom>,<x>) 1004## 1005InstallMethod(ImagesRepresentative,"for recognition mappings", 1006 FamSourceEqFamElm, 1007 [ IsGroupGeneralMapping and 1008 HasRecogDecompinfoHomomorphism,IsMultiplicativeElementWithInverse], 0, 1009function(hom, x) 1010local d; 1011 d:=RecogDecompinfoHomomorphism(hom); 1012 if IsBound(d.isTrivial) then return ();fi; 1013 if IsBound(d.fct) then 1014 return d.fct(x); 1015 else 1016 return CSIElementAct(RecogDecompinfoHomomorphism(hom),x); 1017 fi; 1018end); 1019 1020############################################################################# 1021## 1022#M PreImagesSet(<hom>,<x>) 1023## 1024InstallMethod(PreImagesSet,"for recognition mappings", CollFamRangeEqFamElms, 1025 [ IsGroupGeneralMapping and 1026 HasRecogDecompinfoHomomorphism,IsGroup], 0, 1027function(hom, U) 1028local gens,pre; 1029 gens:=SmallGeneratingSet(U); 1030 pre:=List(gens,x->PreImagesRepresentative(hom,x)); 1031 U:=RecogDecompinfoHomomorphism(hom).LiftSetup; 1032 U:=SubgroupByFittingFreeData(Source(hom),pre,gens,U.pcgs); 1033 return U; 1034end); 1035 1036 1037#TODO: Detect Nonsolvable permuters (via perms) and leave out of pool 1038 1039BasePointsActionsOrbitLengthsStabilizerChain:=function(c) 1040local l,o; 1041 l:=[]; 1042 while c<>false do 1043 o:=c!.orb; 1044 Add(l,[o!.orbit[1],o!.op,Length(o!.orbit)]); 1045 c:=c!.stab; 1046 od; 1047 return l; 1048end; 1049 1050 1051FactorspaceActfun:=function(field,bas) 1052local heads; 1053 bas:=TriangulizedMat(bas); 1054 heads:=HeadsInfoOfSemiEchelonizedMat(bas,Length(bas[1])); 1055 return function(v,g) 1056 local c; 1057 v:=v*g; 1058 v:=SiftedVectorForGaussianRowSpace(field,bas, heads, v ); 1059 c:=PositionNonZero(v); 1060 if c<=Length(v) then 1061 v:=Inverse(v[c])*v; 1062 fi; 1063 MakeImmutable(v); 1064 return v; 1065 end; 1066end; 1067 1068INVTRANSPCACHE:=[]; 1069AsInverseTranspose:=function(x,g) 1070local a,p; 1071 a:=fail; 1072 p:=0; 1073 while a=fail and p<Length(INVTRANSPCACHE) do 1074 p:=p+1; 1075 if IsIdenticalObj(INVTRANSPCACHE[p][1],g) then 1076 a:=INVTRANSPCACHE[p][2]; 1077 if p>50 then 1078 p:=Concatenation([p],[1..p-1],[p+1..Length(INVTRANSPCACHE)]); 1079 INVTRANSPCACHE:=INVTRANSPCACHE{p}; 1080 fi; 1081 fi; 1082 od; 1083 if a=fail then 1084 a:=TransposedMat(Inverse(g)); 1085 if Length(INVTRANSPCACHE)>100 then 1086 # clean out old 1087 INVTRANSPCACHE:=Concatenation([[g,a]],INVTRANSPCACHE{[1..90]}); 1088 else 1089 INVTRANSPCACHE:=Concatenation([[g,a]],INVTRANSPCACHE); 1090 fi; 1091 fi; 1092 return OnRight(x,a); 1093end; 1094 1095 1096ModuleStructureBase:=function(mats) 1097local orbtranslimit,f,total,gens,mo,bas,basn,dims,a,p,vec,orb,t,dict,use,fct, 1098 new,alt,g,pnt,limit,whole,siftchain,sub,b,trymultipleorbits; 1099 1100 orbtranslimit:=function(vec,fct,limit) 1101 local dict,orb,t,pnt,g,img,p,coinc; 1102 dict:=NewDictionary(vec,true,f^Length(vec)); 1103 orb:=[vec]; 1104 AddDictionary(dict,vec,1); 1105 t:=[One(gens[1])]; 1106 pnt:=1; 1107 coinc:=true; 1108 while pnt<=Length(orb) do 1109 for g in gens do 1110 img:=fct(orb[pnt],g); 1111 p:=LookupDictionary(dict,img); 1112 if p=fail then 1113 Add(orb,img); 1114 if Length(orb)>limit then 1115 return fail; 1116 fi; 1117 AddDictionary(dict,img,Length(orb)); 1118 Add(t,t[pnt]*g); 1119 elif coinc then 1120 coinc:=false; 1121 fi; 1122 od; 1123 pnt:=pnt+1; 1124 od; 1125 return [dict,orb,t]; 1126 end; 1127 1128 trymultipleorbits:=function(seeds,fct,limit) 1129 local best,bestl,bestc,i,orb; 1130 best:=fail; 1131 bestl:=0; 1132 bestc:=0; 1133 for i in seeds do 1134 # form orbit with possible initial canonization 1135 orb:=orbtranslimit(fct(i,One(gens[1])),fct,limit); 1136 if orb<>fail and Length(orb[2])>1 then 1137 if Length(orb[2])>bestl then 1138 best:=orb; 1139 bestl:=Length(orb[2]); 1140 bestc:=1; 1141 if bestl>100 then 1142 return best; 1143 fi; 1144 else 1145 bestc:=bestc+1; 1146 # if we got 30 times the same length then use it. 1147 if bestc>30 then 1148 return best; 1149 fi; 1150 fi; 1151 elif bestl>0 then 1152 # we found at least one and others failed 1153 return best; 1154 fi; 1155 od; 1156 return best; 1157 end; 1158 1159 siftchain:=function(use,x) 1160 local i,img,u; 1161 i:=1; 1162 for u in use do 1163 img:=u.fct(u.vec,x); 1164 img:=LookupDictionary(u.dict,img); 1165 if img=fail then 1166 return [fail,i,x]; 1167 fi; 1168 x:=x/u.t[img]; 1169 i:=i+1; 1170 od; 1171 return x; 1172 end; 1173 1174 limit:=1000; 1175 1176 whole:=Group(mats); 1177 f:=FieldOfMatrixList(mats); 1178 total:=1; 1179 gens:=mats; 1180 mo:=GModuleByMats(gens,f); 1181 bas:=MTX.BasesCompositionSeries(mo); 1182 dims:=List(bas,Length); 1183 if ForAll([2..Length(bas)],x->IsSubset(bas[x],bas[x-1])) then 1184 basn:=List([2..Length(bas)],x->Filtered(bas[x],y->not y in bas[x-1])); 1185 bas:=Concatenation(basn); 1186 else 1187 a:=ShallowCopy(bas[2]); # 1 is empty 1188 p:=3; 1189 while p<=Length(bas) do 1190 t:=Difference(bas[p],a); 1191 t:=Difference(t,bas[p-1]); 1192 repeat 1193 vec:=List([1..Length(bas[p])-Length(bas[p-1])],x->Random(t)); 1194 until RankMat(Concatenation(a,vec))=Length(bas[p]); 1195 a:=Concatenation(a,vec); 1196 p:=p+1; 1197 od; 1198 bas:=a; 1199 fi; 1200 use:=[]; 1201 1202 while Length(gens)>0 do 1203 basn:=List(bas,x->OnLines(x,One(gens[1]))); 1204 p:=PositionProperty(basn,x->ForAny(gens,y->OnLines(x,y)<>x)); 1205 if p=fail then 1206 p:=PositionProperty(bas,x->ForAny(gens,y->x*y<>x)); 1207 vec:=bas{[p..Minimum(p+10,Length(bas))]}; 1208 fct:=OnRight; 1209 elif Size(f)^p<limit then 1210 p:=PositionProperty(dims,x->Size(f)^x>limit); 1211 if p=fail then 1212 p:=Length(dims); 1213 else 1214 p:=p-1; 1215 fi; 1216 orb:=OrbitsDomain(Group(gens),NormedRowVectors(VectorSpace(f,bas{[1..dims[p]]},Zero(bas[1]))),OnLines); 1217 orb:=Filtered(orb,x->Length(x)>1); 1218 Sort(orb,function(a,b) return Length(a)>Length(b);end); 1219 vec:=List(orb,x->x[1]); 1220 fct:=OnLines; 1221 else 1222 vec:=bas{[p..Minimum(p+10,Length(bas))]}; 1223 fct:=OnLines; 1224 fi; 1225 1226 # nontrivial orbit 1227 orb:=trymultipleorbits(vec,fct,limit); 1228 if orb=fail then 1229 b:=p; 1230 1231 a:=Reversed(Filtered(dims,x->x<b)); 1232 p:=fail; 1233 while p=fail and a[1]<>0 do 1234 sub:=a[1]; 1235 fct:=FactorspaceActfun(f,bas{[1..sub]}); 1236 basn:=List(bas,x->fct(x,One(gens[1]))); 1237 p:=Filtered([a[1]+1..Minimum(a[1]+10,Length(bas))], 1238 x->ForAny(gens,y->fct(basn[x],y)<>basn[x])); 1239 if p=[] then 1240 p:=Filtered([a[1]+10..Length(bas)], 1241 x->ForAny(gens,y->fct(basn[x],y)<>basn[x])); 1242 fi; 1243 if p=[] then 1244 p:=fail; 1245 else 1246 # try to find one with orbit not just length p but also not too 1247 # large 1248 p:=Reversed(p); 1249 while Length(p)>1 and Size(f)^(p[1]-a[1])/(Size(f)-1)>limit do 1250 p:=p{[2..Length(p)]}; 1251 od; 1252 p:=p[1]; 1253 fi; 1254 a:=a{[2..Length(a)]}; 1255 od; 1256 if p<>fail then 1257 orb:=trymultipleorbits(bas{[p..Minimum(p+10,Length(bas))]},fct,limit); 1258 if orb=fail then 1259 # try dual action 1260 fct:=AsInverseTranspose; 1261 p:=First([b..Length(bas)],x->ForAny(gens,y->fct(bas[x],y)<>bas[x])); 1262 if p<>fail then 1263 p:=[p]; 1264 for orb in [1..5] do 1265 Add(p,p[1]+orb); 1266 Add(p,p[1]-orb); 1267 od; 1268 p:=Filtered(p,x->x>0 and x<Length(bas)); 1269 orb:=trymultipleorbits(bas{p},fct,limit); 1270 fi; 1271 if orb=fail then 1272 Info(InfoFFMat,2,"even too long in dual"); 1273 return rec(points:=[],ops:=[]); 1274 else 1275 Info(InfoFFMat,2,"dual "); 1276 fi; 1277 else 1278 Info(InfoFFMat,2,"factor space "); 1279 fi; 1280 else 1281 return rec(points:=[],ops:=[]); 1282 Error("p=fail -- should not happen"); 1283 fi; 1284 1285 fi; 1286 1287 dict:=orb[1]; 1288 t:=orb[3]; 1289 orb:=orb[2]; 1290 vec:=orb[1]; 1291 Info(InfoFFMat,2,"got Length ",Length(orb),":", 1292 Collected(Factors(total*Length(orb)))); 1293 1294 Add(use,rec(vec:=vec,fct:=fct,orb:=orb,dict:=dict,t:=t,gens:=gens)); 1295 1296 #sift 20 random elements 1297 new:=[]; 1298 p:=Maximum(Minimum(1000,Length(orb)*Length(gens)),10); 1299 alt:=0; 1300 a:=fail; 1301 while a=fail and Length(new)<20 and alt<p do 1302 a:=PseudoRandom(whole); 1303 a:=siftchain(use,a); 1304 alt:=alt+1; 1305 if a[1]<>fail then 1306 if fct(vec,a)<>vec then 1307 Error("not stab"); 1308 fi; 1309 if not IsOne(a) and not a in new then 1310 Add(new,a); 1311 alt:=0; 1312 fi; 1313 a:=fail; # sifted through 1314 fi; 1315 od; 1316 1317 if a<>fail then 1318 gens:=ShallowCopy(use[a[2]].gens); 1319 Add(gens,a[3]); 1320 use:=use{[1..a[2]-1]}; 1321 total:=Product(List(use,x->Length(x.orb))); 1322 Info(InfoFFMat,2,"construction fail on level ",a[2]," of ",Length(use), 1323 " -- redo ",total); 1324 else 1325 total:=total*Length(orb); 1326 gens:=new; 1327 fi; 1328 1329 od; 1330 return rec(points:=List(use,x->x.vec),ops:=List(use,x->x.fct), 1331 order:=total); 1332end; 1333 1334MATGRP_AddGeneratorToStabilizerChain:="2bdefined"; 1335 1336MATGRP_StabilizerChainInner:= 1337 function( gens, size, layer, cand, opt, parentS ) 1338 # Computes a stabilizer chain for the group generated by gens 1339 # with known size size (can be false if not known). This will be 1340 # layer layer in the final stabilizer chain. cand is a (shared) 1341 # record for base point candidates and opt the (shared) option 1342 # record. This is called in StabilizerChain and calls itself. 1343 # It also can be called if a new layer is needed. 1344 local base,gen,S,i,merk,merk2,next,pr,r,stabgens,x; 1345 1346 Info(InfoGenSS,4,"Entering MATGRP_StabilizerChainInner layer=",layer); 1347 next:=rec(point:=cand.points[1],op:=cand.ops[1]); 1348 #next := GENSS_NextBasePoint(gens,cand,opt,parentS); 1349 #cand := next.cand; # This could have changed 1350 1351 S := GENSS_CreateStabChainRecord(parentS,gens,size, 1352 next.point,next.op,cand,opt); 1353 base := S!.base; 1354 1355 Info( InfoGenSS, 3, "Entering orbit enumeration layer ",layer,"..." ); 1356 repeat 1357 Enumerate(S!.orb,opt.OrbitLengthLimit); 1358 if not(IsClosed(S!.orb)) then 1359 if opt.FailInsteadOfError then 1360 return "Orbit too long, increase opt.OrbitLengthLimit"; 1361 else 1362 Error("Orbit too long, increase opt.OrbitLengthLimit"); 1363 fi; 1364 fi; 1365 until IsClosed(S!.orb); 1366 Info(InfoGenSS, 2, "Layer ", layer, ": Orbit length is ", Length(S!.orb)); 1367 1368 if layer > 1 then 1369 parentS!.stab := S; # such that from now on random element 1370 # generation works! 1371 else 1372 if (Length(S!.orb) > 50 or S!.orb!.depth > 5) and 1373 S!.opt.OrbitsWithLog then 1374 Error("ARGH!"); 1375 Info(InfoGenSS, 3, "Trying to make Schreier tree shallower (depth=", 1376 S!.orb!.depth,")..."); 1377 merk := Length(S!.orb!.gens); 1378 merk2 := Length(S!.stronggens); 1379 MakeSchreierTreeShallow(S!.orb); 1380 Append(S!.stronggens,S!.orb!.gens{[merk+1..Length(S!.orb!.gens)]}); 1381 Append(S!.layergens,[merk2+1..Length(S!.stronggens)]); 1382 Info(InfoGenSS, 3, "Depth is now ",S!.orb!.depth); 1383 fi; 1384 fi; 1385 S!.orb!.gensi := List(S!.orb!.gens,x->x^-1); 1386 1387 # as we add normalizing, we are done 1388 return S; 1389 1390 # Are we done? 1391 if size <> false and Length(S!.orb) = size then 1392 S!.proof := true; 1393 Info(InfoGenSS,4,"Leaving MATGRP_StabilizerChainInner layer=",layer); 1394 return S; 1395 fi; 1396 1397 # Now create a few random stabilizer elements: 1398 stabgens := EmptyPlist(opt.RandomStabGens); 1399 for i in [1..opt.RandomStabGens] do 1400 x := GENSS_RandomElementFromAbove(S,layer); 1401 if not(S!.IsOne(x)) then 1402 Add(stabgens,x); 1403 fi; 1404 od; 1405 Info(InfoGenSS,3,"Created ",opt.RandomStabGens, 1406 " random stab els, ", 1407 Length(stabgens)," non-trivial."); 1408 1409 if Length(stabgens) > 0 then # there is a non-trivial stabiliser 1410 Info(InfoGenSS,3,"Found ",Length(stabgens)," non-trivial ones."); 1411 if size <> false then 1412 S!.stab := MATGRP_StabilizerChainInner(stabgens,size/Length(S!.orb), 1413 layer+1,cand,opt,S); 1414 else 1415 S!.stab := MATGRP_StabilizerChainInner(stabgens,false, 1416 layer+1,cand,opt,S); 1417 fi; 1418 if IsString(S!.stab) then return S!.stab; fi; 1419 if opt.ImmediateVerificationElements > 0 then 1420 Info(InfoGenSS,2,"Doing immediate verification in layer ", 1421 S!.layer," (",opt.ImmediateVerificationElements, 1422 " elements)..."); 1423 i := 0; 1424 while i < opt.ImmediateVerificationElements do 1425 i := i + 1; 1426 x := GENSS_RandomElementFromAbove(S,layer); 1427 if MATGRP_AddGeneratorToStabilizerChain(S!.topS,x) then 1428 Info( InfoGenSS, 2, "Immediate verification found error ", 1429 "(layer ",S!.layer,")..." ); 1430 i := 0; 1431 fi; 1432 od; 1433 fi; 1434 1435 S!.proof := S!.stab!.proof; # hand up information 1436 else 1437 # We are not sure that the next stabiliser is trivial, but we believe! 1438 Info(InfoGenSS,3,"Found no non-trivial ones."); 1439 S!.proof := false; 1440 fi; 1441 1442 Info(InfoGenSS,4,"Leaving MATGRP_StabilizerChainInner layer=",layer); 1443 return S; 1444 end; 1445 1446 1447MATGRP_AddGeneratorToStabilizerChain:= 1448 function( S, el,pt,op ) 1449 # Increases the set represented by S by the generator el. 1450 local SS, r, n, pr, i, newstrongnr; 1451 if IsBound(S!.trivialgroup) and S!.trivialgroup then 1452 if S!.IsOne(el) then 1453 return false; 1454 fi; 1455 #SS := StabilizerChain(Group(el),S!.opt); 1456 SS:=StabilizerChain(Group(el),rec(Cand:=rec(ops:=[op],points:=[pt]),Reduced:=false,StrictlyUseCandidates:=true)); 1457 1458 if IsString(SS) then return SS; fi; 1459 for n in NamesOfComponents(SS) do 1460 S!.(n) := SS!.(n); 1461 od; 1462 Unbind(S!.trivialgroup); 1463 return true; 1464 fi; 1465 1466 r := SiftGroupElement( S, el ); 1467 # if this is one then el is already contained in the stabilizer chain 1468 if r.isone then # already in the group! 1469 return false; 1470 fi; 1471 # Now there remain two cases: 1472 # (1) the sift stopped somewhere and we have to add a generator there 1473 # (2) the sift ran all through the chain and the element still was not 1474 # the identity, then we have to prolong the chain 1475 if r.S <> false then # case (1) 1476 SS := r.S; 1477 Info( InfoGenSS, 2, "Adding new generator to stab. chain ", 1478 "in layer ", SS!.layer, " from ",Length(SS!.stronggens) ); 1479 Add(SS!.stronggens,r.rem); 1480 Add(SS!.layergens,Length(SS!.stronggens)); 1481 AddGeneratorsToOrbit(SS!.orb,[r.rem]); 1482 Add(SS!.orb!.gensi,r.rem^-1); 1483 newstrongnr := Length(SS!.stronggens); 1484 Info( InfoGenSS, 4, "Entering orbit enumeration layer ",SS!.layer, 1485 "..." ); 1486 repeat 1487 Enumerate(SS!.orb,S!.opt.OrbitLengthLimit); 1488 if not(IsClosed(SS!.orb)) then 1489 if S!.opt.FailInsteadOfError then 1490 return "Orbit too long, increase S!.opt.OrbitLengthLimit"; 1491 else 1492 Error("Orbit too long, increase S!.opt.OrbitLengthLimit!"); 1493 fi; 1494 fi; 1495 until IsClosed(SS!.orb); 1496 Info( InfoGenSS, 4, "Done orbit enumeration layer ",SS!.layer,"..." ); 1497 SS!.proof := false; 1498 else # case (2) 1499 # Note that we do not create a pr instance here for one 1500 # generator, this will be done later on as needed... 1501 SS := r.preS; 1502 newstrongnr := Length(SS!.stronggens)+1; # r.rem will end up there ! 1503 SS!.stab := MATGRP_StabilizerChainInner([r.rem],false, 1504 SS!.layer+1,rec(points:=[pt],ops:=[op]), SS!.opt, SS ); 1505 if IsString(SS!.stab) then return SS!.stab; fi; 1506 SS := SS!.stab; 1507 fi; 1508 # Now we have added a new generator (or a new layer) at layer SS, 1509 # the new gen came from layer S (we were called here, after all), 1510 # thus we have to check, whether all the orbits between S (inclusively) 1511 # and SS (exclusively) are also closed under the new generator r.rem, 1512 # we add it to all these orbits, thereby also making the Schreier trees 1513 # shallower: 1514 while S!.layer < SS!.layer do 1515 Info(InfoGenSS,2,"Adding new generator to orbit in layer ",S!.layer); 1516 Add(S!.layergens,newstrongnr); 1517 AddGeneratorsToOrbit(S!.orb,[r.rem]); 1518 Add(S!.orb!.gensi,r.rem^-1); 1519 S := S!.stab; 1520 od; 1521 # Finally, we have to add it to the product replacer! 1522 AddGeneratorToProductReplacer(S!.opt!.pr,r.rem); 1523 return true; 1524 end; 1525 1526 1527SolvableBSGS:=function(arg) 1528local CBase, normalizingGenerator,df,ops,firstmoved,i, 1529solvNC,S,pcgs,x,r,c,w,a,bound,U,xp,depths,oldsz,prime,relord,gens,acter,ogens,stabs,n,strongs,stronglevs,laynums,slvec,layerzero,p,laynum,layers,sel,vals,stronglayers,layervecs,slpval,slp,baspts,levp,blocksz,lstrongs,lstrongsinv,bl,opt,check,primes,shortlim,orb,orbs,j,sz,goodbase,CHAINTEST,orblens; 1530 1531 goodbase:=[]; 1532 CHAINTEST:=function(X,str) 1533 while X<>false do 1534 #if IsBound(X!.stronggens) and 1535 # Length(X!.layergens)>Length(Factors(Size(X))) then 1536 # Error("UGH"); 1537 #fi; 1538 if not X!.orb!.orbit[1] in goodbase.points then 1539 Error("new point!"); 1540 fi; 1541 #if Length(X!.orb!.gens)<>Length(Set(X!.orb!.gens)) then 1542 # Error("duplicate!"); 1543 #fi; 1544 #if IsBound(X!.opt) and IsBound(X!.opt.StrictlyUseCandidates) and 1545 # X!.opt.StrictlyUseCandidates=false then Error("eh5!"); fi; 1546 1547 if IsBound(X!.orb!.gens) and Length(X!.orb!.gens)<>Length(Factors(Size(X))) then 1548 Error("length!"); 1549 fi; 1550 1551 if IsBound(X!.orb!.gens) and Length(Difference(X!.orb!.gens,strongs))>0 then 1552 Error("XTRA!"); 1553 fi; 1554 1555 if IsBound(X!.orb!.gensi) and List(X!.orb!.gens,Inverse)<>X!.orb!.gensi then 1556 Error(str,"inverse!"); 1557 fi; 1558 if IsBound(X!.orb!.gensi) and ForAny(X!.orb!.schreiergen,IsInt) and 1559 Maximum(Filtered(X!.orb!.schreiergen,IsInt))>Length(X!.orb!.gensi) 1560 then 1561 Error("length!"); 1562 fi; 1563 X:=X!.stab; 1564 od; 1565 end; 1566 1567 gens:=arg[Minimum(2,Length(arg))]; 1568 if IsGroup(gens) then 1569 a:=gens; 1570 gens:=GeneratorsOfGroup(gens); 1571 else 1572 a:=Group(gens); 1573 fi; 1574 if false and Length(arg)>2 and Length(gens)>5+Length(Factors(arg[3])) then 1575 gens:=List([1..5+Length(Factors(arg[3]))],x->PseudoRandom(a)); 1576 gens:=Set(gens); 1577 a:=Group(gens); 1578 check:=true; 1579 #SetSize(a,sz); 1580 else 1581 check:=false; 1582 fi; 1583 if IsBound(arg[3]) then 1584 sz:=arg[3]; 1585 #SetSize(a,sz); 1586 else 1587 sz:=fail; 1588 fi; 1589 SetIsSolvableGroup(a,true); 1590 1591 if Length(arg)=1 then 1592 acter:=gens; 1593 else 1594 acter:=arg[1]; 1595 if IsGroup(acter) then 1596 acter:=GeneratorsOfGroup(acter); 1597 fi; 1598 fi; 1599 1600 shortlim:=Size(DefaultFieldOfMatrixGroup(a))*3000; 1601 1602 df:=DefaultFieldOfMatrixGroup(a); 1603 if false then 1604 # try vectors from big group 1605 w:=BasisVectorsForMatrixAction(Group(acter)); 1606 w:=ImmutableMatrix(df,w); 1607 if Size(df)=2 then 1608 ops:=List(w,x->OnRight); 1609 else 1610 ops:=List(w,x->OnLines); 1611 fi; 1612 fi; 1613 w:=ModuleStructureBase(gens); 1614 if IsBound(arg[3]) and IsBound(w.order) and w.order<arg[3] then 1615 for i in gens[1] do 1616 Add(w.points,i); 1617 Add(w.ops,OnRight); 1618 od; 1619 fi; 1620 1621 opt:=rec(RandomStabGens:=10, 1622 Cand:=rec(points:=w.points,ops:=w.ops), 1623 #TODO XXX Find better strategy for limit iof field is larger 1624 VeryShortOrbLimit:=shortlim 1625 ); 1626 FUNCSPACEHASH:=[]; # needed in base point selection code 1627 CBase:=StabilizerChain(a,opt); 1628 if sz<>fail and Size(CBase)<sz then 1629 Info(InfoWarning,1,"Wrong size -- redo with `doall' option"); 1630 return fail; 1631 fi; 1632 1633 primes:=Set(Factors(Size(CBase))); 1634 FUNCSPACEHASH:=[]; 1635 w:=BasePointsActionsOrbitLengthsStabilizerChain(CBase); 1636 Info(InfoGenSS,2,"Suggested Base Points",List(w,x->Position(opt.Cand.points,x[1])), 1637 " Lengths ",List(w,x->x[3])); 1638 1639 goodbase:=BaseStabilizerChain(CBase); 1640 w:=1; 1641 opt:=1; 1642 1643 # Dixon's (1968 - the solvable length of a solvable linear group) bounds 1644 if IsPerm(gens[1]) then 1645 a:=Length(MovedPoints(gens)); 1646 bound := Int( LogInt( Maximum(1,a ^ 5), 3 ) / 2 ); 1647 elif IsMatrix(gens[1]) then 1648 a:=Length(gens[1]); # dimension 1649 1650 # we would need proper log 1651 #bound := Int( (8.55*Log(a+1,10)+0.36 ); 1652 # since it does not exist, replace 8.55 by 9 and simplify rounded up. 1653 # this anyhow only matters for catching wrong input, so its harmless 1654 bound := Log((a+1)^9,10)+1; 1655 else 1656 # no bound known 1657 bound:=infinity; 1658 fi; 1659 1660 normalizingGenerator:=function(g) 1661 local oldstrong,oldsz,phq,pbp,pba; 1662 oldsz:=Size(S); 1663 oldstrong:=ShallowCopy(StrongGenerators(S)); 1664 # work around the ommission of prior redundant base points 1665 pba:=Filtered([1..Length(goodbase.points)],x->goodbase.ops[x](goodbase.points[x],g)<>goodbase.points[x]); 1666 pbp:=pba[1]; 1667 1668 #Print("\n"); 1669 #View(S); 1670 #Print("\n use ",pba,"\n"); 1671 1672 MATGRP_AddGeneratorToStabilizerChain(S,g,goodbase.points[pbp],goodbase.ops[pbp]); 1673 1674 while Size(S)=oldsz do 1675 Error("orderr BUG !!!\n"); 1676 MATGRP_AddGeneratorToStabilizerChain(S,g); 1677 od; 1678 return Difference(StrongGenerators(S),oldstrong); 1679 end; 1680 1681 solvNC:=function() 1682 local N,process,layergens,U,g,h,r,V,x,comm,pow,wp,sift; 1683 N:=StabilizerChain(Group(StrongGenerators(S)),rec(Base:=CBase,Size:=Size(S),Reduced:=false,StrictlyUseCandidates:=true)); # really should be copy 1684 1685 # determine relative order 1686 pow:=2; 1687 wp:=w*w; 1688 while not SiftGroupElement(N,wp).isone do 1689 wp:=wp*w; 1690 pow:=pow+1; 1691 od; 1692 if not IsPrimeInt(pow) then 1693 prime:=Factors(pow)[1]; 1694 w:=w^(pow/prime); 1695 else 1696 prime:=pow; 1697 fi; 1698 Info(InfoFFMat,2,"Relative order ",pow," prime ",prime); 1699 1700 process:=[w]; 1701 layergens:=[]; 1702 U:=[]; 1703 for g in process do 1704 sift:=SiftGroupElement(S,g); 1705 if not sift.isone then 1706 Info(InfoFFMat,2,"SS=",Size(S)," ",Length(U)); 1707 1708 for h in layergens do 1709 comm:=Comm(g,h); 1710 #Print(Order(comm),"\n"); 1711 if not SiftGroupElement(N,comm).isone then 1712 # wrong element -- get new generator for layer below 1713 a:=false; 1714 w:=comm; 1715 S:=N; # don't we want to forget what we added? 1716 return false; 1717 fi; 1718 od; 1719 1720 # the sifted element will be as good as generator, but allows for 1721 # cleaner decomposition. 1722 g:=sift.rem; 1723 1724 V:=normalizingGenerator(g); 1725 #g:=V[1]; 1726 if g in U then Error("duplicate");fi; 1727 U:=Concatenation([g],U); 1728 Add(layergens,g); 1729 for x in acter do 1730 Add(process,g^x); 1731 od; 1732 1733 fi; 1734 od; 1735 a:=true; 1736 return U; 1737 1738 end; 1739 1740 S:=StabilizerChain(Group(One(gens[1])),rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true)); 1741 1742 pcgs:=[]; 1743 depths:=[]; 1744 relord:=[]; 1745 xp:=1; 1746 while xp<=Length(gens) do 1747 Info(InfoFFMat,2,"Processing ",xp); 1748 x:=gens[xp]; 1749 # feed through derived series if necessary 1750 while not SiftGroupElement(S,x).isone do 1751 c:=0; 1752 w:=x; 1753 a:=false; 1754 oldsz:=Size(S); 1755 while not a do 1756 c:=c+1; 1757 if c>bound then 1758 # not solvable 1759 Error("not solvable"); 1760 fi; 1761 1762 U:=solvNC(); # will change a,w 1763 1764 od; 1765Info(InfoFFMat,2,"Found Layer ",Size(S)/oldsz); 1766#View(S); 1767Info(InfoFFMat,2,"n"); 1768 if prime^Length(U)<>Size(S)/oldsz then Error("layerlen"); fi; 1769 pcgs:=Concatenation(U,pcgs); 1770 Add(depths,Length(pcgs)); 1771 Append(relord,ListWithIdenticalEntries(Length(U),prime)); 1772 1773 od; 1774 xp:=xp+1; 1775 od; 1776 n:=Length(pcgs); 1777 depths:=Reversed(List(depths,x->n+1-x)); 1778 relord:=Reversed(relord); 1779 Add(depths,n+1); 1780 1781 w:=BasePointsActionsOrbitLengthsStabilizerChain(S); 1782 Info(InfoGenSS,2,"USED Base Points",List(w,x->Position(goodbase.points,x[1])), 1783 " Lengths ",List(w,x->x[3])); 1784 1785 1786 # now build the chain once more with memory so we can decompose 1787 1788 blocksz:=[]; 1789 1790 Info(InfoFFMat,2,"NOW",Size(S)," ",Length(pcgs)); 1791 if Length(Factors(Size(S)))<>Length(pcgs) then Error("redundancies");fi; 1792 1793 # now sift the generators through so that the strong generators *ARE* a pcgs 1794 gens:=Reversed(pcgs); # start with low level gens again 1795 pcgs:=[]; 1796 acter:=[]; 1797 #S:=StabilizerChain(Group(One(gens[1])),rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true)); 1798 #S := GENSS_CreateStabChainRecord(false, 1799# [],1, 1800# goodbase.points[1], 1801# goodbase.ops[1],goodbase, 1802# rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true,InitialHashSize:=10)); 1803 S:=false; # work around special treatment for trivial group. 1804 oldsz:=1; 1805 stabs:=[]; 1806 xp:=1; 1807 strongs:=[]; 1808 stronglevs:=[]; 1809 while xp<=Length(gens) do 1810 1811 1812 if S<>false then 1813 #CHAINTEST(S,"E"); 1814 oldsz:=Size(S); 1815 Info(InfoFFMat,2,"ProcessiNg ",xp," ",Size(S)); 1816 1817 1818#if Length(StrongGenerators(S))>Length(Factors(Size(S))) then Error("more2\n");fi; 1819 x:=gens[xp]; 1820 x:=SiftGroupElement(S,x).rem; 1821 #SiftGroupElement(S,x); 1822 1823 if not IsOne(x) then 1824 normalizingGenerator(x); 1825 fi; 1826 else 1827 x:=gens[1]; 1828 # OrbitsWithLog will add powers as strong generators to make the tree shallower. 1829 # This messes with the correspondence of strong generators and pcgs elements 1830 1831 S:=StabilizerChain(Group(x),rec(Base:=CBase,Reduced:=false,StrictlyUseCandidates:=true, 1832 OrbitsWithLog := false,RandomStabGens:=1,Size:=Order(x))); 1833 1834 # Remove those ~!@#$ extra generators (powers of x) that are added 1835 # to make the tree shallower, but cause problems here. 1836 strongs:=[x]; 1837 1838 a:=S; 1839 1840 while a<>false do 1841 a!.stronggens:=ShallowCopy(strongs); 1842 a!.layergens:=[1]; 1843 if Length(a!.orb!.gens)>0 then 1844 if a!.orb!.gens<>strongs then 1845 a!.orb:=Orb(ShallowCopy(strongs),a!.orb[1],a!.orb!.op,rec(schreier:=true)); 1846 a!.orb!.gensi:=List(a!.orb!.gens,Inverse); 1847 repeat 1848 Enumerate(a!.orb); 1849 until IsClosed(a!.orb); 1850 a!.orb!.depth:=Size(a!.orb)-1; 1851 #if Length(a!.orb)=1 then 1852 # a!.orb!.schreiergen:=ShallowCopy(strongs); 1853 #else 1854 # a!.orb!.schreiergen:=[]; 1855 #fi; 1856 fi; 1857 fi; 1858 a:=a!.stab; 1859 od; 1860 1861 #CHAINTEST(S,"F2"); 1862 1863 strongs:=[]; 1864 1865 fi; 1866 1867 if Size(S)=oldsz then Error("no change!");fi; 1868 1869 Add(pcgs,x); 1870 a:=Set(Filtered(StrongGenerators(S),x->not x in strongs)); 1871 a:=Difference(a,List([0..Order(x)],y->x^y)); 1872 if Length(a)>0 then 1873 # redo 1874 Info(InfoFFMat,1,"nanu-redo"); 1875 #Error("nanu"); 1876 return SolvableBSGS(arg[1],arg[2],arg[3]); 1877 fi; 1878 Append(strongs,[x]); 1879 1880 #CHAINTEST(S,"F"); 1881 1882 1883 Append(stronglevs,ListWithIdenticalEntries(1,n+1-xp)); 1884 #if n+1-xp in depths then 1885 # Add(stabs,StructuralCopy(S)); 1886 #fi; 1887 xp:=xp+1; 1888 od; 1889 1890 pcgs:=Reversed(pcgs); 1891 1892 CBase:=BaseStabilizerChain(S); 1893 w:=BasePointsActionsOrbitLengthsStabilizerChain(S); 1894 Info(InfoGenSS,2,"Used Base Points",List(w,x->Position(goodbase.points,x[1])), 1895 " Lengths ",List(w,x->x[3])); 1896 1897 baspts:=[]; 1898 levp:=[]; 1899 xp:=[1..Length(pcgs)]; 1900 a:=S; 1901 for i in [1..Length(CBase.points)] do 1902 p:=List(xp,x->Position(a!.orb!.orbit,CBase.ops[i](CBase.points[i],pcgs[x]))); 1903 sel:=Filtered([1..Length(xp)],x->p[x]<>1); 1904 baspts{xp{sel}}:=ListWithIdenticalEntries(Length(sel),i); 1905 levp[i]:=xp{sel}; 1906 blocksz{xp{sel}}:=p{sel}-1; 1907 xp:=Difference(xp,xp{sel}); 1908 a:=a!.stab; 1909 od; 1910 1911 slpval:=function(w) 1912 local a,i; 1913 a:=w[2]*vals[w[1]]; 1914 for i in [3,5..Length(w)-1] do 1915 a:=(a+w[i+1]*vals[w[i]]) mod p; 1916 od; 1917 return a; 1918 end; 1919 1920 layerzero:=[]; 1921 laynum:=Length(depths)-1; 1922 layers:=List([1..n],x->First([laynum,laynum-1..1],y->x>=depths[y])); 1923 stronglayers:=layers{stronglevs}; 1924 # get vectors for strong generators (b/c we decompose into these) 1925 slvec:=[]; 1926 for i in [1..laynum] do 1927 p:=relord[depths[i]]; 1928 a:=IdentityMat(depths[i+1]-depths[i]); 1929 layerzero[i]:=Zero(a[1]); 1930 sel:=Positions(stronglayers,i); 1931 slvec{sel}:=a{List(strongs{sel},x->Position(pcgs,x))-depths[i]+1}; 1932 1933# layervecs:=ListWithIdenticalEntries(n,fail); # this will trap use of 1934# # higher gens 1935# layervecs{[depths[i]..depths[i+1]-1]}:=a; 1936# 1937# for j in [depths[i+1]..n] do 1938# layervecs[j]:=Zero(a[1]); 1939# od; 1940# slp:=LinesOfStraightLineProgram(SLPOfElms(strongs{sel})); 1941# vals:=Reversed(layervecs); 1942# for j in slp do 1943# if ForAll(j,IsList) then 1944# #last line 1945# slvec{sel}:=List(j,x->slpval(x)); 1946# elif IsList(j[1]) then 1947# Error("reassign syntax?"); 1948# else 1949# Add(vals,slpval(j)); 1950# fi; 1951# od; 1952 od; 1953 ForgetMemory(strongs); 1954 ForgetMemory(gens); 1955 ForgetMemory(S); 1956 1957 S!.pcgs:=pcgs; 1958 S!.layerzero:=layerzero; 1959 S!.layerprimes:=relord{depths{[1..laynum]}}; 1960 S!.strongvecs:=slvec; 1961 S!.layranges:=List([1..laynum],x->[depths[x]..depths[x+1]-1]); 1962 S!.revranges:=List([1..laynum],x->Reversed([1..Length(S!.layranges[x])])); 1963 1964 # reindexing 1965 a:=S; 1966 while a<>false do 1967 p:=List(a!.orb!.gensi,x->Position(strongs,x^-1)); 1968 if fail in p then 1969 Error("not found!"); 1970 fi; 1971 a!.gensistrongpos:=p; 1972 a!.gensilayers:=stronglayers{p}; 1973 a:=a!.stab; 1974 od; 1975 1976 #TODO: Use chain for decomposition -- at the moment it uses subgroups 1977 a:=pcgs; 1978 pcgs:=PcgsByPcSequenceCons(IsPcgsDefaultRep, 1979 IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs, 1980 FamilyObj(gens[1]),pcgs,[]); 1981 pcgs!.stabilizerChain:=S; 1982 #pcgs!.stabs:=Reversed(stabs); 1983 pcgs!.laynums:=[1..laynum]; 1984 pcgs!.layranges:=S!.layranges; 1985 pcgs!.baspts:=baspts; 1986 pcgs!.blocksz:=blocksz; 1987 pcgs!.levp:=levp; 1988 pcgs!.invpows:= 1989 List([1..Length(pcgs)],x->List([1..relord[x]-1],y->a[x]^(-y))); 1990 SetRelativeOrders(pcgs,relord); 1991 SetIndicesEANormalSteps(pcgs,depths); 1992 1993 return rec( 1994 pcgs:=pcgs, 1995 depths:=depths, 1996 relord:=relord, 1997 baspts:=baspts, 1998 blocksz:=blocksz); 1999end; 2000 2001# old ersion of exponent routines, allowing for arbitrary orbit arrangemnts 2002MatPcgsSift:=function( S, x,l ) 2003local o,p,po,preS,r,v,vecs,prime; 2004 v:=S!.layerzero[l]; 2005 prime:=S!.layerprimes[l]; 2006 vecs:=S!.strongvecs; 2007 preS := false; 2008 while S <> false do 2009 o := S!.orb; 2010 if IsObjWithMemory(x) then 2011 p := o!.op(o[1],x!.el); 2012 else 2013 p := o!.op(o[1],x); 2014 fi; 2015 po := Position(o,p); 2016 if po = fail then # not in current stabilizer 2017 Error("element not in group"); 2018 return rec( v:=v,rem := x, S := S ); 2019 fi; 2020 # Now sift through Schreier tree: 2021 while po > 1 do 2022 r:=o!.schreiergen[po]; 2023 x := x * S!.orb!.gensi[r]; 2024 po := o!.schreierpos[po]; 2025 if S!.gensilayers[r]=l then 2026 # generator on the desired layer -- add up vectors 2027 v:=(v+vecs[S!.gensistrongpos[r]]) mod prime; 2028 fi; 2029 od; 2030 S := S!.stab; 2031 od; 2032 return v; 2033end ; 2034 2035MatPcgsExponentsOld:=function(S,laynums,x) 2036local a,j,i,v; 2037 v:=[]; 2038 for i in [1..Length(laynums)] do 2039 #S:=stabs[laynums[i]]; 2040 a:=MatPcgsSift(S,x,laynums[i]); 2041 Append(v,a); 2042 if i<Length(laynums) then 2043 # need to divide off what is given by the vector 2044 for j in [1..Length(a)] do 2045 if a[j]<>0 then 2046 x:=LeftQuotient(S!.pcgs[S!.layranges[laynums[i]][j]]^a[j],x); 2047 fi; 2048 od; 2049 fi; 2050 od; 2051 return v; 2052end; 2053 2054 2055# new routine, using orbit positions 2056#decomposition as product, but not in canonical order, thus in multi stages 2057#TODO compare cost matrix multiplication vs. collection 2058MatPcgsExponents:=function(arg) 2059local pcgs,laynums,ox,o,p,po,preS,r,isone,ind,i,prd,S,q,rem,bs,pS,x,dep,e,layer,delta, 2060 z,prime,curran,seli,md,pos,sel,deponly; 2061 2062 pcgs:=arg[1]; 2063 laynums:=arg[2]; 2064 ox:=arg[3]; 2065 deponly:=Length(arg)>3 and arg[4]; 2066 2067 dep:=IndicesEANormalSteps(pcgs); 2068 md:=laynums[Length(laynums)]; #the layer depth which we ignore. 2069 2070 z:=ListWithIdenticalEntries(dep[Length(dep)]-1,0); 2071 e:=ShallowCopy(z); 2072 pS:=pcgs!.stabilizerChain; 2073 2074 repeat 2075 # factor the element using orbit positions -- this is correct only for 2076 # the topmost layer used 2077 x:=ox; 2078 S:=pS; 2079 prd:=[]; 2080 pos:=[]; 2081 preS := false; 2082 ind:=1; 2083 layer:=md; # Maximal layer we do 2084 delta:=ShallowCopy(z); # we need to forget all in the previous layer 2085 curran:=pcgs!.layranges[layer]; 2086 prime:=RelativeOrders(pcgs)[curran[1]]; 2087 while S <> false do 2088 sel:=pcgs!.levp[ind]; 2089 bs:=pcgs!.blocksz{sel}; 2090 2091 o := S!.orb; 2092 if IsObjWithMemory(x) then 2093 p := o!.op(o[1],x!.el); 2094 else 2095 p := o!.op(o[1],x); 2096 fi; 2097 po := Position(o,p)-1; 2098 2099 # decompose 2100 for i in [1..Length(sel)] do 2101 seli:=sel[i]; 2102 2103 q:=QuoInt(po,bs[i]); 2104 rem:=po-q*bs[i]; 2105 if q>0 then 2106 if seli<dep[layer] then 2107 # decrease layer in which we decompose 2108 layer:=layer-1; 2109 while seli<dep[layer] do 2110 layer:=layer-1; 2111 od; 2112 delta:=ShallowCopy(z); # we need to forget all in the previous layer 2113 prime:=RelativeOrders(pcgs)[seli]; 2114 curran:=pcgs!.layranges[layer]; 2115 fi; 2116 #Add(prd,[sel[i],q]);; 2117 2118 # remember exponent entry if right layer 2119 if seli<dep[layer+1] then 2120 delta[seli]:=delta[seli]+q mod prime; 2121 fi; 2122 2123 #x:=x/(pcgs[sel[i]]^q); 2124 x:=x*pcgs!.invpows[seli][q]; # use stored inverse powers 2125 2126 fi; 2127 po:=rem; 2128 od; 2129 2130 #if o[1]<>o!.op(o[1],x) then Error("not all off!");fi; 2131 2132 S := S!.stab; 2133 ind:=ind+1; 2134 od; 2135 2136 if delta<>fail then 2137 # now the affected (topmost) layer is correct 2138 e:=e+delta; 2139 2140 # we don't need to correct the last layer we use 2141 if deponly and not IsZero(delta) then 2142 delta:=fail; # pop out. 2143 elif layer<laynums[Length(laynums)] then 2144 # we have d describe a the highest level. So we want x=d*newx 2145 for i in curran do 2146 if delta[i]>0 then 2147 ox:=pcgs!.invpows[i][delta[i]]*ox; 2148 fi; 2149 od; 2150 fi; 2151 2152 fi; 2153 2154 # stop when we've set the lowest layer or if the remainder is the identity 2155 until layer=laynums[Length(laynums)] or delta=fail; 2156 2157 # do we only want some? 2158 if laynums=pcgs!.laynums then 2159 return e; 2160 else 2161 return e{Concatenation(pcgs!.layranges{laynums})}; 2162 fi; 2163 2164end; 2165 2166InstallMethod(ExponentsOfPcElement,"matrix pcgs",IsCollsElms, 2167 [IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,IsMatrix], 2168function(pcgs,x) 2169 return MatPcgsExponents(pcgs,pcgs!.laynums,x); 2170end); 2171 2172InstallOtherMethod(ExponentsOfPcElement,"matrix pcgs,indices",IsCollsElmsX, 2173 [IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,IsMatrix,IsList], 2174function(pcgs,x,inds) 2175local i,j,sel,lay,ip,sl,r,use; 2176 # is it a layer? 2177 for i in pcgs!.laynums do 2178 if inds=pcgs!.layranges[i] then 2179 # only this layer 2180 return MatPcgsExponents(pcgs,[i],x); 2181 fi; 2182 od; 2183 # which layers do we need? 2184 if not IsSortedList(inds) then 2185 # perverse case -- old 2186 return MatPcgsExponents(pcgs,pcgs!.laynums,x){inds}; 2187 fi; 2188 sel:=[]; 2189 lay:=[]; 2190 ip:=1; 2191 sl:=0; 2192 for i in pcgs!.laynums do 2193 r:=pcgs!.layranges[i]; 2194 use:=false; 2195 for j in [1..Length(r)] do 2196 if ip<=Length(inds) and r[j]=inds[ip] then 2197 Add(sel,j+sl); 2198 use:=true; 2199 ip:=ip+1; 2200 fi; 2201 od; 2202 if use then 2203 AddSet(lay,i); 2204 sl:=sl+Length(r); 2205 fi; 2206 od; 2207 return MatPcgsExponents(pcgs,lay,x){sel}; 2208end); 2209 2210InstallMethod(DepthOfPcElement,"matrix pcgs",IsCollsElms, 2211 [IsPcgsMatGroupByStabChainRep and IsPcgs and IsPrimeOrdersPcgs,IsMatrix], 2212function(pcgs,x) 2213local e; 2214 e:=MatPcgsExponents(pcgs,pcgs!.laynums,x,true); 2215 return PositionNonZero(e); 2216end); 2217 2218BindGlobal("TFDepthLeadExp",function(pcgs,x) 2219local p; 2220 x:=MatPcgsExponents(pcgs,pcgs!.laynums,x,true); # first nonzero 2221 p:=PositionNonZero(x); 2222 if p>Length(x) then 2223 # identity 2224 return [p,0]; 2225 else 2226 return [p,x[p]]; 2227 fi; 2228end); 2229 2230InstallMethod(FittingFreeLiftSetup,"fields, using recognition",true, 2231 [IsMatrixGroup],0, 2232function(G) 2233local csi,r,factorhom,sbs,k,pc,hom,rad,it,i,sz,x,stop; 2234 r:=DefaultFieldOfMatrixGroup(G); 2235 if not IsField(r) then 2236 TryNextMethod(); 2237 fi; 2238 2239 r:=RecognizeGroup(G); 2240 SetSize(G,Size(r)); 2241 csi:=GetInformationFromRecog(r); 2242 G!.storedrecog:=csi; 2243 factorhom:=FindAct(csi); 2244 2245 2246 #TODO: Better kernel gens by random selection 2247 sz:=Size(G)/Size(Image(factorhom)); 2248 if Size(Image(factorhom))=1 then 2249 # solvable 2250 k:=GeneratorsOfGroup(G); 2251 else 2252 it:=CoKernelGensIterator(InverseGeneralMapping(factorhom)); 2253 k:=Filtered(List([1..csi.genum],x->CSINiceGens(csi,x)),x->IsOne(ImagesRepresentative(factorhom,x))); 2254 if sz>1 then 2255 stop:=true; 2256 repeat 2257 for i in [1..3*Length(Factors(sz))] do 2258 if not IsDoneIterator(it) then 2259 x:=NextIterator(it); 2260 if not IsOne(x) and not x in k then 2261 Add(k,x); 2262 fi; 2263 fi; 2264 od; 2265 if sz>1 and Length(k)<Length(Factors(sz)) then stop:=false; fi; # work around issue 2266 if ValueOption("doall")=true then stop:=false;fi; 2267 until stop or IsDoneIterator(it); 2268 fi; 2269 fi; 2270 2271 Info(InfoGenSS,3,"|k|=",Length(k)); 2272 # TODO: Proper test 2273 2274 if ForAll(k,IsOne) then 2275 rad:=TrivialSubgroup(G); 2276 sbs:=rec(pcgs:=[],depths:=[1],relord:=[],pcisom:=IsomorphismPcGroup(rad),radical:=rad); 2277 else 2278 sbs:=SolvableBSGS(G,k,sz); 2279 while sbs=fail do 2280 sbs:=SolvableBSGS(G,k,sz:doall); 2281 od; 2282 rad:=SubgroupNC(G,sbs.pcgs); 2283 SetSize(rad,Product(RelativeOrders(sbs.pcgs))); 2284 sbs.radical:=rad; 2285 pc:=PcGroupWithPcgs(sbs.pcgs); 2286 RUN_IN_GGMBI:=true; # hack to skip Nice treatment 2287 hom:=GroupHomomorphismByImagesNC(rad,pc,sbs.pcgs,FamilyPcgs(pc)); 2288 RUN_IN_GGMBI:=false; 2289 SetIsBijective(hom,true); 2290 sbs.pcisom:=hom; 2291 fi; 2292 sbs.csi:=csi; 2293 SetKernelOfMultiplicativeGeneralMapping(factorhom,rad); 2294 sbs.factorhom:=factorhom; 2295 RecogDecompinfoHomomorphism(factorhom)!.LiftSetup:=sbs; 2296 2297 return sbs; 2298end); 2299 2300FFStats:=function(g) 2301local start,f; 2302 start:=Runtime(); 2303 f:=FittingFreeLiftSetup(g); 2304 Print("Time:",Runtime()-start,"\n"); 2305 Print("Factordegree ",NrMovedPoints(Range(f.factorhom)),"\n"); 2306 Print("PcgsDegrees ",Maximum(List( 2307 BasePointsActionsOrbitLengthsStabilizerChain(f.pcgs!.stabilizerChain), 2308 x->x[3])),"\n"); 2309 2310end; 2311 2312# work over residue class rings 2313 2314MyZmodnZObj:=function(a,b) 2315 if IsPrimeInt(b) and b<65536 then 2316 return Z(b)^0*a; 2317 else 2318 return ZmodnZObj(a,b); 2319 fi; 2320end; 2321 2322ReduceModM:=function(a,m) 2323 local b,r,i,j; 2324 if IsList(a) then 2325 if IsList(a[1]) then 2326 # matrix 2327 b:=[]; 2328 for i in a do 2329 r:=[]; 2330 for j in i do 2331 if IsZmodnZObjNonprime(j) then 2332 Add(r,MyZmodnZObj(j![1],m)); 2333 else 2334 Add(r,MyZmodnZObj(j,m)); 2335 fi; 2336 od; 2337 Add(b,r); 2338 od; 2339 if IsPrimeInt(m) then 2340 b:=ImmutableMatrix(m,b); 2341 fi; 2342 return b; 2343 else 2344 # vector 2345 b:=[]; 2346 for j in a do 2347 if IsZmodnZObjNonprime(j) then 2348 Add(b,MyZmodnZObj(j![1],m)); 2349 else 2350 Add(b,MyZmodnZObj(j,m)); 2351 fi; 2352 od; 2353 return b; 2354 fi; 2355 elif IsZmodnZObjNonprime(a) then 2356 a:=a![1]; 2357 fi; 2358 return MyZmodnZObj(a,m); 2359end; 2360 2361ReduceModMFunc:=function(m) 2362 return a->ReduceModM(a,m); 2363end; 2364 2365UnreduceModM:=Error; 2366 2367InstallMethod(FittingFreeLiftSetup,"residue class rings",true, 2368 [IsMatrixGroup],0, 2369function(g) 2370local r,m,f,a,ao,p,i,homs,hom,img,ff,ffp,ffpi,ffppc,ffhoms,ffsubs,d,elmimg, 2371 upperpcgs,upperexp,it,e,moli,pli,j,idx,depths,pcgs,levs,relord,idmat, 2372 fac,idmats,bas,basrep,basrepi,s,triv,addPcElement,procrels,addCleanUpper, 2373 k,l,bl,stack,stacks,gens,gnew,layerlimit,fertig; 2374 2375 triv:=TrivialSubgroup(CyclicGroup(2)); 2376 r:=DefaultFieldOfMatrixGroup(g); 2377 if IsField(r) or 2378 not CategoryCollections(IsZmodnZObjNonprime)(r) then 2379 TryNextMethod(); 2380 fi; 2381 idmat:=IdentityMat(Length(One(g)),1); 2382 2383 # convert to compact type 2384 gens:=GeneratorsOfGroup(g); 2385 if not ForAll(gens,IsZmodnZMat) then 2386 f:=FamilyObj(One(r)); 2387 gens:=List(gens,x->MakeZmodnZMat(f,List(x,r->List(r,Int)))); 2388 gnew:=Group(gens); 2389 if HasSize(g) then SetSize(gnew,Size(g));fi; 2390 else 2391 gnew:=g; 2392 fi; 2393 2394 m:=Size(r); 2395 # the prime power factors occurring 2396 f:=List(Collected(Factors(m)),x->x[1]^x[2]); 2397 Sort(f); 2398 homs:=[]; 2399 ffp:=[]; 2400 ffppc:=[]; 2401 ffhoms:=[]; 2402 moli:=[1]; 2403 pli:=[]; 2404 idx:=[false]; 2405 fac:=[false]; 2406 idmats:=[]; 2407 for i in f do 2408 a:=Factors(i); 2409 p:=a[1]; 2410 if Length(a)>1 then 2411 Add(moli,moli[Length(moli)]*p^2); 2412 Add(fac,1/p); 2413 Add(pli,p); 2414 Add(idx,Length(pli)); 2415 for j in [3..Length(a)] do 2416 Add(moli,moli[Length(moli)]*p); 2417 Add(idx,Length(pli)); 2418 Add(fac,fac[Length(fac)]/p); 2419 od; 2420 fi; 2421 hom:=List(gens,x->ReduceModM(x,p)); 2422 if ForAll(hom,IsOne) then 2423 hom:=GroupHomomorphismByFunction(gnew,Group(()),x->()); 2424 Info(InfoFFMat,2,"alltrivial ",p); 2425 Add(ffp,[]); 2426 Add(homs,hom); 2427 Add(ffhoms,hom); 2428 hom:=GroupHomomorphismByFunction(Group(()),triv,x->One(triv)); 2429 Add(ffppc,hom); 2430 else 2431 img:=Group(hom); 2432 hom:=GroupHomomorphismByFunction(gnew,img,ReduceModMFunc(p),false, 2433 x->UnreduceModM(x,m)); 2434 SetImagesSource(hom,img); 2435 Add(homs,hom); 2436 ff:=FittingFreeLiftSetup(img); 2437 Add(ffp,ff.pcgs); 2438 Add(ffppc,ff.pcisom); 2439 Add(ffhoms,hom*ff.factorhom); 2440 fi; 2441 od; 2442 if Length(ffhoms)=0 then 2443 hom:=GroupHomomorphismByFunction(gnew,Group(()),x->()); 2444 else 2445 d:=List(ffhoms,Image); 2446 d:=DirectProduct(d); 2447 elmimg:=function(x) 2448 local p,i; 2449 p:=One(d); 2450 for i in [1..Length(ffhoms)] do 2451 p:=p*ImagesRepresentative(Embedding(d,i), 2452 ImagesRepresentative(ffhoms[i],x)); 2453 od; 2454 return p; 2455 end; 2456 a:=List(gens,elmimg); 2457 hom:=GroupHomomorphismByFunction(gnew,SubgroupNC(d,a),elmimg); 2458 fi; 2459 2460 # we can't rescue the existing pcgs for the factors as we will not know 2461 # preimages. Compute all anews. 2462 ffpi:=List([1..Length(ffp)],x->[]); 2463 #$ffpi[1]:=ffp[1]; # the first pcgs is guaranteed to be always fully there 2464 2465 ffsubs:=List([1..Length(ffpi)],x-> 2466 TrivialSubgroup(Image(ffppc[x]))); 2467 2468 upperpcgs:=List([1..Length(ffpi)],x->[]); 2469 2470 # Do we have an upper limit for the space on each layer? 2471 layerlimit:=ValueOption("layerlimit"); 2472 if layerlimit=fail then 2473 layerlimit:=Length(One(g))^2; # full matrix space 2474 fi; 2475 2476 it:=CoKernelGensIterator(InverseGeneralMapping(hom)); 2477 bas:=List(moli,x->[]); 2478 basrep:=List(moli,x->[]); 2479 basrepi:=List(moli,x->[]); 2480 2481 addCleanUpper:=function(i,a) 2482 local r,p,s,e; 2483 r:=ImagesRepresentative(homs[i],a); 2484 p:=ImagesRepresentative(ffppc[i],r); 2485 if not p in ffsubs[i] then 2486 # need to extend induced pcgs 2487 s:=CanonicalPcgsByGeneratorsWithImages(ffp[i], 2488 Concatenation(ffpi[i],[r]), 2489 Concatenation(upperpcgs[i],[a])); 2490 ffpi[i]:=s[1]; 2491 upperpcgs[i]:=s[2]; 2492 ffsubs[i]:=ClosureGroup(ffsubs[i],p); 2493 fi; 2494 # divide off the nontrivial part in the quotient 2495 if Length(ffpi[i])>0 then 2496 e:=ExponentsOfPcElement(ffpi[i],r); 2497 e:=LinearCombinationPcgs(upperpcgs[i],e); 2498 a:=LeftQuotient(e,a); 2499 fi; 2500 return a; 2501 end; 2502 2503 addPcElement:=function(a,start) 2504 local i,j,p,e,s,added,bot,tmp; 2505 added:=fail; 2506 for i in [start..Length(moli)] do 2507 bot:=i=Length(moli); # last step -- powers are multiples 2508 p:=pli[idx[i]]; 2509 e:=List(a,x->List(x,y->y![1] mod moli[i])); 2510 e:=e-idmat; 2511 e:=fac[i]*Concatenation(e); 2512 e:=ImmutableVector(GF(p),e*Z(p)^0); 2513 if not IsZero(e) then 2514 if bot and IsBound(bas[i]) and Length(bas[i])=layerlimit then 2515 # Bottom layer is full, nothing else needed to do 2516 a:=fail; 2517 elif IsBound(bas[i]) and Length(bas[i])>0 then 2518 s:=SolutionMat(bas[i],e); 2519 if s=fail then 2520 Add(bas[i],e); 2521 Add(basrep[i],a); 2522 Add(basrepi[i],a^-1); 2523 if added=fail then added:=i;fi; 2524 if not bot then a:=a^p;fi; # no need to do if on bottom 2525 else 2526 s:=List(s,Int); 2527 # avoid inverse by imediately dividing off 2528 for j in [Length(s),Length(s)-1..1] do 2529 if not IsZero(s[j]) then 2530 if not bot then 2531 a:=basrepi[i][j]^s[j]*a; 2532 else 2533 # else case is not needed, as we are on the bottom :-) 2534 a:=fail; 2535# # multiplication is addition of p-residues 2536# tmp:=(basrepi[i][j]![1]*s[j]-s[j]*One(basrepi[i][j]![1])+a![1]); 2537# tmp:=tmp mod Characteristic(a); 2538# tmp:=MakeZmodnZMat(ElementsFamily(ElementsFamily(FamilyObj(a))),tmp); 2539# a:=tmp; 2540 fi; 2541 fi; 2542 od; 2543 fi; 2544 else 2545 bas[i]:=[e]; 2546 basrep[i]:=[a]; 2547 basrepi[i]:=[a^-1]; 2548 if added=fail then added:=i;fi; 2549 if not bot then a:=a^p;fi; # no need to do if on bottom 2550 fi; 2551 fi; 2552 od; 2553 return added; 2554 end; 2555 2556 fertig:=false; 2557 stacks:=[]; 2558 repeat 2559 a:=NextIterator(it); 2560 if not IsOne(a) then 2561 ao:=a; 2562 2563 for i in [1..Length(ffpi)] do 2564 a:=addCleanUpper(i,a); 2565 od; 2566 2567 bl:=List([2..Length(moli)],x->Length(bas[x])); 2568 addPcElement(a,2); 2569 if ForAny([2..Length(moli)],x->Length(bas[x])>bl[x-1]) then 2570 AddSet(stacks,ao); 2571 fi; 2572 2573 fi; 2574 fertig:=Length(moli)>1 and ForAll([2..Length(moli)],x->Length(basrep[x])=layerlimit); 2575 until IsDoneIterator(it) or fertig; 2576 2577 if fertig then stack:=[];fi; 2578 Info(InfoFFMat,2,"layerdimsc:",List(bas,Length)); 2579 2580 stack:=ShallowCopy(stacks); # we'll add to the list 2581 2582 # G-conjugates 2583 for k in stack do 2584 for j in gens do 2585 a:=k^j; 2586 2587 for i in [1..Length(ffpi)] do 2588 a:=addCleanUpper(i,a); 2589 od; 2590 2591 if not a in stacks then 2592 2593 AddSet(stacks,a); 2594 2595 # old base length 2596 bl:=List([2..Length(moli)],x->Length(bas[x])); 2597 addPcElement(a,2); 2598 if ForAny([2..Length(moli)],x->Length(bas[x])>bl[x-1]) then 2599 if not a in stack then Add(stack,a); fi; 2600 fi; 2601 fi; 2602 2603 od; 2604 od; 2605 Info(InfoFFMat,2,"layerdimsb:",List(bas,Length)); 2606 2607 Unbind(stack); 2608 Unbind(stacks); 2609 2610 pcgs:=[]; 2611 relord:=[]; 2612 levs:=[]; 2613 p:=0; 2614 for i in [1..Length(ffp)] do 2615 if Length(upperpcgs[i])>0 then 2616 2617 r:=RelativeOrders(ffpi[i]); 2618 if not fertig then 2619 j:=1; 2620 while j<=Length(upperpcgs[i]) do 2621 a:=upperpcgs[i][j]^r[j]; 2622 for k in [i..Length(ffp)] do 2623 a:=addCleanUpper(k,a); 2624 od; 2625 addPcElement(a,2); 2626 for l in Concatenation(pcgs,upperpcgs[i]) do 2627 a:=upperpcgs[i][j]^l; 2628 for k in [i..Length(ffp)] do 2629 a:=addCleanUpper(k,a); 2630 od; 2631 addPcElement(a,2); 2632 od; 2633 j:=j+1; 2634 od; 2635 fi; 2636 2637 Append(pcgs,upperpcgs[i]); 2638 if not HasIndicesEANormalSteps(ffpi[i]) then 2639 s:=FamilyPcgs(PcGroupWithPcgs(ffpi[i])); 2640 if not IsPcgsElementaryAbelianSeries(s) then Error("not EA pcgs");fi; 2641 s:=IndicesEANormalSteps(s); 2642 else 2643 s:=IndicesEANormalSteps(ffpi[i]); 2644 fi; 2645 s:=s{[1..Length(s)-1]}+p; 2646 Append(levs,s); 2647 Append(relord,RelativeOrders(ffpi[i])); 2648 p:=Length(pcgs); 2649 fi; 2650 od; 2651 Add(levs,Length(pcgs)+1); 2652 2653 Info(InfoFFMat,2,"layerdimsa:",List(bas,Length)); 2654 2655 # conjugation relations in linear bits. Will rarely add new elements (so the 2656 # danger of running through multiple times is minimal). 2657 procrels:=function() 2658 local i,j,k,l,a,b; 2659 b:=true; 2660 for i in [2..Length(moli)] do 2661 # if j>=Length(moli)-1+@ the conjugate is guaranteed the same 2662 for j in [i..Length(moli)-i+1] do 2663 if IsBound(basrep[i]) and IsBound(basrep[j]) then 2664 for k in basrep[i] do 2665 for l in basrep[j] do 2666 a:=addPcElement(l^k,j); 2667 if a<>fail then b:=false;fi; 2668 od; 2669 od; 2670 fi; 2671 od; 2672 od; 2673 for j in [2..Length(moli)] do 2674 if IsBound(basrep[j]) then 2675 for k in pcgs do 2676 for l in basrep[j] do 2677 a:=addPcElement(l^k,j); 2678 if a<>fail then b:=false;fi; 2679 od; 2680 od; 2681 fi; 2682 od; 2683 return b; 2684 end; 2685 2686 repeat until fertig or procrels(); 2687 2688 Info(InfoFFMat,2,"layerdims:",List(bas,Length)); 2689 2690 for i in [2..Length(bas)] do 2691 if IsBound(bas[i]) then 2692 Append(pcgs,basrep[i]); 2693 Append(relord,ListWithIdenticalEntries(Length(basrep[i]),pli[idx[i]])); 2694 fi; 2695 Add(levs,Length(pcgs)+1); 2696 od; 2697 2698 pcgs:=PcgsByPcSequenceCons(IsPcgsDefaultRep, 2699 IsPcgsResidueMatGroupRep and IsPcgs and IsPrimeOrdersPcgs, 2700 FamilyObj(One(g)),pcgs,[]); 2701 2702 # decomposition info for pcgs 2703 r:=rec(pcgs:=pcgs,factorhom:=hom, 2704 homs:=homs,ffp:=ffp,ffpi:=ffpi,ffppc:=ffppc,idmat:=idmat, 2705 upperpcgs:=upperpcgs,moli:=moli,pli:=pli,idx:=idx,fac:=fac, 2706 depths:=levs,bas:=bas,basrep:=basrep); 2707 2708 pcgs!.decompInfo:=r; 2709 SetRelativeOrders(pcgs,relord); 2710 SetIndicesEANormalSteps(pcgs,levs); 2711 2712 s:=SubgroupNC(g,pcgs); 2713 SetSize(s,Product(RelativeOrders(pcgs))); 2714 SetSize(g,Size(Image(hom))*Size(s)); 2715 SetKernelOfMultiplicativeGeneralMapping(hom,s); 2716 2717 r:=rec(pcgs:=pcgs,radical:=s,factorhom:=hom,depths:=levs); 2718 2719 if ValueOption("pcisom")<>false then 2720 i:=List(pcgs,x->ExponentsOfPcElement(pcgs,x^-1)); 2721 p:=PcGroupWithPcgs(pcgs:inversehints:=i); 2722 RUN_IN_GGMBI:=true; # hack to skip Nice treatment 2723 p:=GroupHomomorphismByImagesNC(s,p,pcgs,FamilyPcgs(p)); 2724 RUN_IN_GGMBI:=false; 2725 SetIsBijective(p,true); 2726 r.pcisom:=p; 2727 fi; 2728 2729 SetRecogDecompinfoHomomorphism(hom,rec(fct:=elmimg,LiftSetup:=r)); 2730 return r; 2731 2732end); 2733 2734ExponentsResiduePcgs:=function(r,elm) 2735local e,s,i,p,exp; 2736 exp:=[]; 2737 # upper pcgs part 2738 for i in [1..Length(r.ffpi)] do 2739 if Length(r.ffpi[i])>0 then 2740 s:=ExponentsOfPcElement(r.ffpi[i],ImagesRepresentative(r.homs[i],elm)); 2741 Append(exp,s); 2742 elm:=LeftQuotient(LinearCombinationPcgs(r.upperpcgs[i],s),elm); 2743 fi; 2744 od; 2745 2746 for i in [2..Length(r.moli)] do 2747 p:=r.pli[r.idx[i]]; 2748 e:=List(elm,x->List(x,y->y![1] mod r.moli[i])); 2749 e:=e-r.idmat; 2750 e:=r.fac[i]*Concatenation(e); 2751 e:=e*Z(p)^0; 2752 if not IsZero(e) then 2753 s:=SolutionMat(r.bas[i],e); 2754 if s=fail then 2755 Error("not in span of pcgs"); 2756 else 2757 s:=List(s,Int); 2758 Append(exp,s); 2759 s:=LinearCombinationPcgs(r.basrep[i],s); 2760 elm:=LeftQuotient(s,elm); 2761 fi; 2762 elif IsBound(r.bas[i]) then 2763 Append(exp,ListWithIdenticalEntries(Length(r.bas[i]),0)); 2764 fi; 2765#Print(exp); 2766 od; 2767 2768 return exp; 2769end; 2770 2771InstallMethod(ExponentsOfPcElement,"matrix residue pcgs",IsCollsElms, 2772 [IsPcgsResidueMatGroupRep and IsPcgs and IsPrimeOrdersPcgs,IsMatrix], 2773function(pcgs,x) 2774 return ExponentsResiduePcgs(pcgs!.decompInfo,x); 2775end); 2776 2777############################################################################# 2778## 2779#M RestrictedMapping(<hom>,<U>) 2780## 2781InstallMethod(RestrictedMapping,"for recognition mappings", 2782 CollFamSourceEqFamElms,[IsGroupGeneralMapping and 2783 HasRecogDecompinfoHomomorphism,IsGroup], 2784 # make this ranked higher than even some default subset methods in grp.gi 2785 SUM_FLAGS+100, 2786function(hom, U) 2787local d,gens,imgs,rest; 2788 d:=RecogDecompinfoHomomorphism(hom); 2789 gens:=GeneratorsOfGroup(U); 2790 if IsBound(d.fct) then 2791 imgs:=List(gens,d.fct); 2792 elif IsBound(d.isTrivial) then 2793 imgs:=List(gens,x->()); 2794 else 2795 imgs:=List(gens,x->CSIElementAct(d,x)); 2796 fi; 2797 2798 if IsPermGroup(U) and AssertionLevel()>2 then 2799 rest:=GroupHomomorphismByImages(U,Range(hom),gens,imgs); 2800 else 2801 RUN_IN_GGMBI:=true; # hack to skip Nice treatment 2802 if IsBound(d.fct) then 2803 rest:=GroupHomomorphismByFunction(U,Range(hom),d.fct); 2804 else 2805 rest:=GroupHomomorphismByImagesNC(U,Range(hom),gens,imgs); 2806 fi; 2807 RUN_IN_GGMBI:=false; 2808 fi; 2809 2810 SetRecogDecompinfoHomomorphism(rest,d); 2811 return rest; 2812end); 2813 2814# temporary -- missing in library 2815InstallMethod(MaximalSubgroupClassReps,"TF method",true, 2816 [IsGroup and IsFinite and 2817 HasFittingFreeLiftSetup],OVERRIDENICE,DoMaxesTF); 2818 2819TFSubgroupMembership:=function(g,u,elm) 2820local ffu,ff,x; 2821 ffu:=FittingFreeSubgroupSetup(g,u); 2822 ff:=ffu.parentffs; 2823 x:=ImagesRepresentative(ff.factorhom,elm); 2824 if not x in Image(ffu.rest) then 2825 return false; 2826 fi; 2827 elm:=elm/PreImagesRepresentative(ffu.rest,x); 2828 elm:=ImagesRepresentative(ff.pcisom,elm); 2829 if not IsBound(ffu.pcsub) then 2830 ffu.pcsub:=Subgroup(Image(ff.pcisom), 2831 List(ffu.pcgs,x->ImagesRepresentative(ff.pcisom,x))); 2832 fi; 2833 return elm in ffu.pcsub; 2834end; 2835 2836# careful -- this implicitly assumes we always stay in the parent 2837InstallMethod(\in,"ff subgroup",IsElmsColls, 2838 [IsMultiplicativeElementWithInverse,IsGroup and IsMatrixGroup], 2839 OVERRIDENICE, 2840function(elm,u) 2841 if not HasParentAttr(u) or not HasFittingFreeLiftSetup(ParentAttr(u)) then 2842 TryNextMethod(); 2843 fi; 2844 return TFSubgroupMembership(ParentAttr(u),u,elm); 2845end); 2846 2847# generic method, avoid nice method 2848TFNormalClosure:=function ( G, N ) 2849local gensG, # generators of the group <G> 2850 genG, # one generator of the group <G> 2851 gensN, # generators of the group <N> 2852 genN, # one generator of the group <N> 2853 cnj; # conjugated of a generator of <U> 2854 2855 gensG := GeneratorsOfGroup( G ); 2856 gensN := ShallowCopy( GeneratorsOfGroup( N ) ); 2857 for genN in gensN do 2858 for genG in gensG do 2859 cnj := genN ^ genG; 2860 if not TFSubgroupMembership(G,N,cnj) then 2861 Add( gensN, cnj ); 2862 N := SubgroupNC(G,gensN); 2863 fi; 2864 od; 2865 2866 od; 2867 2868 # return the normal closure 2869 return N; 2870 2871end; 2872 2873