As far as I can see, you are computing Lagrange polynomials.
In the specific case of 3 data points (x_0, y_0), (x_1, y_1), (x_2, y_2) - which in your example are (0, 4), (1, 2), (3, 3), the calculation is quite easy.
f(x) = y_0*l_0(x) = y_0/((x_0-x_1)*(x_0-x_2))*(x^2 + -(x_1+x_2)*x + (x_1*x_2))
The other two polynomials can be computed similarly.
In their sum, you just have to group together the corresponding coefficients, and make the modular arithmetic. (Division can be made with the multiplication of the inverse element, and the inverse element can be easily computed with the help of Fermat's little theorem as a^(p-2) in case of prime modulus.)