Implementing secret reconstruction in Shamir's Secret Sharing I am trying to implement Shamir's secret sharing in C++.
I have got the generation of shares working.

However, I am very confused with the reconstruction of shares. I get the part on how three users can reconstruct the secret. But what if there are more than 3 users needed to reconstruct the secret? How do I implement the formula then?

Any help would be appreciated.
Suppose $s_0, s_1, s_2, \ldots, s_{k-1}$ are 
elements from the finite field you are working in, where $s_0$ is the secret
to be shared, and the $s_i, i &gt; 0$ are randomly chosen nonzero elements of
the field. Then, the polynomial used to construct the shares is
$$S(x) = s_0 + s_1x + s_2x^2+ \cdots + s_{k-1}x^{k-1}$$ and the shares
themselves are $y_i = S(x_i), 1 \leq i \leq n$ where the $x_i$ are
$n$ distinct nonzero elements of the field. Note that $n$ must be at least as
large as $k$, and that the polynomial $S(z)$ can be recovered
from any $k$ of these $n$ shares. This is commonly done using
Lagrange interpolation as has been already pointed out in the other
answers.  Assuming without loss of generality that $y_1, y_2, \ldots, y_k$
are the shares available for reconstruction, $S(x)$ is a weighted
sum (or linear combination) of $k$ different
polynomials $g_i(x)$.  We have
$$\begin{align}
S(x) &amp;= \sum_{i=1}^k y_i \cdot g_i(x)\\
&amp;=\sum_{i=1}^k y_i \prod_{j=1, j\neq i}^k\frac{x-x_j}{x_i-x_j}\\
&amp;= y_1 \frac{(x-x_2)(x-x_3)\cdots(x-x_k)}{(x_1-x_2)(x_1-x_3)\cdots(x_1-x_k)}
+ y_2 \frac{(x-x_1)(x-x_3)\cdots(x-x_k)}{(x_2-x_1)(x_2-x_3)\cdots(x_2-x_k)}\\
&amp;\qquad 
+y_3 \frac{(x-x_1)(x-x_2)(x-x_4)\cdots(x-x_k)}{(x_3-x_1)(x_3-x_2)(x_3-x_4)\cdots(x_3-x_k)}
+\cdots\\
&amp;\qquad 
+y_k \frac{(x-x_1)(x-x_2)\cdots(x-x_{k-1})}{(x_k-x_1)(x_k-x_2)\cdots(x_k-x_{k-1})}
\end{align}$$
Notice that each $y_i$ is being multiplied by the polynomial
$$g_i(x) = \frac{1}{c_i}(x-x_1)(x-x_2)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_k)$$
which has as factors all terms of the form $(x-x_j)$ except $(x-x_i)$, and
the constant $c_i$ is given by
$$c_i = (x_i-x_1)(x_i-x_2)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_k).
\tag{1}$$
In other words, $g_i(x_i)=1$ and $g_i(x_j)=0$ for all $j \neq i$.

So there is a lot of calculation required, and care must be exercised
when writing the computer program that carries out these calculations.
But note that we really don't need to compute the entire polynomial 
$S(z)$  since all we want to find is the secret $s_0 = S(0)$. This
saves us a fair amount of work because we have
$$\begin{align}
S(0) &amp;= \sum_{i=1}^k y_i \cdot g_i(0)\\
&amp;=\sum_{i=1}^k y_i \prod_{j=1, j\neq i}^k\frac{-x_j}{x_i-x_j}\\
&amp;= y_1 \frac{(-x_2)(-x_3)\cdots(-x_k)}{(x_1-x_2)(x_1-x_3)\cdots(x_1-x_k)}
+ y_2 \frac{(-x_1)(-x_3)\cdots(-x_k)}{(x_2-x_1)(x_2-x_3)\cdots(x_2-x_k)}\\
&amp;\qquad 
+y_3 \frac{(-x_1)(x-x_2)(-x_4)\cdots(-x_k)}{(x_3-x_1)(x_3-x_2)(x_3-x_4)\cdots(x_3-x_k)}
+\cdots\\
&amp;\qquad 
+y_k \frac{(-x_1)(-x_2)\cdots(-x_{k-1})}{(x_k-x_1)(x_k-x_2)\cdots(x_k-x_{k-1})}\\
&amp;=  \sum_{i=1}^k y_i \cdot (-1)^k\frac{x_1x_2x_3\cdots x_k}{x_i\cdot c_i}\\
s_0 &amp;= (-1)^k (x_1x_2x_3\cdots x_k) \sum_{i=1}^k \frac{y_i}{x_i\cdot c_i}
\end{align}$$
So now the task seems a little easier. We need to
compute $x_1x_2x_3\cdots x_k$ only once using $k-1$ multiplications,
each $c_i$ has to be calculated separately and needs $k-1$ subtractions
and $k-2$ multiplications. We multiply $c_i$ by $x_i$ and then
do the division $y_i/(x_i\cdot c_i)$.  The hard part in all this is thus writing
the code for computing each $c_i$.  And of course we must remember that
all this is using finite-field arithmetic so that we must use operations
or subroutine calls as described in the answer by Ilmari Karonen.