1#############################################################################
2##
3##  snksetswrsr.g
4##                                recog package
5##                                                            Maska Law
6##                                                 &  Alice C. Niemeyer
7##                                                       & 'Akos Seress
8##
9##
10##  Copyright (C) 2006-2008 by the authors.
11##                             The University of Western Australia.
12##  This file is free software, see license information at the end.
13##
14##
15##  This file provides code for recognising whether a permutation group
16##  on N points is isomorphic to G wr S_r where G is S_n acting on
17##  k-sets and N = Binomial(n,k)^r and kr > 1.
18##
19#############################################################################
20
21
22DeclareInfoClass( "InfoJellyfish" );
23SetInfoLevel( InfoJellyfish, 1 );
24
25
26#############################################################################
27##
28#F  NkrGetParameters( <N> )
29##
30##  Computes all [n,k,r] such that N = Binomial(n,k)^r. Returns those
31##  satisfying k*r > 1 and 2*r*k^2 < n (sufficient condition from
32##  BLS (Lemma 2.10) that ensures that there is a unique minimal and a
33##  unique maximal suborbit of the form which is needed). Note that
34##  this last condition can be changed in the last line of the code.
35
36RECOG.NkrGetParameters := function( N )
37
38    local maxr, list, r, k, n, newN, num;
39
40    list:=[];
41    maxr:=Gcd(List(Collected(FactorsInt(N)),x->x[2]));
42    for r in DivisorsInt(maxr) do
43        newN:=RootInt(N,r);
44        k:=0;
45        repeat
46            k:=k+1;
47            num:=RootInt(newN*Factorial(k),k);
48            for n in [num..num+k-1] do
49                if Binomial(n,k)=newN then
50                    Add(list, [n,k,r]);
51                fi;
52            od;
53        until Binomial(2*k,k)>newN;
54    od;
55    #return list;
56    return Filtered(list,x-> x[2]*x[3]>1 and x[1]>2*x[2]^2*x[3]);
57
58end;
59
60
61#############################################################################
62##
63#F  NkrSchreierTree( <alpha>, <gens>, <gensinv> )
64##
65
66RECOG.NkrSchreierTree := function ( alpha, gens, gensinv )
67
68    local sv, stack, i, j, beta, gamma;
69
70    sv := List( [1 .. LargestMovedPoint(gens)], i -> false );
71    sv[alpha] := One(gens[1]);
72
73    stack := [alpha];
74    i := 1;
75
76    while i <= Length(stack) do
77        beta := stack[i];
78        for j in [1 .. Length(gens)] do
79            gamma := beta^gens[j];
80            if sv[gamma] = false then
81                sv[gamma] := gensinv[j];
82                Add(stack, gamma);
83            fi;
84        od;
85        i := i + 1;
86    od;
87
88    return sv;
89
90end;
91
92
93#############################################################################
94##
95#F  NkrTraceSchreierTree( <beta>, <sv> )
96##
97
98RECOG.NkrTraceSchreierTree := function( beta, sv )
99
100    local o, h;
101
102    if sv[beta] = false then return fail; fi;
103
104    o := One(sv[beta]);
105    h := One(sv[beta]);
106    while sv[beta] <> o do
107        h := h * sv[beta];
108        beta := beta^sv[beta];
109    od;
110
111    return h;
112
113end;
114
115
116#############################################################################
117##
118#F  NkrOrbitsOfStabiliser( <alpha>, <sv>, <grp> )
119##
120##  The following function attempts to construct the stabiliser in <G>
121##  of a point <alpha>.
122
123RECOG.NkrOrbitsOfStabiliser := function ( alpha, sv, grp )
124
125    local  NextGeneratorStabiliser, N, gens, i, orbs;
126
127    NextGeneratorStabiliser := function( alpha, sv, grp )
128        local g;
129        g := PseudoRandom( grp );
130        return g * RECOG.NkrTraceSchreierTree( alpha^g, sv );
131    end;
132
133    N := LargestMovedPoint( grp );
134
135    gens := [];
136    for i in [1 .. Maximum( 4, Int(Log(N,2)/2) )] do
137        Add( gens, NextGeneratorStabiliser( alpha, sv, grp ) );
138    od;
139
140    orbs := OrbitsPerms( gens, Difference( [1..N], [alpha] ) );
141
142    return orbs;
143
144end;
145
146
147#############################################################################
148##
149#F  FirstJellyfish( <alpha>, <svalpha>, <gammaalpha>, <deltaalpha>,
150##                  <G>, <n>, <k>, <r> )
151##
152##  The following function attempts to construct a jellyfish, namely the
153##  set of all points whose i-th component contains the letter l, for
154##  some fixed i in [1..r] and l in [1..n].
155##  A jellyfish has size (n-1 choose k-1)*(n choose k)^(r-1).
156##
157##  The set of points used to construct the jellyfish are returned for
158##  use later as an identifier.
159
160RECOG.FirstJellyfish := function( alpha, svalpha, gammaalpha, deltaalpha,
161                            G, n, k, r )
162
163    local beta, g, deltabeta, gams, jellyfish, I, gamma, deltagamma, idjf;
164
165    # we need a point <beta> in <gammaalpha> : best to choose random
166    # <beta> in <gammaalpha> since <gammaalpha> is small
167    # wlog : assume <alpha> and <beta> differ in component 1 and the
168    # letter x satisfies x in <beta>[1] and not x in <alpha>[1]
169    beta := Random( gammaalpha );
170    g := RECOG.NkrTraceSchreierTree( beta, svalpha );
171    if g = fail then
172        Info(InfoRecog,2,"Schreier tracing failed");
173        return fail;
174    else g := g^-1;
175    fi;
176
177    # use <g> to obtain the longest <G>_<beta>-orbit
178    deltabeta := List( deltaalpha, i -> i^g );
179
180    # <gams> will be those points <gamma> such that
181    # (a)  <gamma>[i] is disjoint from <alpha>[i] for i >= 1;
182    # (b)  <gamma>[1] meet <beta>[1] is x.
183    gams := Difference( deltaalpha, deltabeta );
184
185    # the domain <D> of all points
186    jellyfish := [1..LargestMovedPoint(G)];
187
188    # for each <gamma> in <gams>, we will remove from <D> those points
189    # in deltagamma, leaving the jellyfish on x in component 1
190
191    # we can calculate some small number <I> to determine how many
192    # points in <gams> are enough, up to some acceptable probability
193    I := 4*k*r;
194    # <idjf> will record those <gamma> used to construct the jellyfish
195    idjf := [ ];
196    while I>0 and Size(jellyfish)>Binomial(n-1,k-1)*Binomial(n,k)^(r-1) do
197        I := I - 1;
198        gamma := Random( gams );
199        g := RECOG.NkrTraceSchreierTree( gamma, svalpha );
200        if g = fail then
201            Info(InfoRecog,2,"Schreier tracing failed - 2");
202            return fail;
203        else g := g^-1;
204        fi;
205        deltagamma := List( deltaalpha, i -> i^g );
206        Add( idjf, gamma );
207        jellyfish := Difference( jellyfish, deltagamma );
208    od;
209
210    if I = 0 and Size(jellyfish)>Binomial(n-1,k-1)*Binomial(n,k)^(r-1)
211        then Info(InfoRecog,2,"gammas ran out for ",[n,k,r]);
212        return fail;
213    elif Size(jellyfish) < Binomial(n-1,k-1)*Binomial(n,k)^(r-1)
214        then Info(InfoRecog,2,"too small set remained");
215        return fail;
216    fi;
217
218    return [ jellyfish, Set(idjf) ];
219
220end;
221
222
223#############################################################################
224##
225#F  GetAllJellyfish( <jellyfish>, <idjf>, <G>, <n>, <k>, <r> )
226##
227##  The following function takes a jellyfish <jellyfish> and constructs
228##  n*r jellyfish by building up the <G>-orbit of <jellyfish>, since <G>
229##  is transitive on jellyfish. The order in which the jellyfish are
230##  constructed (as images of some previous jellyfish) determines the
231##  number in [1..n*r] that it represents. This correspondence is recorded
232##  in the table <T>, the rows of which represent the points and contain
233##  the numbers of the jellyfish containing that point.
234##
235##  Each jellyfish is identified by a subset of its tentacles : the size
236##  of these identifying subsets is determined by the number of points
237##  <gamma> in Difference( <deltaalpha>, <deltabeta> ) originally needed
238##  to construct the first jellyfish. These identifying subsets are also
239##  returned for use in constructing more efficiently the permutation
240##  corresponding to a given group element.
241
242RECOG.GetAllJellyfish := function( jellyfish, idjf, G, n, k, r )
243
244    local HaveSeen, T, seen, orb, i, top, cnt, gens, jf, id, im, mm, g;
245
246    # if one jellyfish contains the identifying tentacles of another
247    # jellyfish, then the two jellyfish are equal
248
249    HaveSeen := function( jellyfish )
250        local s;
251        for s in seen do
252            if IsSubset(jellyfish,s) then return true;fi;
253        od;
254        return false;
255    end;
256
257    T := List( [1..LargestMovedPoint(G)], i -> [ ] );
258
259    seen := [ idjf ];
260
261    orb := [ rec(jf:=jellyfish,id:=idjf) ];
262
263    for i in jellyfish do
264        Add( T[i], 1 );
265    od;
266
267    top := 1;
268    cnt := 1;
269
270    gens := GeneratorsOfGroup( G );
271
272    while top > 0 do
273        jf := orb[top].jf;
274        id := orb[top].id;
275        for g in gens do
276            im := OnSets( jf, g );
277            mm := OnSets( id ,g );
278            if not HaveSeen(im) then
279                cnt := cnt + 1;
280                for i in im do
281                    Add( T[i], cnt );
282                od;
283                Add( seen, mm );
284                orb[top] := rec(jf:=im,id:=mm);
285                top := top + 1;
286            fi;
287        od;
288        top := top - 1;
289    od;
290
291    if Length( seen ) <> n*r then return fail; fi;
292
293    return [ T, seen ];
294
295end;
296
297
298#########################################################################
299##
300#F OtherGetAllJellyfish( <jellyfish>, <idjf>, <G>, <n>, <k>, <r> )
301##
302## A variant of GetAllJellyfish.
303##
304
305# FIXME: unused?
306RECOG.OtherGetAllJellyfish:=function ( jellyfish, idjf, G, n, k, r )
307
308    local  HaveSeen, T, seen, orb, i, g, gens, cnt, pnt, jf, id, img, mm;
309
310    HaveSeen := function( jellyfish )
311        local s;
312        for s in seen do
313            if IsSubset(jellyfish,s) then return true;fi;
314        od;
315        return false;
316    end;
317
318    T := List( [1..LargestMovedPoint(G)], i -> [ ] );
319    seen := [ idjf ];
320    orb := [ rec(jf:=jellyfish,id:=idjf) ];
321
322    for i in jellyfish do
323        Add( T[i], 1 );
324    od;
325
326    gens := GeneratorsOfGroup( G );
327    cnt := 1;
328
329    for pnt in orb do
330        jf := pnt.jf;
331        id := pnt.id;
332        for g in gens do
333            img := OnSets( jf, g );
334            mm := OnSets( id ,g );
335            if not HaveSeen(img) then
336                cnt := cnt + 1;
337                for i in img do
338                    Add( T[i], cnt );
339                od;
340                Add( orb, rec(jf:=img,id:=mm) );
341                Add( seen, mm );
342            fi;
343        od;
344    od;
345
346    if Length( seen ) <> n*r then return fail; fi;
347
348    return [ T, seen ];
349
350end;
351
352
353#########################################################################
354##
355#F AllJellyfish( <G> )
356##
357## This is the main function.
358##
359## Note re: the shortest orbit and the longest orbit of <G>_<alpha> :
360## (i) the shortest orbit consists of those points in Omega for which
361## the k-sets in each of r-1 components are identical to the k-sets of
362## <alpha> in those r-1 components and for which the k-set in the
363## remaining component meets the k-set of <alpha> in that component in
364## k-1 letters. It has size (k choose k-1)*(n-k)*r. Notation: gammaalpha.
365## (ii) the longest orbit consists of those points in Omega for which
366## the k-set in each component is disjoint from the k-set in the same
367## component of <alpha>. It has size (n-k choose k)^r. Notation: deltaalpha.
368
369RECOG.AllJellyfish := function( G )
370
371    local N, D, gens, gensinv, i, g, alpha, svalpha, alphaorbs, sizes,
372          s, l, params, p, n, k, r, jellyfishone, alljellyfish;
373
374    N := LargestMovedPoint(G);
375    D := [ 1 .. N ];
376
377    gens := ShallowCopy( GeneratorsOfGroup(G) );
378    gensinv := List( gens, g -> g^-1 );
379    for i in [1 .. LogInt(N,10)] do
380        g := PseudoRandom(G);
381        Add( gens, g );
382        Add( gensinv, g^-1 );
383    od;
384
385    alpha := 1;  # alpha:=Random( [1..LargestMovedPoint(G)] );
386    svalpha := RECOG.NkrSchreierTree( alpha, gens, gensinv );
387    alphaorbs := RECOG.NkrOrbitsOfStabiliser( alpha, svalpha, G );
388
389    sizes := List( alphaorbs, i -> Length(i) );
390    s := Minimum( sizes );
391    l := Maximum( sizes );
392    if Length(Filtered(sizes,i->i=s)) > 1 or
393       Length(Filtered(sizes,i->i=l)) > 1 then
394        Info(InfoRecog,2,"no unique max or min sized G_alpha orbit");
395        return fail;
396    fi;
397    s := alphaorbs[ Position(sizes,s) ]; # shortest <G>_<alpha>-orbit
398    l := alphaorbs[ Position(sizes,l) ]; # longest <G>_<alpha>-orbit
399
400    params := RECOG.NkrGetParameters( N );
401
402    for p in params do
403        n := p[1]; k := p[2]; r := p[3];
404        if Size(s) = Binomial(k,k-1) * (n-k) * r and
405           Size(l) = Binomial(n-k,k)^r then
406            jellyfishone := RECOG.FirstJellyfish(alpha,svalpha,s,l,G,n,k,r);
407            if jellyfishone <> fail then
408                Info(InfoRecog,2,"found jellyfish for ",p);
409                alljellyfish := RECOG.GetAllJellyfish( jellyfishone[1],
410                                 jellyfishone[2], G,n,k,r );
411                if alljellyfish <> fail then
412                    return alljellyfish;
413                fi;
414            fi;
415        else
416            Info(InfoRecog,2,
417        "couldn't find orbits of correct size for n=",n," k=",k," r=",r);
418        fi;
419    od;
420
421    Info(InfoRecog,3,"getting jellyfish failed");
422    return fail;
423
424end;
425
426
427########################################################################
428##
429#F FindImageJellyfish( <g>, <T>, <seen> )
430##
431## This function returns the permutation of [1..n*r] corresponding to <g>.
432##
433## A jellyfish is identified by a subset of its tentacles : the size of
434## these identifying subsets is determined by the number of points
435## <gamma> in Difference( <deltaalpha>, <deltabeta> ) originally needed
436## to construct the first jellyfish. The intersection of the images
437## under <g> of the points in an identifying subset determines the image
438## of the letter of the jellyfish.
439
440RECOG.FindImageJellyfish := function( g, T, seen )
441
442    local nr, h, i, gams, tmp;
443
444    nr := Length( seen );
445    h := List( [1..nr], i -> 0 );
446
447    for i in [1..nr] do
448        gams := OnSets( seen[i], g );
449        if ForAny( gams, t -> not IsBound(T[t])) then
450            return fail;
451        fi;
452        tmp := Intersection( List( gams, t -> T[t] ) );
453        if Length(tmp) = 0 then
454            return fail;
455        fi;
456        h[i] := tmp[1];
457    od;
458
459    if Size(Set(h)) <> nr then
460        return fail;
461    else
462        return PermList(h);
463    fi;
464
465end;
466
467########################################################################
468##
469#F  FindHomMethodsPerm.SnkSetswrSr
470##
471
472RECOG.JellyHomFunc := function(data,el)
473  return RECOG.FindImageJellyfish(el,data.T,data.seen);
474end;
475
476FindHomMethodsPerm.SnkSetswrSr :=
477  function(ri,grp)
478    local res,T,seen,imgens,hom;
479    if not IsPermGroup(grp) then
480        return NeverApplicable;
481    fi;
482    if not IsPrimitive(grp) then
483        return NeverApplicable;
484    fi;
485    RECOG.SetPseudoRandomStamp(grp,"Jellyfish");
486    res := RECOG.AllJellyfish(grp);
487    if res = fail then
488        return TemporaryFailure;
489    fi;
490    ri!.jellyfishinfo := res;
491    T := res[1];
492    seen := res[2];
493    # Build the homomorphism:
494    imgens := List(GeneratorsOfGroup(grp),
495                   x->RECOG.FindImageJellyfish(x,T,seen));
496    hom := GroupHomByFuncWithData(grp,Group(imgens),RECOG.JellyHomFunc,
497                                  rec(T := T,seen := seen));
498    SetHomom(ri,hom);
499
500    return Success;
501  end;
502
503##
504##  This program is free software: you can redistribute it and/or modify
505##  it under the terms of the GNU General Public License as published by
506##  the Free Software Foundation, either version 3 of the License, or
507##  (at your option) any later version.
508##
509##  This program is distributed in the hope that it will be useful,
510##  but WITHOUT ANY WARRANTY; without even the implied warranty of
511##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
512##  GNU General Public License for more details.
513##
514##  You should have received a copy of the GNU General Public License
515##  along with this program.  If not, see <http://www.gnu.org/licenses/>.
516##
517
518