Have a look at:  https://dsp.stackexchange.com/questions/69186#6923

                                      Mhz

        500kHz  Bin width             0.5
         50MHz  Minimum tone           50
         60GHz  Maximum tone        60000  
-5MHz to +5MHz  Tone Range          +/- 5
        100MHz  Mean Tone Spacing     100 
        128GHz  Nyquist Frequency  128000

Translating into bins:

  Minimum bin =     100
  Maximum bin =  120000
  Bin range   =  +/- 10
  Spacing     =     200
  Nyquist bin =  256000
  N bins      =  512000

Have I got that right?

So you have room for about 600 tones, wow.  I really have no experience working with frames this large.  You are going to have an inherent error just with truncation errors, so you will not be able to get away with single precision floats, you will need doubles for sure.  With that many bins, you can almost think of your situation as a continuous case, and should have no problem using complex formulas, rather than the more complicated real ones.

"DFT Bin Value Formulas for Pure Complex Tones"

https://www.dsprelated.com/showarticle/1038.php

Equations (28) and (33) shows that the envelope of the magnitude of 'leakage' drops off like a hyperbola in terms of the difference of your frequency in cycles per frame (f) and the bin index (k)

The magnitude envelope with large N, is then roughly M / ( pi * ( x - f ) )

With spacing of 200 bins between the tones, you shouldn't need to worry about interference nearby peaks at any peak, and that certainly includes your conjugate peak on the other side of the spectrum.

There is a direct tradeoff between the frequency calculation error and the following phase calculation error.  I cover that in this StackExchange answer:

https://dsp.stackexchange.com/questions/61922/is-i...

The advantage of having the more precise real vs complex formulas only comes into play with small N and high SNR.  Still, you might want to take a look at

"A Two Bin Solution" 
https://www.dsprelated.com/showarticle/1284.php

It will give you more precise values, but not worth it if your noise level is too high or with large N.  But after calculating a 51200 bin FFT, the difference is trivial in terms of extra computations.

Doing the complex version of the formulas is certainly simpler calculation wise.  There are a lot of considerations that go into your decision of which approach to use, and you haven't really given any of the relevant specs there, like expected SNR is most important.

I am guessing that simply using Jacobsen's estimator will be sufficiently accurate for you with large N you have, and the complex phase calculation should be fine as well, as you already pointed out.  It is much simpler than the real tone that takes the interference of the conjugate in account.  Three bins for your phase calculation should also be adequate.

It really depends on what the quality of your signal is and how precise results you need.

I started in DSP by trying to develop a music compression algorithm, the first step being removing all the best fitting pure tones possible.  I'll skip the rest of the detail of the following steps.  In that case N is much smaller, and the tones much more closely spaced, so it required doing an iterative process of removing the already calculated tones from the DFT, the re-evaluating the current tone.  For an audio signal composed of a mix of pure tones, you can get the same accuracy as with a single tone calculation.

Hope this helps.

This posting is a continuation from the comments at Phase and Amplitude Calculation for a Pure Complex Tone in a DFT using Multiple Bins.

Thinking about it a little more, since you have such a considerable tone spacing, broadening your stance for the phase and magnitude calculation might not be a bad idea.  However, the same is not true for the frequency calculation.

Looking at the chart in Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT, it clearly shows that you should use the "The Projection of Z and and DZ" method, and if the frequency ends up near the center of two bins, switch to the two bin solution of A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT.