This is an excerpt from my math models textbook. It’s about Lagrange Polynomials which is a technique that lets you fit a polynomial to a set of any number of unique points (x_1,y_1) … (x_n,y_n) so long as all your x-values are different (otherwise it wouldn’t be a function, and couldn’t be a polynomial). The polynomial you’ll calculate will be the unique, lowest degree polynomial that passes through all points.
It’s not too bad, once you consider that everything in each term is a constant value, except for “x” itself.
So the numerator of each term is the product of three linear factors, like (x-4)(x-2)(x-6), which should produce a cubic, like x³ - 12x² + 44x - 48. Then the denominator of each term is a pure constant, so it would be like taking that cubic and dividing it by 4, getting 0.25x³ - 3x² + 11x - 12. Then the yₙ terms are also constants, so no different than doing something like multiplying by 2, getting you something like 0.5x³ - 6x² + 22x - 24, if I take that example a bit too far. And at that point it’s just the sum of four cubics, which will be cubic as long as the x³ terms don’t perfectly cancel out - which I believe would only happen if the four pairs of points used to make the function were all perfectly laying on the same line or parabola.
The construction’s also pretty clever: OP said the point was to fit the function to the four points (x₁, y₁), (x₂, y₂), (x₃, y₃), (x₄, y₄). Let’s say we set x = x₂, then. Because (x-x₂) appears in the numerator of every term but the second term, every term but the second term will have a 0 in its numerator and cancel out - so we only need to consider the second term. Its numerator is then (x₂-x₁)(x₂-x₃)(x₂-x₄) - exactly the same as its denominator. So they both cancel out, leaving only y₂ - meaning we get P₃(x₂) = y₂, as desired.