Quintic spline interpolation

quint.jl

type quint{T}
    dx::T
    x0::T
    Y::Array{T,1}
    B::Array{T,1}
    C::Array{T,1}
    D::Array{T,1}
    E::Array{T,1}
    F::Array{T,1}
    function quint(X::Array{T,1})
        Y = similar(X)
        B = similar(X)
        C = similar(X)
        D = similar(X)
        E = similar(X)
        F = similar(X)
        new(X[2]-X[1], X[1], Y, B, C, D, E, F)
    end
end

quint{T}(X::Array{T,1}) = quint{T}(X)

function qeval{T}(q::quint{T}, x::T)
    P = (x - q.x0)/q.dx
    I = int(P)
    P -= I
    I += 1
    @inbounds return ((((q.F[I]*P+q.E[I])*P+q.D[I])*P+q.C[I])*P+q.B[I])*P+q.Y[I]
end

function qeval{T}(q::quint{T}, X::Array{T,1}, Y::Array{T,1})
    N = length(X)
    for i=1:N
        @inbounds Y[i] = qeval(q, X[i])
    end
end

function qset{T}(q::quint{T}, Y::Array{T,1})
    N = length(Y)
    q.Y[:] = Y
    M = N - 3
    P = 0.0
    Q = 0.0
    R = 0.0
    S = 0.0
    TT = 0.0
    q.D[M+1] = 0.0
    q.D[M+2] = 0.0
    for I=1:M
        @inbounds begin
            U = P*R
            q.B[I] = 1.0/(66.0-U*R-Q)
            R = 26.0 - U
            q.C[I] = R
            q.D[I] = q.Y[I+3] - 3.0*(q.Y[I+2]-q.Y[I+1]) - q.Y[I] - U*S - Q*TT
            Q = P
            P = q.B[I]
            TT = S
            S = q.D[I]
        end
    end
    I = N - 2
    for M=4:N
        I = I - 1
        @inbounds q.D[I] = (q.D[I]-q.D[I+2]-q.D[I+1]*q.C[I])*q.B[I]
    end
    M = N - 1
    Q = 0.0
    R = q.D[1]
    TT = R
    V = R
    for I=2:M
        @inbounds begin
            P = Q
            Q = R
            R = q.D[I]
            S = TT
            TT = P - Q - Q + R
            q.F[I] = TT
            U = 5.0*(-P+Q)
            q.E[I] = U
            q.D[I] = 10.0*(P+Q)
            q.B[I] = 0.5*(q.Y[I+1]-q.Y[I-1]-S-TT) - q.D[I]
            q.C[I] = 0.5*(q.Y[I+1]+q.Y[I-1]+S-TT) - q.Y[I] - U
        end
    end 
    @inbounds begin
        q.F[1] = V
        q.E[1] = 0.0
        q.D[1] = 0.0
        q.C[1] = q.C[2] - 10.0*V
        q.B[1] = q.Y[2] - q.Y[1] - q.C[1] - V
        q.E[N] = 0.0
        q.D[N] = 0.0
        q.C[N] = q.C[N-1] + 10.0*TT
        q.B[N] = q.Y[N] - q.Y[N-1] + q.C[N] - TT
    end
end

tquint.jl

include("quint.jl")
x = linspace(-10.0,10.0,200000)
x2 = x[1:end-1] .+ (x[2]-x[1])/2.0
y = sin(x)
y2 = similar(x2)
q = quint(x)
qset(q, y)
@time qset(q, y)
@time qset(q, y)
@time qset(q, y)
qeval(q, x2, y2)
@time qeval(q, x2, y2)
@time qeval(q, x2, y2)
@time qeval(q, x2, y2)

https://github.com/nagynet200/Archimedes.wiki.git