How to generate all involutive permutations?

I wondered how well Involutions has been implemented, so I tried to reimplement it myself. The following implementation can be up to 15 times faster than Involutions.

Needs["Combinatorica`"];
PackedQ = Developer`PackedArrayQ;
ToPack = Developer`ToPackedArray;

myInvolutions[list_List] := Block[{data, A, n, g},
  n = Length[list];
  A = ToPack[
    UpperTriangularize[{Range[3, n]}[[ConstantArray[1, n - 1]]]]];
  A[[2 ;;]] += ToPack[LowerTriangularize[{Range[2, n - 1]}[[ConstantArray[1, n - 2]]]]];
  g[2, {{}}] = ToPack[{{{1, 2}}}];
  g[n_, data_] := With[{m = Length[data]},
    Join[
     Transpose[{{
        ConstantArray[1, {(n - 1) m}],
        Flatten[Table[ConstantArray[i, {m}], {i, 2, n}]]
        }},
      {2, 3, 1}
      ],
     ArrayReshape[
      A[[1 ;; (n - 1), Flatten[data]]], {(n - 1) m, 
      Quotient[n - 2, 2], 2}
      ], 
    2]
    ];

  data = {{}};
  Join @@ Join[
    {{list}},
    Table[
     data = g[2 i, data];
     getPermutationLists[
      list,
      ArrayReshape[
       Subsets[Range[n], {2 i}][[All, Flatten[data]]],
       {Binomial[n, 2 i] Length[data], Sequence @@ Rest[Dimensions[data]]}
       ]
      ],
     {i, 1, Quotient[n, 2]}]
    ]
  ]

getPermutationLists = Compile[{{ran, _Integer, 1}, {idx, _Integer, 2}},
   Block[{a = ran, i, j, k, x},
    Do[
     i = Compile`GetElement[idx, k, 1];
     j = Compile`GetElement[idx, k, 2];
     x = Compile`GetElement[a, i];
     a[[i]] = Compile`GetElement[a, j];
     a[[j]] = x,
     {k, 1, Length[idx]}
     ];
    a
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

Here is a test:

n = 14;
a = Range[n];
aa = myInvolutions[a]; // AbsoluteTiming // First
bb = Involutions[a]; // AbsoluteTiming // First
Sort[aa] == Sort[bb]

5.63692

83.2192

True

Hard to beat Henrik's answer, but I'd like to show how you could do it by hand, without using a package.

From the decomposition of permutations into cycles, you can see that involutive permutations are the compositions of fixed points and transpositions.

Accordingly, you can build all the fixed points i (for identity) and transpositions t with:

n = 10;
a = Range[n];
t = Subsets[a, {2}];
i = Subsets[a, {1}];

Then, the transpositions will be among

s = Subsets[t, Floor[n/2]];

Unfortunately, s also includes non-disjoint cycles (like {1, 2}, {1, 3}). I don't know how to assemble only disjoint cycles together, but this is an unefficient brute-force approach:

inv = Select[s, Length[Flatten@#] == Length[DeleteDuplicates@Flatten@#] &];

That's a bit stupid because it builds one million subsets for $n=10$, to keep only 1000 in the end. This can be improved. Note however that contrary to your fully brute-force approach, here it's only exhaustive for transpositions, so it's still much much better.

Now, the number of permutations is given by

Length@inv
(* 9496 *)

And you can also rebuild the permutations matrices from inv:

cycles = Cycles /@ inv;
invo = Permute[Range[n], #] & /@ cycles
And @@ Map[PermutationProduct[#, #] == Range[n] &, invo]
(* True *)

The partially obsolete Combinatorica package has a function called Involutions that returns all involutive permutations:

Needs["Combinatorica`"];
a = Range[10];
inv = Involutions[a];
And @@ Map[PermutationProduct[#, #] == a &, inv]

True

Note also the existence of NumberOfInvolutions. Guess what it does...

The code for Involutions is actually available and can be printed as follows:

Needs["Combinatorica`"];
Needs["GeneralUtilities`"];
PrintDefinitions[Involutions]