#  GEMMFIP: Unifying GEMM in BLIS

RuQing G. Xu
Department of Physics
The University of Tokyo
Hongo 7-3-1, Bunkyo-Ward, Tokyo 113-0033, Japan
<EMAIL_ADDRESS>
Field G. Van Zee
Robert A. van de Geijn
Department of Computer Science
&
Oden Institute
The University of Texas at Austin
Austin, TX 78712, U. S. A.
<EMAIL_ADDRESS>
###### Abstract

Matrix libraries often focus on achieving high performance for problems
considered to be either ``small'' or ``large'', as these two scenarios tend to
respond best to different optimization strategies. We propose a unified
technique for implementing matrix operations like general matrix
multiplication (gemm) that can achieve high performance for both small and
large problem sizes. The key is to fuse packing – an operation that copies
data to a contiguous layout in memory and which is critical for large matrix
performance – with the first computational ``pass'' over that data. This
boosts performance across the problem size spectrum. As a result, tuning
general-purpose libraries becomes simpler since it obviates the need to
carefully express and parameterize logic that chooses between a ``small
matrix'' strategy and a ``large matrix'' strategy. A prototype implementation
of the technique built with the BLAS-like Library Instantiation Software
(BLIS) framework is described and performance on a range of architectures is
reported.

##  1 Introduction

The Basic Linear Algebra Subprograms (BLAS) [11, 4, 3] interface has had a
profound impact on scientific software development. It is now also of great
importance to fields like machine learning and data analytics. By coding
applications in terms of the BLAS, portable high performance can be achieved.
For this reason, whenever a new high-performance computer architecture
arrives, the instantiation of this interface is a high priority.

Historically, it was expected that vendors leverage their expertise with the
architecture to create proprietary matrix libraries, with key components coded
in assembly language. IBM's algorithms and architectures approach demonstrated
that by co-designing architectures, compilers, and libraries it was possible
to achieve high performance with implementations coded in a high-level
language (Fortran) [1]. This inspired a number of open-source efforts to
provide portable implementations of the BLAS, including the Automatic Tuned
Linear Algebra Software (ATLAS) [16], the GotoBLAS [6, 7], the OpenBLAS [23]
(a fork of GotoBLAS), and the BLAS-like Library Instantiation Software (BLIS)
[20, 22] upon which this paper implements its approach. An added benefit of
BLIS is that it supports an analytical model for determining blocking
parameters so that autotuning can be avoided [12].

Across all publically available efforts towards matrix libraries, a
fundamental problem that complicates the implementation of matrix-matrix
operations, a. k. a. level-3 BLAS, is that packing to improve data locality,
which is necessary for high performance when targeting large matrix sizes,
actually _impedes_ high performance for smaller matrix sizes. It has been
thought that this is an inherent problem that can only be solved by
implementing separate code paths for small and large matrix sizes and then
selecting one of them based on the problem size characteristics. In this
paper, we provide preliminary evidence that this conventional wisdom _may be
wrong_ : the two code paths can be unified in a way that mostly preserves the
benefits of both. This is achieved by integrating packing – which is optional
in the small code path – more tightly with the computation. Importantly, the
problem sizes where optional packing should be turned on or off can be more
easily justified and encoded.

##  2 Goto's Algorithm and its Instantiation in BLIS

Goto's algorithm [8, 6, 7] was first developed for CPUs with two levels of
cache and continues to be the algorithm that underlies most if not all vendor
and open-source implementations of the level-3 BLAS. This section gives a
high-level description of this algorithm and its instantiation in BLIS.

###  2.1 Goto's algorithm for large matrices

Figure 1:  Goto's algorithm for gemm as five loops around the microkernel.
This diagram, which is often used when explaining the fundamental techniques
that underly the BLIS implementation of gemm, was modified from a similar
image first published by [21] and is used with permission.

Goto's algorithm structures a prototypical gemm, $C:=AB+C$, where $A,B$, and
$C$ are matrices of size $m\times k,k\times n$ and $m\times n$ respectively,
as five loops around the update of a small submatrix of $C$ called the
microtile, as illustrated in Figure 1. We only give the highlights here since
the algorithm and this picture have been explained in many previously-
published papers. In this discussion and the figure, $m_{R}$ and $n_{R}$
denote register blocking sizes, while $m_{C}$, $n_{C}$, and $k_{C}$ denote
cache blocking parameters. At the core is the _microkernel_ , which updates an
$m_{R}\times n_{R}$ microtile of $C$ by multiplying a $m_{R}\times k_{C}$
_micropanel_ of $A$ by a $k_{C}\times n_{R}$ micropanel of $B$. On a typical
architecture, the microtile of $C$ is kept in registers while the micropanels
of $B$ and $A$ are streamed from the L1 and L2 caches, respectively. Blocks of
$A$ and row panels of $B$ are rearranged (packed) at strategic points in the
algorithm to allow memory access with unit stride111This refers to loading
consecutive memory addresses, which allows most CPUs to execute in the fewest
number of cycles. as well as to align micropanels from $A$ and $B$ so that
they can fit into their designated levels of cache as indicated in Figure 2.
We refer to this rearranged storage as ``packed memory'' in contrast to
``unpacked memory'' used to store the original $A$ and $B$ matrices.

Unpacked data can map poorly to cache associativity sets:
---
Packed data maps to sets predictably and equitably:
Figure 2:  A schematic illustrating how the storage format of a micropanel in
memory affects how cache lines map to the associativity sets of a cache. Top:
Unpacked columns with arbitrary leading dimensions tend to cause more cache
lines to be mapped to fewer sets of the cache, leading to inefficient use and
unnecessary evictions. Bottom: The same submatrix stored contiguously causes
cache lines to be mapped across all sets, leading to fewer evictions and, by
proxy, fewer subsequent cache misses. A more rigorous analysis can be found in
[12].

###  2.2 BLIS's refactoring of Goto's algorithm

The BLIS implementation of Goto's algorithm recognizes that as long as the
microkernel is expressed with assembly code222This can take the form of so-
called extended inline assembly code in addition to pure assembly code. Vector
intrinsics _may_ also work, depending on the compiler and instruction set
being emitted., high performance can be achieved even if all remaining parts
of the algorithm above the microkernel (including packing) are written in C.
This reduces how much code must be customized for the gemm operation. It also
allows other matrix-matrix operations supported by the BLAS (level-3 BLAS),
such as Hermitian matrix-matrix multiplication (hemm), Hermitian rank-$k$
update (herk), triangular matrix-matrix multiplication (trmm), and triangular
solve with multiple right-hand sides (trsm), to employ the same microkernel
[20, 22]. This contrasts with the original GotoBLAS implementation, inherited
by OpenBLAS, where the two loops around the microkernel form what we call a
_macrokernel_ that must be customized in assembly code for different gemm-like
operations [7].333When the assembly region extends to encompass the
macrokernel, an OpenBLAS-like implementation may choose to maintain separate
macrokernels for each operation _or_ insert conditional logic into a single
macrokernel that allows the code to handle multiple similar operations (e.g.
syrk and herk). The former case yields more regions of assembly code with only
minor differences between them while the latter results in less assembly code
that is nonetheless more difficult to decipher due to its embedded branching.

An important detail is how to handle situations where the matrix dimensions
are not whole multiples of $m_{R}$, $n_{R}$, and/or $k_{C}$. These so-called
``fringe'' or edge cases. During packing, BLIS pads micropanels with zeroes
when a fringe case is encountered. And since the microkernel author only needs
to target microtiles of one size – $m_{R}\times n_{R}$ – the job of writing
and optimizing microkernels becomes much simpler. The result is an easier-to-
develop and easier-to-maintain code base at the cost of a minor decrease in
performance for certain problem sizes [18].

###  2.3 SUP: Supporting small-ish matrices

More recently, projects and implementations like LIBXSMM [9] and BLASFEO [5]
have sought to improve matrix-matrix performance for small-sized problems.
Typically, these solutions either skip packing or require data to be pre-
packed to exploit the fact that for small problems, matrices $A$ and $B$ have
a chance to fit into the L2 (or even L1) cache in their entirety, thus
avoiding the $\mathcal{O}(mk+nk)$ cost of packing.

BLIS's current approach to cases where at least one matrix dimension is small
is referred to as Skinny/UnPacked (SUP). It combines the following techniques:

  * •

It skips packing.

  * •

Rather than isolating architecture specifics only in the microkernel, it
employs a _millikernel_ that absorbs the first loop around the microkernel
into the kernel primitive in an effort to reduce the frame stack cost (that
is, memory overhead due to subroutine calls). This millikernel is written in a
manner similar to that of a corresponding microkernel (i.e., in assembly
code).

  * •

When a millikernel encounters a fringe case during its last iteration, it
dispatches a helper microkernel that specializes in that size. Alternatively,
a set of fringe cases may be called in sequence to emulate the net effect of a
single, larger microkernel. For example, let us assume the size of the
microtile is $6\times 8$ and the millikernel encounters a fringe that is
$5\times 8$. This could result in a call to a single helper microkernel call
that updates a $5\times 8$ microtile, or it could result in a call to a
$3\times 8$ helper microkernel followed by one that computes the remaining
$2\times 8$ part.

  * •

If a fringe case is encountered in the loop that _surrounds_ the millikernel
(i.e., the 2nd loop around the microkernel), a more specialized millikernel is
dispatched – one that dispatches to a different subset of helper microkernels.
Building on the previous example, if the millikernel is called on an $m\times
5$ submatrix, a helper millikernel that operates on at most 5 columns of $B$
is dispatched, which will then loop over a microkernel-like region of code
targeting $6\times 5$ and eventually, if needed at the fringe, dispatch helper
microkernels whose $m_{R}$ values are less than 6.

