The *Online Encyclopedia of Integer Sequences* has recently celebrated its fiftieth anniversary, amassed over 250 000 sequences of varying importance, and heralded the 75th birthday of its creator Neil Sloane. It is somewhat surprising, therefore, that no-one appears to have investigated the intersection of two particularly natural integer sequences (catalogued as A000040 and A000179 in the OEIS).

In this article, I shall attempt to persuade you that investigating this new sequence of *ménage primes* is a good idea for several reasons. But until then, let us consider the problems that motivated this line of research.

### Vortices and knots

Before the structure of the atom was properly understood, Lord Kelvin conjectured that atoms were vortex rings in a hypothetical ubiquitous fluid known as the *luminiferous aether*. The shape of the vortex ring writhes and undulates, but *up to homotopy* it remains constant throughout its lifetime. Begun was the idea that chemical elements correspond to homotopy classes of embeddings of the circle (more commonly known as *knots*), and classifying them would yield a mathematical explanation for the origin of the periodic table*.

- In reality, the periodic table stems from solving the Schrödinger equation for the potential effected by a point charge (the nucleus).

Realising that knots exhibit ‘unique prime factorisation’, Tait et al began enumerating the prime knots according to the minimal number of crossings when the knot is drawn as a diagram in the plane:

Shown above are all distinct prime knots with eight or fewer crossings. With the exception of the last three, these are all *alternating knots*, meaning that the string passes alternately over and under the crossings.

A knot diagram can be considered to be a 4-regular plane graph, where each vertex is decorated to specify which path goes over and which goes under. We can add further structure by 2-colouring the regions and assigning a direction to the knot. For example, the second knot in the third row yields this diagram:

Each edge can be labelled from 1 to 2*n*, where *n* is the number of crossings, according to the order in which they are traversed. We define the *parity* of an edge to be whether the blue region is to the left or the right of the arrow on the edge. Observe that the parity alternates between successive edges in the traversal, so this corresponds exactly to the parity of the label.

At each vertex, an odd edge and an even edge converge. This yields a bijection between the odd and even edges of the diagram. Moreover, note that in a reduced diagram (with no trivial redundant crossings), consecutive edges never point to the same vertex.

Now draw a complete bipartite graph K(*n*, *n*), where each vertex corresponds to an edge in the diagram. We have a natural matching (the bijection between the odd and even edges of the graph described above) and a natural Hamiltonian cycle (the order of traversal). These share no edges. There are two dual ways to view this as a graph-theoretic problem:

- Counting Hamiltonian cycles of a complete bipartite graph with a perfect matching removed.
- Counting perfect matchings of a complete bipartite graph with a Hamiltonian cycle removed.

This problem was later rediscovered by Edouard Lucas in an entirely different setting.

### Formal combinatorics

At High Table, *n* Fellows and *n* guests (one for each Fellow) are seated around the perimeter of a topological disc. To maximise social interaction, the following rules are observed:

- No two Fellows may be seated adjacent to each other.
- No Fellow may be seated adjacent to his or her guest.

In how many ways is this possible? This is known as the *ménage problem*, due to an outmoded heteronormative statement of the problem involving married couples instead of ordered pairs of the form (Fellow, guest). One such arrangement for *n* = 7 is illustrated below and obtained from the knot diagram:

The crimson lines biject each Fellow (odd number) with his or her guest (even number); these correspond to edges converging to the same vertex in the knot diagram.

We can eliminate much redundancy in our enumeration by exploiting the vast amount of symmetry in the situation:

- Any permutation of the
*n*ordered pairs yields another valid arrangement (order n!). - Interchanging each Fellow simultaneously with his or her guest results in another valid arrangement (order 2).

Removing the symmetry is equivalent to imposing that the Fellows are initially seated in a particular arrangement. Hence, our question reduces to:

Given an initial arrangement of

nFellows occupying alternate seats, in how many ways can we seat their guests in the remainingnseats such that no Fellow is adjacent to his or her guest?

The number M(*n*) is called the *n*th *ménage number*. The first few ménage numbers are:

- M(3) = 1 (the unique arrangement involving the Fellows being diametrically opposite their respective guests)
- M(4) = 2 (two chiral solutions)
- M(5) = 13
- M(6) = 80
- M(7) = 579

An explicit formula is given in the OEIS (sequence A000179).

### Which ménage numbers are prime?

