Counting Moufang Loops
Python 3, 475 410 bytes

Thanks to Mr.Xcoder for saving some bytes!

Use the symmetry of the formula to save 65 bytes. Yes, that's a lot.



from itertools import*
n=int(input())
P=permutations
R=[*range(n)]
u=[]
A=all
S=sorted
for T in P(P(R),n):u+=[T]*(A(A(R==S(x)for x in
t)and any([*x]==S(x)for x in t)and
A(t[z][t[x][t[z][y]]]==t[t[t[z][x]][z]][y]and
t[t[z][x]][t[y][z]]==t[t[z][t[x][y]]][z]for x in R
for y in R for z in R)for t
in(T,[*zip(*T)]))and A(A(1-A(p[T[i][j]]==U[p[i]][p[j]]for i in R
for j in R)for p in P(R))for U in u))
print(len(u))


Try it online!



Some and can be replaced by *, results in lower bytecount but at the cost of significantly slower execution time:

Python 3, ??? bytes

[TODO put code here]

(of course not all * makes the program significantly slower, only some of them are critical)



Ungolfed:

from itertools import *
n = 4 # int(input())
rangeN = list(range(n))

def is_moufang_loop(T):
    A = tuple(zip(*T))
    return all(
        all(sorted(x) == rangeN for x in t)
        and any(list(x) == sorted(x) for x in t)
        and all(
                T[z][T[x][T[z][y]]] == T[T[T[z][x]][z]][y]
            and T[T[z][x]][T[y][z]] == T[T[z][T[x][y]]][z]
            for x in rangeN for y in rangeN for z in rangeN)
        for t in (T, A)
    )

def isomorphic(loop1, loop2):
    for p in permutations(rangeN):
        if all(
            p[loop1[i][j]] == loop2[p[i]][p[j]]
            for i in rangeN
            for j in rangeN
        ): return True
    return False

unique_moufang_loops = []
for x in [
        cayley_table 
        for cayley_table in permutations(permutations(rangeN), n)
        if is_moufang_loop(cayley_table)
]:
    if all(not isomorphic(x, y) for y in unique_moufang_loops):
        unique_moufang_loops.append(x)

print(len(unique_moufang_loops))


Try it online!

No scroll bars...



Explanation:

The program is pretty simple.


Each possible "binary operator" is represented by a Cayley table (0-indexing).
The "identity" property implies that there exists e such that both the e'th row and the e'th column are equal to [0, 1, 2, ..., n-1], which is the same condition as  


  both the array T and its transpose has a row equal to [0, 1, 2, ..., n-1].

The "cancellation" property is equivalent to  


  Every row and every column are a permutation of [0, 1, 2, ..., n-1].



So, the part

all(
        all(sorted(x) == rangeN for x in t) 
        and any(list(x) == sorted(x) for x in t) 
        for t in (T, A))


of the code checks that. (for all row in the array T and its transpose A, it being sorted is equal to rangeN, and there exists a row in both T and A that equals to itself being sorted)

The four conditions of a Moufang loops are checked manually.

z + (x + (z + y)) = ((z + x) + z) + y
((y + z) + x) + z = y + (z + (x + z))
(z + x) + (y + z) = (z + (x + y)) + z
(z + x) + (y + z) = z + ((x + y) + z)


In the code, (a + b) is represented as T[a][b]. (because of the representation as Cayley table). Use Python chaining equality comparison to avoid duplication of (z + x) + (y + z).

However, because the formula is symmetric:

If we switch the operands of + in the first formula, we get the second formula; and if we switch the operands of + in the third formula, we get the fourth formula with x and y swapped place.

Note that the transposition of the Cayley table is equivalent to the binary operators having operands swapped. (x + y -&gt; y + x)

Finally, all the candidates Cayley table are picked from

permutations(permutations(rangeN), n) 


so that each row is a permutation of rangeN (which is [0, 1, 2, ..., n-1]) and there are n distinct rows.