The one counting challenge to rule them all!

Question

Background

Mathematicians are interested in all kinds of algebraic structures, e.g. groups, rings, boolean algebras. Generally, we have several operations (addition, multiplication, meet, join, etc.), and some equational relationships on the operations (e.g. associativity of multiplication). This makes up an algebraic theory. Formally, we want \$k\$ functions, the \$i\$-th of which takes \$s_i\$ arguments. We also want some equalities to hold between the functions. In our case, we only consider equalities between expressions of depth at most 2, e.g. \$f(g(x), h(x, y), h(y, z), f(x, y, w, z))\$. The equalities should hold whatever values you plug in for the variables.

Let's take a quick example. A semigroup is such an algebraic structure. We need to have one function \$f_1\$, with \$s_1=2\$. There is just one equational relationship: \$f_1(f_1(x, y), z) = f_1(x, f_1(y, z))\$.

We are interested in models of the algebraic theories. A model, simply put, is just an implementation of the functions. We choose an underlying set \$M\$, and we implement the required functions \$f_i\$ as taking \$s_i\$ arguments from the set \$M\$, and returning an element of \$M\$, so that the equations hold when you plug in any element of \$M\$ for the variables. For example, if we take \$M = \{0,1,2\}\$, and \$f_1\$ to be the addition modulo 3, we get a model for the semigroup. The order of a model is simply the size of \$M\$. So the model above has order 3.

Suppose that, for a fixed algebraic theory, there are two models, with underlying set \$M\$ and \$N\$, and functions \$f_i\$ and \$g_i\$, respectively. If, after some relabeling of the elements, two models become exactly the same, we say they are isomorphic. In that case, the two models are essentially the same, and usually we don't care about their difference.

More formally, if we can find a mapping \$\sigma: M \to N\$ such that \$g_i(\sigma(m_1), \sigma(m_2), \dots, \sigma(m_{s_i})) = \sigma(f_i(m_1, \dots, m_{s_i}))\$, we say that \$\sigma\$ is a morphism. Furthermore, if there is another morphism \$\pi : N \to M\$, and \$\sigma(\pi(x)) = x\$, \$\pi(\sigma(y)) = y\$ holds for every \$x \in N, y\in M\$, we say that the two models are isomorphic.

Task

Your challenge is a very natural question that arises when studying algebraic theories.

Given an algebraic theory (whose equations involve only expressions of depth at most 2) and a natural number \$n\$, compute the number of distinct models of order \$n\$, up to isomorphism.

You may write a function or a complete program. You can suppose the input is reasonably parsed for you. See the examples. Standard loopholes apply. This is code-golf, so the program with the shortest byte length wins.

Example Cases

Input: The theory of groupoids. There is one function f taking 2 arguments.

f(f(x, y), z) = f(x, f(y, z))

n = 2.

Output: 5.

We'll fully work this out. For the rest of the example cases the basic idea is the same.

Let's take M = {0, 1}, the actual labelling of the elements obviously doesn't affect the outcome. We have four values to decide on, f(0, 0), f(0, 1), f(1, 0), f(1, 1). Let's case split on f(0, 1).

• Case 1: f(0, 1) = 0. So 0 = f(0, 1) = f(f(0, 1), 1) = f(0, f(1, 1)). We further case split on f(1, 1).

  • Case A: f(1, 1) = 0. Then 0 = f(0, f(1, 1)) = f(0, 0). If further, f(1, 0) = 0, then f is a constant function, and obviously satisfies the equation. If f(1, 0) = 1, then 1 = f(1, 0) = f(1, f(1, 1)) = f(f(1, 1), 1) = f(0, 1) = 0, contradiction.
  • Case B: f(1, 1) = 1. Let's consider f(0, 0).
    • Case i: f(0, 0) = 0. Then f(1, 0) = f(f(1, 1), 0) = f(1, f(1, 0)). If f(1, 0) = 1, plugging that in yields 1 = f(1, 1) = 0, contradiction. SO f(1, 0) = 0. So f(x, y) = x and y, which satisfies the equation, as you learned early in discrete mathematics.
    • Case ii: f(0, 0) = 1. Then f(1, 0) = f(f(0, 0), 0) = f(0, f(0, 0)) = f(0, 1) = 0. So the whole thing is just the xnor function.

