the fit is a very sensitive function of $\alpha$, so the grid must be sufficiently dense for things to work. \[In fact, it is even better to work with the orbital periods (rather than the semimajor axes), because of uncertainties in stellar mass.\]

The second group is the orbital eccentricities and inclinations. The coefficient dependence on them is simple: they appear in all powers permitted by symmetries. It turns out to be possible to develop recursive routines to compute the higher powers from the lower ones such that the number of arithmetical operations is minimized.

The third group are the orbital angles, three for each planet. A really efficient way to evaluate the Fourier series for any combination of angles is to make use of complex algebra and symmetries. The algorithm computes a few leading terms and combines them to get the remaining ones almost for free. Therefore, in essence, because the evaluation of the Fourier series for angles is so inexpensive, one can set aside the angle dimensions. This effectively reduces the dimension of the problem.

The last group is the masses of the two planets. The TTVs of the inner transiting planet are nearly independent of its own 