Upper Bounds and Probability Heuristic Estimates of the number of π‘Ÿ-Regular Families of Permutations Pavol KollΓ‘r1,*,† , Martin Mačaj2,† 1 Department of Applied Informatics, Faculty of Mathematics, Physics and Informatics, Comenius University, Bratislava, Slovakia 2 Department of Algebra and Geometry, Faculty of Mathematics, Physics and Informatics, Comenius University, Bratislava, Slovakia Abstract A sure-fire way of answering β€œHow many of these objects are there?” is to generate them all and simply keeping a counter. However, if it so happens that the true count is about 1024 , this way of doing it becomes a β€œsure-fire way of getting nowhere”. This is exactly what happens with so called π‘Ÿ -regular families of permutations, which are a generalisation of the notion of transitivity in groups. They also measure how Cayley-like a vertex-transitive graph is. Therefore, another approach is required for this (and many other) cases. In this article, we shall describe a general-purpose algorithm that utilises complex cyclotomic numbers to give upper bounds to the counts of these families. Additionally, we present a probabilistic algorithm to give heuristic estimates of the true counts of these families. Keywords π‘Ÿ-regular families, Generating functions, Fourier transformations, GAP 1. Introduction The simpler version of 1-regular family (or simply just regular family) was introduced much earlier in [2] As defined in [1], for any set 𝑀 , an β€œπ‘Ÿ-regular family” already. In this paper, the author was looking at the well of permutations is a set β„± βŠ† S𝑀 , such that for every known issue of graph theory, that there are some vertex π‘₯, 𝑦 ∈ 𝑀 , there are exactly π‘Ÿ permutations πœ‘ ∈ β„±, transitive graphs, which are not Cayley graphs, most such that πœ‘(π‘₯) = 𝑦. Sometimes we will take 𝑀 = famously the Petersen graph. Her method involved the {1, 2, 3, . . . , 𝑛}, in which case we will use the standard construction of β€œquasi-Cayley graphs”, which are tied to shorthand [𝑛] := {1, 2, 3, . . . , 𝑛}. The π‘Ÿ-regular family the existence of 1-regular families of automorphisms as as an object is closely related to transitive groups, when follows: β€œa graph Ξ“ is quasi-Cayley if and only if 𝐴𝑒𝑑(Ξ“) their stabiliser is of size π‘Ÿ. Indeed, such a group does contains a 1-regular family as a subset”. This is similar follow the property of π‘Ÿ-regular families and is therefore to what Sabidussi did in [3] with Cayley graphs and their a special case of them. On the other hand, the π‘Ÿ-regular automorphism groups: β€œa graph Ξ“ is Cayley if and only families are combinatorial approximations of these tran- if 𝐴𝑒𝑑(Ξ“) contains a subgroup, which acts regularly on sitive groups, as they themselves need not to follow the 𝑉 (Ξ“)”. group axioms, only the β€œπ‘Ÿ regularity”. With definitions laid out as we have, we can rephrase The authors of [1] investigated β€œhow close to a Cayley the result of Sabidussi as β€œa graph Ξ“ is Cayley if and only graph a given vertex-transitive graph is”, which is a ques- if 𝐴𝑒𝑑(Ξ“) contains a subgroup, which is also a 1-regular tion that has been asked many times in prior research in family (with respect to its action on 𝑉 (Ξ“))”. different variations. As such, there are multiple measures One year after the paper [1], in his bachelor thesis of this β€œcloseness” notion. One such notion is to find [4], the author investigated π‘Ÿ-regular families from a the smallest transitive subgroup of 𝐴𝑒𝑑(Ξ“). A different computational perspective, specifically, by generating notion relies precisely on the aforementioned π‘Ÿ-regular these families with the help of a computer. families, and also gets investigated in [1]. There, one In Section 2, we will summarise KerΓ‘k’s main result seeks small π‘Ÿ-regular families as subsets of the 𝐴𝑒𝑑(Ξ“), from [4], which is relevant to our work, introduce the the group of automorphisms of the given graph Ξ“. The main problem we will be solving, and describe how gen- paper [1] is the first instance we could find which men- erating functions and Fourier transformations can help tions the use of π‘Ÿ-regular families. us get upper bounds on the number of π‘Ÿ-regular families. Continuing that, in Section 3, we talk about some of the ITAT’24: Computational Aspects of Large-Scale Problems in Discrete implementation details of these methods to get said up- Mathematics, September 20–24, 2024, Drienica, Slovakia * per bounds. In Section 4, we briefly describe the method Corresponding author. † [4] used and expand on it with a randomised algorithm, These authors contributed equally. $ pavol.kollar@fmph.uniba.sk (P. KollΓ‘r); which can be used as an estimation heuristic to approx- martin.macaj@fmph.uniba.sk (M. Mačaj) imate the number of π‘Ÿ-regular families on 𝑛 elements. Β© 2024 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). Finally, in Section 5, we summarise our findings. CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) CEUR ceur-ws.org Workshop ISSN 1613-0073 Proceedings Table 1 Generated counts of π‘Ÿ-regular families: π‘Ÿ=0 π‘Ÿ=1 π‘Ÿ=2 π‘Ÿ=3 π‘Ÿ=4 π‘Ÿ=5 π‘Ÿ=6 π‘Ÿ=7 Β·Β·Β· 𝑛=1 1 1 - - - - - - Β·Β·Β· 𝑛=2 1 1 - - - - - - Β·Β·Β· 𝑛=3 1 2 1 - - - - - Β·Β·Β· 𝑛=4 1 24 255 640 255 24 1 - Β·Β·Β· 𝑛=5 1 1344 11073216 NG NG NG NG NG Β·Β·Β· 𝑛=6 1 1128960 NG NG NG NG NG NG Β·Β·Β· .. .. .. .. .. .. .. .. .. .. . . . . . . . . . . 2. Ground Work Given any π‘Ÿ ≀ (𝑛 βˆ’ 1)!, we know that βˆ€β„± ∈ π‘…π‘Ÿπ‘€ : |β„±| = π‘›π‘Ÿ. This facilitates βƒ’the fact βƒ’ that the only 0-regular 2.1. Establishing notions set is the empty set and so ⃒𝑅0𝑀 βƒ’ = 1 for all 𝑀 regardless of its size. Since we have chosen |𝑀 | = 𝑛, we have that Our work picks up where [4] left off, where the author |S𝑀 | = 𝑛!, which means that for all βƒ’π‘Ÿ > βƒ’(𝑛 βˆ’ 1)!, the was generating the π‘Ÿ-regular families with the help of a set π‘…π‘Ÿπ‘€ is empty, so βˆ€π‘Ÿ > (𝑛 βˆ’ 1)! : βƒ’π‘…π‘Ÿπ‘€ βƒ’ = 0. computer. One of his outputs was the Table 1 containing Our main effort for the remainder of this paper is to information on how many distinct π‘Ÿ-regular families on get upper bounds for the total number of π‘Ÿ-regular fam- 𝑛 elements are there for some small π‘Ÿ, 𝑛. Within this ilies for 𝑛 = 5. In Section 4 we will also describe a table, β€œNG” means that for those parameters, the families randomised algorithm that provides a heuristic estimate were only partially generated or not at all. for the number of π‘Ÿ-regular families on 𝑛 elements are For the rest of the paper, we fix 𝑛 ∈ N and also let there in that specific case. We are also going to reference 𝑀 = {π‘š1 , π‘š2 , . . . , π‘šπ‘› } be a fixed set of 𝑛 positive in- the values within the row for 𝑛 = 4 in the Table 1 to tegers with π‘š1 ≀ π‘š2 ≀ . . . ≀ π‘šπ‘› . Borrowing some check the correctness of our methods. notation from [4], instead of writing down permutations of S𝑀 with the cycle notation as is usual in algebra, e.g. (1, 2, 4, 6) or (2, 7, 3)(5, 11), we will write down all the 2.2. Generating functions numbers in one sequence into a list. In this list, the num- To bound the values of the numbers of π‘Ÿ-regular families, ber π‘˜ is in the 𝑝th position in the list precisely when π‘šπ‘ is we will use generating functions and their properties. To mapped to the number π‘˜ by this permutation. Put another be able to use generating functions, we require a weight way, to each permutation πœ‘ ∈ S𝑀 , we associate the list of function on the subsets of S , namely π‘Š : 𝒫(S ) β†’ N. 𝑀 𝑀 numbers 𝐿(πœ‘) = [πœ‘(π‘š1 ), πœ‘(π‘š2 ), πœ‘(π‘š3 ), . . . , πœ‘(π‘šπ‘› )]. Instead of producing a weight for every subset, we assign For 𝑀 = [6] and πœ‘ = (1, 2, 4, 6), this list will thus each permutation a weight with the function 𝑀 : S β†’ 𝑀 be [2, 4, 3, 6, 5, 1] and for 𝑀 = {2, 3, 5, 7, 11, 13} and N and the total weight of a subset of S shall simply be 𝑀 πœ‘ = (2, 7, 3)(5, 11), the list will be [7, 2, 11, 3, 5, 13]. the sum of individual weights. When we talk about the elements of S𝑀 from now on, The weight function of single permutations 𝑀 : S𝑀 β†’ we will mostly talk about the set of lists induced by the N will depend on a chosen list w = [𝑀 , . . . , 𝑀 ] of 1 𝑛 previous construction. weights (and also the chosen set 𝑀 ) and is defined by βˆ‘οΈ€π‘› Observation 1. Within this list association, an π‘Ÿ-regular 𝑀(πœ‘) = 𝑖=1 𝑀𝑖 Β· 𝐿(πœ‘)[𝑖]. In a sense, we really should family of S𝑀 is such a collection of the lists β„± βŠ† S𝑀 , so call this weight function 𝑀w , but we will omit it in this that when the lists in β„± are aligned below one another, section. For example, if 𝑀 = [4] and w = [1, 2, 5, 15], each column contains each π‘š ∈ 𝑀 precisely π‘Ÿ times. for the identity permutation we have 𝑀(𝑒) = 1 Β· 1 + 2 Β· 2 + 5 Β· 3 + 15 Β· 4 = 80. For πœ‘ = (1, 2), we have Later, we will want to talk about individual en- 𝐿(πœ‘) = [2, 1, 3, 4] and 𝑀(πœ‘) = 79, and for πœ‘ = (1, 4, 3), tries of these lists and for that reason, we shall ac- we have 𝐿(πœ‘) = [4, 2, 1, 3] and 𝑀(πœ‘) = 62. cess them via indices written with square brackets, i.e. As mentioned before, the weight function for sets 𝐿(πœ‘)[𝑖] = πœ‘(π‘šπ‘– ). For us, the first entry has index 1, e.g. of permutations π‘Š : 𝒫(S𝑀 ) β†’ N, will be the sum [7, 2, 11, 3, 5, 13][1] = 7, [2, 4, 3, 6, 5, 1][4] = 6. of the individual permutation weights, i.e., π‘Š (𝑋) = βˆ‘οΈ€ 𝑀(πœ‘). For any 𝐴 βŠ† N, we can ask for the set Definition 1. For a given non-negative integer π‘Ÿ, let 𝐢 πœ‘βˆˆπ‘‹ 𝑀 𝐴 = {𝐹 βŠ† S𝑀 |π‘Š (𝐹 ) ∈ 𝐴}. Ideally we want to con- π‘…π‘Ÿ βŠ† 𝒫(S𝑀 )⋃︀be the set of π‘Ÿ-regular families of S𝑀 . struct the function π‘Š (that is, choose the list w as well Also let 𝑅𝑀 = π‘Ÿβ‰₯0 π‘…π‘Ÿπ‘€ . as the set 𝑀 in a clever way) so that there exists an 𝐴 with β„± ∈ π‘…π‘Ÿπ‘€ ⇔ π‘Š (𝐹 ) ∈ 𝐴. However, ensuring both we know βˆοΈ€this generating function 𝑔 is also equal to implications hold while also ensuring that he entries of 𝑔(π‘₯) = πœ‘βˆˆS𝑀 (1 + π‘₯𝑀(πœ‘) ), creating the following es- 𝑀 and w are sufficiently small is highly non-trivial (for sential equation: 𝑀, w big, it is not so difficult to come up with examples), βˆ‘οΈ ∏︁ so we will relax the requirement and demand that only π‘Žπ‘˜ Β· π‘₯π‘˜ = 𝑔(π‘₯) = (1 + π‘₯𝑀(πœ‘) ) (2) β„± ∈ π‘…π‘Ÿπ‘€ β‡’ π‘Š (𝐹 ) ∈ 𝐴 holds. This is why our algo- π‘˜βˆˆN πœ‘βˆˆS𝑀 rithms will, in general, produce upper bounds for those Utilising a standard trick of Fourier transformations, counts, and occasionally produce the exact values. let 𝜁 a primitive 𝑑th root of unity. Then to get the bound- Continuing the simple example above with 𝑀 = [4], βˆ‘οΈ€ ing sum π‘ŸβˆˆN π‘Žπ‘Ÿπ‘‘ , we can average the values of 𝑔(π‘₯) at one can see that for any π‘Ÿ-regular family β„±, each weight is multiplied π‘Ÿ times by each π‘šπ‘– ∈ 𝑀 , so π‘Š (β„±) = all the π‘‘π‘‘β„Ž roots of unity, that is: (1 + 2 + 5 + 15) Β· π‘Ÿ Β· (1 + 2 + 3 + 4) = 230π‘Ÿ. In general 𝑑 βˆ‘οΈ 1 βˆ‘οΈ we have the following lemma. π‘Žπ‘Ÿπ‘‘ = Β· 𝑔(𝜁 𝑖 ) (3) π‘ŸβˆˆN 𝑑 𝑖=1 Lemma 1. Fix w = [𝑀1 , . . . , 𝑀𝑛 ], and letβˆ‘οΈ€β„± βŠ† S𝑀 be an π‘Ÿ-regular family. Then π‘Š (β„±) = π‘Ÿ Β· βˆ‘οΈ€ π‘šπ‘– βˆˆπ‘€ π‘šπ‘– Β· βˆ‘οΈ€ Combining Equation 1 and Equation 3, we get βˆ‘οΈ€π‘€π‘– ∈w 𝑀𝑖 . Therefore, with 𝑑 = π‘šπ‘– βˆˆπ‘€ π‘šπ‘– Β· 𝑑 𝑀 and 𝐴 = {π‘˜ Β· 𝑑|π‘˜ ∈ the condition 1 βˆ‘οΈ βƒ’ βƒ’ 𝑖 N}, βƒ’ 𝑀 βƒ’ βˆ‘οΈ 𝑀𝑖 ∈w ⃒𝑅 βƒ’ ≀ π‘Žπ‘Ÿπ‘‘ = Β· 𝑔(𝜁 𝑖 ) (4) β„± ∈ π‘…π‘Ÿπ‘€ β‡’ π‘Š (β„±) ∈ 𝐴 is satisfied for all π‘Ÿ, in par- π‘ŸβˆˆN 𝑑 𝑖=1 ticular β„± ∈ 𝑅𝑀 β‡’ π‘Š (β„±) ∈ 𝐴. The evaluation of this average requires the evaluations Proof. Let β„± be an π‘Ÿ-regular family. By definition, of 𝑔 on the right hand side. Those values depend on βˆ‘οΈ βˆ‘οΈ βˆ‘οΈ the weight function 𝑀, which in turn depends on the π‘Š (β„±) = 𝑀(πœ‘) = 𝑀𝑖 Β· 𝐿(πœ‘)[𝑖], particular choice of 𝑀 and w. This constitutes the first πœ‘βˆˆβ„± πœ‘βˆˆβ„± 𝑀𝑖 ∈w algorithm, which accepts inputs 𝑀, w and outputs the upper bound from Equation 4, see Section 3.1 for imple- where the 𝐿(πœ‘)[𝑖] are elements of 𝑀 . This is a finite sum mentation. and hence we can swap the order of summations: βˆ‘οΈ βˆ‘οΈ π‘Š (β„±) = . . . = 𝑀𝑖 Β· 𝐿(πœ‘)[𝑖]. 2.3. Generating functions in more than 𝑀𝑖 ∈w πœ‘βˆˆβ„± one variable Because β„± is π‘Ÿ-regular, together over all 𝑖 and all πœ‘, In the previous section, weβˆ‘οΈ€ transformed a permutation 𝐿(πœ‘)[𝑖] equals each π‘šπ‘— ∈ 𝑀 precisely π‘Ÿ times: πœ‘ into its weight 𝑀(πœ‘) = 𝑀𝑖 ∈w 𝑀𝑖 Β· 𝐿(πœ‘)[𝑖] and we βˆ‘οΈ βˆ‘οΈ transformed a permutation set 𝐹 βŠ† S𝑀 into π‘Š (𝐹 ) = π‘Š (β„±) = . . . = 𝑀𝑖 Β· π‘Ÿ Β· π‘šπ‘— . βˆ‘οΈ€ πœ‘βˆˆπΉ 𝑀(πœ‘). Enconding each set of permutations into 𝑀𝑖 ∈w π‘šπ‘— βˆˆπ‘€ a single number caused us to lose a lot of information about said set of permutations. Instead, let us assign To finish the proof of the claim, pull π‘Ÿ out of both sums, the β€œweight” of the permutation to be the permutation since it is independent of the sums, and substitute in the itself in the list notation: 𝑀(πœ‘) = 𝐿(πœ‘). Treating these defined 𝑑: lists of integers as vectors for the purposes of addition and scalar multiplication, we may now define π‘Š (𝐹 ) = βˆ‘οΈ βˆ‘οΈ π‘Š (β„±) = . . . = π‘Ÿ Β· 𝑀𝑖 Β· π‘šπ‘— = π‘Ÿ Β· 𝑑. βˆ‘οΈ€ 𝑀𝑖 ∈w π‘šπ‘— βˆˆπ‘€ πœ‘βˆˆπΉ 𝑀(πœ‘) for any 𝐹 βŠ† S𝑀 . This way, we retain more information about 𝐹 . This means that 𝑑 divides π‘Š (β„±), so π‘Š (β„±) ∈ 𝐴. Lemma 2. βˆ‘οΈ€ Let β„± βŠ† S𝑀 be an π‘Ÿ-regular family. Then Said another way that is more usual in the context of π‘Š (β„±) = π‘ŸΒ· π‘šβˆˆπ‘€ π‘šΒ·[1, βˆ‘οΈ€ 1, . . . , 1], where the vector has generating functions, if π‘Žπ‘˜ = |{𝐹 βŠ† S𝑀 : π‘Š (𝐹 ) = π‘˜}|, |𝑀 | = 𝑛 ones. With d = π‘šβˆˆπ‘€ π‘š Β· [1, 1, . . . , 1] and then the following bound is a direct consequence of the 𝐴 = {π‘˜ Β· d|π‘˜ ∈ N}, the condition β„± ∈ π‘…π‘Ÿπ‘€ β‡’ π‘Š (β„±) ∈ previous lemma: 𝐴 is satisfied for all π‘Ÿ, in particular β„± ∈ 𝑅𝑀 β‡’ π‘Š (β„±) ∈ βƒ’ βƒ’ 𝐴. βƒ’ 𝑀 βƒ’ βˆ‘οΈ ⃒𝑅 βƒ’ ≀ π‘Žπ‘Ÿπ‘‘ (1) π‘ŸβˆˆN Proof. The proof is analogous to the proof of the one- dimensional case from the previous subsection. Each We take 𝑔(π‘₯) to be the generating function for this se- position within the lists gets each number π‘š ∈ 𝑀 pre- quence of numbers, that is, 𝑔(π‘₯) = π‘›βˆˆN π‘Žπ‘› Β· π‘₯𝑛 . By cisely π‘Ÿ times in the expression for π‘Š (β„±), when β„± is an βˆ‘οΈ€ standard theory of generating functions (see e.g. [5]), π‘Ÿ-regular family. To avoid the tiny writing of double indexing in a multi- the polynomial 𝑔. This means that the values 𝑔(𝜁 𝑖 ) are variable coefficient, let I = [𝑖1 , 𝑖2 , . . . , 𝑖𝑛 ]. Denoting by unambiguously determined by the input. 𝑏I = |{𝐹 βŠ† S𝑀 : π‘Š (𝐹 ) = I}|, Lemma 2 tells us that To evaluate these products, we will need a program βƒ’ βƒ’ that can handle addition and multiplication of these kinds βƒ’ 𝑀 βƒ’ βˆ‘οΈ of complex numbers - formally known as cyclotomic num- ⃒𝑅 βƒ’ ≀ π‘π‘ŸΒ·d (5) π‘ŸβˆˆN bers - with exact precision. Fortunately, the GAP system [6] has rich support for arithmetic with complex roots Once again, take 𝐺(x) = 𝐺(π‘₯1 , π‘₯2 , . . . , π‘₯𝑛 ) to be the of unity, cyclotomic numbers, their sums, their products multivariate generating function for the values of 𝑏: (and of much more). This - alongside our previous expe- βˆ‘οΈ rience with the system - is why GAP is our software of 𝐺(x) = 𝑏I Β· π‘₯𝑖11 Β· π‘₯𝑖22 Β· . . . Β· π‘₯𝑖𝑛𝑛 . choice for these computations, the pseudocode for which I∈N𝑛 we present in Algorithm 1. Just like in the one dimensional case, 𝐺 can be written as a product, which we record in the next equation: Algorithm 1 One variable bounding algorithm Input: 𝑀, w βƒ’ βƒ’ Output: Upper bound for the size ⃒𝑅𝑀 βƒ’ βˆ‘οΈ 𝑏I Β· π‘₯𝑖11 Β· π‘₯𝑖22 Β· . . . Β· π‘₯𝑖𝑛𝑛 = I∈N𝑛 π‘π‘œπ‘€π‘’π‘Ÿπ‘ _π‘œπ‘“ _π‘₯ ← [ ] = 𝐺(π‘₯1 , π‘₯2 , . . . , π‘₯𝑛 ) = for π‘π‘’π‘Ÿπ‘š ∈ S𝑀 do ◁ Here, βˆ‘οΈ€π‘π‘’π‘Ÿπ‘š is taken as a list. (6) 𝑛 π‘π‘œπ‘€π‘’π‘Ÿπ‘ _π‘œπ‘“ _π‘₯.append( 𝑛 𝑖=1 π‘π‘’π‘Ÿπ‘š[𝑖] Β· 𝑀𝑖 ) ∏︁ ∏︁ 𝑀(πœ‘)[𝑖] = (1 + π‘₯𝑖 ) end for βˆ‘οΈ€ βˆ‘οΈ€ πœ‘βˆˆS𝑀 𝑖=1 𝑑 ← π‘šπ‘– βˆˆπ‘€ π‘šπ‘– Β· 𝑀𝑖 ∈w 𝑀𝑖 βˆ‘οΈ€ 𝜁 ← 𝐸(𝑑) ◁ Primitive root of unity. Let 𝑠 = π‘šπ‘– βˆˆπ‘€ π‘šπ‘– , which means that d can be written π‘Žπ‘π‘π‘’π‘šπ‘’π‘™π‘Žπ‘‘π‘œπ‘Ÿ ← 0 as [𝑠, 𝑠, . . . , 𝑠]. Here too, we want to use the trick of for 1 ≀ 𝑖 ≀ 𝑑 do Fourier transformations. In this case we can start with π‘π‘Ÿπ‘œπ‘‘π‘’π‘π‘‘ ← 1 the observation that: for π‘π‘œπ‘€π‘’π‘Ÿ ∈ π‘π‘œπ‘€π‘’π‘Ÿπ‘ _π‘œπ‘“ _π‘₯ do 1 βˆ‘οΈ βˆ‘οΈ π‘π‘Ÿπ‘œπ‘‘π‘’π‘π‘‘ ← (1 + 𝜁 π‘–Β·π‘π‘œπ‘€π‘’π‘Ÿ ) Β· π‘π‘Ÿπ‘œπ‘‘π‘’π‘π‘‘ 𝑛 Β· 𝐺(𝜁 𝑖1 , . . . , 𝜁 𝑖𝑛 ) = 𝑏I (7) end for 𝑠 I∈Z𝑛 I≑0 (mod 𝑠) 𝑠 π‘Žπ‘π‘π‘’π‘šπ‘’π‘™π‘Žπ‘‘π‘œπ‘Ÿ ← π‘Žπ‘π‘π‘’π‘šπ‘’π‘™π‘Žπ‘‘π‘œπ‘Ÿ + π‘π‘Ÿπ‘œπ‘‘π‘’π‘π‘‘ Here, the second sum is over all vector-indices, whose end for each coordinate is a multiple of 𝑠. Clearly, the sum on return π‘Žπ‘π‘π‘’π‘šπ‘’π‘™π‘Žπ‘‘π‘œπ‘Ÿ/𝑑 ◁ Always an integer. the right-hand-side includes all coefficients with indices π‘˜ Β· d, therefore we have the following upper bound as a Testing this algorithm on various inputs, we recog- consequence: nised that to get a meaningful decrease of the computed upper bound, we needed to substantially increase the 1 βˆ‘οΈ βƒ’ βƒ’ value of 𝑑. Because 𝑔(𝜁 𝑖 ) = πœ‘βˆˆS𝑀 (1 + 𝜁 𝑖·𝑀(πœ‘) ), the βƒ’ 𝑀 βƒ’ βˆ‘οΈ βˆοΈ€ ⃒𝑅 βƒ’ ≀ π‘π‘ŸΒ·d ≀ 𝑛 Β· 𝐺(𝜁 𝑖1 , . . . , 𝜁 𝑖𝑛 ) (8) 𝑠 π‘ŸβˆˆN I∈Z𝑛 𝑠 computation of each 𝑔(𝜁 𝑖 ) is done in linear time with respect to the order of 𝜁 (which is 𝑑), and therefore the The algorithm implementing this approach and its runtime complexity is 𝑂(𝑛! Β· 𝑑2 ). Coupled with the fact optimisations is described in Subsection 3.3. that we had to substantially increase 𝑑 to get a signifi- cant improvement, we stopped using this exact approach, 3. Algorithmic Upper Bounds even if the further algorithms share the same core idea. Still, we include the description of the algorithm here to illustrate this core idea. Later, in the results, we will 3.1. The basic algorithm highlight where the particular approaches had their com- Let us return to the ideas from Section 2.2 and create an putational limits. algorithm based on them. For the Subsection 3.1, fix w = [𝑀1 , βˆ‘οΈ€ βˆ‘οΈ€ 𝑀 . Furthermore, as before, let . . . , 𝑀𝑛 ] and also 3.2. Long lists of coefficients 𝑑 = π‘šπ‘– βˆˆπ‘€ π‘šπ‘– Β· 𝑀𝑖 ∈w 𝑀𝑖 . At the end of Subsection 2.2, we ended up with the bounding inequality described To make matters for the previous algorithm worse, even in 4. Our algorithm will thus, for every 𝜁 𝑖 , evaluate the if we could reduce the time complexity of the algorithm polynomial 𝑔 at that value and average these results. We down, we still have a memory overflow waiting to hap- remark that by fixing 𝑀, w (which will be the inputs of pen. In fact, one of our intermediate algorithm iterations the algorithm), we are also fixing the value 𝑑 as well as had time complexity just 𝑂(𝑛! Β· 𝑑), because instead of evaluating 𝑑 different values 𝑔(𝜁 𝑖 ), this algorithm ex- in Section 3.1, the evaluations themselves are linear in 𝑠 panded the product form of 𝑔 symbolically. This is be- in terms of time, so the total time complexity is roughly cause what the computation actually needs is the sum of 𝑂(𝑛! Β· 𝑠𝑛+1 ). For just 𝑛 = 5 and 𝑠 = 127, that causes certain polynomial coefficients of 𝑔 - specifically those the algorithm to run for about 15 hours, while giving a whose index is ≑ 0 (mod 𝑑). To perform this expansion weak bound. of 𝑔(π‘₯), whose product form we established in Equation Since the evaluations of 𝐺 themselves cannot be easily (2), we represent each intermediate step as a long list sped up, let us try to decrease the number of evaluations. of coefficients, in which each position is the sum of all We will choose some properties of the set 𝑀 that will coefficients, whose indices have the same remainder lead to optimisation of our algorithm’s speed. With the (mod 𝑑). At the end, the result we seek is then simply claims that follow, we will prove that we can reduce the the 0th coefficient of this list. complexity down to only 𝑂(π‘ π‘›βˆ’1 ). Also, to simplify While this provided a significant time-improvement, the further expressions, we will write 𝐺[𝑖1 , . . . , 𝑖𝑛 ] in- given that we only β€œevaluated” the polynomial once, the stead of 𝐺(𝜁 𝑖1 , . . . , 𝜁 𝑖𝑛 ). Because of the equality in (6), memory requirement was roughly 𝑂(𝑑), i.e., an approxi- it easily follows that: mately linear amount of memory is required for this com- putation. On our laptop, even with tricks of the β€œmeet in Lemma 3. For any [𝑖1 , . . . , 𝑖𝑛 ] ∈ Z𝑛 𝑠 and any πœ™ ∈ S𝑛 the middle” character, choosing 𝑑 β‰ˆ 1.2Β·108 was roughly we have the limit of its memory capabilities. A stronger machine 𝐺[𝑖1 , . . . , 𝑖𝑛 ] = 𝐺[π‘–πœ™(1) , . . . , π‘–πœ™(𝑛) ]. could be employed for the task, but that too does not solve the problem. This is because the memory problem For fixed 𝑛, the next lemma will gain us a much more is, in a sense, unavoidable. Once we start to dabble with significant speed-up: cyclotomic numbers that we need to represent exactly, and we need to have πœ™(𝑑), where πœ™ is the Euler totient Lemma 4. Adding 1 to all the exponents of the roots of function, coefficients to represent those. Choosing 𝑑 to unity does not change the value of 𝐺, i.e.: have many distinct prime factors seemed to yield worse 𝐺[𝑖1 , . . . , 𝑖𝑛 ] = 𝐺[𝑖1 + 1, . . . , 𝑖𝑛 + 1]. upper bounds than when 𝑑 had only a few of those and then πœ™(𝑑) is still 𝑂(𝑑), so we cannot escape the issue, at Proof. Focus on the specific product π‘₯π‘š π‘š2 1 Β·π‘₯2 Β·. . .Β·π‘₯𝑛 1 π‘šπ‘› least not in this way. within 𝐺. Evaluating it at the powers of 𝜁 from the left- hand-side, we of course get another power of 𝜁. The 3.3. Thinking in more dimensions value of that power is 𝑛 As we have argued, the basic algorithms have large mem- βˆ‘οΈ ory requirements, some with bad runtimes in relation to π‘–π‘˜ Β· π‘šπ‘˜ . π‘˜=1 the results we were getting from them. Our third and most recently used method for this problem seeks to im- On the other hand, evaluating that product on the right- prove this via generating polynomials of many variables, hand-side gives us (again omitting the 𝜁 itself and writing as was discussed and introduced in Section 2.3. down just the power): In what follows, 𝑀 is still the set {π‘š1 , π‘š2 , . . . , π‘šπ‘› } βˆ‘οΈ€ 𝑛 𝑛 and let 𝑠 = π‘šπ‘– βˆˆπ‘€ π‘šπ‘– . To simplify the optimisation βˆ‘οΈ βˆ‘οΈ (π‘–π‘˜ + 1) Β· π‘šπ‘˜ = π‘–π‘˜ Β· π‘šπ‘˜ + 𝑠. arguments that follow in the later part of the section, π‘˜=1 π‘˜=1 assume that 𝑠 is a prime number bigger than 𝑛, hence 𝑛 and 𝑠 are coprime. Finally, let 𝜁 be a primitive 𝑠th root These two expressions only differ by the +𝑠 at the end of unity, just as in the Section 2.3. In that section, we and both of these expressions represent powers of 𝜁. arrived at the inequality in (8), which bounded the total Since 𝜁 is a primitive 𝑠th root of unity, the +𝑠 makes number of all π‘Ÿ-regular families on 𝑛 elements by the no difference to the product. expression: As such, this product within 𝐺 has the same value on the inputs we started with and, analogously, so will every 1 βˆ‘οΈ Β· 𝐺(𝜁 𝑖1 , . . . , 𝜁 𝑖𝑛 ). other product in 𝐺. This proves our claim. 𝑠𝑛 I∈Z𝑛 𝑠 We can iteratively repeat this argument to get the This means that all our algorithm needs to do is to following equality as an immediate corollary evaluate the multinomial 𝐺 on all possible vectors of 𝐺[𝑖1 , . . . , 𝑖𝑛 ] = 𝐺[𝑖1 + π‘˜, . . . , 𝑖𝑛 + π‘˜] (9) roots of unity and average the result. However, this is outrageously slow. Notice that for that to happen, we that holds for all 1 ≀ π‘˜ ≀ 𝑠, with π‘˜ = 𝑠 being the trivial require 𝑠𝑛 evaluations of 𝐺. Moreover, as we discussed case, as exponents of 𝜁 operate (mod 𝑠). Picking 𝑠 to be prime is essential for the second big (in the case that 𝑠 is prime and greater than 𝑛) for I = optimisation. Before we get to that optimisation, we men- [0, . . . , 0] we have 𝐡I = 1 and for Ξ“I = {𝑒Γ } we have tion a bit about how GAP internally stores cyclotomic 𝐡I = 𝑠 Β· 𝑛!. numbers. We touched on this earlier in Section 3.2, and A standard way to utilize the Theorem 1 is to choose we will be more specific here. For prime 𝑠 the cyclo- canonical representatives of orbits of Ξ“. Then one needs tomic field Q(𝜁) is an 𝑠 βˆ’ 1 dimensional linear space to derive efficient ways of determining whether given over Q. Instead of β€œthe standard basis” 1, 𝜁, . . . , 𝜁 π‘ βˆ’2 , vector I is the representative of its orbit as well as how GAP uses the basis 𝜁, 𝜁 2 , . . . , 𝜁 π‘ βˆ’1 , that is, each element to compute 𝐡I . One possible way to pick the canonical of Q(𝜁) is expressed as 𝑐1 𝜁 + 𝑐2 𝜁 2 , . . . + π‘π‘ βˆ’1 𝜁 π‘ βˆ’1 and representative is to choose those I’s, which are lexico- it is internally represented as vector [𝑐1 , 𝑐2 . . . , π‘π‘ βˆ’1 ]. In graphically smallest as vectors. For this choice of repre- GAP, it is possible to obtain this vector via the function sentatives, we have the following restrictions: CoeffsCyc. Let us note that any rational number π‘ž is β€’ A canonical representative has non-decreasing represented by the vector [βˆ’π‘ž, βˆ’π‘ž, . . . , βˆ’π‘ž]. elements, because of Lemma 3. Let us also note that this basis is invariant under the Galois group of Q(𝜁) which consists of automorphisms β€’ Thanks to Lemma 4, the representative has sum πœ“π‘˜ , which map 1 ↦→ 1, 𝜁 ↦→ 𝜁 π‘˜ and act on the roots of its coefficients ≑ 0 (mod 𝑠). of unity as a permutation for all 1 ≀ π‘˜ < 𝑠. Note β€’ As a consequence of Theorem 1, the first non-zero that πœ“1 is the identity map. With all the notation and coordinate of a canonical representative is 1. knowledge of this subsection, we are ready to state the However not all such vectors are representatives of next optimisation claim: their orbits and the algorithm needs to be able to de- 𝑛 termine which is β€œthe representative”. For example, if βˆ‘οΈ€π‘ βˆ’1 5.𝑗 Let [𝑖1 , . . . , 𝑖𝑛 ] ∈ Z𝑠 and let 𝐺[𝑖1 , . . . , 𝑖𝑛 ] = Lemma 𝑀 = [1, 2, 3, 4, 7] and 𝑠 = 17 both [0, 1, 2, 3, 11] and 𝑗=1 𝑐𝑗 𝜁 . Then [0, 1, 9, 10, 14] belong to the same orbit of Ξ“. π‘ βˆ’1 The pseudocode for this improved algorithm is in Algo- βˆ‘οΈ 𝐺[π‘˜π‘–1 , . . . , π‘˜π‘–π‘› ] = rithm 2. For full implementation details of this algorithm, π‘˜=1 see our Github repository [7]. π‘ βˆ’1 π‘ βˆ’1 (10) βˆ‘οΈ βˆ‘οΈ πœ“π‘˜ (𝐺[𝑖1 , . . . , 𝑖𝑛 ]) = βˆ’ 𝑐𝑗 π‘˜=1 𝑗=1 4. Heuristic Counting with Backtracking and Probability Proof. Let us note that for any [π‘Ž1 , . . . , π‘Žπ‘› ] ∈ Z𝑛 we have πœ“π‘˜ ((𝜁 𝑖1 )π‘Ž1 Β·. . .Β·(𝜁 𝑖𝑛 )π‘Žπ‘› ) = (𝜁 π‘˜π‘–1 )π‘Ž1 Β·. . .Β·(𝜁 π‘˜π‘–π‘› )π‘Žπ‘› . Let us return to the work of KerΓ‘k [4] and his recursive Therefore, by (6), we have πœ“π‘˜ (𝐺[𝑖1 , . . . , 𝑖𝑛 ]) = algorithm to generate all π‘Ÿ-regular families on 𝑛 elements. 𝐺[π‘˜π‘–1 , . . . , π‘˜π‘–π‘› ], which is the first claimed equation. In what follows, assume 𝑛 = 5 is fixed and, since the As 𝑠 is prime, Galois group of Q(𝜁) acts regularly on βˆ‘οΈ€π‘ βˆ’1 specific choice of 𝑀 does not change the number of π‘Ÿ- , 𝜁 π‘ βˆ’1 }. Therefore {𝜁, . . . (︁ π‘˜=1 πœ“π‘˜ (𝐺[𝑖1 , . . . , 𝑖𝑛 ]) = regular families, let 𝑀 = [5] = {1, 2, 3, 4, 5}. On top of βˆ‘οΈ€π‘ βˆ’1 βˆ‘οΈ€π‘ βˆ’1 )︁ 𝑖 βˆ‘οΈ€π‘ βˆ’1 𝑖 that, partition the set S𝑀 into five sets 𝑆1 , . . . , 𝑆5 , so that 𝑖=1 𝑗=1 𝑐𝑗 𝜁 . Because 𝑖=1 𝜁 = βˆ’1, rear- βˆ‘οΈ€π‘ βˆ’1 the set 𝑆𝑖 contains precisely those permutations πœ‘ ∈ S𝑀 , ranging the last double-sum yields βˆ’ 𝑗=1 𝑐𝑗 . such that πœ‘(1) = 𝑖. We remark that this makes all 𝑆𝑖 Previous lemmas describe how corresponding actions have equal size: |𝑆1 | = . . . = |𝑆5 | = (5 βˆ’ 1)! = 24. of S𝑛 , Z𝑠 , and Z*𝑠 on the set Z𝑛 Secondly, with this partition in place, we can see that 𝑠 modify evaluation of 𝐺(X). Combining these lemmas we obtain information any π‘Ÿ-regular family β„± βŠ† S𝑀 must arise as a union of about the action of group Ξ“ = S𝑛 Γ— Z𝑠 Γ— Z*𝑠 on the set π‘Ÿ-sized subsets 𝑇𝑖 βŠ† 𝑆𝑖 , such that all the other positions Z𝑛 contain each 𝑖 precisely π‘Ÿ times too. 𝑠 and its behaviour with the evaluations of 𝐺(X). That is precisely what KerΓ‘k’s algorithm does - recur- Theorem 1. Let I = [𝑖1 , . . . , 𝑖𝑛 ] ∈ Z𝑛 𝑠 and let sively trying out all subsets 𝑇𝑖 βŠ† 𝑆𝑖 of size π‘Ÿ, exiting 𝐺[𝑖1 , . . . , 𝑖𝑛 ] = 𝑐1 𝜁 + . . . + π‘π‘ βˆ’1 𝜁 π‘ βˆ’1 . Then there exists early if any position overflows , i.e., if any position has a number 𝐡I which depends only on (the conjugacy class more than π‘Ÿ of any number 𝑖 ∈ {1, 2, 3, 4, 5}. Once a of) the stabilizer Ξ“I such that 𝑇𝑖 βŠ† 𝑆𝑖 has been found to not overflow any position βˆ‘οΈ with any number, his algorithm proceeds recursing in a 𝐺[𝑗1 , . . . , 𝑗𝑛 ] = 𝐡I (βˆ’π‘1 βˆ’ . . . βˆ’ π‘π‘ βˆ’1 ) (11) depth-first manner, until all possible combinations of sub- J∈IΞ“ sets of the 𝑆𝑖 ’s have been tested. After all these families have been generated, counting them up is easy. The exact correspondence between 𝐡I and Ξ“I is too This approach is correct and it generates everything, complicated to present here. Let us just point out that but the caveat here is that there are too many possibilities. Algorithm 2 Multivariate bounding Table 2 Input: 𝑀 βƒ’ βƒ’ Assorted results of the heuristic counting algorithm: Output: Upper bound for the size ⃒𝑅𝑀 βƒ’ Sum of all counts exactly Scientific notation 𝑛 ← Length(𝑀 ) 1344808311315380068372992 1.345 Β· 1024 𝑠 ← Sum(𝑀 ) ◁ We could test 𝑠 for primality. 1687144238754283380215130 1.687 Β· 1024 𝐺 ← πΆπ‘œπ‘šπ‘π‘’π‘‘π‘’πΊπ‘’π‘›π‘’π‘Ÿπ‘Žπ‘‘π‘–π‘›π‘”πΉ π‘’π‘›π‘π‘‘π‘–π‘œπ‘›( ) ◁ Helper 4952098573211019632607892 4.952 Β· 1024 function. 3192278092063163925627380 3.192 Β· 1024 π΅π‘–π‘”π‘†π‘’π‘š ← 0 10722414177898381929138566 1.072 Β· 1025 for I ∈ Z𝑛 𝑠 do 5578509063361797937216054 5.579 Β· 1024 if Sum(𝐼) ΜΈ= 0 (mod 𝑠) then 2047532970827157893748616 2.048 Β· 1024 continue ◁ Not a canonical representative. Average of the trials: (rounded) end if 4217826489633026395275233 4.218 Β· 1024 if not(IsSortedList(I)) then continue ◁ Not a canonical representative. end if metric. Still, providing concrete bounds for the error this if First(I, π‘₯ β†’ π‘₯ > 0) ΜΈ= 1 then creates is something we have not done, so it cannot be continue ◁ Not a canonical representative. said that this algorithm is much more than a method end if for obtaining a heuristic estimation of the true counts of 𝐼𝑠𝑀 π‘–π‘›π‘–π‘šπ‘Žπ‘™ ← true π‘Ÿ-regular families on 𝑛 elements. for π‘˜ ∈ Difference(I, [0]) do The trick is to only explore a handful of the branches J ← I/π‘˜ (mod 𝑠) into full depth, thus substantially reducing the number if J < I then of explored nodes. Specifically, we will assign each re- 𝐼𝑠𝑀 π‘–π‘›π‘–π‘šπ‘Žπ‘™ ← false cursion depth a probability, with which we let a branch break ◁ Not a canonical representative. deeper. We keep track, for each of the levels, how many end if branches reached the level, and how many branches were end for let through. At the end, thanks to our β€œall branches look if 𝐼𝑠𝑀 π‘–π‘›π‘–π‘šπ‘Žπ‘™ then the same” assumption, we can just divide the final count 𝐡I ← πΆπ‘œπ‘šπ‘π‘’π‘‘π‘’π΅I (I) ◁ Helper function. by the β€œnumbers of branches let through” and multiply π΅π‘–π‘”π‘†π‘’π‘š ← π΅π‘–π‘”π‘†π‘’π‘š + 𝐡I Β· by the β€œnumbers of branches reaching that level”. Sum(CoeffsCyc(𝐺[I], 𝑠)) The probabilities for recursion depths for 𝑛 = 5 and end if π‘Ÿ ∈ {2, 3, 4, . . . , 12} were determined by experimenta- end for tion, so that the expected number of branches that are return π΅π‘–π‘”π‘†π‘’π‘š/𝑠𝑛 let through on each recursion depth is between 20 and 200. This lets us get a non-trivial sample of the decision tree without making the algorithm absurdly long. The Even for 2-regular families on 5 elements, each 𝑆𝑖 has eleven instances usually run to completion in 7-8 hours, 24 permutations, which means that at each level, we depending on the random branches that are picked by have 24 (οΈ€ )οΈ€ 2 = 276 choices to pick from. Doing that on 5 the algorithm. different levels, the total amount of possibilities is 2765 = We experimented with this algorithm, when 𝑛 = 4, 1601568101376, which is starting to be too much for a where the actual count of all π‘Ÿ-regular families is 1200. laptop. Even for a more powerful machine, it quickly gets On the vast majority of the trials, the algorithm outputs unruly and completely infeasible for families of higher a number between 1150-1250, so we are reasonably con- regularity - for example the naive strategy for 8-regular fident that the algorithm does not give wildly incorrect families on 5 elements has the amount of possibilities on results. The average of the results of the last couple of the order of β‰ˆ 2 Β· 1029 . recorded trials is 1202. At the time of writing, we ran the To simplify matters of memory usage, let us use this algorithm for 𝑛 = 5 seven times, and those results are in algorithm to only count the number of the families with the Table 2. We also include the average of the trials. We this algorithm, not fully generate them all. This means considered excluding the largest and the smallest value that at the bottom level, if we confirm that the chosen set from the average, but because we only have 7 data points is an π‘Ÿ-regular family, we just increment a global counter. until now, we decided against it. Since our tree of possibilities is too big, let us assume For more detailed results, see the table in our Github that β€œit will look about the same in every branch”. repository [7]. We remark that in 𝑛 = 5 case, one does This is not an outrageous assumption as we are study- not need to go past π‘Ÿ = 12, even though π‘Ÿ-regular fami- ing regular sets of permutations, which are quite sym- lies exist up to π‘Ÿ = 24, because for each π‘Ÿ-regular family Algorithm 3 Probabilistic backtracking algorithm π‘π‘œπ‘’π‘›π‘‘π‘’π‘Ÿ ← 0 π‘‘π‘’π‘π‘‘β„Ž_π‘π‘Ÿπ‘œπ‘π‘Žπ‘π‘–π‘™π‘–π‘‘π‘–π‘’π‘  ← [350, 2000, 3000, 200, 1] ◁ Chosen for 𝑛 = 5, π‘Ÿ = 4. π‘‘π‘’π‘π‘‘β„Žπ‘ _𝑙𝑒𝑑_𝑓 π‘œπ‘Ÿπ‘€π‘Žπ‘Ÿπ‘‘, π‘‘π‘’π‘π‘‘β„Žπ‘ _π‘Ÿπ‘’π‘Žπ‘β„Žπ‘’π‘‘ ← [0, 0, 0, 0, 0], [0, 0, 0, 0, 0] procedure RecursiveBacktrack(π‘π‘Žπ‘Ÿπ‘‘π‘–π‘Žπ‘™, π‘‘π‘’π‘π‘‘β„Ž) if π‘‘π‘’π‘π‘‘β„Ž = 𝑛 then ◁ End of recursion and it has been checked already. π‘π‘œπ‘’π‘›π‘‘π‘’π‘Ÿ ← π‘π‘œπ‘’π‘›π‘‘π‘’π‘Ÿ + 1 else for 𝑛𝑒𝑀_π‘π‘’π‘Ÿπ‘šπ‘  ∈ π‘Ÿ-sized subsets of 𝑆𝑖 do ◁ The sets 𝑆𝑖 have been created already. 𝑛𝑒𝑀_π‘π‘Žπ‘Ÿπ‘‘π‘–π‘Žπ‘™ ← 𝑛𝑒𝑀_π‘π‘Žπ‘Ÿπ‘‘π‘–π‘Žπ‘™ βˆͺ 𝑛𝑒𝑀_π‘π‘’π‘Ÿπ‘šπ‘  if βˆƒπ‘₯, 𝑦 ∈ [𝑛], βˆƒπ‘‡ βŠ† 𝑛𝑒𝑀_π‘π‘Žπ‘Ÿπ‘‘π‘–π‘Žπ‘™ : |𝑇 | > π‘Ÿ ∧ βˆ€πœ‘ ∈ 𝑇 : πœ‘(π‘₯) = 𝑦 then ◁ Check overflows. continue ◁ Too many permutations mapping π‘₯ to 𝑦. Try another one. else π‘‘π‘’π‘π‘‘β„Žπ‘ _π‘Ÿπ‘’π‘Žπ‘β„Žπ‘’π‘‘[π‘‘π‘’π‘π‘‘β„Ž] ← π‘‘π‘’π‘π‘‘β„Žπ‘ _π‘Ÿπ‘’π‘Žπ‘β„Žπ‘’π‘‘[π‘‘π‘’π‘π‘‘β„Ž] + 1 if randint(1, π‘‘π‘’π‘π‘‘β„Ž_π‘π‘Ÿπ‘œπ‘π‘Žπ‘π‘–π‘™π‘–π‘‘π‘–π‘’π‘ [π‘‘π‘’π‘π‘‘β„Ž]) = 1 then ◁ Going further deep in the recursion with probability 1/π‘‘π‘’π‘π‘‘β„Ž_π‘π‘Ÿπ‘œπ‘π‘Žπ‘π‘–π‘™π‘–π‘‘π‘–π‘’π‘ [π‘‘π‘’π‘π‘‘β„Ž]. π‘‘π‘’π‘π‘‘β„Žπ‘ _𝑙𝑒𝑑_𝑓 π‘œπ‘Ÿπ‘€π‘Žπ‘Ÿπ‘‘[π‘‘π‘’π‘π‘‘β„Ž] ← π‘‘π‘’π‘π‘‘β„Žπ‘ _𝑙𝑒𝑑_𝑓 π‘œπ‘Ÿπ‘€π‘Žπ‘Ÿπ‘‘[π‘‘π‘’π‘π‘‘β„Ž] + 1 RecursiveBacktrack(𝑛𝑒𝑀_π‘π‘Žπ‘Ÿπ‘‘π‘–π‘Žπ‘™, π‘‘π‘’π‘π‘‘β„Ž + 1) end if end if end for end if end procedure RecursiveBacktrack({}, 0) ◁ To start the recursive algorithm off. for (π‘Ÿπ‘’π‘Žπ‘β„Žπ‘’π‘‘, 𝑙𝑒𝑑_𝑓 π‘œπ‘Ÿπ‘€π‘Žπ‘Ÿπ‘‘) ∈ zip(π‘‘π‘’π‘π‘‘β„Žπ‘ _π‘Ÿπ‘’π‘Žπ‘β„Žπ‘’π‘‘, π‘‘π‘’π‘π‘‘β„Žπ‘ _𝑙𝑒𝑑_𝑓 π‘œπ‘Ÿπ‘€π‘Žπ‘Ÿπ‘‘) do ◁ Python’s zip allows us to iterate through two lists in sync. π‘π‘œπ‘’π‘›π‘‘π‘’π‘Ÿ ← π‘π‘œπ‘’π‘›π‘‘π‘’π‘Ÿ Β· π‘Ÿπ‘’π‘Žπ‘β„Žπ‘’π‘‘/𝑙𝑒𝑑_𝑓 π‘œπ‘Ÿπ‘€π‘Žπ‘Ÿπ‘‘ end for return round(π‘π‘œπ‘’π‘›π‘‘π‘’π‘Ÿ) ◁ To have an integral result. β„± βŠ† S𝑀 , its complement 𝒒 = S𝑀 βˆ– β„± is (𝑛 βˆ’ 1)! βˆ’ π‘Ÿ- are disivible by 5 constitute a trivial (and very excessive!) regular. This means, that after 12, whatever the true bound for the total number of π‘Ÿ-regular families: counts were in the first half will also be in the second (οΈƒ )οΈƒ (οΈƒ )οΈƒ (οΈƒ )οΈƒ half, reflected about the column π‘Ÿ = 12. Hence, the total 120 120 120 βƒ’ βƒ’ βƒ’ 𝑀⃒ ⃒𝑅 βƒ’ ≀ + +. . .+ β‰ˆ 2.658Β·1035 . number of π‘Ÿ-regular families that exist on 5 elements is 0 5 120 then calculated as double of all the numbers, except for the π‘Ÿ = 12 case, that should be counted only once. The best bound the basic algorithm could produce before running out of memory was 10757440577219567348022770930 β‰ˆ 1.076 Β· 1028 . 5. Conclusion After that, we switched to the algorithm using the multivariate generating function and we include a In this paper, we described methods of both getting up- few of the output results in the table too. The lowest per bounds for the number of π‘Ÿ-regular families on 𝑛 upper bound we obtained at the time of writing is elements as well as a method to get a heuristic estima- 4263880475370843510800356 β‰ˆ 4.264Β·1024 . Coupled tion of the true counts. We checked our work for 𝑛 = 4 with the estimation heuristic makes us fairly confident elements, where to get a fast enough version of the al- the upper bounds are approaching the true count. gorithm 2, we did not need to use Theorem 1, so 𝑠 did Let us also contrast our method to the naive approach not need to be prime, just coprime to 𝑛. The bounds here of β€œtest every subset of S and check π‘Ÿ-regularity”. We 5 quickly converge to the true count of 1200 (see Table 3) know there are 2120 subsets of S and let us say we 5 and the estimation heuristic also gets really close to this could check 230 of them in a second (this is approaching number (see the discussion in Section 4). roughly the clock speed of a modern computer and is For the case 𝑛 = 5, one can immediately see that an probably a vast exaggeration of the machine’s capabili- π‘Ÿ-regular family has size divisible by 5. Therefore, the ties) - this approach would thus take 290 seconds, roughly number of subsets of a set with 120 elements, whose sizes 3 Β· 109 ages of the universe. Table 3 Assorted results of the bounding multivariate algorithm for 𝑛 = 4: βˆ‘οΈ€ π‘š1 π‘š2 π‘š3 π‘š4 𝑠 = 𝑖 π‘šπ‘– Time (ms) Bound 1 2 4 8 15 5 5112 1 4 16 64 85 1268 1200 1 6 36 216 259 85609 1200 1 8 64 512 585 1566917 1200 Table 4 Assorted results of the multivariate algorithm for 𝑛 = 5: βˆ‘οΈ€ π‘š1 π‘š2 π‘š3 π‘š4 π‘š5 𝑠 = 𝑖 π‘šπ‘– Time (ms) Bound Trivial bound - 0 265845599161775836797571384326161716 One-dimensional method - ? 22669475893267033898649677785186 1 2 4 8 16 31 34 1439304569993444516046531000316 1 2 5 11 24 43 99 388799463844737009990155623596 Long list of coefficients method - ∼ 1800000 10757440577219567348022770930 1 3 9 27 81 127 4925 5109628694684595790776607756 16 24 36 54 93 223 48255 540465527918830847892764076 1 4 16 64 592 677 4561658 12647340576853621194664376 1 5 27 141 733 907 16534644 6248117909711976319109756 1 5 25 125 883 1039 28784230 5348880561458906796605766 1 5 25 125 1403 1559 ∼ 85384605 4263880475370843510800356 1 7 49 343 2401 2801 ? 3926985392178291058321116 On the other hand, for reasons not mentioned in this 0019 grant. paper, choosing 𝑀 = {1, 25, 625, 15625, βƒ’ βƒ’ 390625} is The second author acknowledges support from the sufficient to get the precise value of ⃒𝑅𝑀 βƒ’. Here, 𝑠 = APVV Research Grant APVV-19-0308 and from the VEGA 406901, and the algorithm evaluates roughly (in fact Research Grants 1/0423/20, 1/0727/22 and 1/0437/23. less than) 𝑠3 points of the multinomial 𝐺, which means roughly 6.737 Β· 1016 evaluations. Pesimistically, let us say that one 𝐺 evaluation takes a minute, which would References still β€œonly take” about 128, 177, 239, 751 years, which [1] R. Jajcay, G. A. Jones, π‘Ÿ-regular families of graph au- is 9.3 ages of the universe. All that is to say that the tomorphisms, European Journal of Combinatorics 79 competition is not even close. (2019) 97–110. doi:https://doi.org/10.1016/ However, the 𝑛 = 6 case is tricky. For one, because the j.ejc.2018.12.002. upper bounding algorithm has time complexity 𝑂(π‘ π‘›βˆ’1 ), [2] G. Gauyacq, On quasi-cayley graphs, Discrete Ap- raising 𝑛 increases the runtime drastically. To get a mean- plied Mathematics 77 (1997) 43–58. doi:https:// ingful bound, one will probably need to substantially im- doi.org/10.1016/S0166-218X(97)00098-X. prove this algorithm further. The estimation algorithm [3] G. Sabidussi, Vertex-transitive graphs, Monatshefte too will not work too well, because the branching factors fΓΌr Mathematik 68 (1964) 426–438. URL: https://api. in the probabilistic backtracking themselves get way too semanticscholar.org/CorpusID:120359810. big, e.g. 120 Β· 1034 . (οΈ€ )οΈ€ 60 β‰ˆ 9.661 [4] F. KerΓ‘k, π‘Ÿ-regular sets of permutations, Bachelor’s thesis, Comenius University, Bratislava, 2020. Acknowledgments [5] H. Wilf, generatingfunctionology: Third Edition, CRC Press, 2005. URL: https://books.google.sk/ The first author acknowledges Funding by the EU books?id=bVFuBwAAQBAJ. NextGenerationEU through the Recovery and Resilience [6] The GAP Group, GAP – Groups, Algo- Plan for Slovakia under the project No. 09I03-03-V02- rithms, and Programming, Version 4.12.2, 00036. The first author also acknowledges support from https://www.gap-system.org, 2022. the Grant of Comenius University No. UK/1114/2024. [7] P. Kollar, R-RegularFamilies, https://github.com/ Thirdly, the first author acknowledges support from the OldAnchovyTopping/R-RegularFamilies, 2024. Ac- grants VEGA Research Grant 1/0437/23 and SK-AT-23- cessed: 2024-06-30.