Now let's look back: we've already got the constant function 0, the boolean and, and xnor. We proceed:

• Case 2: f(0, 1) = 1. We could do the case analysis all over again. But note that everything exactly mirrors case 1, except that 0 and 1 are inverted, and the two arguments of f is swapped. Since swapping the arguments doesn't affect associativity, we immediately get what we wanted: the constant function 1, the boolean nand and xor.

Now we need to make sure they are all non-isomorphic. In fact, the two constant functions are isomorphic, since the relabelling that swaps 0 and 1 converts between them. So we've got 5 different possibilities.

Input: There are three functions f, g, h taking 1, 2 and 0 arguments, respectively. 
The equations are:

g(g(x, y), z) = g(x, g(y, z))
g(x, f(x)) = h()
g(f(x), x) = h()
g(h(), x) = x
g(x, h()) = x

n = 3

Output: 1.

You may assume that the functions are labeled with natural numbers, here we use letters to make it more readable. Note that the expressions on both sides of the equalities should be implemented as trees or nested lists. You can assume that such trees are passed to your function as argument, so you can skip the parsing part. But in case you want to fiddle with eval-like stuff you can also accept string input.

The model of this algebraic theory is exactly the groups. And all groups of order 3 are isomorphic to \$\mathbb Z_3\$, which takes the underlying set to be \$\{0,1,2\}\$, and g is the addition modulo 3, f is negation, h() equals 0. Then the equations translate to:

• \$(x + y) + z = x + (y + z)\$;
• \$x + (-x) = 0\$;
• \$(-x) + x = 0\$;
• \$0 + x = x\$;
• \$x + 0 = x\$.

These are just the familiar properties of addition.

Input: Same as previous, except `n = 6`.

Output: 2.

In fact, for this algebraic theory, the answer should be the very first sequence of OEIS, which demonstrates the importance of such sequences.

---

Input: There are 4 functions f, g, h, i, that takes 2, 2, 0, 1 arguments, repectively. The equations are:

f(f(x, y), z) = f(x, f(y, z))
g(g(x, y), z) = g(x, g(y, z))
f(x, y) = f(y, x)
f(i(x), x) = h()
f(h(), x) = x
g(x, f(y, z)) = f(g(x, y), g(x, z))
g(f(y, z), x) = f(g(y, x), g(z, x))

n = 10.

Output: 4.

This is OEIS A027623.

Input: There is just 1 function f taking 1 argument. The equations are:

f(x) = f(x)
f(f(y)) = f(f(y))

n = 30

Output: 10712682919341.

Note that the equations are redundant. f is simply any function. We draw a directed edge from a to b if f(a) = b, this produces a graph where every vertex's out-degree is 1. Since we consider things up to isomorphism, we need to count such graphs with unlabeled nodes. We first consider a connected component. Since the out-degree is always 1, you can uniquely follow the edges until you hit a previously visited vertex. This produces a cycle (possibly with length 1 or 2). We consider the rest of the vertices in this connected component. They have to be on trees rooted on the cycle. Counting such graphs should be standard practice. The number of connected component is A002861, and the total number is A001372.

Input: 2 functions f, g both taking 2 arguments.

f(x, y) = f(y, x)
g(x, y) = g(y, x)
f(f(x, y), z) = f(x, f(y, z))
g(g(x, y), z) = g(x, g(y, z))
f(x, g(x, y)) = x
g(x, f(x, y)) = x

Output: OEIS A006966.

It is clear that this challenge is a very general generalization of various counting problems.

---

Below are some degenerate corner cases.

Input: There are no functions. There are no equalities. n = 7.

Output: 1.

