take too long to execute or it may not converge at all if the likelihood landscape is too complex. It all depends on how much CPU time it takes to compute transits for one planetary configuration.

I developed a fully numerical code and tested it on mock planetary systems back in 2008. Things were annoyingly slow. So, I explored every avenue to speed things up. For example, do we gain anything if the N-body code is replaced by the Fourier routines that calculate TTVs analytically? If so, they can perhaps be used to determine an approximate solution, and the N-body integrator can take over after that.

It turns out that the 12 variables can be split in four groups. The first and the most difficult to deal with are the semimajor axes of the two planets. The coefficients of the Fourier series depend on them in non-trivial ways. Fortunately, they depend only on the ratio of the semimajor axis, $\alpha = a_1/a_2$, and not on $a_1$ and $a_2$ individually. So, it is possible to precompute all coefficients on a grid in $\alpha$ and devise a fast algorithm that interpolates from the grid to any value of $\alpha$. The likelihood of 