extracting ultra accurate phase and frequency with FFT with real world data
You only have 64 samples per FFT. That means your FFT bins are spaced f_sample/64; if your sample rate is 128·50 Hz, then your frequency resolution, without interpolation, is just 100 Hz. To achieve a frequency resolution of $0.01\,\text{Hz}=\frac{50\,\text{Hz}}{50\cdot 100} = \frac{\frac{f_{sample}}{128}}{5000}= \frac{f_{sample}}{640\cdot5000}$, you'd need an FFT of length 3.2 millions. You can't do that on an M4, obviously, so: the FFT is not a tool applicable to your problem here.

Two approaches:

Either, you first reduce the sampling rate (including sufficient analog anti-aliasing filters) or you decimate your signal (including sufficient digital anti-aliasing filters) to a rate that describes a much smaller bandwidth of potential frequencies that you actually care about. Afterwards, do the FFT. Notice that this always always trades length of observation for precision: it's really that simple, you can't get 1/100 Hz of resolution without 100 s of observation. If things change within that period: bad luck.

Or, you don't do that and instead use parametric estimators for frequency and phase. I kind of feel I recommend that here. Parametric estimators can be as simple as zero-crossing detectors, or as complex as Eigenvalue-based estimators, and what you use depends on a lot of factors that you don't discuss – it's hard to guess what you need.

However, I like to point people at Peter K.'s A Review of the Frequency Estimation and Tracking Problems. Many of these algorithms will deliver phase info. 

Regarding an accuracy requirement of 0.1° or 0.01 Hz in general: You can't have arbitrary precision with limited observation if you have any noise in your signal. Estimation theory doesn't allow for that. I think you'll be running into theoretical limits rather fast, depending on which algorithm you end up (and using the lowest-variance method might not always be the best choice!).