This is just finding sets with 7 elements, and if there are bijections between two sets they count as the same. Of course, there are bijections between every two sets with 7 elements. So the answer is 1.

Input: There are no functions. The only equality is x = y. n = 7.

Output: 0.

The equalities should hold for all x and y. So there can only be at most 1 distinct element, but we require 7.

Input: Same as above, except `n = 1`.

Output: 1.

Input: Any algebraic theory, n = 0.

Output: 1. (Regardless of the theory.)

Answer 1

Haskell, 482 ... 408 402 bytes

• -5 thanks to @ovs
• -18 by inlining functions that are used just once
• -12 by introducing short names of constants and removing parentheses
• -11 by making a function local so it can take less arguments
• -29 thanks to @Laikoni and @Trebor
• -6 bytes by preferring list comprehension to filter, and turn F Int[A] into Q[A]Int for 1 space (I changed the letter to ease the changing of the test cases).

import Data.List
data A=Q[A]Int|V Int
t=map
w=sequence
y(V i)=i
y(Q l _)=maximum$0:t y l
(0!_)_=1
(n!e)q=c[b|m<-w$t(a.(n^))e,let b=[\b->f!!sum(zipWith(*)b(t(n^)r))|f<-m],and[h l==h r|(l,r)<-q,v<-a$1+y l`max`y r,let h(V i)=v!!i;h(Q e f)=b!!f$h<$>e]]where;a i=w$r<$[1..i];d#w=all(\p->or[(t(.t(p!!))d!!f)q/=(t((p!!).)w!!f)q|(f,z)<-zip r e,q<-a z])$permutations r;c[]=0;c x=1+c(filter(head x#)x);r=[0..n-1]

Try it online!

The function (!) accepts input as (n ! arities) equations, where n is the given n, arities is a list of \$s_i\$, as specified in the challenge, equations should be in the predefined datatype A, whose ungolfed version reads

data Expression
  = Func Int [Expression] | Var Int

An expression is either a Var variable, indexed by natural numbers (here Int saves bytes), or Func i [arg1, arg2, arg3, ...], which is just the \$i\$-th function, applied to the list of arguments. I have entered the equations in the footer in TIO, in case you are interested. (Func Int [Expression] becomes Q[A]Int in the golfed code to save one space.)

As mentioned by @fireflame241 it is very hard to get anything interesting done within the time limit. Nonetheless I've written a decent set of tests to check against OEIS.

---

Ungolfed version (still messy but hopefully the comment makes it clear):

import           Data.List

data Expression
  = Func Int [Expression] | Var Int deriving (Show)

tabulate n arity f =
  intercalate
    "\n"
    [ "f" ++ show assv ++ " = " ++ show (f assv)
    | assv <- sequence (replicate arity [0 .. n - 1])
    ]

showModel n arities m =
  intercalate "\n\n" [tabulate n a f | (a, f) <- zip arities m]

-- auxiliary function to find the largest variable number
-- this is less efficient, but saves the hassle of having to bookkeep the used variables
vars (Var i)    = i
vars (Func _ l) = maximum $ 0 : map vars l  -- the 0 is to avoid empty lists

-- evaluates an expression in the model `funcs`
evaluate :: [[Int] -> Int] -> [Int] -> Expression -> Int
evaluate funcs varass (Var i) = varass !! i
evaluate funcs varass (Func f es) = funcs !! f $ map (evaluate funcs varass) es

-- auxiliary function to generate all possible variable assignments
varasss n l r = sequence (replicate (1 + max (vars l) (vars r)) [0 .. n - 1])

-- I think this can be further golfed, but havent come up with good ones yet
checkValid ::
     Int -> [Int] -> [(Expression, Expression)] -> [[Int] -> Int] -> Bool
checkValid n arities equations funcs =
  and
    [ evaluate funcs v l == evaluate funcs v r
    | (l, r) <- equations
    , v <- varasss n l r
    ]

