def spline(points, dx =  1 ):
    x = map(lambda p: float(p[ 0 ]), points)
    y = map(lambda p: float(p[ 1 ]), points)
    N = len(points) -  1 
    h = {}
    l = {}
    delta = {}
    _lambda = {}
    c = {}
    c[ 0 ] =  0 
    c[N] =  0 
    d = {}
    b = {}
    r = []
    for i in range( 1 , N+ 1 ):
        h[i] = x[i] - x[i- 1 ]
        l[i] = (y[i] - y[i- 1 ]) / h[i]
    delta[ 1 ] = - h[ 2 ]/( 2  * (h[ 1 ] + h[ 2 ])) 
    _lambda[ 1 ] =  1 . 5  * (l[ 2 ] - l[ 1 ]) / (h[ 1 ] + h[ 2 ])
    for i in range( 3 , N+ 1 ):
        delta[i- 1 ] = - h[i] / ( 2  *h[i- 1 ] +  2  * h[i] + h[i- 1 ] * delta[i- 2 ])
        _lambda[i- 1 ] = ( 3  * l[i] -  3  * l[i- 1 ] - h[i- 1 ] * _lambda[i- 2 ]) / \
                       ( 2  * h[i- 1 ] +  2  * h[i] + h[i- 1 ] * delta[i- 2 ])
    for i in range(N,  1 , - 1 ): c[i- 1 ] = delta[i- 1 ] * c[i] + _lambda[i- 1 ]
    for i in range( 1 , N+ 1 ):
        d[i] = (c[i] - c[i- 1 ]) / ( 3  * h[i])
        b[i] = l[i] + ( 2  * c[i] * h[i] + h[i] * c[i- 1 ]) /  3 

    i = x[ 0 ]
    while (i <= x[- 1 ]):
        for j in range( 1 , N+ 1 ):
            if (i >= x[j- 1 ]) & (i <= x[j]): break
        f = y[j] + b[j] * (i - x[j]) + c[j] * ((i - x[j]) **  2 ) + d[j] * ((i - x[j]) **  3 )
        print f
        r.append([i, f])        
        i += dx
    return r