While these changes to the conventional BLIS code path are conceptually
simple, the examples above suggest that a nontrivial amount of code is
required to support it. Specifically, if a base-2 ``spanning'' set of kernels
is employed, based on a $6\times 8$ microtile, the following millikernels and
microkernels would be needed:

$6\times 8(\textsc{mm})$ | $6\times 4(\textsc{hm})$ | $6\times 2(\textsc{hm})$ | $6\times 1(\textsc{hm})$
---|---|---|---
$4\times 8(\textsc{m}\mu)$ | $4\times 4(\textsc{h}\mu)$ | $4\times 2(\textsc{h}\mu)$ | $4\times 1(\textsc{h}\mu)$
$2\times 8(\textsc{m}\mu)$ | $2\times 4(\textsc{h}\mu)$ | $2\times 2(\textsc{h}\mu)$ | $2\times 1(\textsc{h}\mu)$
$1\times 8(\textsc{m}\mu)$ | $1\times 4(\textsc{h}\mu)$ | $1\times 2(\textsc{h}\mu)$ | $1\times 1(\textsc{h}\mu)$

Here, mm denotes the main millikernel that is called from the 2nd loop, and
$\textsc{m}\mu$ denotes its helper microkernels. Similarly, the hm labels
denote helper millikernels, each of which calls its own helper microkernels
labeled by $\textsc{h}\mu$.

###  2.4 Combining the two code paths

| DGEMM on Xeon E5-2690 |  | DGEMM on AWS Graviton 3
---|---|---|---
GFLOPS/sec |  |  GFLOPS/sec |
| $m=n=k$ |  | $m=n=k$
Figure 3: Performance of: BLIS's conventional gemm and gemm implemented with
the SUP approach (gemmsup) on the Intel Xeon E5-2690 (left) and AWS Graviton 3
with Arm SVE (right) architectures. Note the different cross-over points for
each system.

Figure 3 reports performance of the conventional BLIS gemm and the gemm
implemented via the SUP approach (gemmsup) on two different architectures.
Clearly, SUP outperforms the conventional algorithm for small problem sizes.
However, as the problem size becomes large, data would unavoidably spill out
of the cache, degrading SUP's performance and creating periodic spikes down on
the Xeon E5-2690 processor. Importantly, to roll out solutions to a matrix
library, the cross-over point must be determined and encoded. The heuristic
for this is complicated by the fact that it is not only a function of the
problem size but also of the row and column strides for matrices $A$ and $B$.
Additionally, there may be an architecture-dependent range of sizes where both
algorithms suffer performance degradation (e.g., around $m=n=k=200$ on the
Xeon E5-2690). In such a region, unaligned SUP may spill data out of the
cache, while the $\mathcal{O}(mk+nk)$ packing cost for the conventional
algorithm is non-negligible. We will soon show that fusing the packing with
computation can address these issues.

##  3 A Unified Approach

It would appear that the packing process inherently imposes unreasonable
overhead for small- to medium-sized matrix cases. We now discuss how that is
not always so.

###  3.1 Interleaving of packing and computing

for  _ $\cdots$ (the 5 and 4 loops around the microkernel proceed as before
but without packing) $\cdots$ _ do

for  _first iteration of 3 loop around the microkernel_ do

for  _the first iteration of 2 loop around the microkernel_ do

Pack the first micropanel of $B_{p,j}\rightarrow\widetilde{B}_{p,j}$;

for  _each iteration of 1 loop around the microkernel_ do

Pack current micropanel of $A_{i,p}\rightarrow\widetilde{A}_{i,p}$;

Call the microkernel with packed micropanels $\widetilde{A}_{i,p}$ and
$\widetilde{B}_{p,j}$;

end for

Upon completion, $A_{i,p}$ is left packed in $\widetilde{A}_{i,p}$;

end for

for  _remaining iterations of 2 loop around microkernel_ do

Pack current micropanel of $B_{p,j}\rightarrow\widetilde{B}_{p,j}$;

for  _each iteration of 1 loop around the microkernel_ do

Call the microkernel with packed micropanels $\widetilde{A}_{i,p}$ and
$\widetilde{B}_{p,j}$;

end for

end for

Upon completion, $B_{p,j}$ is left packed in $\widetilde{B}_{p,j}$;

end for

for  _remaining iterations of 3 loop around microkernel_ do

for  _the first iteration of 2 loop around the microkernel_ do

for  _each iteration of 1 loop around the microkernel_ do

Pack current micropanel of $A_{i,p}\rightarrow\widetilde{A}_{i,p}$;

Call the microkernel with packed micropanels $\widetilde{A}_{i,p}$ and
$\widetilde{B}_{p,j}$;

end for

Upon completion, $A_{i,p}$ is left packed in $\widetilde{A}_{i,p}$;

end for

for  _remaining iterations of 2 loop around microkernel_ do

for  _each iteration of 1 loop around the microkernel_ do

Call the microkernel with packed micropanels $\widetilde{A}_{i,p}$ and
$\widetilde{B}_{p,j}$;

end for

end for

end for

end for

Algorithm 1 Algorithm with interleaved packing and computation.

The basic idea behind unifying is simple: Modify Goto's algorithm so that
elements of the micropanels of $A$ and $B$ are packed _just before_ their
first use by a call to the microkernel. More precisely, consider that Goto's
algorithm has proceeded to the first iteration of the 4th loop, in which a
column-panel of $A$ is multiplied by a row-panel of $B$, except let us assume
that neither is packed:

.

This operation can be implemented by modifying Goto's algorithm as described
in Algorithm 1. It allows the microkernel from the conventional BLIS
implementation to be used without modification.

The benefit of this interleaving is that after packing, micropanels of
$\widetilde{B}_{p,j}$ and $\widetilde{A}_{i,p}$ are still in the L1 cache when
used for the first time by the microkernel, while their eventual migration
back to the L3 and L2 caches, respectively, are (likely) masked by
computation. This can be expected to improve performance in general, but in
particular for small matrices.

###  3.2 FIP: Fused packing in the microkernel

While the simple solution mentioned in the previous subsection should reduce
the net execution time of the microkernel for the first time a packed
micropanel is involved in computation, packed data still moves from registers
to the L1 cache and back. This ``reflowing'' of data adds cycles to the total
execution time, and can sometimes trigger a performance penalty for _read-
after-write_ (RAW) data hazards especially on some x86_64 architectures.

Our solution is to fuse individual packing instructions into the microkernel
itself. We call this technique ``fused-in packing'' (FIP).

Conceptually, this is similar to the previous approach, except that the unit
of data around which the packing and computation are interleaved is reduced
from whole micropanels of $A$ and $B$ to a mere handful of elements. This has
at least two benefits: (1) Data that is loaded into registers in the normal
course of packing is reused immediately for useful computation; (2) The cost
of loading data from the cache into registers is partially hidden by the
computation on previously-loaded values.

Implementing FIP requires instantiating four cases of the microkernel:

  * •

The conventional microkernel, where microtiles from both sides are already
packed.

  * •

A microkernel where $\widetilde{A}_{i,p}$ is already packed but $B_{p,j}$ is
unpacked (and thus needs to be packed).

  * •

A microkernel where $\widetilde{B}_{p,j}$ is already packed but $A_{i,p}$ is
unpacked (and thus needs to be packed).

  * •

A microkernel where both $A_{i,p}$ and $B_{p,j}$ are unpacked.

Each of these might be encountered while computing $C:=AB+C$ depending on the
sizes of the matrices.

###  3.3 Options in fused packings

There are some situations in which there is little benefit from packing
$B_{p,j}$ and/or $A_{i,p}$. Consider a gemm call where $C$ is $m\times n$, $A$
is $m\times k$, and $B$ is $k\times n$. If $n\leq n_{R}$, then the 2nd loop
around the microkernel is only executed once, and hence packing $A_{i,p}$
would yield no benefit. And if $m\leq m_{R}$ then $B_{p,j}$ is similarly not
reused. The resulting spectrum of options can be summarized as

| $n\leq n_{R}$ | $n>n_{R}$
---|---|---
$m\leq m_{R}$ | no packing | pack $A_{i,p}$
$m>m_{R}$ | pack $B_{p,j}$ | pack $A_{i,p}$ and $B_{p,j}$

This provides a decent heuristic from which further tuning may be explored.

Another observation concerns the packing of $A_{i,p}$. Let us assume that $A$
is stored in column-major order. If

${\rm csb}(A)\times k_{C}\leq\rm{(L2~{}cache~{}size~{}in~{}bytes)},$

where ${\rm csb}(A)$ equals the stride in bytes between elements in a row of
$A$, then the unpacked storage of $A_{p}$ will never cause an L2 cache spill –
that is, elements of $A_{p}$ will not be evicted from the L2 cache by other
elements. This constitutes another case where we can skip the packing of $A$.
If $A$ is stored in row-major order, the condition becomes:

$m_{C}\times{\rm rsb}(A)\leq\rm{(L2~{}cache~{}size~{}in~{}bytes)},$

where ${\rm rsb}(A)$ equals the stride in bytes between elements in a column
of $A$. The details behind these inequalities go beyond the scope of this
paper and require an understanding of the results by [12].

| DGEMM on Xeon E5-2690 |  | DGEMM on AWS Graviton 3
---|---|---|---
GFLOPS/sec |  |  GFLOPS/sec |
| $m=n=k$ |  | $m=n=k$
Figure 4: Performance improvements observed on the Intel® Xeon E5-2690 (top)
and the AWS Graviton 3 (bottom) processors from interleaving packing and
computing and optionally fusing their kernels (FIP/gemmfip).

Figure 4 illustrates how interleaving packing and computing as Algorithm 1
benefits the performance of gemm and how deploying the FIP technique further
improves it. Additional experiments are available in Section 4.

###  3.4 Coding effort

⬇

#define cond_inst_false(_1)

#define cond_inst_true(_1) _1

