given points (x0,y0)  (x1,y1)  (x2,y2)  (x3,y3)

    where x0 < x1 < x2 < x3

equations for computing Lagrange polynomial function L(x):

    y0d0 = y0 / ((x0 - x1) * (x0 - x2) * (x0 - x3))
    y1d1 = y1 / ((x1 - x0) * (x1 - x2) * (x1 - x3))
    y2d2 = y2 / ((x2 - x0) * (x2 - x1) * (x2 - x3))
    y3d3 = y3 / ((x3 - x0) * (x3 - x1) * (x3 - x2))

    a3 = y0d0 + y1d1 + y2d2 + y3d3

    a2 = - y0d0*(x1+x2+x3) - y1d1*(x0+x2+x3) -
           y2d2*(x0+x1+x3) - y3d3*(x0+x1+x2)

    a1 = y0d0*(x1*x2+x1*x3+x2*x3) + y1d1*(x0*x2+x0*x3+x2*x3) +
         y2d2*(x0*x1+x0*x3+x1*x3) + y3d3*(x0*x1+x0*x2+x1*x2)

    a0 = - y0d0*(x1*x2*x3) - y1d1*(x0*x2*x3) -
           y2d2*(x0*x1*x3) - y3d3*(x0*x1*x2)

    L(x) = (a3 * x^3) + (a2 * x^2) + (a1 * x) + a0

equations for computing Bezier control points:

    L'(x) = (3 * a3 * x^2) + (2 * a2 * x) + a1

    cx1 = x0 + (1/3 * (x3 - x0))
    cx2 = x3 - (1/3 * (x3 – x0))

    m1 = L'(x0)
    b1 = y0 - (m1 * x0)
    cy1 = (m1 * cx1) + b1

    m2 = L'(x3)
    b2 = y3 - (m2 * x3)
    cy2 = (m2 * cx2) + b2

(x0,y0)  (cx1,cy1)  (cx2,cy2)  (x3,y3)