I believe this can be done on $O(n \log n)$ on a RAM, where $n=|G|$; I don't know a reference so I'll just write it down here (but surely this is not original). Nor do I know if this is the fastest known, but clearly you can't do it faster than $\Omega(n)$, so it's gotta be close.
Let's assume the group elements are denoted in the computer by the integers $1,\dotsc,n$, with $g_i \in G$ being the group element corresponding to the integer $i$. Assume WLOG that $g_1$ is the identity.
1) Find a generating set $\Gamma$ of size $\leq \log_2 n$ in $O(n \log n)$ time:
G = list of group elements
Gamma = []
Gr = new graph with vertex set G and no edges // O(n) steps
found = new array of length n, initialized to all false // O(n) steps
for i = 2 to n: // never need to add identity to generating set
if not found[i] then
Gamma.append(i)
// Now update the graph Gr by adding new edges
// Corresponding to the new generator G[i]=g_i
for each vertex g in Gr:
Gr.addEdge(g,g*G[i])
end for
do BFS on Gr starting from 1, ignoring any vertices already found
(for any vertex j encountered, this sets found[j]=true)
end if
end for
2) Build the following sparse graph (that is, in the edge list representation, not the dense adjacency matrix representation): vertices are the elements of $G$, and there is an (undirected) edge $(g,h)$ if $h = \gamma g\gamma^{-1}$ for some generator $\gamma \in \Gamma$. This is a graph with $n$ vertices and maximum degree $O(\log n)$. Its connected components are the conjugacy classes, and finding them can be done by BFS (say) in time $O(v + e) = O(n + n \log n) = O(n \log n)$.