#define kernel_def(pa, pb) \

\

void fused_kernel_## pa ##_## pb ( /* func params */ ) \

{ \

__asm__ volatile \

( \

/* ... */ \

\

/* Multiply and accumulate: ymm0 holds A, ymm2 holds B and ymm4 holds C
fractions */ \

"vfmadd231pd %ymm0, %ymm2, $ymm4 \n\t" \

\

/* Optionally store ymm2 to the packing space indicated by rdx. */ \

cond_inst_## pb ( "vmovapd %ymm2, 16(%rdx) \n\t" ) \

\

/* ... */ \

) \

}

kernel_def(true, true)

kernel_def(true, false)

kernel_def(false, true)

kernel_def(false, false)

Algorithm 2 A code sample for multi-instantiating SUP-based kernels with
packing instructions conditionally (depending on the values of macro arguments
pa and pb) fused in to handle all the four cases required by the FIP approach.

We implemented FIP kernels based on the microkernel portion of the SUP code
path in BLIS introduced in Section 2.3. Though Sections 3.2 and 3.3 imply that
the approach requires handling all four cases within one set of microkernels,
most of the redundancies between the four cases can be tamed with some careful
refactoring. BLIS kernels are written with inlined assembly code in C.
Meanwhile, fused packing merely requires data from registers to be moved to
memory with hard-codable strides. These two facts allowed us to leverage the C
preprocessor to expand a single code template into the four specializations in
a fashion shown in Algorithm 2.

###  3.5 Multithreading

Figure 5: An illustration of how multithreading can be added to our FIP
approach. Each thread starts from a different unpacked micropanel of $A_{p}$
and packs the micropanel to its designated space in $\widetilde{A}_{i,p}$
while performing the corresponding update of a microtile of $C:C_{i,j}$. Once
all threads have finished packing their
$\left.\frac{m_{C}}{m_{R}}\middle/n_{\mathrm{thr}}\right.$ micropanels, _one_
synchronization occurs so they can use all micropanels in $\widetilde{A}_{p}$.
This collaborative packing and computing only happen when each thread is
working on its first $B_{p,j}$ micropanel within the 2 loop around the
microkernel.

We now take a brief look at the multithreading potential of the FIP technique,
which is left for future research on this approach.

BLIS's refactoring of Goto's algorithm conveniently exposes five loops coded
in C99 in which multithreading can be introduced [14]. In practice,
parallelism is usually gained in the 2nd loop around the microkernel to let
multiple threads operate on the same $m_{C}\times k_{C}$ tile from $A$ or in
the 3rd loop around the microkernel to let multiple threads operate on the
same $k_{C}\times n_{C}$ panel from $B$. On a multi-core chip with non-uniform
memory access (NUMA) architecture, cores usually share L2 caches within each
NUMA, making it reasonable to parallelize over the 2nd loop within and the 3rd
loop across NUMA nodes.

With this insight, if we want to apply our insight to the consumption of tiles
from $A$, it seems appropriate to let each thread work on the unpacked memory
first and collaboratively write the microtiles they have loaded to the packing
space for later consumption. Since threads begin referring to the packed data
only from their second iteration of the 2nd loop around the microkernel, for
$A_{p}$ it is only required that all threads working on the same $m_{C}\times
k_{C}$ tile from $A$ synchronize _once_ after finishing their first iteration
there. This synchronization cost is identical to collaboratively packing
everything from $A$ beforehand.

Furthermore, if one still wants to ensure each tile from unpacked memory is
accessed only once, we can even change the order when each iteration of the
first loop around the microkernel gets executed. Letting each thread enter the
microkernel at a different microtile of $C$ allows the requested $m_{C}\times
k_{C}$ to become readily packed after each thread has finished their first
$\left.\frac{m_{C}}{m_{R}}\middle/n_{\text{thr}}\right.$ iterations, as is
illustrated in Figure 5. Here $n_{\text{thr}}$ is denoted as the number of
threads working on the same $m_{C}\times k_{C}$ tile from $A$ and $m_{C}$ is
supposed to be multiples of $m_{R}n_{\text{thr}}$.

On the $B$ side, in each iteration within the 2nd loop around the microkernel,
the microtile from $B$ is reused from the thread-private L1 cache and no
distinct gain in performance can be expected from reusing micropanels of $B$
from the L3 cache. In addition, thread synchronization across NUMA sockets is
likely to cause nontrivial overhead. This suggests that it is better to let
each thread pack its relevant panels into a separate space despite the modest
storage redundancy. The threads then flow data from their private L1 cache to
the registers without interfering with each other.

##  4 Performance results

We now illustrate the benefits of the discussed technique on a broad range of
architectures, for a single core.

###  4.1 Experimental setup

We implemented gemm with the FIP approach (gemmfip) based on BLIS (Release
0.9.0). Our implementation can be integrated back into BLIS as an alternative
backend for gemm through its _sandbox_ interface. The macrokernel structure of
our implementation is essentially the same as BLIS's refactoring of Goto's
algorithm in Figure 1 with a few additional lines to handle kernel selection
for the four cases mentioned in Section 3.2.

A non-trivial modification here is for x86_64 kernels. As we have already
mentioned in Section 2.3, to minimize costs associated with pushing to the
frame stack, BLIS's SUP kernels include the first loop around the microkernel
(i. e. the millikernel) in the assembly region of the code. In our
implementations, this was also incorporated for x86_64-based architectures
since they have a limited number of general-purpose registers compared to
Arm®, Power ISA, and RISC-V®, making it harder for compilers to transition
between microkernel calls without interfering with memory consistency or
causing frame stack RAW data hazards.

We developed microkernels to support double-precision gemmfip on 3
architectures: Intel AVX2, Arm NEON and Arm SVE. Performance experiments were
performed on five different processors: Intel® Xeon E5-2690 (3500 MHz,
compiler: GCC 10.2), AMD Epyc™ 7R32 with x86_64 AVX2 architecture (3300 MHz,
compiler: GCC 11.3), AWS Graviton 2 with Arm NEON (2500 MHz, compiler: Clang
14.0), AWS Graviton 3 with Arm SVE (2600 MHz, compiler: Clang 14.0), and Apple
M2 with Arm NEON (3490 MHz, compiler: Apple Clang 14.0).

All our experiments used a single core. The uppermost depicted $y$-axis value
represents the single-core theoretical peak performance for each processor
tested.

###  4.2 Evaluation

| Xeon E5-2690 | Epyc 7R32
---|---|---
GFLOPS/sec w/ $\text{LDim}=m$ |  |
GFLOPS/sec w/ $\text{LDim}=2000$ |  |
GFLOPS/sec w/ $\text{LDim}=m$ |  |
GFLOPS/sec w/ $\text{LDim}=2000$ |  |
| $m=n=k$ | $m=n=k$
Figure 6: Performance on the Intel® Xeon E5-2690 (left) and AMD Epyc™ 7R32 with x86_64 AVX2 (right) architectures.  | AWS Graviton 2 | AWS Graviton 3 | Apple M2
---|---|---|---
GFLOPS/sec w/ $\text{LDim}=m$ |  |  |
GFLOPS/sec w/ $\text{LDim}=2000$ |  |  |
GFLOPS/sec w/ $\text{LDim}=m$ |  |  |
GFLOPS/sec w/ $\text{LDim}=2000$ |  |  |
| $m=n=k$ | $m=n=k$ | $m=n=k$
Figure 7:  Performance on various Arm architectures. Left: AWS Graviton 2 with
Arm NEON. Middle: AWS Graviton 3 with Arm SVE architecture. Right: Apple M2
with Arm NEON.

In Figures 7 and 6 (top), we report the performance of the conventional BLIS
code path, its SUP code path, and gemmfip. In the top graphs, the leading
dimension (LDim) – that is, the stride between logically adjacent elements
within in a single row – equals the (row) dimension of the matrices. In
Figures 7 and 6 (center-top), the leading dimension equals $2000$, which can
negatively affect performance, particularly when computing on unpacked data.

We see that the curve for gemmfip uniformly matches or outperforms the
conventional BLIS and SUP paths while smoothly bridging the ``medium-sized''
region where both underperform. This result supports the claim that the new
method can provide a unified approach that yields high performance across the
range of problem sizes without having to deploy heuristics to determine a
crossover point between the paths.

###  4.3 Comparison with other BLAS implementations

We now compare how our implementation of gemmsup performs against other matrix
libraries, including BLIS itself, OpenBLAS, and numerous vendor-specific BLAS
implementations.

On the x86_64 side, Intel's Math Kernel Library (MKL) is akin to the ``gold
standard.'' Since MKL is a closed-source implementation, it is difficult to
discern what allows that code to perform so well. In Figure 6 (center-bottom
and bottom), our implementation handily outperforms OpenBLAS and matches or
surpasses BLIS, especially for medium and even large problem sizes. On AMD's
Zen™-microarchitecture-based processors, BLIS is integrated into the AMD
Optimizing CPU Libraries (AOCL) with vendor-side tunings. AMD also has an
open-source fork of BLIS known as AMD BLIS, which we have built and tested
apart from AOCL. On that processor, gemmfip essentially matches AOCL's
performance (with the exception of a narrow range from about $150$ to $250$
when $\mathrm{LDim}=m$) while handily outperforming BLIS, AMD BLIS, and
OpenBLAS. Meanwhile, MKL yields mediocre (and inconsistent) performance on the
7R32. We imagine that Intel made a conscious decision to throttle the
performance of MKL when running on Epyc hardware in hopes of discouraging
their users from switching to their competitor's products.

