Counting groups of a given size Groups
In abstract algebra, a group is a tuple \$(G,\ast)\$, where \$G\$ is a set and \$\ast\$ is a function \$G\times G\rightarrow G\$ such that the following holds:

For all \$x, y, z\$ in \$G\$, \$(x\ast y)\ast z=x\ast(y\ast z)\$.

There exists an element \$e\$ in \$G\$ such that for all \$x\$ in \$G\$, \$x\ast e=x\$.

For each \$x\$ in \$G\$, there exists an element \$y\$ in \$G\$ such that \$x\ast y=e\$.


The order of a group \$(G,\ast)\$ is defined as the number of elements of \$G\$.
For each strictly positive integer \$n\$, there exists at least one group of order \$n\$. For example, \$(C_n,+_n)\$ is such a group, where \$C_n=\{0,...,n-1\}\$ and \$x+_ny=(x+y)\mod n\$.
Isomorphic groups
Let \$G:=\{1,2\}\$ and define \$\ast\$ by \$x\ast y=(x\times y)\mod3\$. Then \$1\ast1=1=2\ast2\$ and \$1\ast2=2=2\ast1\$.
Likewise, \$0+_20=0=1+_21\$ and \$0+_21=1=1+_20\$.
Although the elements and operations of the groups \$(G,\ast)\$ and \$(C_2,+_2)\$ have different names, the groups share the same structure.
Two groups \$(G_1,\ast_1)\$ and \$(G_2,\ast_2)\$ are said to be isomorphic if there exists a bijection \$\phi:G_1\rightarrow G_2\$ such that \$\phi(x\ast_1y)=\phi(x)\ast_2\phi(y)\$ for all \$x,y\$ in \$G_1\$.
Not all groups of the same order are isomorphic. For example, the Klein group is a group of order 4 that is not isomorphic to \$(C_4,+_4)\$.
Task
Write a program or function that accepts a non-negative integer n as input and prints or returns the number of non-isomorphic groups of order n.
Test cases
Input   Output
0       0
1       1
2       1
3       1
4       2
5       1
6       2
7       1
8       5
9       2
10      2
11      1
12      5
13      1
14      2
15      1
16      14
17      1
18      5
19      1
20      5

(taken from OEIS A000001)
Additional rules

There are no limits on execution time or memory usage.

Built-ins that trivialize this task, such as Mathematica's FiniteGroupCount, are not allowed.

Standard code-golf rules apply.
CJam, 189 187 bytes

This one's gonna be tough to explain... Time complexity is guaranteed to be O(scary).

