Stable and efficient computation of binomial expectations Suppose we want to compute the expected value of some function $f(X)$ where $X \sim \text{Bin}(n,\theta)$.  Taking $\mathbf{f} = (f_0,...,f_n)$ to be the function values over all possible outcomes of the binomial random variable, this expected value can be written as a function:
$$\begin{equation} \begin{aligned}
B(\mathbf{f}, \theta) \equiv \mathbb{E}(f(X)) 
&= \sum_{x=0}^n f_x \cdot \text{Bin}(x|n,\theta) \\[6pt]
&= \sum_{x=0}^n f_x {n \choose x} \theta^x (1-\theta)^{n-x}, \\[6pt]
\end{aligned} \end{equation}$$
where the function arguments are the vector $\mathbf{f} \in \mathbb{R}^{n+1}$ and the parameter $0 \leqslant \theta \leqslant 1$.  (Note that $n$ is determined as one less than the length of the vector $\mathbf{f}$, so it is not a separate argument in the function.)  In mathematical parlance, this function is called the Bézier curve, and the vector $\mathbf{f}$ gives the "control points" of the function.  This object is useful for a range of problems involving binomial random variables.  In particular, computation of this function subsumes the problems of computing the mass and distribution for the binomial distribution (which are easily obtained by using values of zeros and ones in the vector $\mathbf{f}$).
When $n$ is large, direct computation of this quantity is problematic, since the binomial coefficient becomes large, and the power terms become small, leading to arithmetic overflow and underflow problems.  This can potentially be dealt with by a range of methods, such as conversion to a log-scale, using a recursive algorithm, etc.  However, I am not sure what algorithms are actually used in practice in statistical computation, and what algorithm is considered the "gold standard".

My Questions: What are some efficient and stable methods of computation of this function?  What algorithms are used in practice in statistical computing?  Is there any algorithm for this problem that is considered to be the "gold standard"?
 A: 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).
