BGW Multiplication: How does it work?

I could not understand how multiplication is carried out in the BGW protocol for 3 parties. Reference: https://crypto.stanford.edu/pbc/notes/crypto/bgw.html

I've understood that in multiplication, we calculate r(x) = p(x)q(x) but this means that r(x) is a polynomial of degree 2t. How is this polynomial redistributed in the circuit to get an O(n^2) complexity in multiplication over BGW protocol?

Note that p(x) and q(x) are both polynomials of degree t and as I am taking the case of semi-honest adversaries, t < n/2.

We can use the notation of the source document, except I will say we are trying to compute $$ab = c$$ (for $$a,b,c\in\mathbb{F}_p$$), so I can use $$x$$ as a variable when discussing polynomials.

Let $$a_1,\dots, a_n$$ and $$b_1,\dots, b_n$$ be $$t$$-out-of-$$n$$ secret sharings of $$a, b$$. Recall this means that they're computed by taking random $$t$$-degree polynomials: $$A(x) = \sum_{i = 0}^t \alpha_i x^i,\quad B(x) = \sum_{i = 0}^t\beta_i x^i$$ Where $$A(0) = a, B(0) = b$$ (so random subject to this condition). The shares are then created by evaluating the polynomials on the points $$\{1,2,\dots,n\}$$. In particular: $$\forall i \in\{1,\dots,n\} : a_i = A(i),\quad b_i = B(i)$$ Now we want to compute the product. The polynomial $$C(x) = A(x)B(x)$$ has the right constant term (as $$C(0) = A(0)B(0) = ab$$), but is too high of degree (as you mention). Moreover the "shares" of $$C(i)$$ can be computed locally, as $$C(i) = A(i)B(i) = a_ib_i$$.

But $$C(x)$$ is degree $$2t$$ (as you mention). We want to find some other polynomial $$\mathcal{C}(x)$$ such that:

1. $$\mathcal{C}(0) = C(0) = ab$$
2. $$\mathcal{C}(x)$$ is degree $$t$$
3. Computing $$\mathcal{C}(i)$$ isn't "too expensive" (in communication) if you already know $$A(i)$$ and $$B(i)$$

What will we do then? The idea is to utilize that there are two different ways of representing degree $$t$$ polynomials:

1. Through their $$t+1$$ coefficients (this is the "obvious") way
2. Through the evaulation of them on (at least) $$t+1$$ distinct points

Either of these is enough information to uniquely reconstruct the polynomial. The surprising thing is that you can convert from one to the other using a linear operation.

To see how we might establish this, recall that for an $$n\times n$$ matrix $$D$$ and vector $$\vec{v} = (v_1,\dots, v_n)$$, we have that: $$(D\vec{v})_i = \sum_{k = 0}^{n-1}D_{i, k} v_k$$

Note that this is similar to the expression: $$A(x) = \sum_{i = 0}^t \alpha_i x^i$$ If we fix evaluation points $$\{1,\dots, n\}$$, then we have in fact that: $$A(i) = \sum_{k = 0}^t \alpha_k i^k$$ This suggests setting $$D_{i, k} = i^k$$ and $$v_k = \alpha_k$$. This is precisely what we'll do, by defining the Vandermonde matrix (with respect to the aforementioned evaluation points): $$V = \begin{pmatrix} 1^0 & 1^1 & \dots & 1^{n-1}\\ 2^0 & 2^1 & \dots & 2^{n-1} \\ \vdots && \ddots & \vdots\\ n^0 & n^1 & \dots & n^{n-1} \end{pmatrix}$$ Note that: $$V\begin{pmatrix}\alpha_0\\\vdots\\\alpha_{n-1}\end{pmatrix} =\begin{pmatrix} \sum_{k = 0}^{n-1} \alpha_i 1^i\\ \sum_{k = 0}^{n-1} \alpha_i 2^i\\ \vdots\\ \sum_{k = 0}^{n-1} \alpha_i n^i\\ \end{pmatrix} = \begin{pmatrix}A(1)\\ A(2)\\ \vdots\\ A(n) \end{pmatrix} = \begin{pmatrix}a_1\\ a_2\\ \vdots\\ a_n \end{pmatrix}$$ So the Vandermonde Matrix precisely maps the "coefficient representation" of a polynomial to its "evaluation representation". This ends up being extremely closely related to the Fourier transform. The Discrete Fourier Transform can be written as a Vandermonde Matrix, and the Fast Fourier Transform can be explained that it's Vandermonde Matrix is Toeplitz (and in fact circulant), so admits especially efficient representations and matrix/vector multiplications, but this is an ahistorical aside.

So, we have an (invertible) matrix that maps an $$n$$-dimensional "coefficient vector" representation of a polynomial to an $$n$$-dimensional "evaluation vector" representation of a polynomial. For the moment, don't worry about how all of the people move around all the shares --- just make sure you understand how to do the computation.

We start with the "evaluation vector" representation $$C(i) = A(i)B(i)$$ for $$i\in\{1,\dots,n\}$$, which we can write as $$\vec c = (c_1,\dots, c_n)$$. We convert this to the "coefficient vector" representation via $$V^{-1}\vec{c}$$. This gives the coefficients of the polynomial $$C(x) = A(x)B(x)$$ as a vector. While there are $$n$$ coefficients, as discussed before this polynomial (uniquely determined from $$\vec{c}$$) is of degree $$2t$$, so the higher-order coefficients are 0.

We can convert this to a degree $$t$$ polynomial via truncation. Let: $$P = \begin{pmatrix} I_{t+1} & 0\\ 0 & 0\end{pmatrix}$$ Be an $$n\times n$$ block matrix, where $$I_{t+1}$$ is the $$(t+1)\times (t+1)$$ identity matrix. Then $$PV^{-1}\vec{c}$$ will "drop" higher order coefficients, leaving a degree $$t$$ polynomial. Importantly, it doesn't touch the constant term (so $$\mathcal{C}(0) = C(0) = ab$$ is preserved).

All that's left is to convert back to shares, so to convert from the "coefficient representation" to the "evaluation representation", again by using $$V$$. Thus $$VPV^{-1}\vec{c}$$ will output the shares (of a degree $$t$$ polynomial) that you want. Moreover, $$VPV^{-1}$$ can be precomputed by all participants in the protocol (it's just some $$n\times n$$ matrix. I could probably even write it out here, but won't).

This reduces the problem of multiplying shares to the problem of "computing a linear equation" of shares, which your source also discusses (at this link). As this is getting long I'll leave the answer here, but if you don't understand the linear case I encourage you to ask a new question about it.