Although this is my own question, I am going to give an answer showing one possible algorithm, to get the ball rolling. One way to compute the Bézier curve is De Casteljau's algorithm, which is a numerically stable computation method that uses a recursive method related to the recursive property of the binomial distribution. This algorithm can be implemented either on the standard probability scale, or on the log-probability scale.
Recursive characterisation of the Bézier curve: To obtain a recursive characterisation of this function, we will take advantage of the well-known recursive equation:
$$\text{Bin}(x|n, \theta) = (1-\theta) \cdot \text{Bin}(x|n-1, \theta) + \theta \cdot \text{Bin}(x-1|n-1, \theta).$$
Using this recursive equation, for any argument vector $\mathbf{f} = (f_0,...,f_n)$, we have:
$$\begin{equation} \begin{aligned}
B(\mathbf{f}, \theta)
&= \sum_{x=0}^n f_x \cdot \text{Bin}(x|n,\theta) \\[6pt]
&= \sum_{x=0}^n f_x \Big[ (1-\theta) \cdot \text{Bin}(x|n-1, \theta) + \theta \cdot \text{Bin}(x-1|n-1, \theta) \Big] \\[6pt]
&= (1-\theta) \times \sum_{x=0}^{n-1} f_x \cdot \text{Bin}(x|n-1, \theta) + \theta \times \sum_{x=0}^{n-1} f_{x+1} \cdot \text{Bin}(x|n-1, \theta) \\[6pt]
&= (1-\theta) \cdot B(\triangleleft \ \mathbf{f}, \theta) + \theta \cdot B(\triangleright \ \mathbf{f}, \theta), \\[6pt]
\end{aligned} \end{equation}$$
where $\triangleleft \ \mathbf{f} = (f_0,...,f_{n-1})$ and $\triangleright \ \mathbf{f} = (f_1,...,f_{n})$. This gives us a recursive equation for the Bézier curve, where the recursion decomposes the function into a weighted sum of the Bézier curve for smaller control vectors. At each step of the iteration, the length of the control vector is reduced by one element.
De Casteljau's algorithm: This algorithm takes advantage of the above recursive equation for the Bézier curve. To explain the algorithm, we define the operators $\triangleleft$ and $\triangleright$ to truncate the argument vector by one element from the right and left respectively. Now, taking a fixed value of $\theta$ we can define the values:
$$B_{k,x} \equiv B(\triangleleft^{n-k-x} \triangleright^x \mathbf{f},\theta),$$
and we can arrange these values in an $(n+1) \times (n+1)$ matrix $\mathbf{B}$ as follows:
$$\mathbf{B} \equiv \begin{bmatrix}
B_{0,0} & B_{0,1} & B_{0,2} & \cdots & B_{0,n-2} & B_{0,n-1} & B_{0,n} \\
B_{1,0} & B_{1,1} & B_{1,2} & \cdots & B_{1,n-2} & B_{1,n-1} & 0 \\
B_{2,0} & B_{2,1} & B_{2,2} & \cdots & B_{2,n-2} & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
B_{n-2,0} & B_{n-2,1} & B_{n-2,2} & \cdots & 0 & 0 & 0 \\
B_{n-1,0} & B_{n-1,1} & 0 & \cdots & 0 & 0 & 0 \\
B_{n,0} & 0 & 0 & \cdots & 0 & 0 & 0 \\
\end{bmatrix}.$$
The bottom element of this matrix is $B_{n,0} = B(\triangleleft^0 \triangleright^0 \mathbf{f},\theta) = B_\mathbf{f}(\theta)$, which is the function output we want to compute, and the top row consists of values $B_{0,x} = B(\triangleleft^{n-x} \triangleright^x \mathbf{f},\theta) = f_x$, which are the initial "control points" of the function.
De Casteljau's algorithm begins with the control vector at the top line of this matrix, and works downward through the above matrix of values by using the above recursive characterisation of the Bézier curve. In scalar form, the recursive equations for the algorithm are:
$$\begin{equation} \begin{aligned}
B_{0,x} &\equiv f_x
\quad \quad \quad \quad \quad \quad \quad
\quad \quad \quad \quad \quad \quad \quad \text{ }
\text{ for } x = 0,...,n, \\[6pt]
B_{k,x} &\equiv (1-\theta) \cdot B_{k-1,x} + \theta \cdot B_{k-1,x+1}
\quad \quad \quad \text{for } x = 0,...,n-k.
\end{aligned} \end{equation}$$
As can be seen from the matrix, this algorithm involves computation of $n(n+1)/2$ values from the initial set of control points, so it has complexity $\mathcal{O}(n^2)$. Each recursive computation is a simple weighted average of the elements above and above-right, so the computational burden of each step involves two multiplications and one addition (and is therefore relatively small).
The algorithm can be implemented in standard probability scale, so long as the computational environment is sufficient to avoid underflow problems for small probabilities. Alternatively, the computation can be done in log-probability scale to avoid underflow problems (see here for discussion of adding small probabilities in log-scale).