For Arm processors, Arm Performance Libraries (ArmPL) provides a state-of-the-
art matrix library solution. Comparisons of gemmfip against ArmPL, OpenBLAS,
and BLIS are plotted in Figure 7 (center-bottom and bottom). On the AWS
Graviton 2, gemmfip outperforms ArmPL and slightly underperforms OpenBLAS.
This defect is presumably because our microkernels for unpacked memories still
leave room for improvement on this specific microarchitecture, as can be
deduced from the fact that gemmsup performance using a similar microkernel
lags far behind. For Apple's M2 processor, its Accelerate framework was, upon
closer inspection, found to be using a hidden co-processor [2] shared by the
whole chip instead of the Arm NEON pipelines, and the peak performance turned
out to be around 380 GFLOPS/sec for $m=n=k=\left(\text{LDim}\textrm{ of
}A,B,\textrm{and }C\right)\approx 800$, regardless of threading options given.
This difference on the hardware side makes it impossible to measure an
isolated single-core throughput. Therefore, the Accelerate framework is
producing curves in Figure 7 that exceed the $y$-axis limits of the graphs.
Finally, on the Arm SVE architecture of the AWS Graviton 3 processor, gemmfip
yields the highest throughput and the best consistency among all tested
libraries, demonstrating its value in this relatively nascent architecture.

##  5 Conclusion

We have described and demonstrated the benefits of reorganizing the memory
access patterns associated with the Goto algorithm, as refactored by BLIS, so
that small- and large-size problems are handled within a unified code path.
This novel approach integrates the packing instructions for $A$ and/or $B$
into the microkernel, where they may be templatized and conditionally
activated at compile-time. Fusing the packing and computation in this manner
creates a spectrum of algorithms, where the SUP approach (no packing) and the
Goto algorithm (full ``classic'' packing) each appear as special cases. A
runtime logic then tracks whether there is a reuse of $A_{i,p}$ and $B_{p,j}$
micropanels to determine which variant should be deployed to specific problem
sizes. Since the result yields uniformly high performance insensitive to the
switching points in that runtime logic, our unified FIP approach appears to
largely obviate the need for empirical crossover points or other runtime
heuristics.

We have demonstrated the method's benefits by creating implementations for
multiple architectures. Only on the Intel Xeon E5-2690 processor does another
implementation – Intel's MKL – outperform FIP. Given MKL's reputation for
achieving extremely high performance, we believe this clearly illustrates the
importance of the work. We expect these techniques to be adopted by BLIS and
other libraries in the near future, thus further positively impacting the user
community, regardless of which library they use.

The results in this paper suggest future work that will have an extended
impact. The most obvious is that the techniques can be applied to other
precisions and other matrix-matrix multiplications (level-3 BLAS). In
addition, a body of papers shows how BLIS's refactoring of Goto's algorithm
can be used to attain high performance and/or reduce the development effort
for various gemm-like operations:

  * •

In the work by [19], BLIS's 1m method leverages the real-domain microkernel to
implement complex-domain matrix-matrix multiplication operations by cleverly
encoding the definition of complex scalar multiplication within the packing
stage of Goto's algorithm.

  * •

The work by [10] shows how Strassen's algorithm can already attain high
performance for rank-$k$ updates and relatively small matrices.

  * •

[17] give a high-performing implementation for solving the $k$-nearest
neighbor problem by fusing computation into Goto's algorithm.

  * •

In the work by [15], it is reasoned and demonstrated that Goto's algorithm is
but one algorithm in a larger family of algorithms, the Multilevel Optimized
Matrix-matrix Multiplication Sandbox (MOMMS) family. The idea is that as the
speed ratio between CPU arithmetics and memory access becomes worse in the
future, blocking for caches must be modified.

  * •

[13] uses BLIS's refactoring of Goto's algorithm to instead implement tensor
contractions, yielding TBLIS, by recognizing that rearrangement of data to
cast tensor contraction as a matrix multiplication can be incorporated into
packing.

Insights demonstrated in this paper can potentially accord base-performance
benefits for all these algorithms.

Finally, if the strategy in Section 3.5 could bring our insight's unifying
effect to the multithreading regime, perhaps it could extend to and yield
speedup on GPUs whose intermediate-level memory operates in a cache-like way
or as a register-controlled scratchpad memory (SPM)444For GPUs whose L1-level
storage is an SPM with controllable direct memory access (DMA), it is expected
that packing via DMA will provide better performance. as well.

##  6 Code Availability

Our code is developed under a fork of BLIS available at
https://github.com/xrq-phys/blis. It may be enabled as an ``sandbox'' that is
optionally integrated into the library by configuring BLIS as:

./configure -s gemmfip -t none x86_64

for x86_64 microarchitectures or:

./configure -s gemmfip -t none arm64

for Arm hardware.

## Acknowledgements

We thank Prof. S. Todo for supervising RuQing Xu's research and providing
access to various architectures. We also thank members of the BLIS community
for their input.

RuQing Xu is funded by The University of Tokyo's GSGC scholarship. The
researchers at The University of Texas at Austin are funded in part by the
National Science Foundation (Award CSSI-2003921) and gifts from AMD, Arm, and
Oracle.

Any opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the National Science Foundation.

## References

  * [1] R.C. Agarwal, F.G. Gustavson and M. Zubair  ``Exploiting functional parallelism of POWER2 to design high-performance numerical algorithms''  In _IBM Journal of Research and Development_ 38.5, 1994
  * [2] Peter Cawley  ``Apple AMX'', https://github.com/corsix/amx, 2022
  * [3] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling and Iain Duff  ``A Set of Level 3 Basic Linear Algebra Subprograms''  In _ACM Trans. Math. Soft._ 16.1, 1990
  * [4] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling and Richard J. Hanson  ``An Extended Set of FORTRAN Basic Linear Algebra Subprograms''  In _ACM Trans. Math. Soft._ 14.1, 1988
  * [5] Gianluca Frison et al.  ``BLASFEO: Basic Linear Algebra Subroutines for Embedded Optimization''  In _ACM Trans. Math. Softw._ 44.4, 2018  DOI: 10.1145/3210754
  * [6] Kazushige Goto and Robert Geijn  ``Anatomy of High-Performance Matrix Multiplication''  In _ACM Trans. Math. Soft._ 34.3: Article 12, 25 pages, 2008
  * [7] Kazushige Goto and Robert Geijn  ``High-performance implementation of the level-3 BLAS''  In _ACM Trans. Math. Soft._ 35.1, 2008  DOI: http://doi.acm.org/10.1145/1377603.1377607
  * [8] Kazushige Goto and Robert A. Geijn  ``On Reducing TLB Misses in Matrix Multiplication'', 2002
  * [9] Alexander Heinecke, Greg Henry, Maxwell Hutchinson and Hans Pabst  ``LIBXSMM: Accelerating Small Matrix Multiplications by Runtime Code Generation''  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC'16)_ , 2016, pp. 981–991  DOI: 10.1109/SC.2016.83
  * [10] Jianyu Huang, Tyler M. Smith, Greg M. Henry and Robert A. Geijn  ``Strassen's Algorithm Reloaded''  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage Analysis (SC'16)_ , 2016  URL: http://dl.acm.org/citation.cfm?id=3014904.3014983
  * [11] C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. Krogh  ``Basic Linear Algebra Subprograms for Fortran Usage''  In _ACM Trans. Math. Soft._ 5.3, 1979
  * [12] Tze Meng Low, Francisco D. Igual, Tyler M. Smith and Enrique S. Quintana-Ortí  ``Analytical Modeling Is Enough for High-Performance BLIS''  In _ACM Trans. Math. Softw._ 43.2, 2016  DOI: 10.1145/2925987
  * [13] Devin Matthews  ``High-Performance Tensor Contraction without Transposition''  In _SIAM J. Sci. Comput._ 40.1, 2018  DOI: 10.1137/16M108968X
  * [14] Tyler M Smith et al.  ``Anatomy of high-performance many-threaded matrix multiplication''  In _Proceedings of the 28th International Parallel and Distributed Processing Symposium (IPDPS)_ , 2014, pp. 1049–1059
  * [15] Tyler M. Smith and Robert A. Geijn  ``The MOMMS Family of Matrix Multiplication Algorithms''  In _CoRR_ abs/1904.05717, 2019  arXiv: http://arxiv.org/abs/1904.05717
  * [16] R. Clint Whaley, Antoine Petitet and Jack J. Dongarra  ``Automated empirical optimizations of software and the ATLAS project''  In _Parallel Computing_ 27.1–2, 2001
  * [17] Chenhan D. Yu et al.  ``Performance Optimization for the K-Nearest Neighbors Kernel on x86 Architectures''  In _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC'15)_ , 2015  DOI: 10.1145/2807591.2807601
  * [18] Field G. Zee  ``BLIS Performance'', https://github.com/flame/blis/blob/master/docs/Performance.md, 2019
  * [19] Field G. Zee  ``Implementing High-Performance Complex Matrix Multiplication via the 1m Method''  In _SIAM J. Sci. Comput._ 42.5, 2020, pp. C221–C244  URL: https://doi.org/10.1137/19M1282040
  * [20] Field G. Zee and Robert A. Geijn  ``BLIS: A Framework for Rapidly Instantiating BLAS Functionality''  In _ACM Trans. Math. Softw._ 41.3, 2015  DOI: 10.1145/2764454
  * [21] Field G. Zee and Tyler M. Smith  ``Implementing high-performance complex matrix multiplication via the 3m and 4m methods''  In _ACM Trans. Math. Soft._ 44.1, 2017, pp. 7:1–7:36
  * [22] Field G. Zee et al.  ``The BLIS Framework: Experiments in Portability''  In _ACM Trans. Math. Softw._ 42.2, 2016  DOI: 10.1145/2755561
  * [23] Xianyi Zhang  ``OpenBLAS'', http://www.openblas.net, 2011


11institutetext: Univ Paris Diderot, Sorbonne Paris Cité, LIAFA, UMR 7089
CNRS, Paris, France
11email<EMAIL_ADDRESS>
# Implementing Computations
in Automaton (Semi)groups

