1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
#############################################################################
##
#W union.gi Karel Dekimpe
#W Bettina Eick
##
#############################################################################
##
#F GenericDeterminantMat( mat )
##
GenericDeterminantMat := function( mat )
local d, det, sig, i, sub;
# set up
d := Length( mat );
if ForAny( mat, x -> Length(x) <> d ) then return fail; fi;
# the trivial cases
if d = 1 then return mat[1][1]; fi;
if d = 2 then return mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]; fi;
# otherwise use first row and recursion
det := 0;
sig := 1;
for i in [1..d] do
sub := Concatenation( [1..i-1], [i+1..d] );
sub := mat{[2..d]}{sub};
det := det + sig * mat[1][i] * GenericDeterminantMat( sub );
sig := - sig;
od;
return det;
end;
#############################################################################
##
#F NullspaceIntMod( base, vec, p ) . . . . . . . . . . . . b * vec = 0 mod p
##
## computes those elements b in <base> with b * vec = 0 mod p. Returns a
## triangulized basis with respect to <base>.
##
NullspaceIntMod := function( base, vec, p )
local imgs, d, null;
# get images
imgs := List( base, x -> x * vec );
imgs := List( imgs, x -> x mod p );
d := Length( imgs );
if ForAll( imgs, x -> x = 0 ) then return IdentityMat( d ); fi;
# compute kernel of imgs vector - must have rank d-1
Add( imgs, p );
null := NullspaceIntMat( TransposedMat( [imgs] ) );
null := List( null, x -> x{[1..d]} );
null := NormalFormIntMat( null, 0 ).normal;
# return this images nullspace
return Filtered( null, x -> PositionNonZero(x) <= d );
end;
#############################################################################
##
#F FindMaximals( sub ) . . . . . . . .subgroups which are maximal within sub
##
FindMaximals := function( sub )
local new, i, tmp;
new := [];
for i in [1..Length(sub)] do
if Size( sub[i] ) > 1 then
if Length(new) = 0 then
Add( new, sub[i] );
elif not ForAny( new, x -> IsSubset( x, sub[i] ) ) then
tmp := Filtered( new, x -> not IsSubset( sub[i], x ) );
Add( new, sub[i] );
fi;
fi;
od;
return new;
end;
if not IsBound( SizeOfUnion ) then SizeOfUnion := false; fi;
#############################################################################
##
#F SizeOfUnionRec( sub ) -- recursive version
##
SizeOfUnionRec := function( list )
local s, i, int, t, n;
if Length( list ) = 0 then return 1; fi;
s := Size( list[1] );
n := Length( list );
for i in [2..n] do
int := List( list{[1..i-1]}, x -> Intersection( x, list[i]));
t := SizeOfUnion( int );
s := s + Size( list[i] ) - t;
od;
return s;
end;
#############################################################################
##
#F SizeOfUnionTriv -- trivial version
##
SizeOfUnionTriv := function( list )
return Length( Union( List( list, Elements ) ) );
end;
#############################################################################
##
#F SizeOfUnion -- main function
##
SizeOfUnion := function( sub )
local list;
list := FindMaximals( sub );
if Length(list) = 0 then return 1; fi;
if Length(list) = 1 then return Size(list[1]); fi;
if Sum(List(list, Size)) < 2000 then
return SizeOfUnionTriv(list);
fi;
return SizeOfUnionRec(list);
end;
#############################################################################
##
#F SizeOfUnionMod( subs, e ) . . . . . . . . .size of the union of subs mod e
##
## <subs> is a list of bases containing (eZ)^d. Compute the size of the union
## of <subs> in (Z/eZ)^d.
##
## This function needs to be profiled.
##
SizeOfUnionMod := function( subs, e )
local d, F, V, b, news;
# the trivial case
if Length( subs ) = 0 then return 1; fi;
d := Length( subs[1] );
if IsPrimeInt( e ) then
F := GF(e);
V := F^d;
b := BasisVectors( Basis( V ) );
news := List( subs, x -> x * b );
news := List( news, x -> Subspace( V, x ) );
if ForAny( news, x -> Size(x) = e^d ) then return e^d; fi;
# return SizeOfUnion( news );
return Length( Union( List( news, Elements ) ) );
fi;
V := AbelianGroup( List( [1..d], x -> e ) );
b := GeneratorsOfGroup(V);
news := List( subs, x -> List( x, y -> MappedVector( y, b ) ) );
news := List( news, x -> Subgroup( V, x ) );
if ForAny( news, x -> Size(x) = e^d ) then return e^d; fi;
# return SizeOfUnion( news );
return Length( Union( List( news, Elements ) ) );
end;
|