-- generates all models. The outer list comprehension runs over a cartesian product M1 * M2 * ...,
-- where Mi is the (flattened) table of all the possible implementation of the i-th function
-- Then the complicated `sum $ zipWith` stuff computes the index of the flattened table Mi
-- from the un-flattened index
allModels :: Int -> [Int] -> [[[Int] -> Int]]
allModels n arities =
  [ [\b -> f !! (sum $ zipWith (*) b (map (n ^) [0 ..])) | f <- m]
  | m <-
      sequence (map (\i -> sequence (replicate (n ^ i) [0 .. n - 1])) arities)
  ]

-- see if two models are identical, by comparing the values.
ident :: Int -> [Int] -> [[Int] -> Int] -> [[Int] -> Int] -> Bool
ident n arities m1 m2 =
  and
    [ (m1 !! f) assv == (m2 !! f) assv
    | f <- [0 .. length arities - 1]
    , assv <- sequence (replicate (arities !! f) [0 .. n - 1])
    ]

-- see if two models are isomorphic, by testing all permutations.
-- (p !!) computes the permutation
-- the mysterious `((p!!).)` was generated by http://pointfree.io/ from `(\f -> \i -> p !! (f i))`.
-- It is understandable, but I couldnt have come up with that ;)
isom :: Int -> [Int] -> [[Int] -> Int] -> [[Int] -> Int] -> Bool
isom n arities m1 m2 =
  or
    [ ident n arities (map (flip (\i -> ($map (p !!) i))) m1) $map ((p !!) .) m2
    | p <- permutations [0 .. n - 1]
    ]

-- used to partition the valid models into isomophic classes
-- This was taken from one of my projects, which ultimately
-- comes from some StackOverflow code snippet, which Ive lost the source
splitOffFirstGroup :: (a -> a -> Bool) -> [a] -> ([a], [a])
splitOffFirstGroup equal xs@(x:_) = partition (equal x) xs
-- splitOffFirstGroup _     []       = ([],[])  -- Not really necessary

equivalenceClasses _ [] = 0
equivalenceClasses equal xs =
  let (fg, rst) = splitOffFirstGroup equal xs
   in 1 + equivalenceClasses equal rst

-- putting everything together...
answer 0 _ _ = 1  -- corner case
answer n arities equations =
  equivalenceClasses (isom n arities) $
  filter (checkValid n arities equations) $ allModels n arities

Answer 2

Python 3, 393 382 bytes

-11 bytes thanks to @ovs

from itertools import*
T=tuple
P=product
def C(A,e,n):
 G=range(n);t=0;c=[] 
 for m in P(*[[dict(zip(P(*[G]*a),o))for o in P(*[G]*n**a)]for a in A]):
  E=lambda x,f:f[x]if x*0==0else m[x[0]][T(E(s,f)for s in x[1:])]
  if m not in c: c+=[T({T(r[a]for a in k):r[f[k]]for k in f}for f in m)for r in permutations(G)];t+=all(E(q[0],f)==E(q[1],f)for q,N in e for f in P(*[G]*N))
 return t

Try it online!

Takes each equation in the list as a pair of two expressions (LHS and RHS) and the number of free variables. Each expression is a tree representation, where the first non-negative integer of each list represents the index of the function, and the rest of the list are the function's arguments. (Uncomment the second-to-last line of the TIO link to see what the input looks like). The footer performs a bunch of parsing on the input string before passing to the function.