Ines Klimann  Jean Mairesse  Matthieu Picantin

###### Abstract

We consider the growth, order, and finiteness problems for automaton
(semi)groups. We propose new implementations and compare them with the
existing ones. As a result of extensive experimentations, we propose some
conjectures on the order of finite automaton (semi)groups.

###### Keywords:

a

utomaton (semi)groups, growth, order, finiteness, minimization

##  1 Introduction

_Automaton (semi)groups_ — short for semigroups generated by Mealy automata or
groups generated by invertible Mealy automata — were formally introduced a
half century ago (for details, see [10, 7] and references therein). Over the
years, important results have started revealing their full potential. For
instance, the article [9] constructs simple Mealy automata generating infinite
torsion groups and so contributes to the Burnside problem, and, the article
[5] produces Mealy automata generating the first examples of (semi)groups with
intermediate growth and so answers the Milnor problem.

The classical decision problems have been investigated for such (semi)groups.
The word problem is solvable using standard minimization techniques, while the
conjugacy problem is undecidable [16]. Here we concentrate on the problems
related to growth, order, and finiteness.

(-3,-1)(5,4) [1](-4,0.5)a [3](4,0.5)b [2](0,3)c ac32 cb32 c1123 ab1121
[.7]ba112232

(-4,-2.5)(4.5,1) [1](-4,0)a [2](0,0)b [3](4,0)c a132231 b31 ba13 bc22 cb122331

Figure 1: A Mealy automaton and its dual

To illustrate, consider the two Mealy automata of Fig. 1. They are dual, that
is, they can be obtained one from the other by exchanging the roles of
stateset and alphabet. A (semi)group is associated in a natural way with each
automaton (formally defined below). The two Mealy automata of Fig. 1 are
associated with finite (semi)groups. Their orders are respectively: on the
left a semigroup of order 238, on the right a group of order
$\numprint{1494186269970473680896}=2^{64}\cdot 3^{4}\approx 1.5\times
10^{21}$.

Several points are illustrated by this example:

  * •

An automaton and its dual generate (semi)groups which are either both finite
or both infinite (see [12, 2]).

  * •

The order of a finite automaton (semi)group can be amazingly large. It makes a
priori difficult to decide whether an automaton (semi)group is finite or not.
Actually, the decidability of this question is open (see [10, 2]).

  * •

The order of the (semi)groups generated by a Mealy automaton and its dual can
be strikingly different. It suggests to work with both automata together.

The contributions of the present paper are three-fold:

  * •

We propose new implementations (in GAP [8]) of classical algorithms for the
computation of the growth function; the computation of the order (if finite);
the semidecision procedure for the finiteness.

  * •

We compare the new implementations with the existing ones. Indeed, there exist
two GAP packages dedicated to Mealy automata and their associated
(semi)groups: FR by Bartholdi [4] and automgrp by Muntyan and Savchuk [11].

  * •

We realize systematic experimentations on small Mealy automata as well as
randomly chosen large Mealy automata. These serve as testbeds to some
conjectures on the growth types of the associated (semi)groups, as well as on
the order of a (semi)group.

The structure of the paper is the following. In Section 2, we present basic
notions on Mealy automata and automaton (semi)groups. In Section 3, we give
new implementations and compare them with the existing ones. Section 4 is
dedicated to experimentations and to the resulting conjectures.

##  2 Automaton (Semi)groups

###  2.1 Mealy Automaton

If one forgets initial and final states, a (finite, deterministic, and
complete) automaton ${\mathcal{A}}$ is a triple
$\bigl{(}A,\Sigma,\delta=(\delta_{i}:A\rightarrow A)_{i\in\Sigma}\bigr{)}$,
where the _set of states_ $A$ and the _alphabet_ $\Sigma$ are non-empty finite
sets, and where the $\delta_{i}$’s are functions.

A _Mealy automaton_ is a quadruple

$\bigl{(}A,\Sigma,\delta=(\delta_{i}:A\rightarrow
A)_{i\in\Sigma},\rho=(\rho_{x}:\Sigma\rightarrow\Sigma)_{x\in A}\bigr{)}\>,$

such that both $(A,\Sigma,\delta)$ and $(\Sigma,A,\rho)$ are automata. In
other terms, a Mealy automaton is a letter-to-letter transducer with the same
input and output alphabets. The transitions of a Mealy automaton are

$x\xrightarrow{i|\rho_{x}(i)}\delta_{i}(x)\>.$

The graphical representation of a Mealy automaton is standard, see Fig. 1.

The notation $x\xrightarrow{{\mathbf{u}}|{\mathbf{v}}}y$ with
${\mathbf{u}}=u_{1}\cdots u_{n}$, ${\mathbf{v}}=v_{1}\cdots v_{n}$ is a
shorthand for the existence of a path
$x\xrightarrow{u_{1}|v_{1}}x_{1}\xrightarrow{u_{2}|v_{2}}x_{2}\longrightarrow\cdots\longrightarrow
x_{n-1}\xrightarrow{u_{n}|v_{n}}y$ in ${\mathcal{A}}$.

In a Mealy automaton $(A,\Sigma,\delta,\rho)$, the sets $A$ and $\Sigma$ play
dual roles. So we may consider the _dual (Mealy) automaton_ defined by
${\mathfrak{d}}({\mathcal{A}})=(\Sigma,A,\rho,\delta)$, that is:

$i\xrightarrow{x\mid y}j\ \in{\mathfrak{d}}({\mathcal{A}})\quad\iff\quad
x\xrightarrow{i\mid j}y\ \in{\mathcal{A}}\>.$

It is pertinent to consider a Mealy automaton and its dual together, that is
to work with the pair $\\{{\mathcal{A}},{\mathfrak{d}}({\mathcal{A}})\\}$, see
an example in Fig. 1.

