Find the acceleration-minimizing curve between a list of points

NaturalSplines.lua

-- Given a knot sequence p[1<=i<=n], solve for a sequence of cubic
-- splines s[1<=i<=n-1] such that the square of acceleration is minimized.

local gamma = {
	0.50000000000000000, 0.28571428571428571, 0.26923076923076923, 0.26804123711340206,
	0.26795580110497238, 0.26794966691339748, 0.26794922649742166, 0.26794919487697295,
	0.26794919260672685, 0.26794919244373052, 0.26794919243202791, 0.26794919243118770,
	0.26794919243112737, 0.26794919243112304, 0.26794919243112273, 0.26794919243112271,
}
local gammaLimit = gamma[#gamma]


local function naturalSpline(p)
	-- Solve for d[1;;m]
	-- {2 1          } {d[1]  }    {p[2]-p[1]    }
	-- {1 4 1        } {d[2]  }    {p[3]-p[1]    }
	-- {  1 4 1      } {d[3]  }    {p[4]-p[2]    }
	-- {    . . .    } {. . . } = 3{. . .        }
	-- {      1 4 1  } {d[m-2]}    {p[m-1]-p[m-3]}
	-- {        1 4 1} {d[m-1]}    {p[m]-p[m-2]  }
	-- {          1 2} {d[m]  }    {p[m]-p[m-1]  }

	local m = #p - 1

	-- Eliminate the lower diagonal by iterating up degrees of the
	-- 1st Chebyshev polynomial T(1 <= i <= m, 2)
	-- This sequence converges after 16 iterations in double precision.
	
	--[[ Gamma LUT generator
	if n_gamma < m then
		local prev = gamma[n_gamma]
		for i = n_gamma + 1, m do
			prev = 1/(4 - prev)
			gamma[i] = prev
		end
		n_gamma = m
	end
	--]]

	-- Forward eliminate the input matrix
	local prev = 1.5*(p[2] - p[1])
	local delta = {prev}
	for i = 2, m do
		prev = (3*(p[i + 1] - p[i - 1]) - prev)*(gamma[i] or gammaLimit)
		delta[i] = prev
	end

	-- Solve for knot derivatives by back subsituting while feeding
	-- known values into the Hermite solver
	local out = {}
	local p1 = p[m + 1]
	local d1 = (3*(p1 - p[m]) - prev)/(2 - (gamma[m] or gammaLimit))
	for i = m, 1, -1 do
		-- p0, p1: knot positions
		-- d0, d1: knot derivatives
		local p0 = p[i]
		local d0 = delta[i] - (gamma[i] or gammaLimit)*d1
		out[i] = (Hermite :: any)(p0, p1, d0, d1)
		p1 = p0
		d1 = d0
	end

	return out
end


return naturalSpline