Test for group isomorphism and construct automorphism groups

All right, here is a solution: Find group isomorphisms in Mathematica. It's not pretty, but it's practical for groups of order up to about 100. It takes 30 ms to find out that $\text{d8a}\cong\text{d8b}$, and 43 ms to produce the automorphism group on my Mac mini. It finds an isomorphism from $S_5\to S_5$ (order 120) in 7 s. Producing all of $Aut(S_5)$ takes 27 s, since it has to test all possible mappings rather than stopping at the first success. It manages $S_6$ (order 720, $|\text{Aut}(S_6)|=1440$) in a little under 4 hours.

findGroupIsomorphisms[
  group1_,
  group2_,
  max_: ∞
  ] :=
 Module[{ng1, ng2, reverse, g1, g2, order, gs1, go1, gm1, gm2,
   os1, os2, inng1, abelianQ, targets, isomorphisms, w1,
   mt1, φm1, φ, homomorphismQ, t, ts},
  (*
   * Choose the domain group with fewest generators
   *)
  If[GroupOrder[group1] != GroupOrder[group2],
   Return[{}];
   ];
  {ng1, ng2} = Length[GroupGenerators[#]] & /@ {group1, group2};
  reverse = ng2 < ng1;
  If[reverse,
   {g1, g2} = {group2, group1};
   {ng1, ng2} = {ng2, ng1},
   {g1, g2} = {group1, group2}
   ];
  (*
   * Do some quick checks for isomorphism
   *)
  order = GroupOrder[g1];
  {gm1, gm2} = GroupElements /@ {g1, g2};
  {os1, os2} = Map[PermutationOrder, {gm1, gm2}, {2}];
  If[Sort[Tally[os1]] != Sort[Tally[os2]],
   Return[{}]
   ];
  (*
   * Pick possible targets in g2
   *)
  gs1 = GroupGenerators[g1];
  go1 = PermutationOrder /@ gs1;
  targets = Table[
    Pick[gm2, Thread[os2 == n]],
    {n, go1}
    ];
  targets = Tuples[targets];
  (*
   * List the inner automorphisms
   *)
  inng1 = Outer[
    GroupElementPosition[
      g1, #1\[PermutationProduct]#2\[PermutationProduct]\
InversePermutation[#1]] &,
    gm1, gs1
    ];
  inng1 = Union[inng1];
  abelianQ = Length[inng1] == 1;
  (*
   *
  Here I should do something much more efficient for abelian groups...
   *)
  (* Stub *)
  (*
   *
  Test every possible mapping of the generators
   *)
  w1 = GroupElementToWord[g1, #] & /@ gm1;
  mt1 = GroupMultiplicationTable[g1];
  isomorphisms = {};
  While[targets =!= {},
   t = targets[[-1]]; targets = Drop[targets, -1];
   φ = With[{t = t},
     Function[i, GroupElementFromWord[PermutationGroup[t], w1[[i]]]]
     ];
   ts = Map[φ, inng1, {2}];
   targets = Complement[targets, ts];
   If[GroupOrder[PermutationGroup[t]] != order,
    Continue[]
    ];
   φm1 =
    GroupElementFromWord[PermutationGroup[t], #] & /@ w1;
   homomorphismQ =
    Map[φ, mt1, {2}] ==
     Outer[PermutationProduct, φm1, φm1];
   If[homomorphismQ,
    (*
     * We found some isomorphisms!
     *)
    isomorphisms = Join[isomorphisms, ts]
    ];
   If[Length[isomorphisms] >= max,
    Break[]
    ]
   ];
  (*
   * Return the isomorphisms as rule lists
   *)
  isomorphisms =
   Take[isomorphisms, Min[max, Length[isomorphisms]]];
  isomorphisms = {gs1, #} & /@ isomorphisms;
  If[reverse, isomorphisms = Reverse /@ isomorphisms];
  Apply[Rule, Transpose /@ isomorphisms, {2}]
  ]


isomorphicGroupsQ[group1_, group2_] :=
 MatchQ[findGroupIsomorphisms[group1, group2, 1], {_}]

For future searchers who come across this question, there is now a function FindGroupIsomorphism in the Wolfram Function Repository.