qi:N_3&gt;{,aN*]N({{:L;N,X)-e!{X)_@+L@@t}%{X2+&lt;z{_fe=:(:+}%:+!},}%:+}fX{:G;N3m*{_~{G@==}:F~F\1m&gt;~F\F=}%:*},:L,({LX=LX)&gt;1$f{\_@a\a+Ne!\f{\:M;~{M\f=z}2*\Mff==}:|{;}|}\a+}fX]:~$e`{0=1=},,}{!!}?


If you're brave enough, try it online. On my crappy laptop I can get up to 6 with the Java interpreter or 5 in the online interpreter.

Explanation

I don't have a big math background (just finished high school, starting CS at uni next week). So bear with me if I make mistakes, state the obvious, or do things in horribly inefficent ways.

My approach is a brute force, though I tried to make it a little more clever. The main steps are:


Generate all the possible operands ∗ for a group of order n (i.e., enumerate all groups of order n);
Generate all the possible bijections φ between two groups of order n;
Using the results from steps 1 and 2, determine all isomorphisms between two groups of order n;
Using the result from step 3, count the number of groups up to isomorphism.


Before looking at the how each step is done, let's get some trivial code out of the way:

qi:N_             e# Get input as integer, store in N, make a copy
     3&gt;{...}    ? e# If N &gt; 3, do... (see below)
            {!!}  e# Else, push !!N (0 if N=0, 1 otherwise)


The following algorithm doesn't work correctly with n &lt; 4, the cases from 0 to 3 are handled with a double negation.

From now on, the elements of a group will be written as {1, a, b, c, ...}, where 1 is the identity element. In the CJam implementation, the corresponding elements are {0, 1, 2, 3, ...}, where 0 is the identity element.

Let's start from step 1. Writing all possible operators for a group of order n is equivalent to generating all the valid n×n Cayley tables. The first row and column are trivial: they are both {1, a, b, c, ...} (left-to-right, up-to-down).

      e# N is on the stack (duplicated before the if)
,a    e# Generate first row [0 1 2 3 ...] and wrap it in a list
  N*  e# Repeat row N times (placeholders for next rows)
    ] e# Wrap everything in a list
      e# First column will be taken care of later


Knowing that a Cayley table is also a reduced Latin square (due to the cancellation property) allows to generate the possible tables row-by-row. Starting from the second row (index 1), we generate all the unique permutations for that row, keeping the first column fixed to the value of the index.

N({                                 }fX e# For X in [0 ... N-2]:
   {                            }%      e#   For each table in the list:
    :L;                                 e#     Assign the table to L and pop it off the stack
       N,                               e#     Push [0 ... N-1]
         X)                             e#     Push X+1
           -                            e#     Remove X+1 from [0 ... N-1]
            e!                          e#     Generate all the unique permutations of this list
              {         }%              e#     For each permutation:
               X)_                      e#       Push two copies of X+1
                  @+                    e#       Prepend X+1 to the permutation
                    L@@t                e#       Store the permutation at index X+1 in L
                          {...},        e#     Filter permutations (see below)
                                  :+    e#   Concatenate the generated tables to the table list


Not all those permutations are valid, of course: each row and column must contain all the elements exactly one time. A filter block is used for this purpose (a truthy value keeps the permutation, a falsy one removes it):

X2+                 e# Push X+2
   &lt;                e# Slice the permutations to the first X+2 rows
    z               e# Transpose rows and columns
     {        }%    e# For each column:
      _fe=          e#   Count occurences of each element
          :(        e#   Subtract 1 from counts
            :+      e#   Sum counts together
                :+  e# Sum counts from all columns together
                  ! e# Negate count sum:
                    e#   if the sum is 0 (no duplicates) the permutation is kept
                    e#   if the sum is not zero the permutation is filtered away


Note that I am filtering inside the generation loop: this makes the code a bit longer (compared to distinct generation and filtering), but greatly improves performance. The number of permutations of a set of size n is n!, so the shorter solution would require a lot of memory and time.

A list of valid Cayley tables is a great step towards enumerating the operators, but being a 2D structure, it can't check for associativity, which is a 3D property. So next step is filtering out non-associative functions.

{                                 }, e# For each table, keep table if result is true:
 :G;                                 e#   Store table in G, pop it off the stack
    N3m*                             e#   Generate triples [0 ... N-1]^3
        {                     }%     e#   For each triple [a b c]:
         _~                          e#     Make a copy, unwrap top one
           {    }:F                  e#     Define function F(x,y):
            G@==                     e#       x∗y (using table G)
                   ~F                e#     Push a∗(b∗c)
                     \1m&gt;            e#     Rotate triple right by 1
                         ~           e#     Unwrap rotated triple
                          F\F        e#     Push (a∗b)∗c
                             =       e#     Push 1 if a∗(b∗c) == (a∗b)∗c (associative), 0 otherwise
                                :*   e#   Multiply all the results together
                                     e#   1 (true) only if F was associative for every [a b c]


Phew! Lots of work, but now we have enumerated all groups of order n (or better, the operations on it - but the set is fixed, so it's the same thing). Next step: find isomorphisms. An isomorphism is a bijection between two of those groups such that φ(x ∗ y) = φ(x) ∗ φ(y). Generating those bijections in CJam is trivial: Ne! will do it. How can we check them? My solution starts from two copies of the Cayley table for x ∗ y. On one copy, φ is applied to all elements, without touching the order of rows or columns. This generates the table for φ(x ∗ y). On the other one the elements are left as they are, but rows and columns are mapped through φ. That is, the row/column x becomes the row/column φ(x). This generates the table for φ(x) ∗ φ(y). Now that we have the two tables, we just have to compare them: if they are the same, we have found an isomorphism.

Of course, we also need to generate the pairs of groups to test isomorphism on. We need all the 2-combinations of the groups. Looks like CJam has no operator for combinations. We can generate them by taking each group and combining it only with the elements following it in the list. Fun fact: the number of 2-combinations is n × (n - 1) / 2, which is also the sum of the first n - 1 natural numbers. Such numbers are called triangular numbers: try the algorithm on paper, one row per fixed element, and you'll see why.

:L                          e# List of groups is on stack, store in L
  ,(                        e# Push len(L)-1
    {                  }fX  e# For X in [0 ... len(L)-2]:
     LX=                    e#   Push the group L[X]
        LX)&gt;                e#   Push a slice of L excluding the first X+1 elements
            1$              e#   Push a copy of L[X]
              f{...}        e#   Pass each [L[X] Y] combination to ... (see below)
                            e#   The block will give back a list of Y for isomorphic groups
                    \a+     e#   Append L[X] to the isomorphic groups
                          ] e# Wrap everything in a list


The code above fixes the first element of the pair, L[X], and combines it with other groups (let's call each one of those Y). It passes the pair to a isomorphism test block that I'll show in a moment. The block gives back a list of values of Y for which L[X] is isomorphic to Y. Then L[X] is appended to this list. Before understanding why the lists are set up in such a way, let's look at the isomorphism test:

\_@                                      e# Push a copy of Y
   a\a+                                  e# L[X] Y -&gt; [L[X] Y]
       Ne!                               e# Generate all bijective mappings
          \f{                    }       e# For each bijection ([L[X] Y] extra parameter):
             \:M;                        e#   Store the mapping in M, pop it off the stack
                 ~                       e#   [L[X] Y] -&gt; L[X] Y
                  {     }2*              e#   Repeat two times (on Y):
                   M\f=                  e#     Map rows (or transposed columns)
                       z                 e#     Transpose rows and columns
                                         e#     This generates φ(x) ∗ φ(y)
                           \Mff=         e#   Map elements of L[X], generates φ(x ∗ y)
                                =        e#   Push 1 if the tables are equal, 0 otherwise
                                  :|     e#   Push 1 if at least a mapping was isomorphic, 0 otherwise
                                    {;}| e#   If no mapping was isomorphic, pop the copy of Y off the stack


Great, now we have a list of sets like [{L[0], Y1, Y2, ...}, {L[1], Y1, ...}, ...]. The idea here is that, by transitive property, if any two sets have at least one element in common then all the groups in the two sets are isomorphic. They can be aggregated into a single set. As L[X] will never appear in the combinations generated by L[X+...], after aggregating each set of isomorphisms will have one unique element. So, to get the number of isomorphisms, it's sufficient to count how many groups appear exactly once in all sets of isomorphic groups. To do this, I unwrap the sets so they look like [L[0], Y1, Y2, ..., L[1], Y1, ...], sort the list to create clusters of the same group and finally RLE-encode it.

:~            e# Unwrap sets of isomorphic groups
  $           e# Sort list
   e`         e# RLE-encode list
     {    },  e# Filter RLE elements:
      0=      e#   Get number of occurrences
        1=    e#   Keep element if occurrences == 1
            , e# Push length of filtered list
              e# This is the number of groups up to isomorphism


That's all, folks.