Observe that M(4) and M(5) are both prime, which confirms that we have successfully factored out all of the redundant symmetries. Continuing further, one discovers that the next *ménage prime* is M(13) = 775596313, followed a while afterwards by M(53) = 567525075138663383127158192994765939404930668817780601409606090331861.

I decided to do an overnight *Mathematica* run for values of *n* up to circa 10000. In the process, it found a spectacular ** 30281-digit** ménage prime, M(8645). To put this into perspective, this is larger than any prime known before the discovery of 2^132049 − 1 in 1983. This suggests the following questions:

- Are there infinitely many ménage primes?
- If so, what is their asymptotic growth rate?

The first of these is extremely difficult, and the second is at least as hard (although we provide some heuristics in the next section). A more achievable goal is the following:

Find another ménage prime beyond the five already discovered.

### Heuristics

Although the ménage numbers are not exactly equidistributed modulo *p* for all *p* (for example, only 2/9 of them are divisible by 3), they do seem to asymptotically be prime as often as one would expect from a randomly-chosen number of a similar size. That is to say, M(*n*) appears to be prime with probability 1/log(M(*n*)).

Now we need an asymptotic formula for the ménage numbers. Since M(*n*) counts the number of permutations such that both σ(*i*) and σ(*i*) + 1 are derangements, we expect there to be roughly *n*!/e² arrangements. It can be shown that this estimate is asymptotically correct (the ratio tends to 1 as *n* tends to infinity), which by Stirling’s formula is in turn asymptotically equal to:

The sum of the reciprocals of the logarithms of these numbers diverges, which (by Borel-Cantelli) suggests that there are infinitely many ménage primes. However, it diverges sufficiently slowly that we only expect around log(log(*m*)) such primes for *n* < *m*. This means that in the long run, the number of ménage primes below a given value N is expected to be around log(log(invgamma(N))), where invgamma is the inverse Gamma function.

With these heuristic estimates, the expected number of ménage primes M(*n*) with *n* between 10^4 and 10^5 is roughly 0.25. It is entirely plausible that people may attempt to search further than 10^5; increasing this upper bound to 1.5 × 10^7 (corresponding to primes of over 100 million digits, therefore eligible for the 150 000 dollar prize offered by the Electronic Frontier Foundation) gives an expected value closer to 0.64.

Anyway, good night for now — I have a couple of sequences to add to the OEIS (due to appear as A249510 and A249394).

“To put this into perspective, this is larger than any prime known before the discovery of 2^132049 − 1 in 1983.” And yet the number doesn’t rate an entry in Henri & Renaud Lifchitz’s current top 10000 probable prime list! A word of caution to anyone using PrimeQ in Mathematica 10: It’s not as robust as Mathematica 9 in its small-factors subroutine. I’d say it’s a bug.

Indeed, computing power has increased massively! #mooreslaw

I used PrimeQ in Mathematica 8, although I was generating the ménage numbers quite inefficiently in my preliminary overnight run (i.e. individually instead of using the much faster recurrence formula, and not bothering to only test odd-indexed numbers), so I’m sure I can go much further. And thanks for the suggestion, since I’m considering parallelising the task over a cluster of ~ 30 Mathematica-10-enabled CPUs to which I have access.

One can do much better than restrict oneself to odd-indexed n, but it’s a bit of a red herring insofar as the small-factor subroutine should solve these in microseconds — even for very large numbers. I count 4351 M(n) candidates for 8645<n<10^5, based on Mathematica 9 reaching a TimeConstrained two seconds to compute PrimeQ[M(n)] on my Mac. These then are the "difficult to factor" M(n) in that range, the implication being that anything else (being solved in less than two seconds) is necessarily False by virtue of having come across a small factor. While I did that candidate calculation I saved the values of M(n) so that I wouldn't have to recompute them for the full evaluation, which I started five days ago. I've just reached #916, which takes me to 28000.

After checking small factors, I presume PrimeQ attempts a Fermat-bash to ensure that it’s at least pseudoprime before subjecting it to heavy primality testing?

Heuristically one expects 0.129 ménage primes in the interval [28000, 10^5], so it looks like my dignity may survive a few more years, at the very least…

I believe that PrimeQ uses the Miller-Rabin strong pseudoprime test base 2 & base 3, and then the Lucas test.

Pingback: Holyhedron update | Complex Projective 4-Space