fastsquad.lua

local ln   = math.log
local abs  = math.abs
local sin  = math.sin
local cos  = math.cos
local exp  = math.exp
local acos = math.acos
local sqrt = math.sqrt

local function MulLn(w0, x0, y0, z0, w1, x1, y1, z1)
	local w, x, y, z =
		w0*w1 - x0*x1 - y0*y1 - z0*z1,
		x0*w1 + w0*x1 - z0*y1 + y0*z1,
		y0*w1 + z0*x1 + w0*y1 - x0*z1,
		z0*w1 - y0*x1 + x0*y1 + w0*z1
	local normv = x*x + y*y + z*z
	local absq = sqrt(w*w + normv)
	local t = normv < 1e-15 and 0 or acos(abs(w)/absq)/sqrt(normv)
	return ln(absq), x*t, y*t, z*t
end

local function CMap(w0, x0, y0, z0, w1, x1, y1, z1, w2, x2, y2, z2)
	local norm1 = w1*w1 + x1*x1 + y1*y1 + z1*z1
	local iqw, iqx, iqy, iqz = w1/norm1, x1/-norm1, y1/-norm1, z1/-norm1
	local bw0, bx0, by0, bz0 = MulLn(w2, x2, y2, z2, iqw, iqx, iqy, iqz)
	local bw1, bx1, by1, bz1 = MulLn(w0, x0, y0, z0, iqw, iqx, iqy, iqz)
	local mx, my, mz =
		(bx0 + bx1)/-4,
		(by0 + by1)/-4,
		(bz0 + bz1)/-4
	local normv = sqrt(mx*mx + my*my + mz*mz)
	local u = exp((bw0 + bw1)/-4)
	local v = normv > 1e-15 and u*sin(normv)/normv or 0
	local ew, ex, ey, ez = u*cos(normv), v*mx, v*my, v*mz
	return
		ew*w1 - ex*x1 - ey*y1 - ez*z1,
		ex*w1 + ew*x1 - ez*y1 + ey*z1,
		ey*w1 + ez*x1 + ew*y1 - ex*z1,
		ez*w1 - ey*x1 + ex*y1 + ew*z1
end

local function NSlerpParams(p)
	return
		(1.0904 + p*(p*(3.55645 - p*1.43519) - 3.2452)),
		0.848013 + p*(p*0.215638 - 1.06021)
end

local function Squad(q0, q1, q2, q3)
	local q1w, q1x, q1y, q1z = q1[1], q1[2], q1[3], q1[4]
	local q2w, q2x, q2y, q2z = q2[1], q2[2], q2[3], q2[4]

	local p0w, p0x, p0y, p0z = CMap(
		q0[1], q0[2], q0[3], q0[4],
		q1w, q1x, q1y, q1z,
		q2w, q2x, q2y, q2z
	)
	local p1w, p1x, p1y, p1z = CMap(
		q1w, q1x, q1y, q1z,
		q2w, q2x, q2y, q2z,
		q3[1], q3[2], q3[3], q3[4]
	)

	local d0, d1 =
		q1w*q2w + q1x*q2x + q1y*q2y + q1z*q2z,
		p0w*p1w + p0x*p1x + p0y*p1y + p0z*p1z

	local u0, v0 = NSlerpParams(abs(d0))
	local u1, v1 = NSlerpParams(abs(d1))

	d0, d1 = d0 < 0, d1 < 0

	return function(t)
		local awx, awm =
			t*(t - 0.5)*(t - 1),
			(t - 0.5)*(t - 0.5)

		local o = t + awx*(u0*awm + v0)
		local t0, t1 = 1 - o, d0 and -o or o
		local w, x, y, z =
			q1w*t0 + q2w*t1,
			q1x*t0 + q2x*t1,
			q1y*t0 + q2y*t1,
			q1z*t0 + q2z*t1
		local n = sqrt(w*w + x*x + y*y + z*z)
		local c0w, c0x, c0y, c0z = w/n, x/n, y/n, z/n

		o = t + awx*(u1*awm + v1)
		t0, t1 = 1 - o, d1 and -o or o
		w, x, y, z =
			p0w*t0 + p1w*t1,
			p0x*t0 + p1x*t1,
			p0y*t0 + p1y*t1,
			p0z*t0 + p1z*t1
		n = sqrt(w*w + x*x + y*y + z*z)
		local c1w, c1x, c1y, c1z = w/n, x/n, y/n, z/n

		local d2 = c0w*c1w + c0x*c1x + c0y*c1y + c0z*c1z
		local d2a = abs(d2) 
		local h = 2*t*(1 - t)
		o = h + h*(h - 0.5)*(h - 1)*(
			(1.0904 + d2a*(d2a*(3.55645 - d2a*1.43519) - 3.2452))*
			(h - 0.5)*(h - 0.5) + 0.848013 + d2a*(d2a*0.215638 - 1.06021)
		)
		t0, t1 = 1 - o, d2 < 0 and -o or o
		w, x, y, z =
			c0w*t0 + c1w*t1,
			c0x*t0 + c1x*t1,
			c0y*t0 + c1y*t1,
			c0z*t0 + c1z*t1
		n = sqrt(w*w + x*x + y*y + z*z)
		return w/n, x/n, y/n, z/n
	end
end

localized sqrt call is about the same speed as a^b on Lua 5.1.5 and substantially faster than a^b on LuaJIT