Let ${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ and
${\mathcal{B}}=(B,\Sigma,\gamma,\pi)$ be two Mealy automata acting on the same
alphabet; their _product_ ${\mathcal{A}}\times{\mathcal{B}}$ is defined as the
Mealy automaton with stateset $A\times B$, alphabet $\Sigma$, and transitions:

$xy\xrightarrow{i|\pi_{y}(\rho_{x}(i))}\delta_{i}(x)\gamma_{\rho_{x}(i)}(y)\>.$

###  2.2 Generating (Semi)groups

Let ${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ be a Mealy automaton. We view
${\mathcal{A}}$ as an automaton with an input and an output tape, thus
defining mappings from input words over $\Sigma$ to output words over
$\Sigma$. Formally, for $x\in A$, the map
$\rho_{x}:\Sigma^{*}\rightarrow\Sigma^{*}$, extending
$\rho_{x}:\Sigma\rightarrow\Sigma$, is defined by:

$\rho_{x}({\mathbf{u}})={\mathbf{v}}\quad\textrm{if}\quad\exists y,\
x\xrightarrow{{\mathbf{u}}|{\mathbf{v}}}y\>.$

By convention, the image of the empty word is itself. The mapping $\rho_{x}$
is length-preserving and prefix-preserving. It satisfies

$\forall u\in\Sigma,\
\forall{\mathbf{v}}\in\Sigma^{*},\qquad\rho_{x}(u{\mathbf{v}})=\rho_{x}(u)\rho_{\delta_{u}(x)}({\mathbf{v}})\>.$

We say that $\rho_{x}$ is the _production function_ associated with
$({\mathcal{A}},x)$. For ${\mathbf{x}}=x_{1}\cdots x_{n}\in A^{n}$ with $n>0$,
set
$\rho_{\mathbf{x}}:\Sigma^{*}\rightarrow\Sigma^{*},\rho_{\mathbf{x}}=\rho_{x_{n}}\circ\cdots\circ\rho_{x_{1}}\>$.

Denote dually by $\delta_{i}:A^{*}\rightarrow A^{*},i\in\Sigma$, the
production mappings associated with the dual Mealy automaton
${\mathfrak{d}}({\mathcal{A}})$. For ${\mathbf{v}}=v_{1}\cdots
v_{n}\in\Sigma^{n}$ with $n>0$, set $\delta_{\mathbf{v}}:A^{*}\rightarrow
A^{*},\ \delta_{\mathbf{v}}=\delta_{v_{n}}\circ\cdots\circ\delta_{v_{1}}$.

###### Definition 1

Consider a Mealy automaton ${\mathcal{A}}$. The semigroup of mappings from
$\Sigma^{*}$ to $\Sigma^{*}$ generated by $\rho_{x},x\in A$, is called the
_semigroup generated by ${\mathcal{A}}$_ and is denoted by
$\langle{{{\mathcal{A}}}}\rangle_{+}$. A semigroup $G$ is an _automaton
semigroup_ if there exists a Mealy automaton ${\mathcal{A}}$ such that
$G=\langle{{{\mathcal{A}}}}\rangle_{+}$.

A Mealy automaton ${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ is _invertible_ if
all the mappings $\rho_{x}:\Sigma\to\Sigma$ are permutations. Then the
production functions $\rho_{x}:\Sigma^{*}\to\Sigma^{*}$ are invertible.

###### Definition 2

Let ${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ be invertible. The _group generated
by ${\mathcal{A}}$_ is the group generated by the mappings
$\rho_{x}:\Sigma^{*}\to\Sigma^{*}$, $x\in A$. It is denoted by
$\langle{{\mathcal{A}}}\rangle$.

Let ${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ be an invertible Mealy automaton.
Its _inverse_ is the Mealy automaton ${{\mathcal{A}}}^{-1}$ with stateset
$A^{-1}=\\{x^{-1},x\in A\\}$ and set of transitions

$x^{-1}\xrightarrow{j\mid i}y^{-1}\ \in{\mathcal{A}}^{-1}\quad\iff\quad
x\xrightarrow{i\mid j}y\ \in{\mathcal{A}}\>.$

A Mealy automaton is _reversible_ if its dual is invertible. A Mealy automaton
${\mathcal{A}}$ is _bireversible_ if both ${\mathcal{A}}$ and
${{\mathcal{A}}}^{-1}$ are invertible and reversible.

######  Theorem 2.1 ([2, 12, 13])

The (semi)group generated by a Mealy automaton is finite if and only if the
(semi)group generated by its dual is finite.

###  2.3 Minimization and the Word Problem

Let $\mathcal{A}=(A,\Sigma,\delta,\rho)$ be a Mealy automaton. The _Nerode
equivalence on $A$_ is the limit of the sequence of increasingly finer
equivalences $(\equiv_{k})$ recursively defined by:

$\displaystyle\forall x,y\in A,\qquad\qquad x\equiv_{0}y$ $\displaystyle\
\Longleftrightarrow\ \rho_{x}=\rho_{y}\>,$ $\displaystyle\forall k\geqslant
0,\,x\equiv_{k+1}y$ $\displaystyle\ \Longleftrightarrow\
x\equiv_{k}y\quad\text{and}\quad\forall i\in\Sigma,\
\delta_{i}(x)\equiv_{k}\delta_{i}(y)\>.$

Since the set $A$ is finite, this sequence is ultimately constant; moreover if
two consecutive equivalences are equal, the sequence remains constant from
this point. The limit is therefore computable. For every element $x$ in $A$,
we denote by $[x]$ the class of $x$ w.r.t. the Nerode equivalence.

###### Definition 3

Let $\mathcal{A}=(A,\Sigma,\delta,\rho)$ be a Mealy automaton and let $\equiv$
be the Nerode equivalence on $A$. The _minimization_ of $\mathcal{A}$ is the
Mealy automaton
${\mathfrak{m}}(\mathcal{A})=(A/\mathord{\equiv},\Sigma,\tilde{\delta},\tilde{\rho})$,
where for every $(x,i)$ in $A\times\Sigma$,
$\tilde{\delta}_{i}([x])=[\delta_{i}(x)]$ and $\tilde{\rho}_{[x]}=\rho_{x}$.

This definition is consistent with the standard minimization of “deterministic
finite automata” where instead of considering the mappings
$(\rho_{x}:\Sigma\to\Sigma)_{x}$, the computation is initiated by the
separation between terminal and non-terminal states. Using Hopcroft algorithm,
the time complexity of minization is ${\cal O}(\Sigma A\log{A})$, see [1].

By construction, a Mealy automaton and its minimization generate the same
semigroup. Indeed, two states of a Mealy automaton belong to the same class
w.r.t the Nerode equivalence if and only if they represent the same element in
the generated (semi)group.

Consider the _word problem_ :

Input: a Mealy automaton $(A,\Sigma,\delta,\rho)$;
${\mathbf{x}},{\mathbf{y}}\in A^{*}$.

Question:
$(\rho_{{\mathbf{x}}}:\Sigma^{*}\to\Sigma^{*})=(\rho_{{\mathbf{y}}}:\Sigma^{*}\to\Sigma^{*})$?

The word problem is solvable by extending the above minimization procedure. FR
uses this approach, while automgrp uses a method based on the wreath recursion
[7].

##  3 Fully Exploiting the Minimization

Consider the following problems for the (semi)group given by a Mealy
automaton: compute the growth function, compute the order (if finite), detect
the finiteness. The packages FR and automgrp provide implementations of the
three problems. Here we propose new implementations based on a simple idea
which fully uses the automaton structure.

###  3.1 Growth

Consider a Mealy automaton ${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ and an
element ${\mathbf{x}}\in A^{*}$. The _length_ of $\rho_{{\mathbf{x}}}$,
denoted by $|\rho_{{\mathbf{x}}}|$, is defined as follows:

$|\rho_{{\mathbf{x}}}|=\min\\{n\mid\exists{\mathbf{y}}\in
A^{n},\,\rho_{{\mathbf{x}}}=\rho_{{\mathbf{y}}}\\}\>.$

The _growth series_ of ${\mathcal{A}}$ is the formal power series given by

$\sum\limits_{g\in\langle{{{\mathcal{A}}}}\rangle_{+}}t^{|g|}=\sum\limits_{n\in\mathbb{N}}\\#\\{g\in\langle{{{\mathcal{A}}}}\rangle_{+}\,;\,|g|=n\\}\>t^{n}\>.$

In words, the growth series enumerates the semigroup elements according to
their length. This is an instanciation of the notion of spherical growth
series for a finitely generated semigroup. Observe that the series is a
polynomial if and only if the semigroup is finite.

##### Using the Generic Algorithm.

Since the word problem is solvable, it is possible to compute an arbitrary but
finite number of coefficients of the growth series. Indeed for each $n$,
generate the set of elements of length $n$ by multiplying elements of length
$n-1$ with generators and detecting-deleting duplicated elements by solving
the word problem. The functions Growth from automgrp and WordGrowth from FR
both follow this pattern. Therefore the structure of the underlying Mealy
automaton is used only to get a solution to the word problem (in fact, both
Growth and WordGrowth are generic, in the sense that they are applicable for
any (semi)group with an implemented solution to the word problem).

##### New Implementation.

We propose a new implementation based on a simple observation. Knowing the
elements of length $n-1$, Nerode minimization can be used in a global manner
to obtain simultaneously the elements of length $n$. Concretely, with each
integer $n\geq 1$ is associated a new Mealy automaton ${{\mathcal{A}}}_{n}$
defined recursively as follows:

${{\mathcal{A}}}_{n}={\mathfrak{m}}({{\mathcal{A}}}_{n-1}\times{\mathfrak{m}}({\mathcal{A}}))\qquad\hbox{and}\qquad{{\mathcal{A}}}_{1}={\mathfrak{m}}({\mathcal{A}})\>.$

Here, we assume, without real loss of generality, that the identity element is
one of the generators (otherwise simply add a new state to the Mealy automaton
coding the identity). This way, the elements of ${\mathcal{A}}_{n}$ are
exactly the elements of length at most $n$.

⬇

AutomatonGrowth:= function(arg)

local aut, radius, growth, sph, curr, next, r;

aut:=arg[1]; # Mealy automaton

if Length(arg)>1 then radius:=arg[2];

else radius:=infinity;

fi;

r := 0; curr := TrivialMealyMachine([1]);

next := Minimized(aut);

aut := Minimized(next+TrivialMealyMachine(Alphabet(aut)));

sph := aut!.nrstates - 1; # number of non-trivial states

growth := [next!.nrstates-sph];

while sph>0 and r<radius

do Add(growth,sph);

r := r+1; curr := next;

next := Minimized(next*aut);

sph := next!.nrstates-curr!.nrstates;

od;

return growth;

end;

Note that AutomatonGrowth(aut) computes the growth of the semigroup
$\langle{\tt aut}\rangle_{+}$, while AutomatonGrowth(aut+aut^-1) computes the
growth of the group $\langle{\tt aut}\rangle$.

##### Experimental Results.

First we run AutomatonGrowth and FR’s WordGrowth on the Grigorchuk automaton,
a famous Mealy automaton generating an infinite group. For radius 10,
AutomatonGrowth is much faster, 76 ms as opposed to 9 912 ms111All timings
displayed in this paper have been obtained on an Intel Core 2 Duo computer
with clock speed 3,06 GHz.. The explanation is simple: WordGrowth calls the
minimization procedure 57 577 times while AutomatonGrowth calls it only 12
times. Here are the details.

⬇

gap> aut := GrigorchukMachine;; radius:= 10;;

gap> ProfileFunctions([Minimized]);

gap> WordGrowth(SCSemigroupNC(aut), radius); time;

[ 1, 4, 6, 12, 17, 28, 40, 68, 95, 156, 216 ]

9912

gap> DisplayProfile();

count self/ms chld/ms function

57577 7712 0 Minimized

7712 TOTAL

gap> ProfileFunctions([Minimized]);

gap> AutomatonGrowth(aut, radius); time;

[ 1, 4, 6, 12, 17, 28, 40, 68, 95, 156, 216 ]

76

gap> DisplayProfile();

count self/ms chld/ms function

12 72 0 Minimized

72 TOTAL

Now we compare the running times of the implementations for the computation of
the first terms of the growth series for all 335 bireversible 3-letter 3-state
Mealy automata (up to equivalence). In Tab. 1, some computations with FR’s
WordGrowth or with automgrp’s Growth could not be completed in reasonable time
for radius 7.

Table 1: Average time (in ms) radius | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
FR’s WordGrowth |  3.4 |  29.0 |  555.0 |  8 616.5 |  131 091.4 |  2 530 170.3 |  ?
automgrp’s Growth | 0.7 | 2.8 | 16.9 | 158.9 |  1 909.0 |  22 952.8 |  ?
AutomatonGrowth | 0.6 | 1.8 | 5.9 | 28.9 |  187.3 |  1 005.9 |  7 131.4

###  3.2 Order of the (Semi)group

Although the finiteness problem is still open, some semidecision procedures
enable to find the order of an expected finite (semi)group. FR and automgrp
use orthogonal approaches. Our new implementation refines the one of FR and
remains orthogonal to the one of automgrp.

#####  automgrp’s Implementation.

The GAP package automgrp provides the function LevelOfFaithfulAction, which
allows to compute—very efficiently in some cases—the order of the generated
group. The principle is the following. Let
${\mathcal{A}}=(A,\Sigma,\delta,\rho)$ be an invertible Mealy automaton and
let $G_{k}$ be the group generated by the restrictions of the production
functions to $\Sigma^{k}$. If $\\#G_{k}=\\#G_{k+1}$ for some $k$, then
$\langle{{\mathcal{A}}}\rangle$ is finite of order $\\#G_{k}$. This function
can be easily adapted to a non-invertible Mealy automaton.

Observe that LevelOfFaithfulAction cannot be used to compute the growth
series. Indeed at each step a quotient of the (semi)group is computed. On the
other hand LevelOfFaithfulAction is a good bypass strategy for the order
computation. Furthermore, it takes advantage from the special ability of GAP
to manipulate permutation groups.

#####  FR’s Implementation and the New Implementation.

Any algorithm computing the growth series can be used to compute the order of
the generated (semi)group if finite. It suffices to compute the growth series
until finding a coefficient equal to zero. This is the approach followed by
FR. Since we proposed, in the previous section, a new implementation to
compute the growth series, we obtain as a byproduct a new procedure to compute
the order. We call it AutomSGrOrder.

##### Experimental Results.

The orthogonality of the two previous approaches can be simply illustrated by
recalling the introductory example of Fig. 1. Neither FR’s Order nor
AutomSGrOrder are able to compute the order of the large group, while automgrp
via LevelOfFaithfulAction succeeds in only 14 338 ms. Conversely,
AutomSGrOrder computes the order of the small semigroup in 17 ms, while an
adaptation of LevelOfFaithfulAction (to non-invertible Mealy automata) takes 2
193 ms.

###  3.3 Finiteness

There exist several criteria to detect the finiteness of an automaton
(semi)group, see [2, 3, 6, 14, 15, …]. But the decidability of the finiteness
is still an open question. Each procedure to compute the order of a
(semi)group yields a semidecision procedure for the finiteness problem. Both
packages FR and automgrp apply a number of previously known criteria of
(in)finiteness and then intend to conclude by ultimately using an order
computation.

We propose an additional ingredient which uses minimization in a subtle way.
Here, the semigroup to be tested is successively replaced by new ones which
are finite if and only if the original one is finite. It is possible to
incorporate this ingredient to get two new implementations, one in the spirit
of FR and one in the spirit of automgrp. The new implementations are order of
magnitudes better than the old ones. Both are useful since the fastest one
depends on the cases.

####  3.3.1 ${\mathfrak{m}}{\mathfrak{d}}$-reduction of Mealy Automata and
Finiteness

The ${\mathfrak{m}}{\mathfrak{d}}$-reduction was introduced in [2] to give a
sufficient condition of finiteness. The new semidecision procedures start with
this reduction.

###### Definition 4

A pair of dual Mealy automata is _reduced_ if both automata are minimal.
Recall that ${\mathfrak{m}}$ (resp. ${\mathfrak{d}}$) is the operation of
minimization (resp. dualization). The _${\mathfrak{m}}{\mathfrak{d}}$
-reduction_ of a Mealy automaton ${\mathcal{A}}$ consists in minimizing the
automaton or its dual until the resulting pair of dual Mealy automata is
reduced.

The ${\mathfrak{m}}{\mathfrak{d}}$-reduction is well-defined: if both a Mealy
automaton and its dual automaton are non-minimal, the reduction is confluent
[2]. An example of ${\mathfrak{m}}{\mathfrak{d}}$-reduction is given in Fig.
2.

.7 (-3,-19)(26,-5) [a](0,-8)AA [b](6,-8)BB (8,-9)${\mathcal{A}}$ AABB0123
BBAA0321 [.2]AA1032 [.8]BB1032 (10,-8)E2 (13,-8)F2 E2F2d [0](17,-5)A0
[1](23,-5)A1 [3](17,-9)A3 [2](23,-9)A2 [.3]A1A0aabb A0A1ab [.3]A3A2aabb A2A3ab
A0A3ba A2A1ba (20,-11)E3 (20,-13)F3 E3F3m [13](17,-16)A13 [02](23,-16)A02
A13A02aabb A02A13abba (13,-16)E4 (10,-16)F4 dashed E4F4dmdmd solid
[ab](3,-16)Z (7,-17)${\mathfrak{m}}{\mathfrak{d}}^{*}{({\mathcal{A}})}$
[.2]Z01230123

Figure 2: The ${\mathfrak{m}}{\mathfrak{d}}$-reduction of a pair of dual Mealy
automata

The sequence of minimization-dualization can be arbitrarily long: the
minimization of a Mealy automaton with a minimal dual can make the dual
automaton non-minimal.

If ${\mathcal{A}}$ is a Mealy automaton, we denote by
${\mathfrak{m}}{\mathfrak{d}}^{*}{({\mathcal{A}})}$ the corresponding Mealy
automaton after ${\mathfrak{m}}{\mathfrak{d}}$-reduction.

######  Theorem 3.1 ([2])

A Mealy automaton ${\mathcal{A}}$ generates a finite (semi)group if and only
if ${\mathfrak{m}}{\mathfrak{d}}^{*}{({\mathcal{A}})}$ generates a finite
(semi)group.

This is the starting point of the new implementations. We use an additional
fact. We can prune a Mealy automaton by deleting the states which are not
accessible from a cycle. This does not change the finiteness or infiniteness
of the generated (semi)group [3].

####  3.3.2 The New Implementations

The design of procedure IsFinite1 is consistent with the one of
AutomatonGrowth. Hence IsFinite1 is much closer to FR than to automgrp. Here
we propose a version that works with the automaton and its dual in parallel.

⬇

IsFinite1 := function (aut, limit)

local radius, dual, curr1, next1, curr2, next2;

radius := 0;

aut := MDReduced(Prune(aut)); dual := DualMachine(aut);

curr1 := MealyMachine([[1]],[()]); curr2 := curr1;

next1 := aut; next2 := dual;

while curr2!.nrstates<>next2!.nrstates and radius<limit

do radius := radius + 1; curr1 := next1;

next1 := Minimized(next1*aut);

if curr1!.nrstates<>next1!.nrstates

then curr2 := next2;

next2 := Minimized(next2*dual);

else return true;

fi;

od;

if curr2!.nrstates = next2!.nrstates then return true; fi;

return fail;

end;

The procedure IsFinite2 is a refinement of automgrp’s LevelOfFaithfulAction:
the minimization is called on the dual and can be enhanced again to work in
parallel on the Mealy automaton and its dual.

⬇

IsFinite2 := function(aut,limit)

local f1, f2, next, cs, ns, lev;

aut := MDReduced(Prune(aut));

if IsInvertible(aut) then f1:=Group; f2:=PermList;

else f1:=Semigroup; f2:=Transformation;

fi;

lev := 0; cs := 1; ns := Size(f1(List(aut!.output,f2)));

aut := DualMachine(aut); next := aut;

while cs<ns and lev<limit

do lev := lev+1; cs := ns; next := Minimized(next*aut);

ns := Size(f1(List(DualMachine(next)!.output,f2)));

od;

if cs=ns then return true; else return fail; fi;

end;

##### Experimental Results.

Tab. 2 presents the average time to detect finiteness of (semi)groups
generated by $p$-letter $q$-state invertible or reversible Mealy automata with
$p+q\in\\{5,6\\}$. To get a fair comparison of the implementations, what is
given is the minimum of the running times for an automaton and its dual (see
Theorem 2.1).

Table 2: Average time (in ms) to detect finiteness of (semi)groups

2- 3- | 2- 4- | 3- 3-
---|---|---
FR | aut | Fin1 | Fin2 | FR | aut | Fin1 | Fin2 | FR | aut | Fin1 | Fin2
0.68 | 0.81 | 0.49 | 0.49 | 36.36 | 1.79 | 0.52 | 0.62 |  1 342.12 | 3.78 | 0.61 | 0.70

FR: FR’s IsFinite; aut: automgrp’s IsFinite; Fin1: IsFinite1; Fin2: IsFinite2

##  4 Conjectures

The efficiency of the new implementations enables to carry out extensive
experimentations. We propose several conjectures supported by these
experiments.

Recall the example given in the introduction. The (semi)groups generated by
the Mealy automaton and its dual were strikingly different, with a very large
one and a rather small one. This seems to be a general fact that we can state
as an informal conjecture:

_Whenever a Mealy automaton generates a finite (semi)group which is very large
with respect to the number of states and letters of the automaton, then its
dual generates a small one._

_Observation:_ Any pair of finite (semi)groups can be generated by a pair of
dual Mealy automata, see [2, Prop. 9]. The standard construction leads to
automata whose sizes are related to the orders of the (semi)groups. Therefore
it does not contradict the informal conjecture.

$\\#\langle{{{\mathcal{A}}}}\rangle_{+}$$\\#\langle{{{\mathfrak{d}}({{\mathcal{A}}})}}\rangle_{+}$4
000$10^{2}$$10^{4}$$10^{6}$$10^{8}$$10^{10}$$10^{12}$$10^{14}$$10^{16}$$10^{18}$$10^{20}$$10^{22}$
Figure 3: Size of $\langle{{{\mathcal{A}}}}\rangle_{+}$ vs. size of
$\langle{{{\mathfrak{d}}({{\mathcal{A}}})}}\rangle_{+}$

Fig. 3 illustrates this informal conjecture: for ${\mathcal{A}}$ covering the
set of all 3-letter 3-state invertible Mealy automata, the endpoints of each
segment represent respectively the order of
$\langle{{{\mathcal{A}}}}\rangle_{+}$ and of
$\langle{{{\mathfrak{d}}({{\mathcal{A}}})}}\rangle_{+}$, for all pairs
detected as being finite.

To assess finiteness, the procedures IsFinite1 and IsFinite2 have been used.
If the tested Mealy automaton and its dual were both found to have more than
4000 elements, the procedures were stopped, and the (semi)groups were supposed
to be infinite. Based on the informal conjecture, we believe to have captured
all finite groups. If true:

  * •

There are 14 089 Mealy automata generating finite (semi)groups among the 233
339 invertible or reversible 3-letter 3-state Mealy automata;

  * •

The group generated by Fig. 1-right is the largest finite group.

.4 2.2 (-4.5,-1.3)(1.5,.5) (-.62,.78)A1 [x](.22,.97)AQ (-1,0)A2 (0,0)A0
(.9,-.43)A5(.9,.43)A6 (-.62,-.78)A3 (.22,-.97)A4 A1A2$\footnotesize 2$
A2A3$\footnotesize 2$ A3A4$\footnotesize 2$ A4A5$\footnotesize 2$
A6AQ$\footnotesize 2$ [.3]A1A0$\footnotesize 1$ 135A1 75A0 AQA1 A2A0 A3A0 A4A0
A5A0 A4A0 A5A0 A6A0 dotted A5A6 (-3.3,-.2)$\rho_{x}=(1,2,\dots,p)$
(-3.3,-.6)$\forall y\neq x,\,\rho_{y}=(1,3,\dots,p)$

4.1 among invertible automata: ${\mathcal{M}}_{p,q}$

.4 2.2 (-1.5,-1.3)(5,.5) (-.62,.78)A1 [x](.22,.97)AQ (-1,0)A2 (0,0)A0
(.9,-.43)A5(.9,.43)A6 (-.62,-.78)A3 (.22,-.97)A4 A1A2 A2A3 A3A4 A5A6 A6AQ
[.3]A6A0$\footnotesize 1$ 75A0 AQA1 dotted A4A5 (3.2,-.2)$\rho_{x}=(1,2)$
(3.2,-.6)$\forall z\neq x,\,\rho_{z}=()$

4.2 among 2-letter invertible automata: ${\mathcal{M}}_{2,q}$

.4 2.2 (-5,-1.3)(1.5,1.1) [y](-0.82,.28)A1 [¯x](.22,.28)AQ AQA1 250A1
[.4]290AQ (.4,-.45)1,2 (.4,-.7)(plus $p$ if even)
(-3.05,.2)$\rho_{\bar{x}}=t(1,2,\dots,p)t^{-1}$
(-3.2,-.2)$\rho_{y}=(1,3,\dots,p)$
(-3.2,-.75)$t=\begin{cases}()&\hbox{for~{}$p$ even}\\\
(p,\frac{p+1}{2})&\hbox{for~{}$p$ odd}\end{cases}$

4.3 among 2-state invertible automata: ${\mathcal{M}}_{p,2}$

.4 2.2 (-1.5,-1.3)(5,1.1) [x](.22,.97)AQ [y](-.62,.78)A1 (-1,0)A2
(.9,-.43)A5(.9,.43)A6 (-.62,-.78)A3 (.22,-.97)A4 [.7]A1A2 [.3]A2A3 A3A4 A4A5
A6AQ AQA1 dotted A5A6 (3.2,.2)$\rho_{x}=(1,2,\dots,p)$
(3.2,-.2)$\rho_{y}=(1,3,\dots,p)$ (3.2,-.6)$\forall
z\not\in\\{x,y\\},\,\rho_{z}=(\,)$

4.4 among bireversible automata: ${\mathcal{B}}_{p,q}$

Figure 4: Automata conjectured to generate the largest finite automaton groups

Our next conjectures are concerned with the largest finite groups that can be
generated by automata of a given size.

Consider the family of $p$-letter $q$-state Mealy automata
$({\mathcal{M}}_{p,q})_{p+q>5}$ displayed on Fig. 4.1 for $p>2$ and $q>2$,
while the specializations for $p=2$ and $q=2$ are displayed on Fig. 4.2 and
Fig. 4.3. The example of Fig. 1-right is ${\mathcal{M}}_{3,3}$.

###### Conjecture 1

The group $\langle{{\mathcal{M}}_{p,q}}\rangle$ is finite. Every $p$-letter
$q$-state invertible Mealy automaton generates a group which is either
infinite or has an order smaller than
$\\#\langle{{\mathcal{M}}_{p,q}}\rangle$.

If true, Conjecture 1 implies the decidability of the finiteness problem for
automaton groups. Without entering into the details of the experimentations,
we consider that Conj. 1 is reasonably well supported for $p+q<9$. As for
actually computing $\\#\langle{{\mathcal{M}}_{p,q}}\rangle$, here are the only
cases with $q>2$ for which we succeeded:

$\displaystyle\forall q,\,4\leq q\leq
8,\qquad\\#\langle{{\mathcal{M}}_{2,q}}\rangle=2^{2^{q-1}+\frac{(q-2)(q-1)}{2}-2}\>,\qquad\qquad$
$\displaystyle\\#\langle{{\mathcal{M}}_{3,3}}\rangle=2^{64}\cdot
3^{4},\qquad\\#\langle{{\mathcal{M}}_{3,4}}\rangle=2^{325}\cdot
3^{13},\qquad\\#\langle{{\mathcal{M}}_{4,3}}\rangle=2^{288}\cdot 3^{422}\>.$

These groups are indeed huge. Incidentally, the finiteness of
$\langle{{\mathcal{M}}_{p,q}}\rangle$ is checked for $p+q<11$ and the informal
conjecture is supported further by computing the order of the much smaller
semigroups generated by the duals:

$\\#\langle{{\mathfrak{d}}({{\mathcal{M}}_{p,q}})}\rangle_{+}$ |  2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---
2 | - | - | 219 | 1 759 | 13 135 | 94 143 | 656 831
3 | - | 238 | 1 552 | 8 140 | 37 786 | 162 202 | $\cdots$
4 | 89 | 1 381 | 12 309 | 87 125 | 543 061 | $\cdots$ | $\cdots$
5 | 131 | 6 056 | 67 906 | 602 656 | $\cdots$ | $\cdots$ | $\cdots$
6 | 337 | 22 399 | 302 011 | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$
7 | 351 | 74 194 | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$

Experimentally, the finite groups generated by bireversible Mealy automata
seem to be much smaller. Consider the family of bireversible automata
$({\mathcal{B}}_{p,q})_{p,q}$ of Fig. 4.4. The group
$\langle{{\mathcal{B}}_{p,q}}\rangle$ is isomorphic to
${\mathfrak{S}}_{p}^{q}$, while the group
$\langle{{\mathfrak{d}}({{\mathcal{B}}_{p,q}})}\rangle$ is isomorphic to
$\mathbb{Z}_{q}$. Again, the following is reasonably well supported for
$p+q<9$:

###### Conjecture 2

Every $p$-letter $q$-state bireversible Mealy automaton generates a group
which is either infinite or has an order smaller than
$\\#\langle{{\mathcal{B}}_{p,q}}\rangle=p!^{q}$.

Our last conjecture is of a different nature and deals with the structure of
infinite automaton semigroups.

###### Conjecture 3

Every $2$-state reversible Mealy automaton generates a semigroup which is
either finite or free of rank 2.

The conjecture has been tested and seems correct for reversible 2-state Mealy
automata up to 6 letters. In the experiments, a semigroup generated by a
$p$-letter automaton is conjectured to be free if its growth series coincides
with $(2t)^{n}$ up to radius $p^{2}/2$ and if its dual generates a seemingly
infinite group.

## References

  * [1] A.V. Aho, J.E. Hopcroft, and J.D. Ullman.  The Design and Analysis of Computer Algorithms.  Addison-Wesley, 1974.
  * [2] A. Akhavi, I. Klimann, S. Lombardy, J. Mairesse, and M. Picantin.  On the finiteness problem for automaton (semi)groups.  Internat. J. Algebra Comput., (accepted), 2012.  arXiv:cs.FL/1105.4725.
  * [3] A. S. Antonenko.  On transition functions of Mealy automata of finite growth.  Matematychni Studii., 29(1):3–17, 2008.
  * [4] L. Bartholdi.  ${\sf FR}$ Functionally recursive groups — a GAP package, v.1.2.4.2, 2011.
  * [5] L. Bartholdi, I. I. Reznykov, and V. I. Sushchanskiĭ.  The smallest Mealy automaton of intermediate growth.  J. Algebra, 295(2):387–414, 2006.
  * [6] I. V. Bondarenko, N. V. Bondarenko, S. N. Sidki, and F. R. Zapata.  On the conjugacy problem for finite-state automorphisms of regular rooted trees.  arXiv:math.GR/1011.2227.
  * [7] A. J. Cain.  Automaton semigroups.  Theor. Comput. Sci., 410:5022–5038, 2009.
  * [8] The GAP Group.  GAP – Groups, Algorithms, and Programming, v.4.4.12, 2008.
  * [9] R. I. Grigorchuk.  On Burnside’s problem on periodic groups.  Funktsional. Anal. i Prilozhen., 14(1):53–54, 1980.
  * [10] R. I. Grigorchuk, V. V. Nekrashevich, and V. I. Sushchanskiĭ.  Automata, dynamical systems, and groups.  Tr. Mat. Inst. Steklova, 231:134–214, 2000.
  * [11] Y. Muntyan and D. Savchuk.  ${\sf automgrp}$ Automata Groups — a GAP package, v.1.1.4.1, 2008.
  * [12] V. Nekrashevych.  Self-similar groups, volume 117 of Mathematical Surveys and Monographs.  American Mathematical Society, Providence, RI, 2005.
  * [13] D. M Savchuk and Y. Vorobets.  Automata generating free products of groups of order 2.  J. Algebra, 336(1):53–66, 2011.
  * [14] S. N. Sidki.  Automorphisms of one-rooted trees: growth, circuit structure, and acyclicity.  J. Math. Sci. (New York), 100(1):1925–1943, 2000.  Algebra, 12.
  * [15] P. V. Silva and B. Steinberg.  On a class of automata groups generalizing lamplighter groups.  Internat. J. Algebra Comput., 15(5-6):1213–1234, 2005.
  * [16] Z. Šuniḱ and E. Ventura.  The conjugacy problem is not solvable in automaton groups.  2010\.  arXiv:math.GR/1010.1993.