def computeC2Spline(nodes):
    '''Compute a cubic C2-spline by the list of 3 nodes.
    The spline consists of 2 cubic polynomials p0(x) and p1(x).
    p0(x) = p00 + p01*x + p02*x^2 + p03*x^3
    p1(x) = p10 + p11*x + p12*x^2 + p13*x^3
    The function returns the list of quadruples:
    [(p00, p01, p02, p03), (p10, p11, p12, p13)]
    The coefficients of polynomials are computed as the solution
    of the linear system. Let (x0, y0), (x1, y1), (x2, y2) are
    the nodes of interpolation. Then the system of linear equations is

    p0(x0) = y0
    p0''(x0) = 0

    p0(x1) = y1
    p1(x1) = y1
    p0'(x1) = p1'(x1)
    p0''(x1) = p1''(x1)

    p1(x2) = y2
    p1''(x2) = 0
    '''

    assert(len(nodes) == 3)
    n = len(nodes)
    a = []      # The matrix is to be composed as a list of rows
    matrixLine = [0.]*(4*(n-1))         # The matrix row
    b = [0.]*(4*(n-1))
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    x0 = nodes[0][0]; y0 = nodes[0][1]
    # Equation 0: 
    #     p0(x0) = y0
    matrixLine[0] = 1
    matrixLine[1] = x0
    matrixLine[2] = x0^2
    matrixLine[3] = x0^3
    a.append(matrixLine.copy())
    b[0] = y0
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    # Equation 1:
    #    p0''(x0) = 0
    #    2*p02 + 6*p03*x0 = 0
    matrixLine[2] = 2
    matrixLine[3] = 6*x0
    a.append(matrixLine.copy())
    b[1] = 0
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    x1 = nodes[1][0]; y1 = nodes[1][1]
    # Equation 2:
    #     p0(x1) = y1
    matrixLine[0] = 1
    matrixLine[1] = x1
    matrixLine[2] = x1^2
    matrixLine[3] = x1^3
    a.append(matrixLine.copy())
    b[2] = y1
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    # Equation 3:
    #     p1(x1) = y1
    matrixLine[4] = 1
    matrixLine[5] = x1
    matrixLine[6] = x1^2
    matrixLine[7] = x1^3
    a.append(matrixLine.copy())
    b[3] = y1
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    # Equation 4:
    #     p0'(x1) - p1'(x1) = 0
    matrixLine[1] = 1
    matrixLine[2] = 2*x1
    matrixLine[3] = 3*x1^2
    matrixLine[5] = -1
    matrixLine[6] = -2*x1
    matrixLine[7] = -3*x1^2
    a.append(matrixLine.copy())
    b[4] = 0
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    # Equation 5:
    #     p0''(x1) - p1''(x1) = 0
    matrixLine[2] = 2
    matrixLine[3] = 6*x1
    matrixLine[6] = -2
    matrixLine[7] = -6*x1
    a.append(matrixLine.copy())
    b[5] = 0
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    x2 = nodes[2][0]; y2 = nodes[2][1]
    # Equation 6:
    #     p1(x2) = y2
    matrixLine[4] = 1
    matrixLine[5] = x2
    matrixLine[6] = x2^2
    matrixLine[7] = x2^3
    a.append(matrixLine.copy())
    b[6] = y2
    
    for i in range(len(matrixLine)):    # Clear the matrix row
        matrixLine[i] = 0

    # Equation 7:
    #     p1''(x2) = 0
    matrixLine[6] = 2
    matrixLine[7] = 6*x2
    a.append(matrixLine.copy())
    b[7] = 0

    m = matrix(a)
    bb = vector(b)
    print("Matrix of the system:")
    print(m)
    print("Column of free terms:", bb)
    
    p = m.solve_right(vector(b))
    print("Solution of the system:", p)
    s = [p[0:4], p[4:8]]

    print("Resulting spline:", s)
    return s

def splineValue(nodes, spline, x):
    '''The value of the spline in the point x'''
    i = 0
    n = len(nodes)
    # x_i <= x < x_i+1
    while i < n-1 and x >= nodes[i+1][0]:
        i += 1
    assert(
        (i == 0 and x < nodes[1]) or
        (i == n-1 and x >= nodes[n-1]) or
        nodes[i][0] <= x < nodes[i+1][0]
    )  
    p = spline[i]
    return p[0] + p[1]*x + p[2]*x^2 + p[3]*x^3

def splineC2(nodes):
    '''Plot the C2-spline constructed by nodes.
    This is the restricted version that works only in the case
    len(nodes) == 3'''

    assert(len(nodes) == 3)

    s = computeC2Spline(nodes)
    var("x")
    p0 = s[0][0] + s[0][1]*x + s[0][2]*x^2 + s[0][3]*x^3    
    shape = plot(p0, (x, nodes[0][0], nodes[1][0]))

    p1 = s[1][0] + s[1][1]*x + s[1][2]*x^2 + s[1][3]*x^3
    shape += plot(p1, (x, nodes[1][0], nodes[2][0]))

    # Draw nodes
    for node in nodes:
        shape += point(node, size=40, color="red")

    return shape