Approximating a function's derivative in machine precision Let's say I have a very complex analytical function $f$, in which finding its derivative at point $Q$, i.e. $f'(Q)$ is not practical. So we resort to finding the derivative numerically without analytically differentiate the function at point $Q$ (note that $Q$ is only a point, not an array).
Method 1: I have tried the $d$ order of finite difference. This doesn't give me a machine precision result when the order is too high or singularity behaviour occurs when steps are too small.
Method 2: Interpolation method. This doesn't give machine precision accuracy neither when differentiate the respected interpolant at point $Q$. i.e. $f'(Q) ≅ \sum_i L_i'(Q)f(x_i)$ where $L_i(Q)$ is the interpolant with basis $i$ and $x_i$ is the nodal point.
Method 3: Spectral method such as Legendre expansion and the like. No machine precision.
Method 4: least square polynomial fit. Just no.
Method 5: minimax polynomial approximation. Seems alright with around 13 decimal places accuracy in matlab.
Any better method? Suggestion? Alternatives?
I will be much appreciated. Thank you all!

Obtaining numerical derivative up to machine precision using real arithmetics only is hard due to roundoff errors. 
If derivative is estimated using $d$-point differentiation method, then expected error using optimal step lenght is of order $\epsilon^\frac{d}{d+1}$ times term depending on derivative of order $d+1$.
