Hamilton cycle decompositions of the complete graph I'm looking for the number of Hamilton cycle decompositions of the labelled complete graph $K_n$ for small $n$.  From such a decomposition, we can construct a special type of Latin square (called a row-Hamiltonian Latin square).
Edit:  Clearly, we require $n$ to be odd.  To ensure that each Hamilton cycle decomposition is counted once, we only include the $n$-cycle permutations $\alpha$ of $\{1,2,\ldots,n\}$ that have $\alpha(1)<\alpha^{-1}(1)$.  We also write the decomposition $\alpha\beta\ldots$ such that $\alpha(1)<\beta(1)<\cdots$.
The count for $n=3$ is $1$ counting (123).  The count for $n=5$ is $6$, counting the following:  $(12345)(13524)$, $(12354)(13425)$, $(12453)(14325)$, $(12435)(13254)$, $(12543)(14235)$ and $(12534)(13245)$.  Assuming my code is correct, the count for $n=7$ is $960$.
 A: In Two-factorizations of complete graphs it is stated that $K_9$ has 122 non-isomorphic Hamiltonian decompositions, and the corresponding number for $K_{11}$ is 3140 (EDIT: the actual figure is much more than this - see comment). I don't think they know any other values. (Sloane's database does not have any sequences with these numbers in.)
Now you are interested in the labeled case, which may be easier. However I have not been able to find anything (on Google).
A: Just reporting that I wrote another algorithm for this and found the following values:
3 1
5 6
7 960
9 40037760

I ran this through the superseeker on Sloane and it came up with nothing (so perhaps nobody has counted these before).
Here's my code below (it uses GAP).  We generate a (n-1) x n Latin rectangle where each row is an n-cycle and the i-th and (i+(n-1)/2)-th rows are inverses.
EnumerateHamiltonDecompositionsBacktrackingAlgorithm:=function(n,L,step)
  local i,j,k,count,A;
  i:=Int((step-1)/n)+1;
  j:=(step-1) mod n+1;
  count:=0;

  if(n mod 2=0 or n<3) then return fail; fi;
  if(j=1) then A:=[Minimum(Filtered([2..n],i->ForAll([1..n-1],t->L[t][1]<>i)))]; else A:=Filtered([1..n],s->ForAll([1..n-1],t->L[t][j]<>s) and ForAll([1..n],t->L[i][t]<>s)); fi;
  for k in A do
    L[i][j]:=k;
    L[i+(n-1)/2][k]:=j;
    if((j=n and CycleLengths(PermList(L[i]),[1..n])=[n]) or j<n) then
      if(i=(n-1)/2 and j=n) then count:=count+1;
      else count:=count+EnumerateHamiltonDecompositionsBacktrackingAlgorithm(n,L,step+1); fi;
    fi;
    L[i][j]:=0;
    L[i+(n-1)/2][k]:=0;
  od;
  return count;
end;;

EnumerateHamiltonDecompositions:=function(n)
  local L;
  if(n mod 2=0 or n<3) then return fail; fi;
  if(n=3) then return 1; fi;
  L:=List([1..n-1],i->List([1..n],j->0));
  L[1]:=List([1..n],i->i mod n+1);
  L[1+(n-1)/2]:=ListPerm(Inverse(PermList(List([1..n],i->i mod n+1))));
  return Factorial(n-2)*EnumerateHamiltonDecompositionsBacktrackingAlgorithm(n,L,n+1);
end;;

The extra data point comes from assuming that (12..n) is one of the cycles, then multiplying the result by (n-2)!.  This is legitimate since each decomposition contains a unique cycle with the edge 12, and by permuting the remaining n-2 edges, we generate a unique decomposition with the cycle (12..n).  There are no automorphisms under this group action, so each orbit has cardinality (n-2)!.
A: Gah, I commented but my answer was wrong. I don't have a copy of Mathematica available, but here's (I think) a description of how to compute small cases in Mathematica.
There's a package (Combinatorica) with a function called HamiltonianCycle[graph, All] that returns a list of all the directed Hamiltonian cycles beginning and ending at a single node (as lists). Set the graph to be CirculantGraph[n, {2, 3, ..., (n-1)/2}] and compute this list. This is the graph resulting after we remove the first Hamiltonian cycle.
Now if you're doing Mathematica, it counts directed cycles, so we only want to consider half the lists. Throw out every cycle where the second element is larger than the second-to-last element (the first and last elements are both 1). (N.B. I originally described this step incorrectly, whence the comments.) Create this sublist, which we'll call hamcyc, and then compute
partitions := Subsets[hamcyc, (n-1)/2].
This is a 3D array. Count the number of elements (2D arrays) in this such that every pair of distinct integers in {1, ..., n} is contained as adjacent elements in exactly one of the lists in this 2D array. (Not sure how to do this, but this is the only thing I don't know how to do.)
Multiply this count by n!/(n-1) to get the number of partitions into Hamiltonian cycles.
