Multivariate multipoint evaluation is the problem of evaluating a multivariate polynomial, given as a coefficient vector, simultaneously at multiple evaluation points. In this work, we show that there exists a deterministic algorithm for multivariate multipoint evaluation over any finite field $\mathbb{F}$ that outputs the evaluations of an $m$-variate polynomial of degree less than $d$ in each variable at $N$ points in time
\[
(d^m+N)^{1+o(1)}\cdot poly(m,d,\log|\F|)
\]
for all $m\in\mathbb{N}$ and all sufficiently large $d\in \mathbb{N}$.

A previous work of Kedlaya and Umans (FOCS 2008, SICOMP 2011) achieved the same time complexity when the number of variables $m$ is at most $d^{o(1)}$ and had left the problem of removing this condition as an open problem. A recent work of Bhargava, Ghosh, Kumar and Mohapatra (STOC 2022) answered this question when the underlying field is not too large and has characteristic less than $d^{o(1)}$. In this work, we remove this constraint on the number of variables over all finite fields, thereby answering the question of Kedlaya and Umans over all finite fields.

Our algorithm relies on a non-trivial combination of ideas from three seemingly different previously known algorithms for multivariate multipoint evaluation, namely the algorithms of Kedlaya and Umans, that of Bj\"orklund, Kaski and Williams (IPEC 2017, Algorithmica 2019), and that of Bhargava, Ghosh, Kumar and Mohapatra, together with a result of Bombieri and Vinogradov from analytic number theory about the distribution of primes in an arithmetic progression.

We also present a second algorithm for multivariate multipoint evaluation that is completely elementary and in particular, avoids the use of the Bombieri--Vinogradov Theorem. However, it requires a mild assumption that the field size is bounded by an exponential-tower in $d$ of bounded height. More specifically, our second algorithm solves the multivariate multipoint evaluation problem over a finite field $\mathbb{F}$ in time
\[
(d^m+N)^{1+o(1)}\cdot poly(m,d,\log |\mathbb{F}|)\]
for all $m\in \mathbb{N}$ and all sufficiently large $d\in \mathbb{N}$, provided that the size of the finite field $\mathbb{F}$ is at most $(\exp(\exp(\exp(\cdots (\exp(d)))))$, where the height of this tower of exponentials is fixed.