This generates all possible mappings of inputs to outputs for each function and tests if each one works. Ideally implemented, this general algorithm would have worst-case time complexity something like \$O(\frac{n^{\sum A}LEn^L}{n!})\$ (very slow), where \$n\$ is the order the given, \$\sum A\$ is the sum of the arities of the functions, \$L\$ is the average length of each equation (the whole equation is evaluated), and \$E\$ is the number of equations. There are \$n^{\sum A}\$ possible models, but sets of \$n!\$ on average are mutually isomorphic. Each class of models gets tested once (Ideally, we could ignore all but \$1\$ of each class, but it is not trivial to skip over isomorphic models. Hence, there is probably a full \$n^{\sum A}\$ term added for my implementation, and even higher since it checks for dict membership in a list of dicts). In the worst case, all of these \$n^{\sum A}/n!\$ models are valid, but they have to be verified by checking all \$E\$ equations and all \$n^L\$ possible inputs, each of which take (on average across the equations) \$L\$ time to evaluate. In the average case, I"m guessing only the log of this number has to be checked, so the average case would be roughly \$O(\frac{n^{\sum A}L\log(En^L)}{n!})=O(\frac{n^{\sum A}L^2\log(E)}{n!})\$. In the best case, for example with the equation x=y, each model gets invalidated with the first or second check (0=0 ... 0=1), so best case is \$O(\frac{n^{\sum A}}{n!})\$.

With these unwieldy time complexities, I have had to (greatly) shrink some of the orders input in order to run on TIO. The OEIS sequences that list millions of elements were generated with some specialized algorithm or formula. However, this has to run in the general case, so it can't easily make logical deductions that a human can (SAT is NP-hard, so I doubt the time complexity can be much improved).

Commented golfed code

from itertools import*

T=tuple
P=product

def C(A,e,n):
    G=range(n);
    t=0;# t: total number of models found so far
    c=[]    # c: list of isomorphic models of those already considered
            # A set would be better for time, but it takes a few too many bytes
            # to stringify dicts in a consistent way
    for m in P(   # m is a model
        *[
            [
                dict(zip(
                    P(*[G]*a), # inputs
                    o
                ))
                for o in P(*[G]*n**a) # possible outputs
            ]
            for a in A    # Why does this have to be enclosed in a list?
        ]   # P: choose one mapping for each function
    ):
        # I place the definition of E here so it can take advantage of scope to access m,n
        E=lambda x,f:f[x]if x*0==0else m[x[0]][T(E(s,f)for s in x[1:])] # E: evaluate an eXpression given free variables f
        # haven't already considered an isomorphic model:
        if m not in c:\
            c+=[    # append all isomorphic models
                T(
                    {
                        T(r[a]for a in k):r[f[k]]  # apply mapping
                        for k in f     # to keys (inputs) and values (outputs)
                    }
                    for f in m
                )
                for r in permutations(G) # r is a remapping, e.g. [0,2,1] on n=3 swaps 2 and 1
            ];\
            t+=all( # is the model valid?
                E(q[0],f)==E(q[1],f)        # LHS and RHS are equal
                    for q,N in e            # when given each possible list of free variables
                    for f in P(*[G]*N)
            )
    return t

General idea

from itertools import*

def execute(expr,free_vars,model,n):
    try:
        func_id, *params = expr
        params = tuple(execute(subexpr,free_vars,model,n) for subexpr in params)
        return model[func_id][params]
    except:
        # expr is an int ==> a variable
        return free_vars[expr]

def get_all_models():
    return product(
        *[
            [
                {i:o for i,o in zip(
                    product(*[range(n)]*arity), # inputs
                    output
                )}
                for output in product(*[range(n)]*n**arity)
            ] # possible outputs
            for arity in arities    # Why does this have to be enclosed in a list?
        ]
    )

def relabellings(model,n):
    return [
        tuple(
            {
                tuple(relabelling[a] for a in k): relabelling[e]
                for k,e in func.items()
            }
            for func in model
        )
        for relabelling in permutations(range(n))
    ]

def model_works(model, equations, n):
    return all(
        execute(eq[0],free_vars,model,n) == execute(eq[1],free_vars,model,n)
            for eq, num_free_vars in equations
            for free_vars in product(*[range(n)]*num_free_vars)
    )

def count(arities, equations, n):
    total=0
    considered = []
    for model in get_all_models():
        if model in considered:
            continue
        for r in relabellings(model, n):
            considered.append(r)
        if model_works(model, equations, n):
            total += 1
    return total