#############################################################################
##
## This file is part of recog, a package for the GAP computer algebra system
## which provides a collection of methods for the constructive recognition
## of groups.
##
## This files's authors include Stephen Howe, Maska Law, Max Neunhöffer,
## Ákos Seress.
##
## Copyright of recog belongs to its developers whose names are too numerous
## to list here. Please refer to the COPYRIGHT file for details.
##
## SPDX-License-Identifier: GPL-3.0-or-later
##
##
## This file contains an implementation of an algorithm for recognising
## An or Sn acting as matrix groups on their fully deleted permutation
## modules. The particular algorithm is described in the paper [BLGN+05].
##
#############################################################################
## General utility functions
RECOG.LogRat := function(rat, base)
# LogRat(<rat>, <base>) is an upper bound
# on log <rat> to base <base>.
return (LogInt(NumeratorRat(rat), base) + 1)
- LogInt(DenominatorRat(rat), base);
end;
RECOG.IterateWithRandomElements := function(func, times, group)
# IterateWithRandomElements(<func>, <times>, <group>) randomly
# selects up to <times> random elements (from <group>) to find
# an element such that <func> (with input being this random element)
# does not return 'fail'. The function <func> must take a group
# element as its only argument.
local itr, # number of random elements selected
re, # a random element
ret; # the result of func(re)
itr := 1;
while itr <= times do
# select random element
re := PseudoRandom(group);
# call func with random element
ret := func(re);
# if function successful
if ret <> fail then
return ret;
fi;
# else try again
itr := itr + 1;
od;
return fail;
end;
## Functions for efficient matrix operations
RECOG.Commute := function(g, h);
# Commute(g,h) returns true if and only if
# gh = hg
ConvertToMatrixRep(g);
ConvertToMatrixRep(h);
return g*h = h*g;
end;
RECOG.Conj := function(A, B)
# Conj(A, B) returns the conjugation of
# A by B (where A and B are two matrices)
# That is Conj(A, B) = B^-1*A*B
ConvertToMatrixRep(B);
return A^B;
end;
RECOG.ConjInv := function(A, B)
# ConjInv(A, B) returns the conjugation of
# A by B^-1 (where A and B are two matrices)
# That is ConjInv(A, B) = B*A*B^-1
ConvertToMatrixRep(B);
return B*A*B^-1;
end;
## Finding a 3-cycle or double transposition
RECOG.LinearFactors := function(poly, field)
# LinearFactors(<poly>, <field>) returns a list with two
# elements when given a univariate polynomial <poly> with
# coefficients from the finite field <field>. The first
# element is a list of field elements [c1,...,cn] where
# t - ci divides <poly> for i = 1,...,n. The second element
# is a function from the multiplicative group of <field> to
# the integers such that: for every non-zero field element c,
# the image of c is the multiplicity of t - c in <poly>.
local linearFactors, # [c1,...,cn] s.t. t - ci is a factor
ix, # list index
z, # primitive root of field, i.e., z = Z(q)
multList, # a list of multiplicities where t - Z(q)^(i-1)
# has multiplicity multList[i]
Multiplicity; # function mapping non-zero field elem. c to
# mult. of t - c
# find the linear factors of poly (as polynomials),
# so linearFactors = [[ai*t + bi, mi]] where ai, bi are
# field elements and mi are the multiplicities
linearFactors := Filtered(Collected(Factors(PolynomialRing(field), poly)),
f -> DegreeOfLaurentPolynomial(f[1]) = 1);
# put linearFactors into form [[ci, mi]] where mi is
# multiplicity of t - ci
Apply(linearFactors,
f -> [-Iterated(CoefficientsOfUnivariatePolynomial(f[1]), \/),
f[2]]);
# create the function Multiplicity
z := PrimitiveRoot(field); # find Z(q) where q is size of field
# create list s.t. the ith entry is the multiplicity
# of t - Z(q)^(i-1) in poly
multList := [];
for ix in [1..Length(linearFactors)] do
multList[LogFFE(linearFactors[ix][1], z) + 1] := linearFactors[ix][2];
od;
Multiplicity := function(ffe)
local ffeIndex;
# calc index for ffe
ffeIndex := LogFFE(ffe, z) + 1;
# if t - Z(q)^i is not a factor of poly
# then multList[i] is not bound
if IsBound(multList[ffeIndex]) then
return multList[ffeIndex];
else
return 0;
fi;
end;
return [List(linearFactors, f -> f[1]), Multiplicity];
end;
RECOG.IsPreDoubleTransposition := function(charPoly, t, linearFactors,
fdpmDim, nAmbig, field)
# IsPreDoubleTransposition(<charPoly>, <t>, <linearFactors>,
# <fdpmDim>, <nAmbig>, <field>)
# determines whether some scalar multiple of the matrix described
# by <charPoly> (with indeterminate <t>) represents a
# pre-double-transposition. If so the scalar multiple is returned,
# otherwise 'false' is returned.
# <linearFactors> - describes the linear factors of <charPoly>
# <fdpmDim> - the dimension of the matrices described
# <nAmbig> - true when there are two possible values for n
# <field> - the field of <charPoly>
# This function works for any field with odd characteristic.
# Note: requires n >= 8 or n = 7 and p <> 7, <=> fdpmDim >= 6
local b, c, # field elements
tpList, # list of elements c s.t. t^2 - c^2 | charPoly
r; # a prime integer
if Length(linearFactors[1]) < 2 then
# we must have t + c and t - c as linear factors
# for some c
return false;
else
tpList := [];
# we expect tpList to only have two elements
# at the end of this loop
for c in linearFactors[1] do
# we add to tpList in (+/-) pairs
# so check c is not -d for some d in tpList
if c in tpList then
continue;
elif linearFactors[2](-c) > 0 then
# t^2 - c^2 divides charPoly
if not IsEmpty(tpList) then
# tpList would have more than 2 elements
return false;
else
Add(tpList, c);
Add(tpList, -c);
fi;
fi;
od;
fi;
# determine whether tpList is of the form { c, -c } where
# t + c has multiplicity 2
if IsEmpty(tpList) then
return false;
elif linearFactors[2](tpList[1]) = 2 then
# (t - tpList[1]) has multiplicity 2
c := -tpList[1];
elif linearFactors[2](tpList[2]) = 2 then
# (t - tpList[2]) has multiplicity 2
c := -tpList[2];
else
return false;
fi;
# determine the scalar b associated with the matrix
# whose characteristic polynomial we have
if linearFactors[2](c) <> 2 then
b := c;
# some extra cases: Assuming s = 2 and 2 <= r < n/2 in proof,
# so there may be a (unique) odd prime r >= n/2
if IsPrime(fdpmDim - 2)
and charPoly = (t^2 - b^2)*(t^(fdpmDim - 2) + b^(fdpmDim - 2)) then
# corresponds to permutation with cycle type 1^1 2^1 (n-3)^1
return false;
elif nAmbig
and fdpmDim - 4 <> 2 and IsPrime(fdpmDim - 4) # fdpmDim-4 is odd prime
and charPoly = (t - b)^3*(t+b)*(t^(fdpmDim - 4) + b^(fdpmDim - 4)) then
# corresponds to permutation with cycle type 2^3 (n-6)^1
return false;
fi;
elif fdpmDim mod 2 = 0 then # also mult (t - c) = 2
# get coeffs of t^(fdpmDim - 1), t^(fdpmDim - 2), and t^(fdpmDim - 3)
# in charPoly
c := CoefficientsOfUnivariatePolynomial(charPoly){
[fdpmDim - 3 .. fdpmDim - 1] + 1};
# Check case 3 (ii) and first part of 3 (iii) together (they are
# identical in terms of fdpmDim = n - delta). That is,
# check t^2-c[3]^2 | charPoly and
# charPoly = (t-c[3])(t+c[3])^2(t^(fdpmDim-3)-c[3]^(fdpmDim - 3))
if c[3] in tpList
and charPoly = (t - c[3]) * (t + c[3])^2 *
(t^(fdpmDim - 3) - c[3]^(fdpmDim - 3)) then
# since we know the form of the characteristic polynomial we don't
# need to check t^r + b^r does not divide charPoly for r etc.
return c[3];
elif nAmbig then # do the rest of 3 (iii) only if delta could be 2
# expect coeffs of t^(fdpmDim-1) and t^(fdpmDim - 2) are
# 2b and b^2 respectively with b in tpList
if c[3]/2 in tpList and c[2] = (c[3]/2)^2 then
# if c[1] <> 0 then we must have cycle type 2^2 3^1 (n-7)^1
# if c[1] = 0 then we just set b
b := c[3]/2;
if not IsZero(c[1])
and charPoly = (t + b)^2 * (t^3 - b^3) *
(t^(fdpmDim - 5) - b^(fdpmDim - 5)) then
# we can return b since we have check the form
# of charPoly
return b;
fi;
else
return false; # c[3]/2 is not in tpList
fi;
else
return false; # n not ambiguous
fi;
else
return false; # no special cases for this value of fdpmDim
fi; # we have now either returned or found a possible scalar b
# check there is no prime r = 2 or 5 <= r < n/2 with
# t^r + b^r dividing charPoly
r := 2;
while r < QuoInt(fdpmDim + 2, 2) do
if r <> Characteristic(field)
and IsZero(QuotientRemainder(charPoly, t^r + b^r)[2]) then
return false;
fi;
r := NextPrimeInt(r);
od;
return b;
end;
RECOG.IsPre3Cycle :=
function(charPoly, t, linearFactors, fdpmDim, nAmbig, field)
# IsPre3Cycle(<charPoly>, <t>, <linearFactors>, <fdpmDim>, <nAmbig>, <field>)
# determines whether some scalar multiple of the matrix described
# by <charPoly> (with indeterminate <t>) represents a
# pre-3-cycle. If so the scalar multiple is returned,
# otherwise 'false' is returned.
# <linearFactors> - describes the linear factors of <charPoly>
# <fdpmDim> - the dimension of the matrices described
# <nAmbig> - true when there are two possible values for n
# <field> - the field of <charPoly>
# Note: requires n >= 5 or n = 5 and p <> 5, <=> fdpmDim >= 4
local tcList, # list of field elems c s.t. t^3 - c^3 | charPoly
lMult, # a list of multiplicities
r, # a prime
b, c, # field elements
y; # field element of order 3
if Size(field) mod 3 = 2 then
# find all c s.t. t^3 - c^3 | charPoly
tcList := Filtered(linearFactors[1],
c -> IsZero(QuotientRemainder(charPoly,
t^3 - c^3)[2]));
if Length(tcList) = 1 then
# check whether tcList = {c} where t^2+ct+c^2 has mult 1
if IsZero(QuotientRemainder(charPoly,
(t^2 + tcList[1]*t + tcList[1]^2)^2)[2]) then
# t^2+ct+c^2 has mult > 1
return false;
else
# set scalar
b := tcList[1];
fi;
elif IsEmpty(tcList) and nAmbig and (fdpmDim + 2) mod 3 <> 0 then
if Characteristic(field) <> 2 then # if char of field is odd
# coeff of t^(fdpmDim - 1) in charPoly is 2b
c := CoefficientsOfUnivariatePolynomial(charPoly)[(fdpmDim - 1)+1]/2;
else
# coeff of t^(fdpmDim - 2) in charPoly is b^2
c := CoefficientsOfUnivariatePolynomial(charPoly)[(fdpmDim - 2)+1];
c := SquareRoots(field, c)[1];
fi;
if charPoly = (t^2 + c*t + c^2)*(t^(fdpmDim-1) - c^(fdpmDim-1))
/ (t - c) then
# return the scalar since we know the form of charPoly
return c;
else
return false;
fi;
else # Length(tcList) >= 2
return false;
fi;
else # q mod 3 = 1 since we assume q mod 3 <> 0
# calc field element of order 3
y := PrimitiveRoot(field)^((Size(field)-1)/3);
# find all elements c such that t-c | charPoly and t-cy^i | charPoly
# for i = 1 or 2
tcList := Filtered(linearFactors[1],
c -> linearFactors[2](c*y) > 0
or linearFactors[2](c*y^2) > 0);
if Length(tcList) = 2
and Order(tcList[1]*tcList[2]^-1) = 3
and linearFactors[2](tcList[1]) = 1
and linearFactors[2](tcList[2]) = 1 then
# set the scalar
b := tcList[1]^2*tcList[2]^-1;
elif Length(tcList) = 3 and ForAll(tcList, c -> y*c in tcList) then
# check whether y*tcList = tcList (setwise)
# find elements of tcList that have mult > 1
lMult := Filtered(tcList, c -> linearFactors[2](c) > 1);
if Length(lMult) = 1 then
# set the scalar
b := lMult[1];
elif IsEmpty(lMult) then
# get coeffs of t^(fdpmDim-1), t, and t^0
c := CoefficientsOfUnivariatePolynomial(charPoly){[1,2,(fdpmDim-1)+1]};
# check for pre-3-cycles with cycle type 3^1 (n-3)^1, if delta = 1,
# or cycle type 1^1 3^1 (n-4)^1, if delta = 2. The tests are
# identical when written in terms of n - delta.
if (fdpmDim + 1) mod 3 <> 0
and c[3] in tcList
and charPoly = (t^2 + c[3]*t + c[3]^2)
*(t^(fdpmDim - 2) - c[3]^(fdpmDim - 2)) then
# return the scalar since we know the form of charPoly
return c[3];
elif nAmbig # this case only occurs when delta = 2
and Characteristic(field) >= 5
and c[3]/2 in tcList
and c[1] = -(c[3]/2)^fdpmDim # check coeff of t^0
and c[2] = -2*(c[3]/2)^(fdpmDim-1) then # check coeff of t
# set the scalar
b := c[3]/2;
else
return false;
fi;
else # Length(lMult) > 1
return false;
fi;
else # Length(tcList) <> 2 or 3
return false;
fi;
fi; # end of q mod 3 = 1 case
# check there are no primes r (<= n/3, r <> p) such that
# t^2r + b^r*t^r + b^2r | charPoly
r := 2;
while r <= QuoInt(fdpmDim + 2, 3) do
if r <> Characteristic(field)
and IsZero(QuotientRemainder(charPoly,
t^(2*r) + b^r*t^r + b^(2*r))[2]) then
return false;
fi;
r := NextPrimeInt(r);
od;
return b;
end;
RECOG.Construct3Cycle := function(mat, fdpm)
# Construct3Cycle(<mat>, <fdpm>) attempts
# to construct a matrix representing a 3-cycle from
# <mat>. If it is successful the new matrix is
# returned in a record, otherwise fail is returned.
# The record has two components; the first is
# 'matrix' and stores the matrix, the second is 'order'
# and stores the order of the matrix i.e. 3.
local charPoly, # the characteristic polynomial of mat
t, # the indeterminate of charPoly
linearFactors, # the linear factors of charPoly
b, # a field element
order; # order of the matrix
# calculate the characteristic polynomial and its linear factors
charPoly := CharacteristicPolynomial(mat);
t := IndeterminateOfLaurentPolynomial(charPoly);
# calculate the linear factors of charPoly
linearFactors := RECOG.LinearFactors(charPoly, fdpm.field);
b := RECOG.IsPre3Cycle(charPoly, t, linearFactors,
fdpm.dim, fdpm.nAmbig, fdpm.field);
if b <> false then
mat := b^-1*mat;
# check the order
order := Order(mat);
if order < (fdpm.dim + 2)^(18*RECOG.LogRat(fdpm.dim + 2, 2))
and order mod 3 = 0 then
ConvertToMatrixRep(mat);
mat := mat^(order / 3);
return rec( matrix := mat, order := 3 );
fi;
fi;
return fail;
end;
RECOG.ConstructDoubleTransposition := function(mat, fdpm)
# ConstructDoubleTransposition(<mat>, <fdpm>)
# attempts to construct a matrix representing a double
# transposition from <mat>. If it is successful the new
# matrix is returned in a record, otherwise 'fail' is
# returned. The record has two components; the first is
# 'matrix' and stores the matrix, the second is 'order'
# and stores the order of the matrix i.e. 2.
local charPoly, # the characteristic polynomial of mat
t, # the indeterminate of charPoly
linearFactors, # the linear factors of charPoly
b, # a field element
order; # order of the matrix
# calculate the characteristic polynomial and its linear factors
charPoly := CharacteristicPolynomial(mat);
t := IndeterminateOfLaurentPolynomial(charPoly);
# calculate the linear factors of charPoly
linearFactors := RECOG.LinearFactors(charPoly, fdpm.field);
b := RECOG.IsPreDoubleTransposition(charPoly, t, linearFactors,
fdpm.dim, fdpm.nAmbig, fdpm.field);
if b <> false then # mat is a pre-double-transposition
mat := b^-1*mat;
# check the order
order := Order(mat);
if order < (fdpm.dim + 2)^(18*RECOG.LogRat(fdpm.dim + 2, 2))
and order mod 2 = 0 then
ConvertToMatrixRep(mat);
mat := mat^(order / 2);
return rec( matrix := mat, order := 2 );
fi;
fi;
return fail;
end;
## Functions for finding the first element of our linked basis
RECOG.DoubleAndShrink := function(g, group, n)
# DoubleAndShrink(<g>, <group>, <n>) returns a conjugate
# h of <g> such that <g> and h have only 1 moved
# point in common where <group> is some matrix representation of
# a finite symmetric group of order <n>, and <g> moves between
# 3 and sqrt(<n>) points.
local gs, # list of generated elements
rs, # list of random elements, their inverses,
# and the conjugations
ri, # the current random element
gi, # the current generated element
temp, # temp storage for some element
gconjs, # gs[j-1]^s
gconjrs, # gs[j-1]^(r[j-1]s) respectively
i, j, # counters
s, # conjugating element, g^s is returned
maxItr; # maximum number of random element selections
# initialise variables
maxItr := RECOG.LogRat(n, 2);
gs := [g];
gi := g;
# get random element and calculate its inverse
# and gi^ri
ri := PseudoRandom(group);
ConvertToMatrixRep(ri);
temp := ri^-1;
rs := [[ri, temp, temp*gi*ri]];
i := 1;
# find elements that do not commute
while RECOG.Commute(gi, rs[i][3]) do
if i >= maxItr then
return fail;
fi;
# generate next element
gi := gi * rs[i][3];
Add(gs, gi);
# get next random element
ri := PseudoRandom(group);
temp := ri^-1;
Add(rs, [ri, temp, temp * gi * ri]);
i := i + 1;
od;
# go back through the generated elements to
# get a conjugating element s such that g
# and g^s do not commute
s := rs[i][1];
for j in [i,i-1..2] do
# calc gs[j-1]^s
gconjs := RECOG.Conj(gs[j-1], s);
# determine new s
if RECOG.Commute(gs[j-1], gconjs) then
# calc gs[j-1]^(rs[j-1]*s)
gconjrs := RECOG.Conj(rs[j-1][3], s);
if not RECOG.Commute(gs[j-1], gconjrs) then
s := rs[j-1][1] * s;
elif not RECOG.Commute(rs[j-1][3], gconjs) then
s := s * rs[j-1][2];
else
# s := s^(rs[j-1][1]^-1)
s := rs[j-1][1] * s * rs[j-1][2];
fi;
# else
# s stays the same
fi;
od;
return RECOG.Conj(g, s);
end;
RECOG.FixedPointSubspace := function(mat, k)
# FixedPointSubspace(<mat>, <k>) for a matrix <mat> with
# order <k> returns [F, W] where
# - F is a list of vectors that form a basis
# for the fixed point subspace of <mat>, and
# - W is a list of vectors that form a basis for
# the complement W(<mat>) of F such that for all w
# in W(<mat>), w+w<mat>+w<mat>^2+...+w<mat>^(k-1) = 0
# This function uses the row echelon form of a matrix in order
# to compute the relevent bases. Requires k mod p <> 0 where
# p is the characteristic of the field.
local moveLT, trans;
# linear transformation mapping elements of the
# vector space to elements of W(<mat>), takes
# a matrix as its argument
moveLT := function(v)
local sum, i;
# calculate (((v*mat + v)*mat + v)*mat + v)...)*mat + v
sum := v;
for i in [1..k-1] do
sum := sum * mat + v;
od;
return k*v - sum;
end;
# SemiEchelonMatTransformation(m)
# returns rec ( heads, vectors, coeffs, relations )
# where the rank of m is the length of coeffs
# and coeffs*m are the row vectors of m in row echelon form
# In particular [ coeffs ] is the echelonizing matrix
# [ relations ]
# v is a fixed point of mat iff v is in the null space of Xmat
# compute the transformation taking 'mat - One(mat)' to row echelon form
trans := SemiEchelonMatTransformation(mat - One(mat));
return [trans.relations, moveLT(trans.coeffs)];
end;
RECOG.SubspaceIntersection := function(basisA, basisB)
# SubspaceIntersection(basisA, basisB) returns a vector in the intersection
# of the subspaces <basisA[1], basisA[2]> and <basisB[1], basisB[2]>
# We assume that basisA and basisB are both lists of two linearly
# independent vectors.
local mat, # matrix formed from rows of basisA and basisB
row, # a row vector in the intersection
nullv; # basis for nullspace of mat
# check we are given lists of length 2
if Length(basisA) <> 2 or Length(basisB) <> 2 then
# this function only works for intersection of subspaces
# of dimension 2
return fail;
fi;
# create the matrix whose rows are basisA[1], basisA[2],
# basisB[1], and basisB[2]
mat := [ShallowCopy(basisA[1]),
ShallowCopy(basisA[2]),
ShallowCopy(basisB[1]),
ShallowCopy(basisB[2])];
# find the null space of the matrix
nullv := SemiEchelonMatTransformation(mat).relations;
if Length(nullv) <> 1 or ( IsZero(nullv[1][1])
and IsZero(nullv[1][2])) then
# the two subspaces don't intersect in a one dimensional
# subspace
return fail;
else
# calculate vector in subspace intersection
ConvertToVectorRep(mat[1]);
MultRowVector(mat[1], nullv[1][1]);
AddRowVector(mat[1], mat[2], nullv[1][2]);
return mat[1];
fi;
end;
RECOG.FindBasisElement := function(g, group, epsilon, fdpm)
# FindBasisElement(<g>, <group>, <epsilon>, <fdpm>) returns,
# with probability 1 - <epsilon>, a vector representing a basis
# element of the fully deleted permutation module, that is,
# a vector representing b(e_i - e_j) + (W \cap E) for some
# non-zero field element b and integers i and j. A matrix
# representing a 3-cycle or a double transposition is passed
# in the record <g>. The record has two components; the first is
# 'matrix' and stores the matrix, the second is 'order' and stores
# the order of the matrix (from which we can determine whether
# the permutation type the matrix represents).
local itr, h, gh, temp;
itr := 1;
# Check the bounds
while itr < RECOG.LogRat(epsilon^-1, 2)*7 do # 1/log(10/9) < 7
# get random conjugate
h := RECOG.DoubleAndShrink(g.matrix, group, fdpm.dim + 2);
# check if it has exactly one moved point in common
if h <> fail then
gh := g.matrix * h;
if ( (g.order = 2 and Order(gh) = 6)
or (g.order = 3 and Order(gh) = 5)) then
# return intersection of the totally moved subspaces
ConvertToMatrixRep(h);
temp := h^-1;
return RECOG.SubspaceIntersection(
RECOG.FixedPointSubspace(g.matrix, g.order)[2],
RECOG.FixedPointSubspace(temp * gh, g.order)[2]);
fi;
fi;
itr := itr + 1;
od;
return fail;
end;
## Functions for finding a linked sequence of vectors
RECOG.Power := function(xs, k)
# Power(<xs>, <k>) computes <xs>[1]^<k> in
# log k time where <xs> is the list
# a list [x, x^2, x^4, ...., x^(2^m)] for some
# matrix x. If k >= 2^(m+1) then the list will
# be updated
local q, r, prod, xn;
q := k;
xn := 1; # list index x[xn] = x^(2^(xn-1))
prod := One(xs[1]); # partial product which will finish as x^k
while q > 0 do
# extend the list if we have to
if xn > Length(xs) then # we have xn = Length(xs) + 1
ConvertToMatrixRep(xs[xn - 1]);
Add(xs, xs[xn - 1]^2);
fi;
# if this bit in the binary expansion of k is 1
if RemInt(q, 2) = 1 then
# update the product
prod := prod * xs[xn];
fi;
xn := xn + 1;
# move to the next bit
q := QuoInt(q, 2);
od;
return prod;
end;
RECOG.VectorImages := function(v, xs)
# VectorImages(<v>, <xs>) returns the list
# [<v>, <v>*x, ..., <v>*x^(2^(k+1)-1)]
# where xs = [x,x^2,x^4,...,x^(2^k)].
# The lists are computed in O(logm) time (wrt matrix mult).
local k, # counter
vxs; # the list [v,vx,...,vx^2logm]
# calc v,vx,...,vx^2logm
vxs := [v];
for k in [1..Length(xs)] do
Append(vxs, vxs * xs[k]);
od;
return vxs;
end;
RECOG.VectorImagesUnder := function(v, h, k)
# VectorImagesUnder(<v>,<h>, <k>) returns
# the list [<v>, <v>*h, ..., <v>*<h>^(2^(LogInt(k,2)+1)-1)]
local hs, p;
hs := [h];
for p in [1..LogInt(k, 2)] do
ConvertToMatrixRep(hs[p]);
Add(hs, hs[p]^2);
od;
return RECOG.VectorImages(v, hs);
end;
RECOG.MoreBasisVectors := function(x, g, v, fdpm)
# MoreBasisVectors(<x>, <g>, <v>, <fdpm>), where
# <g> is a record and <g>.matrix represents a 3-cycle
# or double transposition, <g>.order is
# the order of <g>.matrix, and <x> is a random group element;
# returns a set of linked basis vectors of prime length r where
# 0.6n+0.4 < r < 0.95n - 0.85
local i, # the number of vectors (obtained from v) we consider
vs, # the list [vg^i | i = 0, 1, 2, if g reps. a 3 cycle
# i = 0, if g reps. a dt ]
xs, # the list [x,x^2,x^4,...,x^(2^k)]
# where 2^k <= n < 2^(k+1)
vBases, # a list of Basis objects for the vector spaces
# spanned by the vectors in vs
vImgs, # images of vectors in vs under x, ..., x^(n-1)
rs, # [r(vg^i,x) | i = 0, 1, 2, if g reps. a 3 cycle
# i = 0, if g reps. a dt ]
R, # returns r(v, x)
lList, # { l : 1 <= l < r, vx^l <> vx^lg }
mat, # temporary matrix for calc. vImgsG
vImgsG, # [ vg, vxg, vx^2g, ..., vx^(n-1)g ]
u, # a vector
d, # u = dv or u = dvg
h, # a matrix rep. a (large) prime cycle
k; # loop counter
R := function(vBasis, vxs)
# R(vBasis, vxs) returns r(v, x) where
# vBasis is a Basis object for the vector space
# spanned by v and vxs is the list of images of v
# under x, ..., x^(n-1)
local r;
# find least +ve integer k such that vx^k in
# the vector space spanned by v
# Note: opened up the bounds a bit so it works for
# n < 13
r := NextPrimeInt(Int(6*(fdpm.dim + 1)/10 + 4/10) - 1);
# we can just check primes since, if r = r(v, x), and
# vx^s \in <v> then r divides s. If we are given a
# scalar multiple of the identity matrix then
# we will fail when considering the list
# { l : 1 <= l < r, vx^l <> vx^lg } as v <> vg
while r < (19*(fdpm.dim + 2) - 17)/20 do
# check if vx^r = dv for some field element d
if Coefficients(vBasis, vxs[r+1]) <> fail then
# check r is prime and in the right range
return r;
fi;
r := NextPrimeInt(r);
od;
return fail;
end;
# Determine how many vectors (all obtained from v) that
# we will consider. This number depends on whether we
# have a 3-cycle or double transposition
if g.order = 3 then
i := 3;
else
i := 1;
fi;
# Calculate the vectors we will consider
ConvertToMatrixRep(g.matrix);
vs := [];
for k in [1..i] do
Add(vs, v);
v := v*g.matrix;
od;
# Construct the matrices [x,x^2,x^4,...,x^(2^k)]
# where 2^k <= n < 2^(k+1) in log n time
xs := [x];
for k in [1..Int(LogInt(fdpm.dim + 2, 2))] do
ConvertToMatrixRep(xs[k]);
Add(xs, xs[k]^2);
od;
# Calculate Basis objects for the vector spaces
# spanned by the vectors of vs. We can use these
# Basis objects to determine whether a vector lies
# in the corresponding vector space (which is need
# for calculating r(x,vg^i)).
vBases := List(vs, b -> BasisNC(VectorSpace(fdpm.field, [b]), [b]));
# Construct the images of the vectors in vs under
# x, ..., x^(n-1) (in log n time). Then we will have
# the vectors vg^ix, ..., vg^ix^(n-1) for i = 0, 1, 2,
# if g rep. a 3-cycle and i = 0, otherwise
ConvertToMatrixRep(vs);
vImgs := TransposedMat(RECOG.VectorImages(vs, xs));
# Calculate r(vg^i, x) for i = 0,1,2 if
# g is a 3-cycle and, i = 0, if g is a double
# transposition
rs := ListN(vBases, vImgs, R);
# for all vg^i such that r(vg^i, x) = r <> fail we check
# whether the list { l | 1 <= l < r, vx^l <> vx^lg } = 2.
# If the list does have length two then we can compute
# a linked sequence of vectors using vg^i
for k in [1..i] do
# check that r(vg^k, x) succeeded
if rs[k] = fail then
continue;
fi;
# calculate vx^lg for l = 1,..., r-1
mat := vImgs[k]{[1..rs[k]]};
vImgsG := mat * g.matrix;
# calculate { l | 1 <= l < r, vx^l <> vx^lg }
lList := Filtered([1..rs[k]-1], l -> vImgs[k][l+1] <> vImgsG[l+1]);
if Length(lList) = 2 then
# try to construct linked sequence of vectors
u := vImgs[k][lList[1]+1] - vImgsG[lList[1]+1];
if i = 3 then # if g is a 3-cycle
# check whether u in span of vs[k]
d := Coefficients(vBases[k], u);
if d <> fail then
h := -d[1]^-1*RECOG.Power(xs, lList[1]);
return RECOG.VectorImagesUnder(-vs[k], h, rs[k]-1){[1..rs[k]-1]};
fi;
# if u not in span of vs[k] then check
# whether u in span of vs[k]*g = vs[(k + 1) mod 3 + 1]
d := Coefficients(vBases[k^(1,2,3)], u);
if d <> fail then
h := d[1]^-1*RECOG.Power(xs, lList[1]);
return RECOG.VectorImagesUnder(vs[k], h, rs[k]-1){[1..rs[k]-1]};
fi;
# if neither then fail
return fail;
else # if g is a double transposition
# check whether u in span of vs[k]
d := Coefficients(vBases[k], u);
if d <> fail then
h := -d[1]^-1*RECOG.Power(xs, lList[1]);
return RECOG.VectorImagesUnder(vs[k], h, rs[k]-1){[1..rs[k]-1]};
fi;
# otherwise fail
return fail;
fi;
fi;
od;
return fail;
end;
## Functions for finding a basis from the linked sequence of
## vectors found by MoreBasisVectors
RECOG.RowSpaceBasis := function(mat)
# RowSpaceGenerators(<mat>) returns the rows of
# <mat> that form a basis for the row space
# of <mat>. This code borrows from the
# implementation of SemiEchcelonMatDestructive.
local basis, # vectors of mat forming a basis
# for the row space of mat
nzheads, # columns with leading 1 in the
# row echelon form of mat
vectors, # vectors spanning subspace of row
# space of mat
row, # current row of mat
nrows, # number of rows of mat
ncols, # number of columns of mat
x, # a field element
r, c; # row and column indexes
# calc size of matrix
nrows := Length(mat);
ncols := Length(mat[1]);
# list of the non-zero heads, that is, nzheads[i]
# is the first non-element of vectors[i]
nzheads := [];
# vectors contains a linearly independent
# set of vectors spanning a subset of the row-space
vectors := [];
# basis contains vectors that span the same subspace
# spanned by vectors, but the vectors in basis are
# from mat
basis := [];
for r in [1..nrows] do
row := ShallowCopy(mat[r]);
# reduce row using vectors
for c in [1..Length(nzheads)] do
x := row[nzheads[c]];
# if row does not already have a zero in this
# column
if not IsZero(x) then
AddRowVector(row, vectors[c], -x);
fi;
od;
# check row is not the zero row
c := PositionNonZero(row);
if c <= ncols then
# mult row by row[c]^-1 so leading coeff. is 1
MultRowVector(row, row[c]^-1);
Add(vectors, row);
Add(nzheads, c);
# add mat[r] to basis
Add(basis, mat[r]);
fi;
od;
return basis;
end;
RECOG.ExtendToBasisEchelon := function(partial, d, field)
# ExtendToBasisEchelon(<partial>, <d>, <field>)
# returns a basis for <field>^<d> such that <partial>
# is a prefix.
local vectors, v;
# copy partial
vectors := ShallowCopy(partial);
# append identity matrix
Append(vectors, IdentityMat(d, field));
# return row space basis
return RECOG.RowSpaceBasis(vectors);
end;
RECOG.IsIntervalVector := function(v, s)
# IsIntervalVector(<v>, <s>)
# returns [b, supp] if <v> is zero except
# for a uninterrupted sequence of b's within the
# first s positions. If v is not an interval
# vector then fail is returned
local col, b, i, j;
# find first non-zero entry in v
col := 1;
while col <= s and IsZero(v[col]) do
col := col + 1;
od;
# check first non-zero entry is within first s coords
if col > s then
return fail;
fi;
# else set i and b
i := col;
b := v[i];
# find first entry not equal to b
while col <= s and v[col] = b do
col := col + 1;
od;
# set j
j := col;
# check the rest of the entries are zero
while col <= Length(v) do
if not IsZero(v[col]) then
return fail;
fi;
col := col + 1;
od;
return [b, [i, j]];
end;
RECOG.FindBasis := function(x, vs, us, fdpm)
# FindBasis(<x>, <vs>, <us>, <fdpm>) extends the list
# of vectors in <vs> (and in <us>) to a linked basis
# for the fully deleted permutation module, described
# by <fdpm>, when given a random group element <x>.
# The lists are updated as a side effect. If a basis
# is found during the execution of this procedure
# then the basis is returned, otherwise 'fail' is
# returned.
local s, # the length of our linked sequence
P, # change of basis matrix
Pinv, # the inverse of P, if it exists
vsx, # the vectors vis*x w.r.t. basis vs
vsh, # vs*h
inti, # the interval of vi*x
ints, # interval of a sum of vi*x
i, j, # integers
ih, # the image of i under h where x = bh
vsum, # sum of vectors
vsumh, # sum of vectors in vsh
temp, # a row vector
mone, # -One(field)
b; # non-zero field element
s := Length(vs); # Length(us) should be s + 1
# calc vi*x w.r.t. basis vi: for finding interval vectors
P := RECOG.ExtendToBasisEchelon(vs, fdpm.dim, fdpm.field);
ConvertToMatrixRep(P);
Pinv := P^-1;
if Pinv = fail then
return fail;
else
vsx := P * x * Pinv;
ConvertToMatrixRep(vsx);
fi;
# Find scalar b and i,j such that x = bh and j = i^h
# find i <= s such that vix is an interval vector in <vs>
i := 1;
inti := RECOG.IsIntervalVector(vsx[i], s);
while inti = fail and i < s do
i := i + 1;
inti := RECOG.IsIntervalVector(vsx[i], s);
od;
# check whether we found an i
if inti <> fail then # we have; now we find the image of i
ih := fail; # we have not found the image of i yet
vsum := ShallowCopy(vsx[i]);
ConvertToVectorRep(vsum);
# find j s.t. vi*x + ... + vj*x is an interval vector in
# span of vs
for j in [i+1..s] do
AddRowVector(vsum, vsx[j]);
# find interval (if any) of vsum
ints := RECOG.IsIntervalVector(vsum, s);
if ints <> fail then
# we can find the image of i under h and
# the scalar b s.t. x = bh for h in H0
ih := Intersection(inti[2], ints[2]);
if Length(ih) <> 1 then
return fail; # won't happen if input group is Sn or An
else
ih := ih[1];
fi;
if inti[2][1] = ih then
b := inti[1];
else
b := -inti[1];
fi;
break; # leave the for loop
fi;
od;
# if we haven't found i^h try another way
if ih = fail then
vsum := ShallowCopy(vsx[i]);
ConvertToVectorRep(vsum);
# find j s.t. vj*x + ... + vi*x is an interval vector
for j in [i-1,i-2..1] do
AddRowVector(vsum, vsx[j]);
# find interval (if any) of vsum
ints := RECOG.IsIntervalVector(vsum, s);
if ints <> fail then # we can find (i+1)^h and so i^h and b
ih := Difference(inti[2], ints[2]);
if Length(ih) <> 1 then
return fail; # only happens when group not rep. of Sn or An
else
ih := ih[1];
fi;
if inti[2][1] = ih then
b := inti[1];
else
b := -inti[1];
fi;
break; # leave the for loop
fi;
od;
# check we have found ih
if ih = fail then
return fail; # if we have found it now we never will
fi;
fi;
else # look for i < s such that vi*x + v(i+1)*x is an interval vector
i := 1;
temp := ShallowCopy(vsx[i]);
ConvertToVectorRep(temp);
AddRowVector(temp, vsx[i+1]);
inti := RECOG.IsIntervalVector(temp, s);
while inti = fail and i < s - 1 do
i := i + 1;
temp := ShallowCopy(vsx[i]);
ConvertToVectorRep(temp);
AddRowVector(temp, vsx[i+1]);
inti := RECOG.IsIntervalVector(temp, s);
od;
# check whether we found an i
if inti <> fail then # we have; now we find the image of i
ih := fail; # we have not found the image of i yet
# we only try this for i <= s - 2
if i = s - 1 then
vsum := ShallowCopy(vsx[i]);
ConvertToVectorRep(vsum);
AddRowVector(vsum, vsx[i+1]);
# find j s.t. vi*x + v(i+1)^x + ... + vj*x is an interval vector
for j in [i+2..s] do
AddRowVector(vsum, vsx[j]);
# find interval (if any) of vsum
ints := RECOG.IsIntervalVector(vsum, s);
if ints <> fail then # we can find i^h and b
ih := Intersection(inti[2], ints[2]);
if Length(ih) <> 1 then
return fail; # only when group not rep. of Sn or An
else
ih := ih[1];
fi;
if inti[2][1] = ih then
b := inti[1];
else
b := -inti[1];
fi;
break; # leave the for loop
fi;
od;
fi;
# if we haven't found ih try another way
if ih = fail then
vsum := ShallowCopy(vsx[i]);
ConvertToVectorRep(vsum);
AddRowVector(vsum, vsx[i+1]);
# find j s.t. vj*x + ... + vi*x + v(i+1)*x is an interval vector
for j in [i-1,i-2..1] do
AddRowVector(vsum, vsx[j]);
# find interval (if any) of vsum
ints := RECOG.IsIntervalVector(vsum, s);
if ints <> fail then # we can find (i+1)^h and so i^h and b
ih := Difference(inti[2], ints[2]);
if Length(ih) <> 1 then
return fail; # only when group not rep. of Sn or An
else
ih := ih[1];
fi;
if inti[2][1] = ih then
b := inti[1];
else
b := -inti[1];
fi;
break; # leave the for loop
fi;
od;
# check we have found ih
if ih = fail then
return fail; # if we have found it now we never will
fi;
fi;
else
return fail;
fi;
fi;
# found ih and b
# calc -1 in field, for use with AddRowVector
mone := -One(fdpm.field);
# calc vi*h w.r.t. basis vs
vsx := b^-1 * vsx;
# calc vi*h not w.r.t. basis vs
vsh := vs * b^-1*x;
ConvertToMatrixRep(vsh);
# extend linked sequence
vsum := ShallowCopy(vsx[i]);
ConvertToVectorRep(vsum);
vsumh := ShallowCopy(vsh[i]);
ConvertToVectorRep(vsumh);
for j in [i+1..s] do
# calc vi*h + v(i+1)*h + ... + vj*h w.r.t to basis vs
AddRowVector(vsum, vsx[j]);
# calc vi*h + v(i+1)*h + ... + vj*h w.r.t. to standard basis
AddRowVector(vsumh, vsh[j]);
if RECOG.IsIntervalVector(vsum,s) = fail then
# vsum is not an interval vector so vsumh is an interval
# vector not in the span of vs
temp := ShallowCopy(vsumh);
ConvertToVectorRep(temp);
AddRowVector(temp, us[ih]);
Add(us, temp); # us[s+2] := vi*h + ... + vj*h + us[ih]
temp := ShallowCopy(us[s+2]);
ConvertToVectorRep(temp);
AddRowVector(temp, us[s+1], mone);
Add(vs, temp); # vs[s+1] := us[s+2] - us[s+1]
s := s + 1;
# if vs and us span the fully deleted permutation module
if s = fdpm.dim then
return [vs, us{[2..s+1]}];
fi;
fi;
od;
if i <> 1 then
# vsum = v1*h + ... + v(i-1)*h
vsum := ShallowCopy(vsx[1]);
ConvertToVectorRep(vsum);
for temp in [2..i-1] do
AddRowVector(vsum, vsx[temp]);
od;
vsumh := ShallowCopy(vsh[1]);
ConvertToVectorRep(vsumh);
for temp in [2..i-1] do
AddRowVector(vsumh, vsh[temp]);
od;
j := 0;
repeat
if RECOG.IsIntervalVector(vsum, s) = fail then
# vsum is an interval vector not in the span of vs
temp := ShallowCopy(us[ih]);
ConvertToVectorRep(temp);
AddRowVector(temp, vsumh, mone);
Add(us, temp);
temp := ShallowCopy(us[s+2]);
ConvertToVectorRep(temp);
AddRowVector(temp, us[s+1], mone);
Add(vs, temp);
s := s + 1;
if s = fdpm.dim then
# return vs and us if they span the fully deleted permutation
# module. Note: us[1] = u1 = 0
ConvertToMatrixRep(vs);
return [vs, us{[2..s+1]}];
fi;
fi;
j := j + 1;
# calc vj*h + ... + v(i-1)*h
AddRowVector(vsum, vsx[j], mone);
AddRowVector(vsumh, vsh[j], mone);
until j = i - 1;
fi;
# we have not found a basis for the fully deleted permutation
# module yet
return fail;
end;
## Functions for finding the permutation and scalar represented
## by a matrix given the vectors u_2,...,u_(n-delta+1)
RECOG.PositionsNot := function(list, elem)
# PositionsNot(<list>, <elem>) returns a
# list of all the positions in <list> whose
# value is not <elem>.
local pos, posList, zero;
# find position of first element <> elem
posList := [];
pos := PositionNot(list, elem);
# find subsequent positions
while pos <= Length(list) do
Add(posList, pos);
pos := PositionNot(list, elem, pos);
od;
return posList;
end;
RECOG.ExpectedForm := function(row, fdpmDim, nAmbig)
# ExpectedForm(<row>, <fdpmDim>, <nAmbig>) checks
# that <row> is of the form b*(u_ix - u_1x)
# for some b and {ix, 1x}; if so either [b, [ix, 1x]] is
# returned or [-b, [1x, ix]] is returned.
local posNonZero, # the positions of row whose values
# are non-zero
posNotB, # the positions of row whose values
# are not equal to b
b; # a (non-zero) field element
# find the non-zero positions of row
posNonZero := RECOG.PositionsNot(row, Zero(row[1]));
if Length(posNonZero) = 1 then
# if only one position is non-zero then
# 1x = 1 or ix = 1.
return [row[posNonZero[1]], [posNonZero[1] + 1, 1]];
elif Length(posNonZero) = 2 then
# check row = [..b..-b..] for some b
if row[posNonZero[1]] = -row[posNonZero[2]] then
return [row[posNonZero[1]], posNonZero + 1];
else
return fail;
fi;
elif Length(posNonZero) >= fdpmDim - 1 and nAmbig then # only if delta = 2
# determine +/- b: find two entries that are the same
if row[1] = row[2] or row[1] = row[3] then
b := row[1];
elif row[2] = row[3] then
b := row[2];
else # we have three different field elems when there
# can be only two
return fail;
fi;
# find positions in row that are not b
posNotB := RECOG.PositionsNot(row, b);
if IsEmpty(posNotB) then
# {ix,1x} = {1,n} so b(u_ix - u_1x) = bu_ix
return [b, [1, fdpmDim + 2]];
elif Length(posNotB) = 1 then
# {ix, 1x} = {j,n} with j <> 1 and
# row[k] = b for k <> j - 1 and row[j-1] = 2b
if row[posNotB[1]] = 2*b then
return [b, [posNotB[1]+1, fdpmDim + 2]];
fi;
else # there are three different field elements in row
return fail;
fi;
else # wrong number positions with non-zero entries
return fail;
fi;
end;
RECOG.FindPermutation := function(rs, fdpm)
# FindPermutation(<rs>, <fdpmDim>) returns the scalar
# and permutation associated with a given matrix. If
# the matrix is b*x (for some scalar b) then <rs> contains
# the row vectors ui*b*x (where ui = e1 - ei) written
# with respect to the vectors ui. If the matrix is not
# a representation of a permutation then fail is
# returned.
local perm, # a permutation
permList, # the image list of the perm represented
r, # a row
oneIm, # the image of 1 under the permutation
imgsA, # image of ui under bx
imgsB, # image of uj under bx
FindB, # function to determine b from image of ui
b, # the scalar that the matrix represents
pts; # a list of images of points under a permutation
FindB := function(pmb, imgs);
# input is b, [u_ix, 1x] or
# -b and [1_x, u_ix]
if oneIm = imgs[1] then
return -pmb;
else
return pmb;
fi;
end;
# determine +/- b and [1x,2x]
imgsA := RECOG.ExpectedForm(rs[1], fdpm.dim, fdpm.nAmbig);
# determine +/- b and [1x,3x]
imgsB := RECOG.ExpectedForm(rs[2], fdpm.dim, fdpm.nAmbig);
if imgsA <> fail and imgsB <> fail then
# determine 1x
oneIm := Intersection(imgsA[2], imgsB[2]);
if Length(oneIm) = 1 then
oneIm := oneIm[1];
else
return fail;
fi;
permList := [oneIm];
# determine b
b := CallFuncList(FindB, imgsA);
else
return fail;
fi;
# build up permList by looking at rs[i] = u_(i+1)x (w.r.t to basis
# us) for i = 2,...,n-delta
for r in rs do
# check whether rs[i] in is an expected form
imgsB := RECOG.ExpectedForm(r, fdpm.dim, fdpm.nAmbig);
if imgsB = fail then
return fail;
else # rs[i] is in an expected form so imgsB[2] = [1x, (i+1)x]
pts := Difference(imgsB[2], [oneIm]);
# check (i+1)x has been defined uniquely and b is correct
if Length(pts) = 1 and b = CallFuncList(FindB, imgsB) then
Append(permList, pts);
else
return fail;
fi;
fi;
od;
# convert permList into permutation
# if delta = 2 then n will appear in permList but Length(permList) = n - 1
# so we need to append the missing point to permList
Append(permList, Difference([1..Maximum(permList)], permList));
perm := PermListList([1..Length(permList)], permList);
if perm = fail then
return fail;
fi;
return [b, perm];
end;
## The main function that puts all the above together to
## recognise our matrix group:
RECOG.RecogniseFDPM := function(group, field, eps)
# RecogniseFDPM(<group>, <field>, <eps>) is an example
# of using the algorithm described in ' ' to recognise
# a matrix group representation of An or Sn acting on
# their fully deleted permutation module. It actually computes
# data to apply RECOG.FindPermutation.
local fdpm, # a record containing information on
# the fully deleted permutation module
g, # group element representing 3-cycle or double trans
v, # a vector from a linked basis
k, # a counter
vs, # a linked sequence of vectors
us, # the alterate basis
uvs, # output of CompleteBasis
cob, # change of basis matrix, to the basis u_i
# for i = 2,...,n-delta+1
cobi, # inverse matrix
mgens; # the matrix generators
# construct matrix group
mgens := GeneratorsOfGroup(group);
# information on the fully deleted permutation module
fdpm := rec( field := field,
fieldChar := Characteristic(field),
dim := DimensionOfMatrixGroup(group),
nAmbig := (DimensionOfMatrixGroup(group) + 2) mod
Characteristic(field) = 0
);
# check the dimensions allow for the possibility that group is a
# representation of An or Sn on its fully deleted permutation module
if (fdpm.dim + 1) mod fdpm.fieldChar = 0 then
return fail;
fi;
Info(InfoFDPM, 2, "Searching for 3-cycle or double-transposition.");
# find 3-cycle or double transposition
if fdpm.fieldChar = 3 then
# check fdpm.dim is largest for IsPreDoubleTransposition to
# work
if fdpm.dim < 6 then
ErrorNoReturn( "This function only works for matrix groups of dimension",
" at least 6 over a field of this characteristic.");
return fail;
fi;
# search for double transpositions
g := RECOG.IterateWithRandomElements(
m -> RECOG.ConstructDoubleTransposition(m, fdpm),
RECOG.LogRat(eps^-1, 2)*(RootInt(fdpm.dim + 2) + 1)*41,
# 1/0.0249 < 41
group);
else
# check fdpm.dim is largest for IsPre3Cycle to work
if fdpm.dim < 4 then
ErrorNoReturn( "This function only works for matrix groups of dimension",
" at least 4 over a field of this characteristic.");
return fail;
fi;
# search for a 3-cycle
g := RECOG.IterateWithRandomElements(
m -> RECOG.Construct3Cycle(m, fdpm),
RECOG.LogRat(eps^-1, 2)*(RootInt(fdpm.dim+2,3) + 1)*15,
# 1/0.07 < 15
group);
fi;
# Note: when the field characteristic is >= 5 we could search
# for both 3-cycles and double transpositions. However we would still
# require about as many random element selections as if we searched
# only for 3-cycles and, if we found a double transposition, we would
# require many more random element selections for MoreBasisVectors.
# This is especially important when group is not a representation of
# An or Sn.
# check we found a 3-cycle or double transposition
if g = fail then
Info(InfoFDPM, 2, "Could not find 3-cycle or double-transposition.");
return fail;
fi;
if g.order = 3 then
Info(InfoFDPM, 2, "Found 3-cycle.");
else
Info(InfoFDPM, 2, "Found double-transposition.");
fi;
# find first basis vector
v := RECOG.FindBasisElement(g, group, eps, fdpm);
if v = fail then
Info(InfoFDPM, 2, "Could not find basis element.");
return fail;
fi;
Info(InfoFDPM, 2, "Found first basis element.");
# extend basis
if g.order = 3 then
vs := RECOG.IterateWithRandomElements(
m -> RECOG.MoreBasisVectors(m, g, v, fdpm),
34*RECOG.LogRat(eps^-1, 2)*RECOG.LogRat(fdpm.dim + 2, 2),
# 1/0.03 < 34
group);
else
vs := RECOG.IterateWithRandomElements(
m -> RECOG.MoreBasisVectors(m, g, v, fdpm),
2000*RECOG.LogRat(eps^-1, 2)*RECOG.LogRat(fdpm.dim + 2, 2),
group);
fi;
if vs = fail then
Info(InfoFDPM, 2, "Could not find linked sequence of vectors.");
return fail;
fi;
Info(InfoFDPM, 2, "Found linked sequence of vectors.");
# complete basis
# compute start of alternate basis
us := [Zero(vs[1])];
for k in [1..Length(vs)] do
Add(us, us[k] + vs[k]);
od;
ConvertToMatrixRep(us);
uvs := RECOG.IterateWithRandomElements(
m -> RECOG.FindBasis(m, vs, us, fdpm),
RECOG.LogRat(fdpm.dim + 2, 2)*( RECOG.LogRat(eps^-1, 2)
+ RECOG.LogRat(fdpm.dim + 2, 2)),
group);
if uvs = fail then
Info(InfoFDPM, 2, "Could not complete basis.");
return fail;
fi;
Info(InfoFDPM, 2, "Found change of basis matrix.");
# determine permutations and scalars: uvs[2] contains the
# change of basis matrix that takes group elements into
# our prefered form
cob := uvs[2];
ConvertToMatrixRep(cob);
cobi := cob^-1;
if cobi = fail then
return fail;
fi;
return rec( cob := cob, cobi := cobi, fdpm := fdpm );
end;
