I am looking numerically compute some $$\partial_{x_0} \partial_{x_1} \dots \partial_{x_{n-1}}f$$ where $f = f(x_0, \dots, x_{n-1})$.
My question is to know if this can be done, at a given precision, with a polynomial numbers of terms to sum (and not exponentially many).
Indeed, using naively, $$\partial_{x_0}f \approx \frac{f(x_0 + \epsilon, ...)-f(x_0, ...)}{\epsilon}$$ iteratively, e.g. for the derivative in the second variable
$$\partial_{x_0}\partial_{x_1}f \approx \frac{f(x_0 + \epsilon, x_1 + \epsilon, ...)-f(x_0 + \epsilon, x_1, ...)-f(x_0, x_1 + \epsilon, ...) + f(x_0, x_1, ...) }{\epsilon^2}$$ we see that there is an exponential increase of terms in the numerator as we derive again and again. I would like to know if there is a method to evaluate $\partial_{x_0} \partial_{x_1} \dots \partial_{x_{n-1}}f$ where the number of terms to be summed increases at most polynomially in $n$, for a given $\epsilon$.