1#############################################################################
2##
3##  This file is part of GAP, a system for computational discrete algebra.
4##  This file's authors include Bettina Eick and 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
12#############################################################################
13##
14#F  CollectedWordSQ( C, u, v )
15##
16##  The tail of  a conjugate  i^j  (i>j) or a   power i^p (i=j) is  stored at
17##  posiition (i^2-i)/2+j
18##
19InstallGlobalFunction( CollectedWordSQ, function( C, u, v )
20    local   w, p, c, m, g, n, i, j, x, mx, l1, l2, l;
21
22    # convert lists in to word/module pair
23    if IsList(v)  then
24        v := rec( word := v,  tail := [] );
25    fi;
26    if IsList(u)  then
27        u := rec( word := u,  tail := [] );
28    fi;
29
30    # if <v> is trivial  return <u>
31    if 0 = Length(v.word) and 0 = Length(v.tail)  then
32        return u;
33    fi;
34
35    # if <u> is trivial  return <v>
36    if 0 = Length(u.word) and 0 = Length(u.tail)  then
37        return v;
38    fi;
39
40    # if <v> has trivial word but a nontrivial tail add tails
41    if 0 = Length(v.word)  then
42        u := ShallowCopy(u);
43        for i  in [ 1 .. Length(v.tail) ]  do
44            if IsBound(v.tail[i])  then
45                if IsBound(u.tail[i])  then
46                    u.tail[i] := u.tail[i] + v.tail[i];
47                else
48                    u.tail[i] := v.tail[i];
49                fi;
50            fi;
51        od;
52        return u;
53    fi;
54
55    # unpack <u> into <x>
56    x := C.list;
57    n := Length(x);
58    for i  in [ 1 .. n ]  do
59        x[i] := 0;
60    od;
61    for i  in [ 1, 3 .. Length(u.word)-1 ]  do
62        x[u.word[i]] := u.word[i+1];
63    od;
64
65    # <mx> contains the tail of <x>
66    mx := ShallowCopy(u.tail);
67
68    # get stacks
69    w := C.wstack;
70    p := C.pstack;
71    c := C.cstack;
72    m := C.mstack;
73
74    # put <v> onto the stack
75    w[1] := v.word;
76    p[1] := 1;
77    c[1] := 1;
78    m[1] := ShallowCopy(v.tail);
79
80    # run until the stack is empty
81    l := 1;
82    while 0 < l  do
83
84        # remove next generator from stack
85        g := w[l][p[l]];
86
87        # apply generator to <mx>
88        for i  in [ 1 .. Length(mx) ]  do
89            if IsBound(mx[i])  then
90
91                # we use the transposed for technical reasons
92                mx[i] := C.module[g] * mx[i];
93            fi;
94        od;
95
96        # raise current exponent
97        c[l] := c[l]+1;
98
99        # if exponent is too big
100        if w[l][p[l]+1] < c[l]  then
101
102            # reset exponent
103            c[l] := 1;
104
105            # move position
106            p[l] := p[l] + 2;
107
108            # if position is too big
109            if Length(w[l]) < p[l]  then
110
111                # modify tail (add both tails)
112                l1 := Length(mx);
113                l2 := Length(m[l]);
114                for i  in [ 1 .. Minimum(l1,l2) ]  do
115                    if IsBound(mx[i])  then
116                        if IsBound(m[l][i])  then
117                            mx[i] := mx[i]+m[l][i];
118                        fi;
119                    elif IsBound(m[l][i])  then
120                        mx[i] := m[l][i];
121                    fi;
122                od;
123                if l1 < l2  then
124                    for i  in [ l1+1 .. l2 ]  do
125                        if IsBound(m[l][i])  then
126                            mx[i] := m[l][i];
127                        fi;
128                    od;
129                fi;
130
131                # and unbind word
132                m[l] := 0;
133                l := l - 1;
134            fi;
135        fi;
136
137        # now move generator to correct position
138        for i  in [ n, n-1 .. g+1 ]  do
139            if x[i] <> 0  then
140                l := l+1;
141                w[l] := C.relators[i][g];
142                c[l] := 1;
143                p[l] := 1;
144                l1   := (i^2-i)/2+g;
145                if not l1 in C.avoid  then
146                    if IsBound(mx[l1])  then
147                        mx[l1] := C.mone + mx[l1];
148                    else
149                        mx[l1] := C.mone;
150                    fi;
151                fi;
152                m[l] := mx;
153                mx := [];
154                l2 := [];
155                if not l1 in C.avoid  then
156                    l2[l1] := C.mone;
157                fi;
158                for j  in [ 2 .. x[i] ]  do
159                    l := l+1;
160                    w[l] := C.relators[i][g];
161                    c[l] := 1;
162                    p[l] := 1;
163                    m[l] := l2;
164                od;
165                x[i] := 0;
166            fi;
167        od;
168
169        # raise exponent
170        x[g] := x[g] + 1;
171
172        # and check for overflow
173        if x[g] = C.orders[g]  then
174            x[g] := 0;
175            if C.relators[g][g] <> 0  then
176                l1 := C.relators[g][g];
177                for i  in [ 1, 3 .. Length(l1)-1 ]  do
178                    x[l1[i]] := l1[i+1];
179                od;
180            fi;
181            l1 := (g^2+g)/2;
182            if not l1 in C.avoid  then
183                if IsBound(mx[l1])  then
184                    mx[l1] := C.mone + mx[l1];
185                else
186                    mx[l1] := C.mone;
187                fi;
188            fi;
189        fi;
190    od;
191
192    # and return result
193    w := [];
194    for i  in [ 1 .. Length(x) ]  do
195        if x[i] <> 0  then
196            Add( w, i );
197            Add( w, x[i] );
198        fi;
199    od;
200    return rec( word := w,  tail := mx );
201end );
202
203#############################################################################
204##
205#F  CollectorSQ( G, M, isSplit )
206##
207InstallGlobalFunction( CollectorSQ, function( G, M, isSplit )
208    local  r,  pcgs,  o,  i,  j,  word,  k, Gcoll;
209
210    # convert word into gen/exp form
211    word := function( pcgs, w )
212        local   r,  l,  k;
213        if OneOfPcgs( pcgs ) = w  then
214            r := 0;
215        else
216            l := ExponentsOfPcElement( pcgs, w );
217            r := [];
218            for k  in [ 1 .. Length(l) ]  do
219                if l[k] <> 0  then
220                    Add( r, k );
221                    Add( r, l[k] );
222                fi;
223            od;
224        fi;
225        return r;
226    end;
227
228    # convert relators into list of lists
229    if IsPcgs( G ) then
230        pcgs := G;
231    else
232        pcgs := Pcgs( G );
233    fi;
234    r := [];
235    o := RelativeOrders( pcgs );
236    for i  in [ 1 .. Length(pcgs) ]  do
237        r[i] := [];
238        for j  in [ 1 .. i-1 ]  do
239            r[i][j] := word( pcgs, pcgs[i]^pcgs[j]);
240        od;
241        r[i][i] := word( pcgs, pcgs[i]^o[i] );
242    od;
243
244    # create collector for G
245    Gcoll := rec( );
246    Gcoll.relators := r;
247    Gcoll.orders := o;
248
249    # create stacks
250    Gcoll.wstack := [];
251    Gcoll.estack := [];
252    Gcoll.pstack := [];
253    Gcoll.cstack := [];
254    Gcoll.mstack := [];
255
256    # create collector list
257    Gcoll.list := List( pcgs, x -> 0 );
258
259    # in case we are not interested in the module
260    if IsBool( M ) then return Gcoll; fi;
261
262    # copy collector and add module generators
263    r := ShallowCopy(Gcoll);
264
265    # create module gens (the transposed is a technical detail)
266    r.module:=List(M.generators,i->ImmutableMatrix(M.field,TransposedMat(i)));
267    r.mone  :=ImmutableMatrix(M.field,(IdentityMat(M.dimension, M.field)));
268    r.mzero :=ImmutableMatrix(M.field,Zero(r.mone));
269
270    # add avoid
271    r.avoid := [];
272    if isSplit  then
273        k := Characteristic( M.field );
274        for i  in [ 1 .. Length(r.orders) ]  do
275            for j  in [ 1 .. i ]  do
276                # we can avoid
277                if Order(pcgs[i]) mod k <>0 and Order(pcgs[j]) mod k <>0 then
278                # was: r.orders[i] <> k and r.orders[j] <> k  then
279                    AddSet( r.avoid, (i^2-i)/2 + j );
280                fi;
281            od;
282        od;
283    fi;
284
285    # and return collector
286    return r;
287end );
288
289#############################################################################
290##
291#F  AddEquationsSQ( eq, t1, t2 )
292##
293InstallGlobalFunction( AddEquationsSQ, function( eq, t1, t2 )
294    local   i,  j,  l,  v,  w,  x,  n,  c;
295
296    # if <t1> = <t2>  return
297    if t1 = t2  then return; fi;
298
299    # compute <t1> - <t2>
300    t1 := ShallowCopy(t1);
301    for i  in [ 1 .. Length(t2) ]  do
302        if IsBound(t2[i])  then
303            if IsBound(t1[i])  then
304                t1[i] := t1[i] - t2[i];
305            else
306                t1[i] := -t2[i];
307            fi;
308        fi;
309    od;
310
311    # make lines
312    l := List( eq.vzero,  x -> [] );
313    v := [];
314    for i  in [ 1 .. Length(t1) ]  do
315        if IsBound(t1[i])  then
316            for j  in [ 1 .. eq.dimension ]  do
317                if t1[i][j] <> eq.vzero then
318                    l[j][i] := ShallowCopy(t1[i][j]);
319                    AddCoeffs( l[j][i], v );
320                    ShrinkRowVector(l[j][i]);
321                fi;
322            od;
323        fi;
324    od;
325
326    # and reduce lines
327    n := eq.dimension;
328    for j  in [ 1 .. n ]  do
329        x := l[j];
330        v := Length(x);
331        if 0 < v  then w := (v-1)*n + Length(x[v]);  fi;
332        while 0 < v and IsBound(eq.system[w])  do
333            c := -x[v][Length(x[v])];
334            for i  in eq.spos[w]  do
335                if IsBound(x[i])  then
336                    x[i] := ShallowCopy( x[i] );
337                    AddCoeffs( x[i], eq.system[w][i], c );
338                    ShrinkRowVector(x[i]);
339                    if 0 = Length(x[i])  then
340                        Unbind(x[i]);
341                    fi;
342                else
343                    x[i] := c * eq.system[w][i];
344                fi;
345            od;
346            v := Length(x);
347            if 0 < v  then w := (v-1)*n + Length(x[v]);  fi;
348        od;
349        if 0 < v  then
350            eq.system[w] := x * (1/x[v][Length(x[v])]);
351            eq.spos[w]   := Filtered( [1..eq.nrels], t -> IsBound(x[t]) );
352        fi;
353    od;
354end );
355
356#############################################################################
357##
358#F  SolutionSQ( C, eq )
359##
360InstallGlobalFunction( SolutionSQ, function( C, eq )
361    local   x,  e,  d,  t,  j,  v,  i,  n,  p,  w;
362
363    # construct null vector
364    n := [];
365    for i  in [ 1 .. eq.nrels-Length(C.avoid) ]  do
366        Append( n, eq.vzero );
367    od;
368
369    # generated position
370    C.unavoidable := [];
371    j := 1;
372    for i  in [ 1 .. eq.nrels ]  do
373        if not i in C.avoid  then
374            C.unavoidable[i] := j;
375            j := j+1;
376        fi;
377    od;
378
379    # blow up vectors
380    t := [];
381    w := [];
382    for j  in [ 1 .. Length(eq.system) ]  do
383        if IsBound(eq.system[j])  then
384            v := ShallowCopy(n);
385            for i  in eq.spos[j]  do
386                if not i in C.avoid  then
387                    p := eq.dimension*(C.unavoidable[i]-1);
388                    v{[p+1..p+Length(eq.system[j][i])]} := eq.system[j][i];
389                fi;
390            od;
391            ShrinkRowVector(v);
392            t[Length(v)] := v;
393            AddSet( w, Length(v) );
394        fi;
395    od;
396
397    # normalize system
398    v := 0*eq.vzero[1];
399    for i  in w  do
400        for j  in w  do
401            if j > i  then
402                p := t[j][i];
403                if p <> v  then
404                    t[j] := ShallowCopy( t[j] );
405                    AddCoeffs( t[j], t[i], -p );
406                    ShrinkRowVector(t[j]);
407                fi;
408            fi;
409        od;
410    od;
411
412    # compute homogeneous solution
413    d := Difference( [ 1 .. (eq.nrels-Length(C.avoid))*eq.dimension ],  w );
414    v := [];
415    e := eq.vzero[1]^0;
416    for i  in d  do
417        x := ShallowCopy(n);
418        x[i] := e;
419        for j  in w  do
420            if j >= i  then
421                x[j] := -t[j][i];
422            fi;
423        od;
424        Add( v, x );
425    od;
426    if 0 = Length(C.avoid)  then
427        return v;
428    fi;
429
430    # construct null vector
431    n := [];
432    for i  in [ 1 .. eq.nrels ]  do
433        Append( n, eq.vzero );
434    od;
435
436    # construct a blow up matrix
437    i := [];
438    for j  in [ 1 .. eq.nrels ]  do
439        if not j in C.avoid  then
440            Append( i, (j-1)*eq.dimension + [ 1 .. eq.dimension ] );
441        fi;
442    od;
443
444    # blowup the vectors
445    e := [];
446    for x  in v  do
447        d := ShallowCopy(n);
448        d{i} := x;
449        Add( e, d );
450    od;
451
452    # and return
453    return e;
454end );
455
456#############################################################################
457##
458#F  TwoCocyclesSQ( C, G, M )
459##
460InstallGlobalFunction( TwoCocyclesSQ, function( C, G, M )
461    local   pairs, i,  j,  k, w1, w2, eq, p, n;
462
463    # get number of generators
464    n := Length(Pcgs(G));
465
466    # collect equations in <eq>
467    eq := rec( vzero     := C.mzero[1],
468               mzero     := C.mzero,
469               dimension := Length(C.mzero),
470               nrels     := (n^2+n)/2,
471               spos      := [],
472               system    := [] );
473
474    # precalculate (ij) for i > j
475    pairs := List( [1..n], x -> [] );
476    for i  in [ 2 .. n ]  do
477        for j  in [ 1 .. i-1 ]  do
478            pairs[i][j] := CollectedWordSQ( C, [i,1], [j,1] );
479        od;
480    od;
481
482    # consistency 1:  k(ji) = (kj)i
483    for i  in [ n, n-1 .. 1 ]  do
484        for j  in [ n, n-1 .. i+1 ]  do
485            for k  in [ n, n-1 .. j+1 ]  do
486                w1 := CollectedWordSQ( C, [k,1], pairs[j][i] );
487                w2 := CollectedWordSQ( C, pairs[k][j], [i,1] );
488                if w1.word <> w2.word  then
489                    Error( "k(ji) <> (kj)i" );
490                else
491                    AddEquationsSQ( eq, w1.tail, w2.tail );
492                fi;
493            od;
494        od;
495    od;
496
497    # consistency 2: j^(p-1) (ji) = j^p i
498    for i  in [ n, n-1 .. 1 ]  do
499        for j  in [ n, n-1 .. i+1 ]  do
500            p  := C.orders[j];
501            w1 := CollectedWordSQ( C, [j,p-1],
502                    CollectedWordSQ( C, [j,1], [i,1] ) );
503            w2 := CollectedWordSQ( C, CollectedWordSQ( C, [j,p-1], [j,1] ),
504                    [i,1] );
505            if w1.word <> w2.word  then
506                Error( "j^(p-1) (ji) <> j^p i" );
507            else
508                AddEquationsSQ( eq, w1.tail, w2.tail );
509            fi;
510        od;
511    od;
512
513    # consistency 3: k (i i^(p-1)) = (ki) i^p-1
514    for i  in [ n, n-1 .. 1 ]  do
515        p := C.orders[i];
516        for k  in [ n, n-1 .. i+1 ]  do
517            w1 := CollectedWordSQ( C, [k,1],
518                    CollectedWordSQ( C, [i,1], [i,p-1] ) );
519            w2 := CollectedWordSQ( C, CollectedWordSQ( C, [k,1], [i,1] ),
520                    [i,p-1] );
521            if w1.word <> w2.word  then
522                Error( "k i^p <> (ki) i^(p-1)" );
523            else
524                AddEquationsSQ( eq, w1.tail, w2.tail );
525            fi;
526        od;
527    od;
528
529    # consistency 4: (i i^(p-1)) i = i (i^(p-1) i)
530    for i  in [ n, n-1 .. 1 ]  do
531        p  := C.orders[i];
532        w1 := CollectedWordSQ( C, CollectedWordSQ( C, [i,1], [i,p-1] ),
533                [i,1] );
534        w2 := CollectedWordSQ( C, [i,1],
535                CollectedWordSQ( C, [i,p-1], [i,1] ) );
536        if w1.word <> w2.word  then
537            Error( "i i^p-1 <> i^p" );
538        else
539            AddEquationsSQ( eq, w1.tail, w2.tail );
540        fi;
541    od;
542
543    # and return solution
544    return SolutionSQ( C, eq );
545end );
546
547#############################################################################
548##
549#F  TwoCoboundariesSQ( C, G, M )
550##
551InstallGlobalFunction( TwoCoboundariesSQ, function( C, G, M )
552    local   n,  R,  MI,  j,  i,  x,  m,  e,  k,  r,  d;
553
554    # start with zero matrix
555    n := Length(Pcgs( G ));
556    R := [];
557    r := n*(n+1)/2;
558    for i  in [ 1 .. n ]  do
559        R[i] := [];
560        for j in [ 1 .. r ] do
561            R[i][j] := C.mzero;
562        od;
563    od;
564
565    # compute inverse generators
566    M  := M.generators;
567    MI := List( M, x -> x^-1 );
568    d  := Length(M[1]);
569
570    # loop over all relators
571    for j  in [ 1 .. n ]  do
572        for i in  [ j .. n ]  do
573            x := (i^2-i)/2 + j;
574
575            # power relator
576            if i = j  then
577                m := C.mone;
578                for e  in [ 1 .. C.orders[i] ]  do
579                    R[i][x] := R[i][x] - m;  m := M[i] * m;
580                od;
581
582            # conjugate
583            else
584                R[i][x] := R[i][x] - M[j];
585                R[j][x] := R[j][x] + MI[j]*M[i]*M[j] - C.mone;
586            fi;
587
588            # compute fox derivatives
589            m := C.mone;
590            r := C.relators[i][j];
591            if r <> 0  then
592                for k  in [ Length(r)-1, Length(r)-3 .. 1 ]  do
593                    for e  in [ 1 .. r[k+1] ]  do
594                        R[r[k]][x] := R[r[k]][x] + m;
595                        m := M[r[k]] * m;
596                    od;
597                od;
598            fi;
599        od;
600    od;
601
602    # make one list
603    m := [];
604    r := n*(n+1)/2;
605    for i  in [ 1 .. n ]  do
606        for k  in [ 1 .. d ]  do
607            e := [];
608            for j  in [ 1 .. r ]  do
609                Append( e, R[i][j][k] );
610            od;
611            Add( m, e );
612        od;
613    od;
614
615    # compute a base for <m>
616    return BaseMat(m);
617end );
618
619#############################################################################
620##
621#F  TwoCohomologySQ( C, G, M )
622##
623InstallGlobalFunction( TwoCohomologySQ, function( C, G, M )
624    local cc, cb;
625    cc := TwoCocyclesSQ( C, G, M );
626    if Length( cc ) > 0 then
627        cb := TwoCoboundariesSQ( C, G, M );
628        if Length( C.avoid ) > 0 then
629            cb := SumIntersectionMat( cc, cb )[2];
630        fi;
631        if Length( cb ) > 0 then
632            cc := BaseSteinitzVectors( cc, cb ).factorspace;
633        fi;
634    fi;
635    return cc;
636end );
637
638#############################################################################
639##
640#M  TwoCocycles( G, M )
641##
642InstallMethod( TwoCocycles,
643    "generic method for pc groups",
644    true,
645    [ IsPcGroup, IsObject ],
646    0,
647
648function( G, M )
649    local C;
650    C := CollectorSQ( G, M, false );
651    return TwoCocyclesSQ( C, G, M );
652end );
653
654#############################################################################
655##
656#M  TwoCoboundaries( G, M )
657##
658InstallMethod( TwoCoboundaries,
659    "generic method for pc groups",
660    true,
661    [ IsPcGroup, IsObject ],
662    0,
663
664function( G, M )
665    local C;
666    C := CollectorSQ( G, M, false );
667    return TwoCoboundariesSQ( C, G, M );
668end );
669
670# non-solvable cohomology based on rewriting
671
672# return isomorphism G-fp and fp->mon, such that presentation of monoid is
673# confluent (wrt wreath order). Returns list [fphom,monhom,ordering]
674BindGlobal("ConfluentMonoidPresentationForGroup",function(G)
675local iso,fp,n,dec,homs,mos,i,j,ffp,imo,m,k,gens,fm,mgens,rules,
676      loff,off,monreps,left,right,fmgens,r,diff,monreal,nums,reduce,hom,dept;
677  if IsSymmetricGroup(G) then
678    i:=SymmetricGroup(SymmetricDegree(G));
679    iso:=CheapIsomSymAlt(G,i)*IsomorphismFpGroup(i);
680    fp:=Range(iso);
681    hom:=IsomorphismFpMonoid(fp);
682    m:=Range(hom);
683    fm:=FreeMonoidOfFpMonoid(m);
684    k:=KnuthBendixRewritingSystem(m);
685    MakeConfluent(k);
686    rules:=Rules(k);
687    dept:=fail;
688  else
689    iso:=IsomorphismFpGroupByChiefSeries(G:rewrite);
690
691    fp:=Range(iso);
692    gens:=GeneratorsOfGroup(fp);
693    n:=Length(gens);
694    dec:=iso!.decompinfo;
695
696    fmgens:=[];
697    mgens:=[];
698    for i in gens do
699      Add(fmgens,i);
700      Add(fmgens,i^-1);
701      Add(mgens,String(i));
702      Add(mgens,String(i^-1));
703    od;
704    nums:=List(fmgens,x->LetterRepAssocWord(UnderlyingElement(x))[1]);
705    fm:=FreeMonoid(mgens);
706    mgens:=GeneratorsOfMonoid(fm);
707    rules:=[];
708    reduce:=function(w)
709    local red,i,p;
710      w:=LetterRepAssocWord(w);
711      repeat
712        i:=1;
713        red:=false;
714        while i<=Length(rules) and red=false do
715          p:=PositionSublist(w,LetterRepAssocWord(rules[i][1]));
716          if p<>fail then
717            #Print("Apply ",rules[i],p,w,"\n");
718            w:=Concatenation(w{[1..p-1]},LetterRepAssocWord(rules[i][2]),
719              w{[p+Length(rules[i][1])..Length(w)]});
720            red:=true;
721          else
722            i:=i+1;
723          fi;
724        od;
725      until red=false;
726      return AssocWordByLetterRep(FamilyObj(One(fm)),w);
727    end;
728
729
730    homs:=ShallowCopy(dec.homs);
731    mos:=[];
732    off:=Length(mgens);
733    dept:=[];
734    # go up so we may reduce tails
735    for i in [Length(homs),Length(homs)-1..1] do
736      Add(dept,off);
737      if IsPcgs(homs[i]) then
738        ffp:=AbelianGroup(IsFpGroup,RelativeOrders(homs[i]));
739      else
740        ffp:=Range(dec.homs[i]);
741      fi;
742      imo:=IsomorphismFpMonoid(ffp);
743      Add(mos,imo);
744      m:=Range(imo);
745      loff:=off-Length(GeneratorsOfMonoid(m));
746      monreps:=fmgens{[loff+1..off]};
747      monreal:=mgens{[loff+1..off]};
748      if IsBound(m!.rewritingSystem) then
749        k:=m!.rewritingSystem;
750      else
751        k:=KnuthBendixRewritingSystem(m);
752      fi;
753      MakeConfluent(k);
754      # convert rules
755      for r in Rules(k) do
756        left:=MappedWord(r[1],FreeGeneratorsOfFpMonoid(m),monreps);
757        right:=MappedWord(r[2],FreeGeneratorsOfFpMonoid(m),monreps);
758        diff:=LeftQuotient(PreImagesRepresentative(iso,right),
759                PreImagesRepresentative(iso,left));
760        diff:=ImagesRepresentative(iso,diff);
761
762        left:=MappedWord(r[1],FreeGeneratorsOfFpMonoid(m),monreal);
763        right:=MappedWord(r[2],FreeGeneratorsOfFpMonoid(m),monreal);
764        if not IsOne(diff) then
765          right:=right*Product(List(LetterRepAssocWord(UnderlyingElement(diff)),
766            x->mgens[Position(nums,x)]));
767        fi;
768        right:=reduce(right); # monoid word might change
769        Add(rules,[left,right]);
770      od;
771      for j in [loff+1..off] do
772        # if the generator gets reduced away, won't need to use it
773        if reduce(mgens[j])=mgens[j] then
774          for k in [off+1..Length(mgens)] do
775            if reduce(mgens[k])=mgens[k] then
776              right:=fmgens[j]^-1*fmgens[k]*fmgens[j];
777              #collect
778              right:=ImagesRepresentative(iso,PreImagesRepresentative(iso,right));
779              right:=Product(List(LetterRepAssocWord(UnderlyingElement(right)),
780                x->mgens[Position(nums,x)]));
781              right:=reduce(mgens[j]*right);
782              Add(rules,[mgens[k]*mgens[j],right]);
783            fi;
784          od;
785        fi;
786      od;
787      #if i<Length(homs) then Error("ZU");fi;
788      off:=loff;
789    od;
790    Add(dept,off);
791    # calculate levels for ordering
792    dept:=dept+1;
793    dept:=List([1..Length(mgens)],
794      x->PositionProperty(dept,y->x>=y)-1);
795
796    if ForAny(rules,x->x[2]<>reduce(x[2])) then Error("irreduced right");fi;
797
798    # inverses are true inverses, also for extension
799    for i in [1..Length(gens)] do
800      left:=mgens[2*i-1]*mgens[2*i];
801      left:=reduce(left);
802      if left<>One(fm) then Add(rules,[left,One(fm)]); fi;
803      left:=mgens[2*i]*mgens[2*i-1];
804      left:=reduce(left);
805      if left<>One(fm) then Add(rules,[left,One(fm)]); fi;
806    od;
807  fi;
808
809  # finally create
810  m:=FactorFreeMonoidByRelations(fm,rules);
811  mgens:=GeneratorsOfMonoid(m);
812
813  hom:=MagmaIsomorphismByFunctionsNC(fp,m,
814        function(w)
815          local l,i;
816          l:=[];
817          for i in LetterRepAssocWord(UnderlyingElement(w)) do
818            if i>0 then Add(l,2*i-1);
819            else Add(l,-2*i);fi;
820          od;
821          return ElementOfFpMonoid(FamilyObj(One(m)),
822                  AssocWordByLetterRep(FamilyObj(One(fm)),l));
823        end,
824        function(w)
825          local g,i,x;
826          g:=[];
827          for i in LetterRepAssocWord(UnderlyingElement(w)) do
828            if IsOddInt(i) then x:=(i+1)/2;
829            else x:=-i/2;fi;
830            # word must be freely cancelled
831            if Length(g)>0 and x=-g[Length(g)] then
832              Unbind(g[Length(g)]);
833            else Add(g,x); fi;
834          od;
835          return ElementOfFpGroup(FamilyObj(One(fp)),
836                  AssocWordByLetterRep(FamilyObj(One(FreeGroupOfFpGroup(fp))),g));
837        end);
838
839  hom!.type:=1;
840  if not HasIsomorphismFpMonoid(G) then
841    SetIsomorphismFpMonoid(G,hom);
842  fi;
843  if dept=fail then
844    return [iso,hom,k!.ordering];
845  else
846    return [iso,hom,WreathProductOrdering(fm,dept)];
847  fi;
848end);
849
850#############################################################################
851##
852#M  TwoCohomology( G, M )
853##
854InstallMethod( TwoCohomology,
855    "generic method for pc groups",
856    true,
857    [ IsPcGroup, IsObject ],
858    0,
859
860function( G, M )
861    local C, d, z, co, cb, pr;
862    C := CollectorSQ( G, M, false );
863    d := Length( C.orders );
864    d := d * (d+1) / 2;
865    z := Flat( List( [1..d], x -> C.mzero[1] ) );
866    co := TwoCocyclesSQ( C, G, M );
867    co := VectorSpace( M.field, co, z );
868    cb := TwoCoboundariesSQ( C, G, M );
869    cb := SubspaceNC( co, cb );
870    pr := FpGroupPcGroupSQ( G );
871    return rec( group := G,
872                module := M,
873                collector := C,
874                isPcCohomology:=true,
875                cohom :=
876                NaturalHomomorphismBySubspaceOntoFullRowSpace(co,cb),
877                presentation := FpGroupPcGroupSQ( G ) );
878end );
879
880# code for generic 2-cohomology (solvable not required)
881
882InstallMethod( TwoCohomologyGeneric,"generic, using rewriting system",true,
883  [IsGroup and IsFinite,IsObject],0,
884function(G,mo)
885local field,fp,fpg,gens,hom,mats,fm,mon,kb,tzrules,dim,rules,eqs,i,j,k,l,o,l1,
886      len1,l2,m,start,formalinverse,hastail,one,zero,new,v1,v2,collectail,
887      findtail,colltz,mapped,mapped2,onemat,zerovec,dict,max,mal,s,p,
888      c,nvars,htpos,zeroq,r,ogens,bds,model,q,pre,pcgs,miso,ker,solvec,rulpos;
889
890
891  # collect the word in factor group
892  colltz:=function(a)
893  local i,j,s,mm,p;
894
895    # collect from left
896    i:=1;
897    while i<=Length(a) do
898
899      # does a rule apply at position i?
900      j:=0;
901      s:=0;
902      mm:=Minimum(mal,Length(a)-i+1);
903      while j<mm do
904        s:=s*max+a[i+j];
905        p:=LookupDictionary(dict,s);
906        if p<>fail then break; fi;
907        j:=j+1;
908      od;
909
910      if p<>fail then
911        a:=Concatenation(a{[1..i-1]},tzrules[p][2],
912          a{[i+Length(tzrules[p][1])..Length(a)]});
913        i:=Maximum(0,i-mal); # earliest which could be affected
914      fi;
915      i:=i+1;
916    od;
917
918    return a;
919  end;
920
921  # matrix corresponding to monoid word
922  mapped:=function(list)
923  local a,i;
924    a:=onemat;
925    for i in list do
926      a:=a*mats[i];
927    od;
928    return a;
929  end;
930
931  # normalform word and collect the tails
932  collectail:=function(wrd)
933  local v,tail,i,j,s,p,mm;
934    v:=List(rules,x->zero);
935
936    # collect from left
937    i:=1;
938    while i<=Length(wrd) do
939
940      # does a rule apply at position i?
941      j:=0;
942      s:=0;
943      mm:=Minimum(mal,Length(wrd)-i+1);
944      while j<mm do
945        s:=s*max+wrd[i+j];
946        p:=LookupDictionary(dict,s);
947        if p<>fail and rulpos[p]<>fail then break; fi;
948        j:=j+1;
949      od;
950
951      if p<>fail and rulpos[p]<>fail then
952        p:=rulpos[p];
953        tail:=wrd{[i+Length(rules[p][1])..Length(wrd)]};
954        wrd:=Concatenation(wrd{[1..i-1]},rules[p][2],tail);
955#Print("Apply ",p,"@",i,":",wrd,"\n");
956        if p in hastail then v[p]:=v[p]+mapped(tail); fi;
957        i:=Maximum(0,i-mal); # earliest which could be affected
958      fi;
959      i:=i+1;
960    od;
961
962    return [wrd,v];
963  end;
964
965  field:=mo.field;
966
967  ogens:=GeneratorsOfGroup(G);
968
969  if false then
970    #  old general KB code, left for debugging
971    fp:=IsomorphismFpGroup(G);
972    fpg:=Range(fp);
973    fm:=IsomorphismFpMonoid(fpg);
974    mon:=Range(fm);
975
976    if IsBound(mon!.confl) then
977      tzrules:=mon!.confl;
978    else
979      kb:=KnuthBendixRewritingSystem(mon);
980      MakeConfluent(kb);
981      tzrules:=kb!.tzrules;
982      mon!.confl:=tzrules;
983    fi;
984  else
985    # new approach with RWS from chief series
986    mon:=ConfluentMonoidPresentationForGroup(G);
987    fp:=mon[1];
988    fpg:=Range(fp);
989    fm:=mon[2];
990    mon:=Range(fm);
991    tzrules:=List(RelationsOfFpMonoid(mon),x->List(x,LetterRepAssocWord));
992  fi;
993
994  # build data structure to find rule applicable at given position. Assumes
995  # that rule set is reduced.
996  max:=Maximum(Union(List(tzrules,x->x[1])))+1;
997  mal:=Maximum(List(tzrules,x->Length(x[1])));
998  dict:=NewDictionary(max,Integers,true);
999  for i in [1..mal] do
1000    p:=Filtered([1..Length(tzrules)],x->Length(tzrules[x][1])=i);
1001    for j in p do
1002      s:=0;
1003      for k in [1..i] do
1004        s:=s*max+tzrules[j][1][k];
1005      od;
1006      AddDictionary(dict,s,j);
1007    od;
1008  od;
1009
1010  gens:=List(GeneratorsOfGroup(FamilyObj(fpg)!.wholeGroup),
1011    x->PreImagesRepresentative(fp,x));
1012
1013  hom:=GroupHomomorphismByImagesNC(G,Group(mo.generators),GeneratorsOfGroup(G),mo.generators);
1014  mo:=GModuleByMats(List(gens,x->ImagesRepresentative(hom,x)),mo.field); # new gens
1015
1016  #rules:=ShallowCopy(kb!.tzrules);
1017  #hastail:=Filtered([1..Length(rules)],x->Length(rules[x][1])<>2 or
1018  #  Length(rules[x][2])>0 or formalinverse[rules[x][1][1]]<>rules[x][1][2]);
1019  #IsSet(hastail); # quick membership test
1020
1021  l1:=GeneratorsOfGroup(fpg);
1022  l1:=Concatenation(l1,List(l1,Inverse));
1023  formalinverse:=[];
1024  for i in l1 do
1025    j:=LetterRepAssocWord(UnderlyingElement(ImagesRepresentative(fm,i)));
1026    o:=LetterRepAssocWord(UnderlyingElement(ImagesRepresentative(fm,i^-1)));
1027    if Length(j)<>1 or Length(o)<>1 then Error("length!");fi;
1028    formalinverse[j[1]]:=o[1];
1029  od;
1030
1031  # rules that describe formal inverses, or delete generators of order 2 get no
1032  # tails.
1033  hastail:=[];
1034  rules:=[];
1035  for r in tzrules do
1036    if Length(r[1])>=2 then
1037      Add(rules,r);
1038      if Length(r[1])>2 or
1039        (Length(r[1])=2 and (Length(r[2])>0 or formalinverse[r[1][1]]<>r[1][2]))
1040        then
1041          AddSet(hastail,Length(rules));
1042      fi;
1043    elif Length(r[1])>1 then
1044      if Length(r[2])=0 then Error("generator is trivial");fi;
1045      if Length(r[2])<>1 or formalinverse[r[1][1]]<>r[2][1] then
1046        Add(rules,r);
1047        AddSet(hastail,Length(rules));
1048      else
1049#Print("Not use: ",r,"\n");
1050        # Do not use these rules for overlaps
1051        if r[2][1]>r[1][1] then
1052          Error("code assumes that larger number gets reduced to smaller");
1053        fi;
1054        if ForAny(rules,x->r[1][1] in x[2] or (x<>r and r[1][1] in x[1])) then
1055          Error("rules are not reduced");
1056        fi;
1057      fi;
1058    else
1059      # Length of r[1] is 1. That is this generator is not used.
1060      # check that it is really just an inverse that goes away, otherwise
1061      # awkward.
1062      if formalinverse[r[1][1]]>r[1][1] then
1063        Error("generator vanishes?");
1064      fi;
1065    fi;
1066  od;
1067
1068  rulpos:=List(tzrules,x->PositionProperty(rules,y->y[1]=x[1]));
1069
1070  htpos:=List([1..Length(rules)],x->Position(hastail,x));
1071
1072  model:=ValueOption("model");
1073  if model<>fail then
1074    q:=GQuotients(model,G)[1];
1075    pre:=List(gens,x->PreImagesRepresentative(q,x));
1076    ker:=KernelOfMultiplicativeGeneralMapping(q);
1077    pcgs:=Pcgs(ker);
1078    l1:=GModuleByMats(LinearActionLayer(Group(pre),pcgs),mo.field);
1079    MTX.IsIrreducible(mo);
1080    miso:=MTX.Isomorphism(mo,l1);
1081    new:=List(miso,x->PcElementByExponents(pcgs,x));
1082    pcgs:=PcgsByPcSequence(FamilyObj(One(model)),new);
1083    # now calculate the vector
1084    solvec:=[];
1085    one:=One(model);
1086    m:=GroupGeneralMappingByImagesNC(fpg,model,GeneratorsOfGroup(fpg),pre);
1087    mats:=List(GeneratorsOfMonoid(mon),
1088      x->ImagesRepresentative(m,
1089      PreImagesRepresentative(fm,x))); #Elements for monoid generators
1090    pre:=mats;
1091    onemat:=One(G);
1092    for i in [1..Length(hastail)] do
1093      r:=rules[hastail[i]];
1094      m:=LeftQuotient(mapped(r[2]),mapped(r[1]));
1095      m:=ExponentsOfPcElement(pcgs,m);
1096      solvec:=Concatenation(solvec,m*One(field));
1097    od;
1098  fi;
1099
1100  onemat:=One(mo.generators[1]);
1101  zerovec:=Zero(onemat[1]);
1102
1103  mats:=List(GeneratorsOfMonoid(mon),
1104    x->ImagesRepresentative(hom,PreImagesRepresentative(fp,
1105    PreImagesRepresentative(fm,x)))); # matrices for monoid generators
1106  one:=One(mats[1]);
1107  zero:=zerovec;
1108  dim:=Length(zero);
1109  nvars:=dim*Length(hastail); #Number of variables
1110
1111  eqs:=[];
1112  zeroq:=ImmutableVector(field,ListWithIdenticalEntries(nvars,Zero(field)));
1113  for i in [1..Length(rules)] do
1114    l1:=rules[i][1];
1115    len1:=Length(l1);
1116    for j in [1..Length(rules)] do
1117      l2:=rules[j][1];
1118      m:=Minimum(len1,Length(l2));
1119      for o in [1..m-1] do # possible overlap Length
1120        start:=len1-o;
1121        if ForAll([1..o],k->l1[start+k]=l2[k]) then
1122          #Print("Overlap ",l1," ",l2," ",o,"\n");
1123
1124          # apply l1 first
1125          new:=Concatenation(rules[i][2],l2{[o+1..Length(l2)]});
1126          c:=collectail(new);
1127          v1:=c[2];c:=c[1];
1128          if i in hastail then
1129            v1[i]:=v1[i]+mapped(l2{[o+1..Length(l2)]});
1130          fi;
1131
1132          # apply l2 first
1133          new:=Concatenation(l1{[1..len1-o]},rules[j][2]);
1134          v2:=collectail(new);
1135          if c<>v2[1] then Error("different in factor");
1136          #else Print("Both reduce to ",c,"\n");
1137          fi;
1138          v2:=v2[2];
1139          if j in hastail then
1140            v2[j]:=v2[j]+one;
1141          fi;
1142          if v1<>v2 then # If entries stay zero they are identical (and thus test cheaply)
1143            c:=List([1..dim],x->ShallowCopy(zeroq));
1144            for k in hastail do
1145              if v1[k]<>v2[k] then
1146                new:=TransposedMat(v1[k]-v2[k]);
1147                r:=(htpos[k]-1)*dim+[1..dim];
1148                for l in [1..dim] do
1149                  if not IsZero(new[l]) then
1150                    c[l]{r}:=new[l];
1151                  fi;
1152                od;
1153              fi;
1154            od;
1155            for k in c do
1156              if not IsZero(k) then
1157                if model<>fail and not IsZero(solvec*k) then
1158                  Error("model does not fit");
1159                fi;
1160                AddSet(eqs,ImmutableVector(field,k));
1161              fi;
1162            od;
1163          fi;
1164        fi;
1165      od;
1166    od;
1167  od;
1168  eqs:=Filtered(TriangulizedMat(eqs),x->not IsZero(x));
1169  eqs:=NullspaceMat(TransposedMat(eqs)); # basis of cocycles
1170
1171  # Now get Coboundaries
1172
1173  # different collection function: Sum up changes for generator i
1174  findtail:=function(wrd,nums,vecs)
1175  local a,i,j,p,v;
1176    a:=zerovec;
1177    for i in [1..Length(wrd)] do
1178      p:=Position(nums,wrd[i]);
1179      if p<>fail then
1180        v:=vecs[p];
1181        for j in [i+1..Length(wrd)] do
1182          v:=v*mats[wrd[j]];
1183        od;
1184        a:=a+v;
1185      fi;
1186    od;
1187    return a;
1188  end;
1189
1190  bds:=[];
1191  for i in [1..Length(mats)] do
1192    if formalinverse[i]>i then
1193      c:=[i,formalinverse[i]];
1194      for k in one do
1195        r:=[k,-k*mats[formalinverse[i]]];
1196
1197        new:=[];
1198        for j in hastail do
1199          v1:=findtail(rules[j][1],c,r);
1200          v2:=findtail(rules[j][2],c,r);
1201          Append(new,v1-v2);
1202        od;
1203        new:=ImmutableVector(field,new);
1204        Assert(0,SolutionMat(eqs,new)<>fail);
1205        Add(bds,new);
1206      od;
1207    fi;
1208  od;
1209  bds:=Filtered(TriangulizedMat(bds),x->not IsZero(x));
1210  bds:=List(bds,Immutable);
1211
1212  if gens<>GeneratorsOfGroup(G) then
1213    G:=GroupWithGenerators(gens);
1214  fi;
1215  r:=rec(group:=G,module:=mo,cocycles:=eqs,coboundaries:=bds,zero:=zeroq,
1216         prime:=Size(field));
1217  new:=List(BaseSteinitzVectors(eqs,bds).factorspace,x->ImmutableVector(field,x));
1218  r.cohomology:=new;
1219
1220  new:=[];
1221  one:=One(FreeGroupOfFpGroup(fpg));
1222
1223  k:=List(GeneratorsOfMonoid(mon),
1224    x->UnderlyingElement(PreImagesRepresentative(fm,x)));
1225  # matrix corresponding to monoid word
1226  mapped2:=function(list)
1227  local a,i;
1228    a:=one;
1229    for i in list do
1230      a:=a*k[i];
1231    od;
1232    return a;
1233  end;
1234
1235  for i in hastail do
1236    # rule that would get the tail
1237    Add(new,LeftQuotient(mapped2(rules[i][1]),mapped2(rules[i][2])));
1238  od;
1239  r.presentation:=rec(group:=FreeGroupOfFpGroup(fpg),relators:=new,
1240    # position of relators with tails in tzrules
1241    monrulpos:=List(rules{hastail},x->Position(tzrules,x)),
1242    prewords:=List(ogens,x->UnderlyingElement(ImagesRepresentative(fp,x))));
1243
1244  # normalform word and collect the tails
1245  r.tailforword:=function(wrd,zy)
1246  local v,i,j,s,p,mm,w,tail;
1247    v:=zerovec;
1248
1249    # collect from left
1250    i:=1;
1251    while i<=Length(wrd) do
1252
1253      # does a rule apply at position i?
1254      j:=0;
1255      s:=0;
1256      mm:=Minimum(mal,Length(wrd)-i+1);
1257      while j<mm do
1258        s:=s*max+wrd[i+j];
1259        p:=LookupDictionary(dict,s);
1260        if p<>fail and rulpos[p]<>fail then break; fi;
1261        j:=j+1;
1262      od;
1263
1264      if p<>fail and rulpos[p]<>fail then
1265        p:=rulpos[p];
1266        tail:=wrd{[i+Length(rules[p][1])..Length(wrd)]};
1267        wrd:=Concatenation(wrd{[1..i-1]},rules[p][2],tail);
1268
1269        if p in hastail then
1270          w:=zy{dim*(htpos[p]-1)+[1..dim]};
1271          for j in tail do
1272            w:=w*mats[j];
1273          od;
1274          v:=v+w;
1275        fi;
1276
1277        i:=Maximum(0,i-mal); # earliest which could be affected
1278      fi;
1279      i:=i+1;
1280    od;
1281    return [wrd,v];
1282  end;
1283
1284  r.fphom:=fp;
1285  r.monhom:=fm;
1286  r.colltz:=colltz;
1287
1288  # inverses of generators
1289  r.myinvers:=function(wrd,zy)
1290    return [colltz(formalinverse{Reversed(wrd)}),
1291    -r.tailforword(Concatenation(wrd,colltz(formalinverse{Reversed(wrd)})),zy)
1292      [2]];
1293  end;
1294
1295  r.pairact:=function(zy,pair)
1296  local autom,mat,i,imagemonwords,imgwrd,left,right,extim,prdout,myinvers,v;
1297
1298    autom:=InverseGeneralMapping(pair[1]);
1299
1300    # cache words for images of monoid generators
1301    if not IsBound(autom!.imagemonwords) then autom!.imagemonwords:=[]; fi;
1302    imagemonwords:=autom!.imagemonwords;
1303    imgwrd:=function(nr)
1304    local a;
1305      if not IsBound(imagemonwords[nr]) then
1306        # apply automorphism
1307        a:=PreImagesRepresentative(fm,GeneratorsOfMonoid(mon)[nr]);
1308        a:=PreImagesRepresentative(fp,a);
1309        a:=ImagesRepresentative(autom,a);
1310        a:=ImagesRepresentative(fp,a);
1311        a:=ImagesRepresentative(fm,a);
1312        a:=LetterRepAssocWord(UnderlyingElement(a));
1313        a:=colltz(a);
1314        imagemonwords[nr]:=a;
1315      fi;
1316      return imagemonwords[nr];
1317    end;
1318
1319    # normalform word and collect the tails
1320    collectail:=function(wrd)
1321    local v,i,j,s,p,mm,w,tail;
1322      v:=zerovec;
1323
1324      # collect from left
1325      i:=1;
1326      while i<=Length(wrd) do
1327
1328        # does a rule apply at position i?
1329        j:=0;
1330        s:=0;
1331        mm:=Minimum(mal,Length(wrd)-i+1);
1332        while j<mm do
1333          s:=s*max+wrd[i+j];
1334          p:=LookupDictionary(dict,s);
1335          if p<>fail and rulpos[p]<>fail then break; fi;
1336          j:=j+1;
1337        od;
1338
1339        if p<>fail and rulpos[p]<>fail then
1340          p:=rulpos[p];
1341          tail:=wrd{[i+Length(rules[p][1])..Length(wrd)]};
1342          wrd:=Concatenation(wrd{[1..i-1]},rules[p][2],tail);
1343
1344          if p in hastail then
1345            w:=zy{dim*(htpos[p]-1)+[1..dim]};
1346            for j in tail do
1347              w:=w*mats[j];
1348            od;
1349            v:=v+w;
1350          fi;
1351
1352          i:=Maximum(0,i-mal); # earliest which could be affected
1353        fi;
1354        i:=i+1;
1355      od;
1356      return [wrd,v];
1357    end;
1358
1359    prdout:=function(l)
1360    local w,v,j;
1361      w:=[];
1362      v:=zerovec;
1363      for j in l do
1364        v:=v*mapped(j[1]); # move right
1365        w:=collectail(Concatenation(w,j[1]));
1366        v:=v+w[2]+j[2];
1367        w:=w[1];
1368      od;
1369      return [w,v];
1370    end;
1371
1372    mat:=pair[2];
1373    new:=[];
1374
1375    # inverses of generators
1376    myinvers:=wrd->[colltz(formalinverse{Reversed(wrd)}),
1377      -collectail(Concatenation(wrd,colltz(formalinverse{Reversed(wrd)})))[2]];
1378
1379    # images of generators
1380    extim:=List([1..Length(mats)/2],x->imgwrd(2*x-1));
1381
1382    # collect inverses and interleave generators/inverses
1383    extim:=Concatenation(List(extim,x->[[x,zerovec],myinvers(x)]));
1384
1385    for i in [1..Length(hastail)] do
1386      # evaluate rules in images
1387      left:=prdout(extim{rules[hastail[i]][1]});
1388
1389      # invert left
1390      v:=myinvers(left[1]);
1391      left:=[v[1],v[2]-left[2]/mapped(left[1])];
1392
1393      right:=prdout(extim{rules[hastail[i]][2]});
1394
1395      # quotient left^-1*right
1396      left:=prdout(Concatenation([left],[right]));
1397      if left[1]<>[] then Error("not rule");fi;
1398
1399      Append(new,-left[2]*mat);
1400    od;
1401
1402# debugging code: Construct group and work there
1403#    new2:=[];
1404#    ext:=FpGroupCocycle(r,zy,true);
1405#    hom:=IsomorphismPermGroup(ext);
1406#    ext:=Image(hom);
1407#    epcgs:=PcgsByPcSequence(FamilyObj(One(ext)),
1408#      GeneratorsOfGroup(ext){[Length(mats)/2+1..Length(GeneratorsOfGroup(ext))]});
1409#
1410#    exta:=Concatenation(List([1..Length(mats)/2],x->[ext.(x),ext.(x)^-1]));
1411#    exte:=exta;
1412#
1413#    myprd:=function(wrd)
1414#    local i,w;
1415#      w:=One(ext);
1416#      for i in wrd do
1417#        w:=w*exte[i];
1418#      od;
1419#      return w;
1420#    end;
1421#
1422#    v:=List([1..Length(mats)/2],x->myprd(imgwrd(2*x-1))); # new generators
1423#    v:=Concatenation(List([1..Length(mats)/2],x->[v[x],v[x]^-1]));
1424#    exte:=v;
1425#
1426#    for i in [1..Length(hastail)] do
1427#      left:=rules[hastail[i]][1];
1428#      right:=rules[hastail[i]][2];
1429#      left:=myprd(left)^-1*myprd(right);
1430#      v:=-ExponentsOfPcElement(epcgs,left)*One(mo.field);
1431#
1432#      Append(new2,v*mat);
1433#    od;
1434#    if new2<>new then Error("wrong ",pair);fi;
1435#    if IsOne(pair[1]) and IsOne(pair[2]) then
1436#      if new<>zy then Error("one");else Print("onegood\n");fi;
1437#    fi;
1438
1439    return ImmutableVector(mo.field,new);
1440  end;
1441
1442  return r;
1443end);
1444
1445MatricesStabilizerOneDim:=function(field,mats)
1446local e,one,m,r,a,c,is,i;
1447    e:=[];
1448    one:=One(mats[1]);
1449    for m in mats do
1450      r:=Set(RootsOfUPol(field,CharacteristicPolynomial(m)));
1451      a:=List(r,x->NullspaceMat(m-x*one));
1452      if Length(a)=0 then return false;fi;
1453      Add(e,a);
1454    od;
1455    for c in Cartesian(e) do
1456      is:=c[1];
1457      for i in [2..Length(c)] do
1458        is:=SumIntersectionMat(is,c[i])[2];
1459      od;
1460      if Length(is)>0 then
1461        return is;
1462      fi;
1463    od;
1464  return false;
1465end;
1466
1467BindGlobal("WreathElm",function(b,l,m)
1468local n,ran,r,d,p,i,j;
1469  n:=Length(l);
1470  ran:=[1..b];
1471  r:=0;
1472  d:=[];
1473  p:=[];
1474  # base bit
1475  for i in [1..n] do
1476    for j in ran do
1477      p[r+j]:=r+j^l[i];
1478    od;
1479    Add(d,ran+r);
1480    r:=r+b;
1481  od;
1482  # permuter bit
1483  p:=PermList(p)/PermList(Concatenation(Permuted(d,m)));
1484  return p;
1485end);
1486
1487BindGlobal("PermrepSemidirectModule",function(G,module)
1488local hom,mats,mode,m,min,i,j,mo,bas,a,l,ugens,gi,r,cy,act,k,it,p;
1489  if not MTX.IsIrreducible(module) then Error("reducible");fi;
1490  p:=Size(module.field);
1491  if not IsPrime(p) then Error("must be over prime field");fi;
1492  k:=Length(module.generators);
1493  hom:=GroupHomomorphismByImagesNC(G,Group(module.generators),
1494    GeneratorsOfGroup(G),module.generators);
1495  # we allow do go immediately to normal subgroup of index up to 4.
1496  # This reduces search space
1497  it:=DescSubgroupIterator(G:skip:=LogInt(Size(G),2));
1498  repeat
1499    m:=NextIterator(it);
1500
1501    if Index(G,m)*p>p^module.dimension then
1502      Info(InfoExtReps,2,"Index reached ",Index(G,m),
1503        ", write out module action");
1504      # alternative, boring version
1505      mo:=AbelianGroup(ListWithIdenticalEntries(module.dimension,p));
1506      bas:=Pcgs(mo);
1507      act:=List(module.generators,m->
1508        GroupHomomorphismByImagesNC(mo,mo,bas,
1509          List(m,x->PcElementByExponents(bas,x)):noassert));
1510      Assert(3,ForAll(act,IsBijective));
1511      a:=Group(act);
1512      SetIsGroupOfAutomorphismsFiniteGroup(a,true);
1513      gi:=SemidirectProduct(G,GroupHomomorphismByImages(G,a,
1514        GeneratorsOfGroup(G),act),mo);
1515      r:=rec(group:=gi,
1516                  ggens:=List(GeneratorsOfGroup(G),x->
1517                    ImagesRepresentative(Embedding(gi,1),x)),
1518                  #vector:=ImagesRepresentative(Embedding(gi,2),bas[1]),
1519        basis:=List(bas,x->ImagesRepresentative(Embedding(gi,2),x)));
1520      Assert(0,Size(gi)=Size(G)*p^module.dimension);
1521      return r;
1522
1523    elif Size(Core(G,m))>1 then
1524      Info(InfoExtReps,4,"Index ",Index(G,m)," has nontrivial core");
1525    else
1526      Info(InfoExtReps,2,"Trying index ",Index(G,m));
1527
1528      mats:=List(GeneratorsOfGroup(m),
1529        x->TransposedMat(ImagesRepresentative(hom,x^-1)));
1530      a:=MatricesStabilizerOneDim(module.field,mats);
1531      if a<>false then
1532        # quotient module
1533        # basis: supplemental vector and submodule basis
1534        bas:=BaseSteinitzVectors(IdentityMat(module.dimension,GF(p)),
1535          NullspaceMat(TransposedMat(a{[1]})));
1536        bas:=Concatenation(bas.factorspace,bas.subspace);
1537
1538        # assume we have a generating set for the SDP consisting of the
1539        # complement gens, and one element of the module.
1540        cy:=CyclicGroup(IsPermGroup,p).1; # p-cycle
1541        ugens:=[];
1542
1543        # also transversal chosen in complement
1544        r:=RightTransversal(G,m);
1545        act:=ActionHomomorphism(G,r,OnRight);
1546        act:=List(GeneratorsOfGroup(G),x->ImagesRepresentative(act,x));
1547        #l:=ListWithIdenticalEntries(Index(G,m),());
1548        for gi in [1..k] do
1549          l:=[];
1550          for j in [1..Length(r)] do
1551            # matrix for Schreier gen
1552            a:=ImagesRepresentative(hom,
1553              r[j]*G.(gi)/r[PositionCanonical(r,r[j]*G.(gi))]);
1554            # how does it act on bas[1]-span, get factor
1555            a:=Int(SolutionMat(bas,bas[1]*a)[1]);
1556            # write down permutation that acts as mult by a
1557            if IsZero(a) then
1558              l[j]:=();
1559            else
1560              l[j]:=MappingPermListList([1..p],Cycle(cy^a,1));
1561            fi;
1562          od;
1563          Add(ugens,WreathElm(p,l,act[gi]) );
1564        od;
1565
1566        # module generator
1567        for j in [1..Length(r)] do
1568          #r[j]*g/r[j^g]); Note j^g=j here as kernel of permrep
1569          l[j]:=
1570            cy^Int(SolutionMat(bas,bas[1]/ImagesRepresentative(hom,r[j]))[1]);
1571        od;
1572        Add(ugens,WreathElm(p,l,()) );
1573        gi:=Group(ugens);
1574        Assert(0,Size(gi)=p^module.dimension*Size(G));
1575        if ValueOption("cheap")<>true then
1576          a:=SmallerDegreePermutationRepresentation(gi:cheap);
1577          if NrMovedPoints(Range(a))<NrMovedPoints(gi) then
1578            gi:=Image(a,gi);
1579            ugens:=List(ugens,x->ImagesRepresentative(a,x));
1580          fi;
1581        fi;
1582
1583        r:=rec(group:=gi,ggens:=ugens{[1..k]});
1584
1585        # compute basis
1586        bas:=bas{[1]};
1587        gi:=ugens{[1..k]};
1588        ugens:=ugens{[k+1]};
1589        i:=1;
1590        while Length(bas)<module.dimension do
1591          for j in [1..k] do
1592            a:=bas[i]*module.generators[j];
1593            if SolutionMat(bas,a)=fail then
1594              Add(bas,a);
1595              Add(ugens,ugens[i]^gi[j]);
1596            fi;
1597          od;
1598          i:=i+1;
1599        od;
1600        # convert back to standard basis
1601        r.basis:=List(Inverse(bas),x->LinearCombinationPcgs(ugens,x));
1602
1603        return r;
1604      fi;
1605    fi;
1606  until false;
1607end);
1608
1609BindGlobal("RandomSubgroupNotIncluding",function(g,n,time)
1610local best,i,u,v,start,cnt,runtime;
1611  runtime:=GET_TIMER_FROM_ReproducibleBehaviour();
1612  start:=runtime();
1613  best:=TrivialSubgroup(g);
1614  cnt:=0;
1615  while runtime()-start<time or cnt<100 do
1616    cnt:=cnt+1;
1617    u:=TrivialSubgroup(g);
1618    repeat
1619      v:=u;
1620      u:=ClosureGroup(u,Random(g));
1621    until IsSubset(u,n);
1622    if IndexNC(g,v)<IndexNC(g,best) then best:=v;fi;
1623  od;
1624  return best;
1625end);
1626
1627InstallGlobalFunction(FpGroupCocycle,function(arg)
1628local r,z,ogens,n,gens,str,dim,i,j,f,rels,new,quot,g,p,collect,m,e,fp,old,sim,
1629      it,hom,trysy,prime,mindeg,fps,ei,mgens,mwrd,nn,newfree,mfpi,mmats,sub,
1630      tab,tab0,evalprod,gensmrep,invsmrep,zerob,step;
1631
1632  # function to evaluate product (as integer list) in gens (and their
1633  # inverses invs) with corresponding action mats
1634  evalprod:=function(w,gens,invs,mats)
1635  local new,i,a;
1636    new:=[[],zerob];
1637    for i in w do
1638      if i>0 then
1639        collect:=r.tailforword(Concatenation(new[1],gens[i][1]),z);
1640        new:=[collect[1],collect[2]+new[2]*mats[i]+gens[i][2]];
1641      else
1642        collect:=r.tailforword(Concatenation(new[1],invs[-i][1]),z);
1643        new:=[collect[1],collect[2]+new[2]/mats[-i]+invs[-i][2]];
1644      fi;
1645    od;
1646    return new;
1647
1648  end;
1649
1650
1651
1652  r:=arg[1];
1653  z:=arg[2];
1654  ogens:=GeneratorsOfGroup(r.presentation.group);
1655
1656  str:=List(ogens,String);
1657  dim:=r.module.dimension;
1658  zerob:=ImmutableVector(r.module.field,
1659    ListWithIdenticalEntries(dim,Zero(r.module.field)));
1660  n:=Length(ogens);
1661
1662  gensmrep:=List(GeneratorsOfGroup(r.presentation.group),
1663   x->[LetterRepAssocWord(UnderlyingElement(ImagesRepresentative(r.monhom,
1664     ElementOfFpGroup(FamilyObj(One(Range(r.fphom))),x)))),zerob]);
1665  # module generators
1666  Append(gensmrep,List(IdentityMat(r.module.dimension,r.module.field),
1667    x->[[],x]));
1668  invsmrep:=List(gensmrep,x->r.myinvers(x[1],z));
1669
1670  for i in [1..dim] do
1671    Add(str,Concatenation("m",String(i)));
1672  od;
1673  f:=FreeGroup(str);
1674  gens:=GeneratorsOfGroup(f);
1675  rels:=[];
1676  for i in [1..Length(r.presentation.relators)] do
1677    new:=MappedWord(r.presentation.relators[i],ogens,gens{[1..n]});
1678    for j in [1..dim] do
1679      new:=new*gens[n+j]^Int(z[(i-1)*dim+j]);
1680    od;
1681    Add(rels,new);
1682  od;
1683  for i in [n+1..Length(gens)] do
1684    Add(rels,gens[i]^r.prime);
1685    for j in [i+1..Length(gens)] do
1686      Add(rels,Comm(gens[i],gens[j]));
1687    od;
1688  od;
1689  for i in [1..n] do
1690    for j in [1..dim] do
1691      Add(rels,gens[n+j]^gens[i]/
1692        LinearCombinationPcgs(gens{[n+1..Length(gens)]},r.module.generators[i][j]));
1693    od;
1694  od;
1695  fp:=f/rels;
1696  prime:=Size(r.module.field);
1697  SetSize(fp,Size(r.group)*prime^r.module.dimension);
1698
1699  if Length(arg)>2 and arg[3]=true then
1700    if IsZero(z) and MTX.IsIrreducible(r.module) then
1701      # make SDP directly
1702      m:=PermrepSemidirectModule(r.group,r.module:cheap);
1703      p:=m.group;
1704      # test is cheap here, thus no NC
1705      new:=GroupHomomorphismByImages(fp,p,GeneratorsOfGroup(fp),
1706        Concatenation(m.ggens,m.basis));
1707    else
1708
1709      #sim:=IsomorphismSimplifiedFpGroup(fp);
1710      sim:=IdentityMapping(fp);
1711      fps:=Image(sim,fp);
1712
1713      g:=r.group;
1714      quot:=InverseGeneralMapping(sim)
1715        *GroupHomomorphismByImages(fp,g,GeneratorsOfGroup(fp),
1716        Concatenation(GeneratorsOfGroup(g),
1717          ListWithIdenticalEntries(r.module.dimension,One(g))));
1718      hom:=GroupHomomorphismByImages(r.group,Group(r.module.generators),
1719        GeneratorsOfGroup(r.group),r.module.generators);
1720      p:=Image(quot);
1721      trysy:=Maximum(1000,50*IndexNC(p,SylowSubgroup(p,prime)));
1722      # allow to enforce test coverage
1723      if ValueOption("forcetest")=true then trysy:=2;fi;
1724
1725      mindeg:=2;
1726      if IsPermGroup(p) then
1727        mindeg:=Minimum(List(Orbits(p,MovedPoints(p)),Length));
1728      fi;
1729      while Size(p)<Size(fp) do
1730        # we allow do go immediately to normal subgroup of index up to 4.
1731        # This reduces search space
1732        it:=DescSubgroupIterator(p:skip:=LogInt(Size(p),2));
1733        repeat
1734          m:=NextIterator(it);
1735          e:=fail;
1736          if Index(p,m)>=mindeg and (hom=false or Size(m)=1 or
1737            false<>MatricesStabilizerOneDim(r.module.field,
1738              List(GeneratorsOfGroup(m),
1739              x->TransposedMat(ImagesRepresentative(hom,x))^-1))) then
1740            Info(InfoExtReps,2,"Attempt index ",Index(p,m));
1741            if hom=false then Info(InfoExtReps,2,"hom is false");fi;
1742
1743            if hom<>false and Index(p,m)>
1744              # up to index
1745              50
1746              # the rewriting seems to be sufficiently spiffy that we don't
1747              # need to worry about this more involved process.
1748              then
1749
1750              # Rewriting produces a bad presentation. Rather rebuild a new
1751              # fp group using rewriting rules, finds its abelianization,
1752              # and then lift that quotient
1753              sub:=PreImage(quot,m);
1754              tab0:=AugmentedCosetTableInWholeGroup(sub);
1755
1756              # primary generators
1757              mwrd:=List(tab0!.primaryGeneratorWords,
1758                x->ElementOfFpGroup(FamilyObj(One(fps)),x));
1759              Info(InfoExtReps,4,"mwrd=",mwrd);
1760              mgens:=List(mwrd,x->ImagesRepresentative(quot,x));
1761              nn:=Length(mgens);
1762              if IsPermGroup(m) then
1763                e:=SmallerDegreePermutationRepresentation(m);
1764                mfpi:=IsomorphismFpGroupByGenerators(Image(e,m),
1765                  List(mgens,x->ImagesRepresentative(e,x)):cheap);
1766                mfpi:=e*mfpi;
1767              else
1768                mfpi:=IsomorphismFpGroupByGenerators(m,mgens:cheap);
1769              fi;
1770              mmats:=List(mgens,x->ImagesRepresentative(hom,x));
1771
1772              i:=Concatenation(r.module.generators,
1773                ListWithIdenticalEntries(r.module.dimension,
1774                  IdentityMat(r.module.dimension,r.module.field)));
1775              e:=List(mwrd,
1776                x->evalprod(LetterRepAssocWord(UnderlyingElement(x)),
1777                gensmrep,invsmrep,i));
1778              ei:=List(e,function(x)
1779                  local i;
1780                    i:=r.myinvers(x[1],z);
1781                    return [i[1],i[2]-x[2]];
1782                  end);
1783
1784              newfree:=FreeGroup(nn+dim);
1785              gens:=GeneratorsOfGroup(newfree);
1786              rels:=[];
1787
1788              # module relations
1789              for i in [1..dim] do
1790                Add(rels,gens[nn+i]^prime);
1791                for j in [1..i-1] do
1792                  Add(rels,Comm(gens[nn+i],gens[nn+j]));
1793                od;
1794                for j in [1..Length(mgens)] do
1795                  Add(rels,gens[nn+i]^gens[j]/
1796                     LinearCombinationPcgs(gens{[nn+1..nn+dim]},mmats[j][i]));
1797                od;
1798              od;
1799
1800              #extended presentation
1801
1802              for i in RelatorsOfFpGroup(Range(mfpi)) do
1803                i:=LetterRepAssocWord(i);
1804                str:=evalprod(i,e,ei,mmats);
1805                if Length(str[1])>0 then Error("inconsistent");fi;
1806                Add(rels,AssocWordByLetterRep(FamilyObj(One(newfree)),i)
1807                         /LinearCombinationPcgs(gens{[nn+1..nn+dim]},str[2]));
1808              od;
1809              newfree:=newfree/rels; # new extension
1810              Assert(2,
1811                AbelianInvariants(newfree)=AbelianInvariants(PreImage(quot,m)));
1812              mfpi:=GroupHomomorphismByImagesNC(newfree,m,
1813                GeneratorsOfGroup(newfree),Concatenation(mgens,
1814                  ListWithIdenticalEntries(dim,One(m))));
1815
1816              # first just try a bit, and see whether this gets all (e.g. if
1817              # module is irreducible).
1818              e:=LargerQuotientBySubgroupAbelianization(mfpi,m:cheap);
1819              if e<>fail then
1820                step:=0;
1821                while step<=1 do
1822                  # Now write down the combined representation in wreath
1823
1824                  e:=DefiningQuotientHomomorphism(e);
1825                  # define map on subgroup.
1826                  tab:=CopiedAugmentedCosetTable(tab0);
1827                  tab.primaryImages:=Immutable(
1828                    List(GeneratorsOfGroup(newfree){[1..nn]},
1829                      x->ImagesRepresentative(e,x)));
1830                  TrySecondaryImages(tab);
1831
1832                  i:=GroupHomomorphismByImagesNC(sub,Range(e),mwrd,
1833                    List(GeneratorsOfGroup(newfree){[1..nn]},
1834                      x->ImagesRepresentative(e,x)):noassert);
1835                  SetCosetTableFpHom(i,tab);
1836                  e:=PreImage(i,TrivialSubgroup(Range(e)));
1837                  e:=Intersection(e,
1838                      KernelOfMultiplicativeGeneralMapping(quot));
1839
1840                  # this check is very cheap, in comparison. So just be safe
1841                  Assert(0,ForAll(RelatorsOfFpGroup(fps),
1842                    x->IsOne(MappedWord(x,FreeGeneratorsOfFpGroup(fps),
1843                      GeneratorsOfGroup(e!.quot)))));
1844
1845                  if step=0 then
1846                    i:=GroupHomomorphismByImagesNC(e!.quot,p,
1847                      GeneratorsOfGroup(e!.quot),
1848                      GeneratorsOfGroup(p));
1849                    j:=PreImage(i,m);
1850                    if AbelianInvariants(j)=AbelianInvariants(newfree) then
1851                      step:=2;
1852                      Info(InfoExtReps,2,"Small bit did good");
1853                    else
1854                      Info(InfoExtReps,2,"Need expensive version");
1855                      e:=LargerQuotientBySubgroupAbelianization(mfpi,m:
1856                        cheap:=false);
1857                      e:=Intersection(e,
1858                        KernelOfMultiplicativeGeneralMapping(quot));
1859                    fi;
1860                  fi;
1861
1862
1863                  step:=step+1;
1864                od;
1865
1866              fi;
1867
1868            else
1869              e:=LargerQuotientBySubgroupAbelianization(quot,m);
1870              if e<>fail then
1871                e:=Intersection(e,KernelOfMultiplicativeGeneralMapping(quot));
1872              fi;
1873            fi;
1874
1875
1876            if e<>fail then
1877              # can we do better degree -- greedy block reduction?
1878              nn:=e!.quot;
1879              if IsTransitive(nn,MovedPoints(nn)) then
1880                repeat
1881                  ei:=ShallowCopy(RepresentativesMinimalBlocks(nn,
1882                    MovedPoints(nn)));
1883                  SortBy(ei,x->-Length(x)); # long ones first
1884                  j:=false;
1885                  i:=1;
1886                  while i<=Length(ei) and j=false do
1887                    str:=Stabilizer(nn,ei[i],OnSets);
1888                    if Size(Core(nn,str))=1 then
1889                      j:=ei[i];
1890                    fi;
1891                    i:=i+1;
1892                  od;
1893                  if j<>false then
1894                    Info(InfoExtReps,4,"deg. improved by blocks ",Size(j));
1895                    # action on blocks
1896                    i:=ActionHomomorphism(nn,Orbit(nn,j,OnSets),OnSets);
1897                    j:=Size(nn);
1898                    # make new group to not cahe anything about old.
1899                    nn:=Group(List(GeneratorsOfGroup(nn),
1900                      x->ImagesRepresentative(i,x)),());
1901                    SetSize(nn,j);
1902                  fi;
1903                until j=false;
1904
1905              fi;
1906              if not IsIdenticalObj(nn,e!.quot) then
1907                Info(InfoExtReps,2,"Degree improved by factor ",
1908                  NrMovedPoints(e!.quot)/NrMovedPoints(nn));
1909
1910                e:=SubgroupOfWholeGroupByQuotientSubgroup(FamilyObj(fps),nn,
1911                  Stabilizer(nn,1));
1912              fi;
1913            elif e=fail and IndexNC(p,m)>trysy then
1914              trysy:=Size(fp); # never try again
1915              # can the Sylow subgroup get us something?
1916              m:=SylowSubgroup(p,prime);
1917              e:=PreImage(quot,m);
1918              Info(InfoExtReps,2,"Sylow test for index ",IndexNC(p,m));
1919              i:=IsomorphismFpGroup(e:cheap); # only one TzGo
1920              e:=EpimorphismPGroup(Range(i),prime,PClassPGroup(m)+1);
1921              e:=i*e; # map onto pgroup
1922              j:=KernelOfMultiplicativeGeneralMapping(
1923                  InverseGeneralMapping(e)*quot);
1924              i:=RandomSubgroupNotIncluding(Range(e),j,20000); # 20 seconds
1925              Info(InfoExtReps,2,"Sylow found ",IndexNC(p,m)," * ",
1926                IndexNC(Range(e),i));
1927              if IndexNC(Range(e),i)*IndexNC(p,m)<
1928                  # consider permdegree up to
1929                  100000
1930                  # as manageable
1931                then
1932                e:=PreImage(e,i);
1933                #e:=KernelOfMultiplicativeGeneralMapping(
1934                #  DefiningQuotientHomomorphism(i));
1935              else
1936                e:=fail; # not good
1937              fi;
1938            fi;
1939
1940          else Info(InfoExtReps,4,"Don't index ",Index(p,m));
1941          fi;
1942
1943        until e<>fail;
1944        i:=p;
1945
1946        quot:=DefiningQuotientHomomorphism(Intersection(e,
1947          KernelOfMultiplicativeGeneralMapping(quot)));
1948        p:=Image(quot);
1949        Info(InfoExtReps,1,"index ",Index(i,m)," increases factor by ",
1950             Size(p)/Size(i)," at degree ",NrMovedPoints(p));
1951        hom:=false; # we don't have hom cheaply any longer as group changed.
1952        # this is not an issue if module is irreducible
1953      od;
1954      quot:=sim*quot;
1955      new:=GroupHomomorphismByImages(fp,p,GeneratorsOfGroup(fp),
1956        List(GeneratorsOfGroup(fp),x->ImagesRepresentative(quot,x)));
1957
1958    fi;
1959    new:=new*SmallerDegreePermutationRepresentation(p:cheap);
1960    SetIsomorphismPermGroup(fp,new);
1961  fi;
1962
1963  return fp;
1964end);
1965
1966InstallGlobalFunction("CompatiblePairOrbitRepsGeneric",function(cp,coh)
1967local bas,ran,mats,o;
1968  if Length(coh.cohomology)=0 then return [coh.zero];fi;
1969  if Length(coh.cohomology)=1 and Size(coh.module.field)=2 then
1970    # 1-dim space over GF(2)
1971    return [coh.zero,coh.cohomology[1]];
1972  fi;
1973  # make basis
1974  bas:=Concatenation(coh.coboundaries,coh.cohomology);
1975  ran:=[Length(coh.coboundaries)+1..Length(bas)];
1976  mats:=List(GeneratorsOfGroup(cp),g->List(coh.cohomology,
1977    v->SolutionMat(bas,coh.pairact(v,g)){ran}));
1978  o:=Orbits(Group(mats),coh.module.field^Length(coh.cohomology));
1979  o:=List(o,x->x[1]*coh.cohomology);
1980  return o;
1981end);
1982