Initial commit of the actual code.

LICENSE.md

@@ -1,6 +1,12 @@
-The MIT License (MIT)
+Licenses
+======================
+CPML is Copyright (c) 2014 Colby Klein <shakesoda@gmail.com>.
+Portions of mat4.lua are from LuaMatrix, (c) 2010 Michael Lutz. MIT.
+Code in simplex.lua is (c) 2011 Stefan Gustavson. MIT.
+Code in quat.lua is from Andrew Stacey and covered under the CC0 license.
 
-Copyright (c) 2014 Colby Klein <shakesoda@gmail.com>
+The MIT License (MIT)
+======================
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
@@ -19,4 +25,3 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.
-

cpml/constants.lua

@@ -0,0 +1,6 @@
+local constants = {}
+
+-- same as C's FLT_EPSILON
+constants.FLT_EPSILON = 1.19209290e-07
+
+return constants

cpml/intersect.lua

@@ -0,0 +1,89 @@
+local current_folder = (...):gsub('%.[^%.]+$', '')
+local vec3 = require(current_folder .. "vec3")
+local constants = require(current_folder .. "constants")
+
+-- http://www.lighthouse3d.com/tutorials/maths/ray-triangle-intersection/
+function cpml.ray_triangle(ray, triangle)
+	assert(ray.point ~= nil)
+	assert(ray.direction ~= nil)
+	assert(#triangle == 3)
+
+	local p, d = ray.point, ray.direction
+
+	local h, s, q = vec3(), vec3(), vec3()
+	local a, f, u, v
+
+	local e1 = triangle[2] - triangle[1]
+	local e2 = triangle[3] - triangle[1]
+
+	h = d:clone():cross(e2)
+
+	a = (e1*h) -- dot product
+
+	if a > -0.00001 and a < 0.00001 then
+		return false
+	end
+
+	f = 1/a
+	s = p - triangle[1]
+	u = f * (s*h)
+
+	if u < 0 or u > 1 then
+		return false
+	end
+
+	q = s:clone():cross(e1)
+	v = f * (d*q)
+
+	if v < 0 or u + v > 1 then
+		return false
+	end
+
+	-- at this stage we can compute t to find out where
+	-- the intersection point is on the line
+	t = f * (e2*q)
+
+	if t > constants.FLT_EPSILON then
+		return p + t * d -- we've got a hit!
+	else
+		return false -- the line intersects, but it's behind the point
+	end
+end
+
+-- Algorithm is ported from the C algorithm of 
+-- Paul Bourke at http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline3d/
+-- Archive.org am hero \o/
+function cpml.line_line(p1, p2, p3, p4)
+	local epsilon = constants.FLT_EPSILON
+	local resultSegmentPoint1 = vec3(0,0,0)
+	local resultSegmentPoint2 = vec3(0,0,0)
+
+	local p13 = p1 - p3
+	local p43 = p4 - p3
+	local p21 = p2 - p1
+
+	if p43:len2() < epsilon then return false end
+	if p21:len2() < epsilon then return false end
+
+	local d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z
+	local d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z
+	local d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z
+	local d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z
+	local d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z
+
+	local denom = d2121 * d4343 - d4321 * d4321
+	if math.abs(denom) < epsilon then return false end
+	local numer = d1343 * d4321 - d1321 * d4343
+
+	local mua = numer / denom
+	local mub = (d1343 + d4321 * (mua)) / d4343
+
+	resultSegmentPoint1.x = p1.x + mua * p21.x
+	resultSegmentPoint1.y = p1.y + mua * p21.y
+	resultSegmentPoint1.z = p1.z + mua * p21.z
+	resultSegmentPoint2.x = p3.x + mub * p43.x
+	resultSegmentPoint2.y = p3.y + mub * p43.y
+	resultSegmentPoint2.z = p3.z + mub * p43.z
+
+	return true, resultSegmentPoint1, resultSegmentPoint2
+end

cpml/mat4.lua

@@ -0,0 +1,456 @@
+-- double 4x4, 1-based, column major
+-- local matrix = {}
+
+local current_folder = (...):gsub('%.[^%.]+$', '')
+local constants = require(current_folder .. "constants")
+
+local mat4 = {}
+mat4.__index = mat4
+setmetatable(mat4, mat4)
+
+-- from https://github.com/davidm/lua-matrix/blob/master/lua/matrix.lua
+-- Multiply two matrices; m1 columns must be equal to m2 rows
+-- e.g. #m1[1] == #m2
+local function matrix_mult_nxn(m1, m2)
+	local mtx = {}
+	for i = 1, #m1 do
+		mtx[i] = {}
+		for j = 1, #m2[1] do
+			local num = m1[i][1] * m2[1][j]
+			for n = 2, #m1[1] do
+				num = num + m1[i][n] * m2[n][j]
+			end
+			mtx[i][j] = num
+		end
+	end
+	return mtx
+end
+
+function mat4:__call(v)
+	local m = {
+		1, 0, 0, 0,
+		0, 1, 0, 0,
+		0, 0, 1, 0,
+		0, 0, 0, 1
+	}
+	if type(v) == "table" and #v == 16 then
+		for i=1,16 do
+			m[i] = v[i]
+		end
+	elseif type(v) == "table" and #v == 9 then
+		m[1], m[2], m[3] = v[1], v[2], v[3]
+		m[5], m[6], m[7] = v[4], v[5], v[6]
+		m[9], m[10], m[11] = v[7], v[8], v[9]
+		m[16] = 1
+	elseif type(v) == "table" and type(v[1]) == "table" then
+		local idx = 1
+		for i=1, 4 do
+			for j=1, 4 do
+				m[idx] = v[i][j]
+				idx = idx + 1
+			end
+		end
+	end
+
+	-- Look in mat4 for metamethods
+	setmetatable(m, mat4)
+
+	return m
+end
+
+function mat4:__eq(b)
+	local abs = math.abs
+	for i=1, 16 do
+		if abs(self[i] - b[i]) > constants.FLT_EPSILON then
+			return false
+		end
+	end
+	return true
+end
+
+function mat4:__tostring()
+	local str = "[ "
+	for i, v in ipairs(self) do
+		str = str .. string.format("%2.5f", v)
+		if i < #self then
+			str = str .. ", "
+		end
+	end
+	str = str .. " ]"
+	return str
+end
+
+function mat4:ortho(left, right, top, bottom, near, far)
+	local out = mat4()
+	out[1] = 2 / (right - left)
+	out[6] = 2 / (top - bottom)
+	out[11] = -2 / (far - near)
+	out[13] = -((right + left) / (right - left))
+	out[14] = -((top + bottom) / (top - bottom))
+	out[15] = -((far + near) / (far - near))
+	out[16] = 1
+	return out
+end
+
+function mat4:persp(fovy, aspect, near, far)
+	local out = mat4()
+	local f = math.tan(fovy / 2)
+	out[1] = 1 / (aspect * f)
+	out[6] = 1 / f
+	out[11] = -((far + near) / (far - near))
+	out[12] = -1
+	out[15] = -((2 * far * near) / (far - near))
+	return out
+end
+
+function mat4:translate(x, y, z)
+	if type(x) == "table" then
+		x, y, z = unpack(x)
+	end
+	local m = {
+		1, 0, 0, 0,
+		0, 1, 0, 0,
+		0, 0, 1, 0,
+		x, y, z, 1
+	}
+	return mat4(m) * mat4(self)
+end
+
+function mat4:scale(x, y, z)
+	if type(x) == "table" then
+		x, y, z = unpack(x)
+	end
+	local m = {
+		x, 0, 0, 0,
+		0, y, 0, 0,
+		0, 0, z, 0,
+		0, 0, 0, 1
+	}
+	return mat4(m) * mat4(self)
+end
+
+local function len(v)
+	return math.sqrt(v[1] * v[1] + v[2] * v[2] + v[3] * v[3])
+end
+
+function mat4:rotate(angle, axis)
+	if type(angle) == "table" then
+		angle, axis = angle:to_axis_angle()
+	end
+	local l = len(axis)
+	if l == 0 then
+		return self
+	end
+	local x, y, z = axis[1] / l, axis[2] / l, axis[3] / l
+	local c = math.cos(angle)
+	local s = math.sin(angle)
+	local m = {
+		x*x*(1-c)+c, y*x*(1-c)+z*s, x*z*(1-c)-y*s, 0,
+		x*y*(1-c)-z*s, y*y*(1-c)+c, y*z*(1-c)+x*s, 0,
+		x*z*(1-c)+y*s, y*z*(1-c)-x*s, z*z*(1-c)+c, 0,
+		0, 0, 0, 1,
+	}
+	return mat4(m) * mat4(self)
+end
+
+-- Set mat4 to identity mat4. Tested OK
+function mat4:identity()
+	local out = mat4()
+	for i=1, 16, 5 do
+		out[i] = 1
+	end
+	return out
+end
+
+function mat4:clone()
+	return mat4(self)
+end
+
+-- Inverse of matrix. Tested OK
+function mat4:invert()
+	local out = mat4()
+
+	out[1] =  self[6]  * self[11] * self[16] -
+		self[6]  * self[12] * self[15] -
+		self[10] * self[7]  * self[16] +
+		self[10] * self[8]  * self[15] +
+		self[14] * self[7]  * self[12] -
+		self[14] * self[8]  * self[11]
+
+	out[5] = -self[5]  * self[11] * self[16] +
+		self[5]  * self[12] * self[15] +
+		self[9]  * self[7]  * self[16] -
+		self[9]  * self[8]  * self[15] -
+		self[13] * self[7]  * self[12] +
+		self[13] * self[8]  * self[11]
+
+	out[9] =  self[5]  * self[10] * self[16] -
+		self[5]  * self[12] * self[14] -
+		self[9]  * self[6]  * self[16] +
+		self[9]  * self[8]  * self[14] +
+		self[13] * self[6]  * self[12] -
+		self[13] * self[8]  * self[10]
+
+	out[13] = -self[5]  * self[10] * self[15] +
+		self[5]  * self[11] * self[14] +
+		self[9]  * self[6]  * self[15] -
+		self[9]  * self[7]  * self[14] -
+		self[13] * self[6]  * self[11] +
+		self[13] * self[7]  * self[10]
+
+	out[2] = -self[2]  * self[11] * self[16] +
+		self[2]  * self[12] * self[15] +
+		self[10] * self[3]  * self[16] -
+		self[10] * self[4]  * self[15] -
+		self[14] * self[3]  * self[12] +
+		self[14] * self[4]  * self[11]
+
+	out[6] =  self[1]  * self[11] * self[16] -
+		self[1]  * self[12] * self[15] -
+		self[9]  * self[3] * self[16] +
+		self[9]  * self[4] * self[15] +
+		self[13] * self[3] * self[12] -
+		self[13] * self[4] * self[11]
+
+	out[10] = -self[1]  * self[10] * self[16] +
+		self[1]  * self[12] * self[14] +
+		self[9]  * self[2]  * self[16] -
+		self[9]  * self[4]  * self[14] -
+		self[13] * self[2]  * self[12] +
+		self[13] * self[4]  * self[10]
+
+	out[14] = self[1]  * self[10] * self[15] -
+		self[1]  * self[11] * self[14] -
+		self[9]  * self[2] * self[15] +
+		self[9]  * self[3] * self[14] +
+		self[13] * self[2] * self[11] -
+		self[13] * self[3] * self[10]
+
+	out[3] = self[2]  * self[7] * self[16] -
+		self[2]  * self[8] * self[15] -
+		self[6]  * self[3] * self[16] +
+		self[6]  * self[4] * self[15] +
+		self[14] * self[3] * self[8] -
+		self[14] * self[4] * self[7]
+
+	out[7] = -self[1]  * self[7] * self[16] +
+		self[1]  * self[8] * self[15] +
+		self[5]  * self[3] * self[16] -
+		self[5]  * self[4] * self[15] -
+		self[13] * self[3] * self[8] +
+		self[13] * self[4] * self[7]
+
+	out[11] = self[1]  * self[6] * self[16] -
+		self[1]  * self[8] * self[14] -
+		self[5]  * self[2] * self[16] +
+		self[5]  * self[4] * self[14] +
+		self[13] * self[2] * self[8] -
+		self[13] * self[4] * self[6]
+
+	out[15] = -self[1]  * self[6] * self[15] +
+		self[1]  * self[7] * self[14] +
+		self[5]  * self[2] * self[15] -
+		self[5]  * self[3] * self[14] -
+		self[13] * self[2] * self[7] +
+		self[13] * self[3] * self[6]
+
+	out[4] = -self[2]  * self[7] * self[12] +
+		self[2]  * self[8] * self[11] +
+		self[6]  * self[3] * self[12] -
+		self[6]  * self[4] * self[11] -
+		self[10] * self[3] * self[8] +
+		self[10] * self[4] * self[7]
+
+	out[8] = self[1] * self[7] * self[12] -
+		self[1] * self[8] * self[11] -
+		self[5] * self[3] * self[12] +
+		self[5] * self[4] * self[11] +
+		self[9] * self[3] * self[8] -
+		self[9] * self[4] * self[7]
+
+	out[12] = -self[1] * self[6] * self[12] +
+		self[1] * self[8] * self[10] +
+		self[5] * self[2] * self[12] -
+		self[5] * self[4] * self[10] -
+		self[9] * self[2] * self[8] +
+		self[9] * self[4] * self[6]
+
+	out[16] = self[1] * self[6] * self[11] -
+		self[1] * self[7] * self[10] -
+		self[5] * self[2] * self[11] +
+		self[5] * self[3] * self[10] +
+		self[9] * self[2] * self[7] -
+		self[9] * self[3] * self[6]
+
+	local det = self[1] * out[1] + self[2] * out[5] + self[3] * out[9] + self[4] * out[13]
+
+	if det == 0 then return self end
+
+	det = 1.0 / det
+
+	for i = 1, 16 do
+		out[i] = out[i] * det
+	end
+
+	return out
+end
+
+-- THREE broken projections! ARGH!
+function mat4.project(objx, objy, objz, model, view, projection, viewport)
+	local position = projection * view * model * { objx, objy, objz, 1 }
+
+	-- perspective division
+	position[4] = 1 / position[4]
+	position[1] = position[1] * position[4]
+	position[2] = position[2] * position[4]
+	position[3] = position[3] * position[4]
+
+	local out = {}
+	out[1] = ((position[1] + 1) * 0.5) * viewport.w + viewport.x
+	out[2] = ((position[2] + 1) * 0.5) * viewport.h + viewport.y
+	out[3] = ((position[3] + 1) * 0.5) * viewport.f + viewport.n
+	return out
+end
+
+function mat4.unproject(winx, winy, winz, modelview, projection, viewport)
+	local m = (projection * modelview):inverse()
+	if not m then
+		return false
+	end
+
+	-- Transformation of normalized coordinates between -1 and 1
+	local view = {}
+	view[1] = (winx - viewport[1]) / viewport[3] * 2 - 1
+	view[2] = (winy - viewport[2]) / viewport[4] * 2 - 1
+	view[3] = 2 * winz - 1
+	view[4] = 1
+
+	local out = m * view
+
+	if out[4] == 0 then
+		return false
+	end
+
+	out[4] = 1 / out[4]
+
+	objectCoordinate = {}
+	objectCoordinate[1] = out[1] * out[4]
+	objectCoordinate[2] = out[2] * out[4]
+	objectCoordinate[3] = out[3] * out[4]
+
+	return objectCoordinate
+end
+
+function mat4:look_at(eye, center, up)
+	local forward = (center - eye):normalize()
+	local side = forward:cross(up):normalize()
+	local new_up = side:cross(forward):normalize()
+
+	local view = mat4()
+	view[1]	= side.x
+	view[5]	= side.y
+	view[9]	= side.z
+
+	view[2]	= new_up.x
+	view[6]	= new_up.y
+	view[10]= new_up.z
+
+	view[3]	= -forward.x
+	view[7]	= -forward.y
+	view[11]= -forward.z
+
+	view[16]= 1
+
+	-- Fix 1u offset
+	local new_eye = eye + forward
+
+	return mat4():identity():translate(-new_eye.x, -new_eye.y, -new_eye.z) * view
+end
+
+
+function mat4:transpose()
+	local m = {
+		self[1], self[5], self[9], self[13],
+		self[2], self[6], self[10], self[14],
+		self[3], self[7], self[11], self[15],
+		self[4], self[8], self[12], self[16]
+	}
+	return mat4(m)
+end
+
+function mat4:__unm()
+	return self:invert()
+end
+
+-- Multiply mat4 by a mat4. Tested OK
+function mat4:__mul(m)
+	if #m == 4 then
+		-- should, hopefully, return a 4x1 matrix
+		-- i.e. a vec4
+		local tmp = matrix_mult_nxn(self:to_vec4s(), { {m[1]}, {m[2]}, {m[3]}, {m[4]} })
+		local v = {}
+		for i=1, 4 do
+			v[i] = tmp[i][1]
+		end
+		return v
+	end
+
+	-- local tmp = matrix_mult_nxn(self:to_vec4s(), m:to_vec4s())
+	-- return mat4(tmp)
+	local out = mat4()
+	for i=0, 12, 4 do
+		for j=1, 4 do
+			-- if type(m[j]) == "table" then
+			-- 	error("Whoa, this is broken!")
+			-- end
+			out[i+j] = m[j] * self[i+1] + m[j+4] * self[i+2] + m[j+8] * self[i+3] + m[j+12] * self[i+4]
+		end
+	end
+	return out
+end
+
+-- Multiply mat4 by a mat4. Tested OK
+function mat4:debug_mul(m)
+	if #m == 4 then
+		-- should, hopefully, return a 4x1 matrix
+		-- i.e. a vec4
+		local tmp = matrix_mult_nxn(self:to_vec4s(), { {m[1]}, {m[2]}, {m[3]}, {m[4]} })
+		local v = {}
+		for i=1, 4 do
+			v[i] = tmp[i][1]
+		end
+		return v
+	end
+
+	-- local tmp = matrix_mult_nxn(self:to_vec4s(), m:to_vec4s())
+	-- return mat4(tmp)
+
+	console.i("beginning mult")
+	console.i("input a: " .. tostring(self))
+	console.i("input b: " .. tostring(m))
+	local out = mat4()
+	local step = 1
+	for i=0, 12, 4 do
+		for j=1, 4 do
+			-- if type(m[j]) == "table" then
+			-- 	error("Whoa, this is broken!")
+			-- end
+			out[i+j] = m[j] * self[i+1] + m[j+4] * self[i+2] + m[j+8] * self[i+3] + m[j+12] * self[i+4]
+			console.i(string.format("step %d %s", step, out))
+			step = step + 1
+		end
+	end
+	console.i(out)
+	return out
+end
+
+function mat4:to_vec4s()
+	return {
+		{ self[1], self[2], self[3], self[4] },
+		{ self[5], self[6], self[7], self[8] },
+		{ self[9], self[10], self[11], self[12] },
+		{ self[13], self[14], self[15], self[16] }
+	}
+end
+
+return mat4

cpml/quat.lua

@@ -0,0 +1,500 @@
+-- quaternions
+-- Author: Andrew Stacey
+-- Website: http://www.math.ntnu.no/~stacey/HowDidIDoThat/iPad/Codea.html
+-- Licence: CC0 http://wiki.creativecommons.org/CC0
+
+--[[
+This is a class for handling quaternion numbers.  It was originally
+designed as a way of encoding rotations of 3 dimensional space.
+--]]
+
+local Class = require "libs.hump.class"
+local quaternion = Class {}
+local FLT_EPSILON = 1.19209290e-07
+
+--[[
+A quaternion can either be specified by giving the four coordinates as
+real numbers or by giving the scalar part and the vector part.
+--]]
+
+
+function quaternion:init(...)
+	-- copy
+	local arg = {...}
+	if #arg == 1 and type(arg[1]) == "table" then
+		self.a = arg[1].a
+		self.b = arg[1].b
+		self.c = arg[1].c
+		self.d = arg[1].d
+	-- four numbers
+	elseif #arg == 4 then
+		self.a = arg[1]
+		self.b = arg[2]
+		self.c = arg[3]
+		self.d = arg[4]
+	-- real number plus vector
+	elseif #arg == 2 then
+		self.a = arg[1]
+		self.b = arg[2].x
+		self.c = arg[2].y
+		self.d = arg[2].z
+	else
+		error("Incorrect number of arguments to quaternion")
+	end
+end
+
+
+function quaternion:to_axis_angle()
+	local tmp = self
+	if tmp.a > 1 then
+		tmp = tmp:normalize()
+	end
+	local angle = 2 * math.acos(tmp.a)
+	local s = math.sqrt(1-tmp.a*tmp.a)
+	local x, y, z
+	if s < FLT_EPSILON then
+		x = tmp.b
+		y = tmp.c
+		z = tmp.d
+	else
+		x = tmp.b / s -- normalize axis
+		y = tmp.c / s
+		z = tmp.d / s
+	end
+	return angle, { x, y, z }
+end
+
+
+--[[
+Test if we are zero.
+--]]
+
+function quaternion:is_zero()
+	-- are we the zero vector
+	if self.a ~= 0 or self.b ~= 0 or self.c ~= 0 or self.d ~= 0 then
+		return false
+	end
+	return true
+end
+
+--[[
+Test if we are real.
+--]]
+
+function quaternion:is_real()
+	-- are we the zero vector
+	if self.b ~= 0 or self.c ~= 0 or self.d ~= 0 then
+		return false
+	end
+	return true
+end
+
+--[[
+Test if the real part is zero.
+--]]
+
+function quaternion:is_imaginary()
+	-- are we the zero vector
+	if self.a ~= 0 then
+		return false
+	end
+	return true
+end
+
+--[[
+Test for equality.
+--]]
+
+function quaternion:is_eq(q)
+	if self.a ~= q.a or self.b ~= q.b or self.c ~= q.c or self.d ~= q.d then
+		return false
+	end
+	return true
+end
+
+--[[
+Defines the "==" shortcut.
+--]]
+
+function quaternion:__eq(q)
+	return self:is_eq(q)
+end
+
+--[[
+The inner product of two quaternions.
+--]]
+
+function quaternion:dot(q)
+	return self.a * q.a + self.b * q.b + self.c * q.c + self.d * q.d
+end
+
+--[[
+Makes "q .. p" return the inner product.
+
+Probably a bad choice and likely to be removed in future versions.
+--]]
+
+function quaternion:__concat(q)
+	return self:dot(q)
+end
+
+--[[
+Length of a quaternion.
+--]]
+
+function quaternion:len()
+	return math.sqrt(self:lensq())
+end
+
+--[[
+Often enough to know the length squared, which is quicker.
+--]]
+
+function quaternion:lensq()
+	return self.a * self.a + self.b * self.b + self.c * self.c + self.d * self.d
+end
+
+--[[
+Normalize a quaternion to have length 1, if possible.
+--]]
+
+function quaternion:normalize()
+	if self:is_zero() then
+		error("Unable to normalize a zero-length quaternion")
+		return false
+	end
+	local l = 1/self:len()
+	return self:scale(l)
+end
+
+--[[
+Scale the quaternion.
+--]]
+
+function quaternion:scale(l)
+	return quaternion(self.a * l,self.b * l,self.c * l, self.d * l)
+end
+
+--[[
+Add two quaternions.  Or add a real number to a quaternion.
+--]]
+
+function quaternion:add(q)
+	if type(q) == "number" then
+		return quaternion(self.a + q, self.b, self.c, self.d)
+	else
+		return quaternion(self.a + q.a, self.b + q.b, self.c + q.c, self.d + q.d)
+	end
+end
+
+--[[
+q + p
+--]]
+
+function quaternion:__add(q)
+	return self:add(q)
+end
+
+--[[
+Subtraction
+--]]
+
+function quaternion:subtract(q)
+	return quaternion(self.a - q.a, self.b - q.b, self.c - q.c, self.d - q.d)
+end
+
+--[[
+q - p
+--]]
+
+function quaternion:__sub(q)
+	return self:subtract(q)
+end
+
+--[[
+Negation (-q)
+--]]
+
+function quaternion:__unm()
+	return self:scale(-1)
+end
+
+--[[
+Length (#q)
+--]]
+
+function quaternion:__len()
+	return self:len()
+end
+
+--[[
+Multiply the current quaternion on the right.
+
+Corresponds to composition of rotations.
+--]]
+
+function quaternion:multiplyRight(q)
+	local a,b,c,d
+	a = self.a * q.a - self.b * q.b - self.c * q.c - self.d * q.d
+	b = self.a * q.b + self.b * q.a + self.c * q.d - self.d * q.c
+	c = self.a * q.c - self.b * q.d + self.c * q.a + self.d * q.b
+	d = self.a * q.d + self.b * q.c - self.c * q.b + self.d * q.a
+	return quaternion(a,b,c,d)
+end
+
+--[[
+q * p
+--]]
+
+function quaternion:__mul(q)
+	if type(q) == "number" then
+		return self:scale(q)
+	elseif type(q) == "table" then
+		return self:multiplyRight(q)
+	end
+end
+
+--[[
+Multiply the current quaternion on the left.
+
+Corresponds to composition of rotations.
+--]]
+
+function quaternion:multiplyLeft(q)
+	return q:multiplyRight(self)
+end
+
+--[[
+Conjugation (corresponds to inverting a rotation).
+--]]
+
+function quaternion:conjugate()
+	return quaternion(self.a, - self.b, - self.c, - self.d)
+end
+
+function quaternion:co()
+	return self:conjugate()
+end
+
+--[[
+Reciprocal: 1/q
+--]]
+
+function quaternion:reciprocal()
+	if self.is_zero() then
+		error("Cannot reciprocate a zero quaternion")
+		return false
+	end
+	local q = self:conjugate()
+	local l = self:lensq()
+	q = q:scale(1/l)
+	return q
+end
+
+--[[
+Integral powers.
+--]]
+
+function quaternion:power(n)
+	if n ~= math.floor(n) then
+		error("Only able to do integer powers")
+		return false
+	end
+	if n == 0 then
+		return quaternion(1,0,0,0)
+	elseif n > 0 then
+		return self:multiplyRight(self:power(n-1))
+	elseif n < 0 then
+		return self:reciprocal():power(-n)
+	end
+end
+
+--[[
+q^n
+
+This is overloaded so that a non-number exponent returns the
+conjugate.  This means that one can write things like q^* or q^"" to
+get the conjugate of a quaternion.
+--]]
+
+function quaternion:__pow(n)
+	if type(n) == "number" then
+		return self:power(n)
+	else
+		return self:conjugate()
+	end
+end
+
+--[[
+Division: q/p
+--]]
+
+function quaternion:__div(q)
+	if type(q) == "number" then
+		return self:scale(1/q)
+	elseif type(q) == "table" then
+		return self:multiplyRight(q:reciprocal())
+	end
+end
+
+--[[
+Returns the real part.
+--]]
+
+function quaternion:real()
+	return self.a
+end
+
+--[[
+Returns the vector (imaginary) part as a Vec3 object.
+--]]
+
+function quaternion:vector()
+	return Vec3(self.b, self.c, self.d)
+end
+
+--[[
+Represents a quaternion as a string.
+--]]
+
+function quaternion:__tostring()
+	local s
+	local im ={{self.b,"i"},{self.c,"j"},{self.d,"k"}}
+	if self.a ~= 0 then
+		s = self.a
+	end
+	for k,v in pairs(im) do
+	if v[1] ~= 0 then
+		if s then 
+			if v[1] > 0 then
+				if v[1] == 1 then
+					s = s .. " + " .. v[2]
+				else
+					s = s .. " + " .. v[1] .. v[2]
+				end
+			else
+				if v[1] == -1 then
+					s = s .. " - " .. v[2]
+				else
+					s = s .. " - " .. (-v[1]) .. v[2]
+				end
+				end
+		else
+			if v[1] == 1 then
+				s = v[2]
+			elseif v[1] == - 1 then
+				s = "-" .. v[2]
+			else
+				s = v[1] .. v[2]
+			end
+		end
+	end
+	end
+	if s then
+		return s
+	else
+		return "0"
+	end
+end
+
+--[[
+(Not a class function)
+
+Returns a quaternion corresponding to the current gravitational vector
+so that after applying the corresponding rotation, the y-axis points
+in the gravitational direction and the x-axis is in the plane of the
+iPad screen.
+
+When we have access to the compass, the x-axis behaviour might change.
+--]]
+
+--[[
+function quaternion.gravity()
+	local gxy, gy, gygxy, a, b, c, d
+	if Gravity.x == 0 and Gravity.y == 0 then
+		return quaternion(1,0,0,0)
+	else
+		gy = - Gravity.y
+		gxy = math.sqrt(math.pow(Gravity.x,2) + math.pow(Gravity.y,2))
+		gygxy = gy/gxy
+		a = math.sqrt(1 + gxy - gygxy - gy)/2
+		b = math.sqrt(1 - gxy - gygxy + gy)/2
+		c = math.sqrt(1 - gxy + gygxy - gy)/2
+		d = math.sqrt(1 + gxy + gygxy + gy)/2
+		if Gravity.y > 0 then
+			a = a
+			b = b
+		end
+		if Gravity.z < 0 then
+			b = - b
+			c = - c
+		end
+		if Gravity.x > 0 then
+			c = - c
+			d = - d
+		end
+		return quaternion(a,b,c,d)
+	end
+end
+--]]
+
+--[[
+Converts a rotation to a quaternion.  The first argument is the angle
+to rotate, the rest must specify an axis, either as a Vec3 object or
+as three numbers.
+--]]
+
+function quaternion.rotation(a,...)
+	local q,c,s
+	q = quaternion(0,...)
+	q = q:normalize()
+	c = math.cos(a/2)
+	s = math.sin(a/2)
+	q = q:scale(s)
+	q = q:add(c)
+	return q
+end
+
+function quaternion:to_euler()
+	local sqw = self.a*self.a
+	local sqx = self.b*self.b
+	local sqy = self.c*self.c
+	local sqz = self.d*self.d
+
+	 -- if normalised is one, otherwise is correction factor
+	local unit = sqx + sqy + sqz + sqw
+	local test = self.b*self.c + self.d*self.a
+
+	local pitch, yaw, roll
+
+	 -- singularity at north pole
+	if test > 0.499*unit then
+		yaw = 2 * math.atan2(self.b,self.a)
+		pitch = math.pi/2
+		roll = 0
+		return pitch, yaw, roll
+	end
+
+	 -- singularity at south pole
+	if test < -0.499*unit then
+		yaw = -2 * math.atan2(self.b,self.a)
+		pitch = -math.pi/2
+		roll = 0
+		return pitch, yaw, roll
+	end
+	yaw = math.atan2(2*self.c*self.a-2*self.b*self.d , sqx - sqy - sqz + sqw)
+	pitch = math.asin(2*test/unit)
+	roll = math.atan2(2*self.b*self.a-2*self.c*self.d , -sqx + sqy - sqz + sqw)
+
+	return pitch, roll, yaw
+end
+
+--[[
+The unit quaternion.
+--]]
+
+function quaternion.unit()
+	return quaternion(1,0,0,0)
+end
+
+return quaternion

cpml/simplex.lua

@@ -0,0 +1,349 @@
+--
+-- Based on code in "Simplex noise demystified", by Stefan Gustavson
+-- www.itn.liu.se/~stegu/simplexnoise/simplexnoise.pdf
+--
+-- Thanks to Mike Pall for some cleanup and improvements (and for LuaJIT!)
+--
+-- Permission is hereby granted, free of charge, to any person obtaining
+-- a copy of this software and associated documentation files (the
+-- "Software"), to deal in the Software without restriction, including
+-- without limitation the rights to use, copy, modify, merge, publish,
+-- distribute, sublicense, and/or sell copies of the Software, and to
+-- permit persons to whom the Software is furnished to do so, subject to
+-- the following conditions:
+--
+-- The above copyright notice and this permission notice shall be
+-- included in all copies or substantial portions of the Software.
+--
+-- THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+-- EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+-- MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+-- IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+-- CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+-- TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+-- SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+--
+-- [ MIT license: http://www.opensource.org/licenses/mit-license.php ]
+--
+
+-- Modules --
+local bit = require("bit")
+local ffi = require("ffi")
+local math = require("math")
+
+-- Imports --
+local band = bit.band
+local bor = bit.bor
+local bxor = bit.bxor
+local floor = math.floor
+local lshift = bit.lshift
+local max = math.max
+local rshift = bit.rshift
+
+-- Module table --
+local M = {}
+
+-- Permutation of 0-255, replicated to allow easy indexing with sums of two bytes --
+local Perms = ffi.new("uint8_t[512]", {
+	151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7, 225,
+	140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148,
+	247, 120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32,
+	57, 177, 33, 88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68,	175,
+	74, 165, 71, 134, 139, 48, 27, 166, 77, 146, 158, 231, 83, 111,	229, 122,
+	60, 211, 133, 230, 220, 105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54,
+	65, 25, 63, 161, 1, 216, 80, 73, 209, 76, 132, 187, 208, 89, 18, 169,
+	200, 196, 135, 130, 116, 188, 159, 86, 164, 100, 109, 198, 173, 186, 3, 64,
+	52, 217, 226, 250, 124, 123, 5, 202, 38, 147, 118, 126, 255, 82, 85, 212,
+	207, 206, 59, 227, 47, 16, 58, 17, 182, 189, 28, 42, 223, 183, 170, 213,
+	119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101, 155, 167, 43, 172, 9,
+	129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232, 178, 185, 112, 104,
+	218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12, 191, 179, 162, 241,
+	81,	51, 145, 235, 249, 14, 239,	107, 49, 192, 214, 31, 181, 199, 106, 157,
+	184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236, 205, 93,
+	222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180
+})
+
+-- The above, mod 12 for each element --
+local Perms12 = ffi.new("uint8_t[512]")
+
+for i = 0, 255 do
+	local x = Perms[i] % 12
+
+	Perms[i + 256], Perms12[i], Perms12[i + 256] = Perms[i], x, x
+end
+
+-- Gradients for 2D, 3D case --
+local Grads3 = ffi.new("const double[12][3]",
+	{ 1, 1, 0 }, { -1, 1, 0 }, { 1, -1, 0 }, { -1, -1, 0 },
+	{ 1, 0, 1 }, { -1, 0, 1 }, { 1, 0, -1 }, { -1, 0, -1 },
+	{ 0, 1, 1 }, { 0, -1, 1 }, { 0, 1, -1 }, { 0, -1, -1 }
+)
+
+do
+	-- 2D weight contribution
+	local function GetN (bx, by, x, y)
+		local t = .5 - x * x - y * y
+		local index = Perms12[bx + Perms[by]]
+
+		return max(0, (t * t) * (t * t)) * (Grads3[index][0] * x + Grads3[index][1] * y)
+	end
+
+	---
+	-- @param x
+	-- @param y
+	-- @return Noise value in the range [-1, +1]
+	function M.Simplex2D (x, y)
+		--[[
+			2D skew factors:
+			F = (math.sqrt(3) - 1) / 2
+			G = (3 - math.sqrt(3)) / 6
+			G2 = 2 * G - 1
+		]]
+
+		-- Skew the input space to determine which simplex cell we are in.
+		local s = (x + y) * 0.366025403 -- F
+		local ix, iy = floor(x + s), floor(y + s)
+
+		-- Unskew the cell origin back to (x, y) space.
+		local t = (ix + iy) * 0.211324865 -- G
+		local x0 = x + t - ix
+		local y0 = y + t - iy
+
+		-- Calculate the contribution from the two fixed corners.
+		-- A step of (1,0) in (i,j) means a step of (1-G,-G) in (x,y), and
+		-- A step of (0,1) in (i,j) means a step of (-G,1-G) in (x,y).
+		ix, iy = band(ix, 255), band(iy, 255)
+
+		local n0 = GetN(ix, iy, x0, y0)
+		local n2 = GetN(ix + 1, iy + 1, x0 - 0.577350270, y0 - 0.577350270) -- G2
+
+		--[[
+			Determine other corner based on simplex (equilateral triangle) we are in:
+			if x0 > y0 then
+				ix, x1 = ix + 1, x1 - 1
+			else
+				iy, y1 = iy + 1, y1 - 1
+			end
+		]]
+		local xi = rshift(floor(y0 - x0), 31) -- x0 >= y0
+		local n1 = GetN(ix + xi, iy + (1 - xi), x0 + 0.211324865 - xi, y0 - 0.788675135 + xi) -- x0 + G - xi, y0 + G - (1 - xi)
+
+		-- Add contributions from each corner to get the final noise value.
+		-- The result is scaled to return values in the interval [-1,1].
+		return 70 * (n0 + n1 + n2)
+	end
+end
+
+do
+	-- 3D weight contribution
+	local function GetN (ix, iy, iz, x, y, z)
+		local t = .6 - x * x - y * y - z * z
+		local index = Perms12[ix + Perms[iy + Perms[iz]]]
+
+		return max(0, (t * t) * (t * t)) * (Grads3[index][0] * x + Grads3[index][1] * y + Grads3[index][2] * z)
+	end
+
+	---
+	-- @param x
+	-- @param y
+	-- @param z
+	-- @return Noise value in the range [-1, +1]
+	function M.Simplex3D (x, y, z)
+		--[[
+			3D skew factors:
+			F = 1 / 3
+			G = 1 / 6
+			G2 = 2 * G
+			G3 = 3 * G - 1
+		]]
+
+		-- Skew the input space to determine which simplex cell we are in.
+		local s = (x + y + z) * 0.333333333 -- F
+		local ix, iy, iz = floor(x + s), floor(y + s), floor(z + s)
+
+		-- Unskew the cell origin back to (x, y, z) space.
+		local t = (ix + iy + iz) * 0.166666667 -- G
+		local x0 = x + t - ix
+		local y0 = y + t - iy
+		local z0 = z + t - iz
+
+		-- Calculate the contribution from the two fixed corners.
+		-- A step of (1,0,0) in (i,j,k) means a step of (1-G,-G,-G) in (x,y,z);
+		-- a step of (0,1,0) in (i,j,k) means a step of (-G,1-G,-G) in (x,y,z);
+		-- a step of (0,0,1) in (i,j,k) means a step of (-G,-G,1-G) in (x,y,z).
+		ix, iy, iz = band(ix, 255), band(iy, 255), band(iz, 255)
+
+		local n0 = GetN(ix, iy, iz, x0, y0, z0)
+		local n3 = GetN(ix + 1, iy + 1, iz + 1, x0 - 0.5, y0 - 0.5, z0 - 0.5) -- G3
+
+		--[[
+			Determine other corners based on simplex (skewed tetrahedron) we are in:
+			local ix2, iy2, iz2 = ix, iy, iz
+
+			if x0 >= y0 then
+				ix2, x2 = ix + 1, x2 - 1
+
+				if y0 >= z0 then -- X Y Z
+					ix, iy2, x1, y2 = ix + 1, iy + 1, x1 - 1, y2 - 1
+				elseif x0 >= z0 then -- X Z Y
+					ix, iz2, x1, z2 = ix + 1, iz + 1, x1 - 1, z2 - 1
+				else -- Z X Y
+					iz, iz2, z1, z2 = iz + 1, iz + 1, z1 - 1, z2 - 1
+				end
+			else
+				iy2, y2 = iy + 1, y2 - 1
+
+				if y0 < z0 then -- Z Y X
+					iz, iz2, z1, z2 = iz + 1, iz + 1, z1 - 1, z2 - 1
+				elseif x0 < z0 then -- Y Z X
+					iy, iz2, y1, z2 = iy + 1, iz + 1, y1 - 1, z2 - 1
+				else -- Y X Z
+					iy, ix2, y1, x2 = iy + 1, ix + 1, y1 - 1, x2 - 1
+				end
+			end		
+		]]
+		local yx = rshift(floor(y0 - x0), 31) -- x0 >= y0
+		local zy = rshift(floor(z0 - y0), 31) -- y0 >= z0
+		local zx = rshift(floor(z0 - x0), 31) -- x0 >= z0
+
+		local i1 = band(yx, bor(zy, zx)) -- x >= y and (y >= z or x >= z)
+		local j1 = band(1 - yx, zy) -- x < y and y >= z
+		local k1 = band(1 - zy, 1 - band(yx, zx)) -- y < z and not (x >= y and x >= z)
+
+		local i2 = bor(yx, band(zy, zx)) -- x >= z or (y >= z and x >= z)
+		local j2 = bor(1 - yx, zy) -- x < y or y >= z
+		local k2 = bxor(yx, zy) -- (x >= y and y < z) xor (x < y and y >= z)
+
+		local n1 = GetN(ix + i1, iy + j1, iz + k1, x0 + 0.166666667 - i1, y0 + 0.166666667 - j1, z0 + 0.166666667 - k1) -- G
+		local n2 = GetN(ix + i2, iy + j2, iz + k2, x0 + 0.333333333 - i2, y0 + 0.333333333 - j2, z0 + 0.333333333 - k2) -- G2
+
+		-- Add contributions from each corner to get the final noise value.
+		-- The result is scaled to stay just inside [-1,1]
+		return 32 * (n0 + n1 + n2 + n3)
+	end
+end
+
+do
+	-- Gradients for 4D case --
+	local Grads4 = ffi.new("const double[32][4]",
+		{ 0, 1, 1, 1 }, { 0, 1, 1, -1 }, { 0, 1, -1, 1 }, { 0, 1, -1, -1 },
+		{ 0, -1, 1, 1 }, { 0, -1, 1, -1 }, { 0, -1, -1, 1 }, { 0, -1, -1, -1 },
+		{ 1, 0, 1, 1 }, { 1, 0, 1, -1 }, { 1, 0, -1, 1 }, { 1, 0, -1, -1 },
+		{ -1, 0, 1, 1 }, { -1, 0, 1, -1 }, { -1, 0, -1, 1 }, { -1, 0, -1, -1 },
+		{ 1, 1, 0, 1 }, { 1, 1, 0, -1 }, { 1, -1, 0, 1 }, { 1, -1, 0, -1 },
+		{ -1, 1, 0, 1 }, { -1, 1, 0, -1 }, { -1, -1, 0, 1 }, { -1, -1, 0, -1 },
+		{ 1, 1, 1, 0 }, { 1, 1, -1, 0 }, { 1, -1, 1, 0 }, { 1, -1, -1, 0 },
+		{ -1, 1, 1, 0 }, { -1, 1, -1, 0 }, { -1, -1, 1, 0 }, { -1, -1, -1, 0 }
+	)
+
+	-- 4D weight contribution
+	local function GetN (ix, iy, iz, iw, x, y, z, w)
+		local t = .6 - x * x - y * y - z * z - w * w
+		local index = band(Perms[ix + Perms[iy + Perms[iz + Perms[iw]]]], 0x1F)
+
+		return max(0, (t * t) * (t * t)) * (Grads4[index][0] * x + Grads4[index][1] * y + Grads4[index][2] * z + Grads4[index][3] * w)
+	end
+
+	-- A lookup table to traverse the simplex around a given point in 4D.
+	-- Details can be found where this table is used, in the 4D noise method.
+	local Simplex = ffi.new("uint8_t[64][4]",
+		{ 0, 1, 2, 3 }, { 0, 1, 3, 2 }, {}, { 0, 2, 3, 1 }, {}, {}, {}, { 1, 2, 3 },
+		{ 0, 2, 1, 3 }, {}, { 0, 3, 1, 2 }, { 0, 3, 2, 1 }, {}, {}, {}, { 1, 3, 2 },
+		{}, {}, {}, {}, {}, {}, {}, {},
+		{ 1, 2, 0, 3 }, {}, { 1, 3, 0, 2 }, {}, {}, {}, { 2, 3, 0, 1 }, { 2, 3, 1 },
+		{ 1, 0, 2, 3 }, { 1, 0, 3, 2 }, {}, {}, {}, { 2, 0, 3, 1 }, {}, { 2, 1, 3 },
+		{}, {}, {}, {}, {}, {}, {}, {},
+		{ 2, 0, 1, 3 }, {}, {}, {}, { 3, 0, 1, 2 }, { 3, 0, 2, 1 }, {}, { 3, 1, 2 },
+		{ 2, 1, 0, 3 }, {}, {}, {}, { 3, 1, 0, 2 }, {}, { 3, 2, 0, 1 }, { 3, 2, 1 }
+	)
+
+	-- Convert the above indices to masks that can be shifted / anded into offsets --
+	for i = 0, 63 do
+		Simplex[i][0] = lshift(1, Simplex[i][0]) - 1
+		Simplex[i][1] = lshift(1, Simplex[i][1]) - 1
+		Simplex[i][2] = lshift(1, Simplex[i][2]) - 1
+		Simplex[i][3] = lshift(1, Simplex[i][3]) - 1
+	end
+
+	---
+	-- @param x
+	-- @param y
+	-- @param z
+	-- @param w
+	-- @return Noise value in the range [-1, +1]
+	function M.Simplex4D (x, y, z, w)
+		--[[
+			4D skew factors:
+			F = (math.sqrt(5) - 1) / 4 
+			G = (5 - math.sqrt(5)) / 20
+			G2 = 2 * G
+			G3 = 3 * G
+			G4 = 4 * G - 1
+		]]
+
+		-- Skew the input space to determine which simplex cell we are in.
+		local s = (x + y + z + w) * 0.309016994 -- F
+		local ix, iy, iz, iw = floor(x + s), floor(y + s), floor(z + s), floor(w + s)
+
+		-- Unskew the cell origin back to (x, y, z) space.
+		local t = (ix + iy + iz + iw) * 0.138196601 -- G
+		local x0 = x + t - ix
+		local y0 = y + t - iy
+		local z0 = z + t - iz
+		local w0 = w + t - iw
+
+		-- For the 4D case, the simplex is a 4D shape I won't even try to describe.
+		-- To find out which of the 24 possible simplices we're in, we need to
+		-- determine the magnitude ordering of x0, y0, z0 and w0.
+		-- The method below is a good way of finding the ordering of x,y,z,w and
+		-- then find the correct traversal order for the simplex we’re in.
+		-- First, six pair-wise comparisons are performed between each possible pair
+		-- of the four coordinates, and the results are used to add up binary bits
+		-- for an integer index.
+		local c1 = band(rshift(floor(y0 - x0), 26), 32)
+		local c2 = band(rshift(floor(z0 - x0), 27), 16)
+		local c3 = band(rshift(floor(z0 - y0), 28), 8)
+		local c4 = band(rshift(floor(w0 - x0), 29), 4)
+		local c5 = band(rshift(floor(w0 - y0), 30), 2)
+		local c6 = rshift(floor(w0 - z0), 31)
+
+		-- Simplex[c] is a 4-vector with the numbers 0, 1, 2 and 3 in some order.
+		-- Many values of c will never occur, since e.g. x>y>z>w makes x<z, y<w and x<w
+		-- impossible. Only the 24 indices which have non-zero entries make any sense.
+		-- We use a thresholding to set the coordinates in turn from the largest magnitude.
+		local c = c1 + c2 + c3 + c4 + c5 + c6
+
+		-- The number 3 (i.e. bit 2) in the "simplex" array is at the position of the largest coordinate.
+		local i1 = rshift(Simplex[c][0], 2)
+		local j1 = rshift(Simplex[c][1], 2)
+		local k1 = rshift(Simplex[c][2], 2)
+		local l1 = rshift(Simplex[c][3], 2)
+
+		-- The number 2 (i.e. bit 1) in the "simplex" array is at the second largest coordinate.
+		local i2 = band(rshift(Simplex[c][0], 1), 1)
+		local j2 = band(rshift(Simplex[c][1], 1), 1)
+		local k2 = band(rshift(Simplex[c][2], 1), 1)
+		local l2 = band(rshift(Simplex[c][3], 1), 1)
+
+		-- The number 1 (i.e. bit 0) in the "simplex" array is at the second smallest coordinate.
+		local i3 = band(Simplex[c][0], 1)
+		local j3 = band(Simplex[c][1], 1)
+		local k3 = band(Simplex[c][2], 1)
+		local l3 = band(Simplex[c][3], 1)
+
+		-- Work out the hashed gradient indices of the five simplex corners
+		-- Sum up and scale the result to cover the range [-1,1]
+		ix, iy, iz, iw = band(ix, 255), band(iy, 255), band(iz, 255), band(iw, 255)
+
+		local n0 = GetN(ix, iy, iz, iw, x0, y0, z0, w0)
+		local n1 = GetN(ix + i1, iy + j1, iz + k1, iw + l1, x0 + 0.138196601 - i1, y0 + 0.138196601 - j1, z0 + 0.138196601 - k1, w0 + 0.138196601 - l1) -- G
+		local n2 = GetN(ix + i2, iy + j2, iz + k2, iw + l2, x0 + 0.276393202 - i2, y0 + 0.276393202 - j2, z0 + 0.276393202 - k2, w0 + 0.276393202 - l2) -- G2
+		local n3 = GetN(ix + i3, iy + j3, iz + k3, iw + l3, x0 + 0.414589803 - i3, y0 + 0.414589803 - j3, z0 + 0.414589803 - k3, w0 + 0.414589803 - l3) -- G3
+		local n4 = GetN(ix + 1, iy + 1, iz + 1, iw + 1, x0 - 0.447213595, y0 - 0.447213595, z0 - 0.447213595, w0 - 0.447213595) -- G4
+
+		return 27 * (n0 + n1 + n2 + n3 + n4)
+	end
+end
+
+-- Export the module.
+return M

cpml/vec3.lua

@@ -0,0 +1,228 @@
+--[[
+Copyright (c) 2010-2013 Matthias Richter
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+Except as contained in this notice, the name(s) of the above copyright holders
+shall not be used in advertising or otherwise to promote the sale, use or
+other dealings in this Software without prior written authorization.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+]]--
+
+-- Modified to include 3D capabilities by Bill Shillito, April 2014
+-- Various bug fixes by Colby Klein, October 2014
+
+local assert = assert
+local sqrt, cos, sin, atan2, acos = math.sqrt, math.cos, math.sin, math.atan2, math.acos
+
+local vector = {}
+vector.__index = vector
+
+local function new(x,y,z)
+	return setmetatable({x = x or 0, y = y or 0, z = z or 0}, vector)
+end
+local zero = new(0,0,0)
+
+local function isvector(v)
+	return type(v) == 'table' and type(v.x) == 'number' and type(v.y) == 'number' and type(v.z) == 'number'
+end
+
+function vector:clone()
+	return new(self.x, self.y, self.z)
+end
+
+function vector:unpack()
+	return self.x, self.y, self.z
+end
+
+function vector:__tostring()
+	return "("..tonumber(self.x)..","..tonumber(self.y)..","..tonumber(self.z)..")"
+end
+
+function vector.__unm(a)
+	return new(-a.x, -a.y, -a.z)
+end
+
+function vector.__add(a,b)
+	assert(isvector(a) and isvector(b), "Add: wrong argument types (<vector> expected)")
+	return new(a.x+b.x, a.y+b.y, a.z+b.z)
+end
+
+function vector.__sub(a,b)
+	assert(isvector(a) and isvector(b), "Sub: wrong argument types (<vector> expected)")
+	return new(a.x-b.x, a.y-b.y, a.z-b.z)
+end
+
+function vector.__mul(a,b)
+	if type(a) == "number" then
+		return new(a*b.x, a*b.y, a*b.z)
+	elseif type(b) == "number" then
+		return new(b*a.x, b*a.y, b*a.z)
+	else
+		-- This is the dot product.
+		assert(isvector(a) and isvector(b), "Mul: wrong argument types (<vector> or <number> expected)")
+		return a.x*b.x + a.y*b.y + a.z*b.z
+	end
+end
+
+function vector.__div(a,b)
+	assert(isvector(a) and type(b) == "number", "wrong argument types (expected <vector> / <number>)")
+	return new(a.x / b, a.y / b, a.z / b)
+end
+
+function vector.__eq(a,b)
+	return a.x == b.x and a.y == b.y and a.z == b.z
+end
+
+function vector.__lt(a,b)
+	-- This is a lexicographical order.
+	return a.x < b.x or (a.x == b.x and a.y < b.y) or (a.x == b.x and a.y == b.y and a.z < b.z)
+end
+
+function vector.__le(a,b)
+	-- This is a lexicographical order.
+	return a.x <= b.x and a.y <= b.y and a.z <= b.z
+end
+
+function vector.permul(a,b)
+	assert(isvector(a) and isvector(b), "permul: wrong argument types (<vector> expected)")
+	return new(a.x*b.x, a.y*b.y, a.z*b*z)
+end
+
+function vector:tuple()
+	return self.x, self.y, self.z
+end
+
+function vector:len2()
+	return self.x * self.x + self.y * self.y + self.z * self.z
+end
+
+function vector:len()
+	return sqrt(self.x * self.x + self.y * self.y + self.z * self.z)
+end
+
+function vector.dist(a, b)
+	assert(isvector(a) and isvector(b), "dist: wrong argument types (<vector> expected)")
+	local dx = a.x - b.x
+	local dy = a.y - b.y
+	local dz = a.z - b.z
+	return sqrt(dx * dx + dy * dy + dz * dz)
+end
+
+function vector.dist2(a, b)
+	assert(isvector(a) and isvector(b), "dist: wrong argument types (<vector> expected)")
+	local dx = a.x - b.x
+	local dy = a.y - b.y
+	local dz = a.z - b.z
+	return (dx * dx + dy * dy + dz * dz)
+end
+
+function vector:normalize_inplace()
+	local l = self:len()
+	if l > 0 then
+		self.x, self.y, self.z = self.x / l, self.y / l, self.z / l
+	end
+	return self
+end
+
+function vector:normalize()
+	return self:clone():normalize_inplace()
+end
+
+function vector:rotated(phi, axis)
+	if axis == nil then return self end
+
+	local u = axis:normalize() or Vector(0,0,1) -- default is to rotate in the xy plane
+	local c, s = cos(phi), sin(phi)
+
+	-- Calculate generalized rotation matrix
+	local m1 = new((c + u.x * u.x * (1-c)),       (u.x * u.y * (1-c) - u.z * s), (u.x * u.z * (1-c) + u.y * s))
+	local m2 = new((u.y * u.x * (1-c) + u.z * s), (c + u.y * u.y * (1-c)),       (u.y * u.z * (1-c) - u.x * s))
+	local m3 = new((u.z * u.x * (1-c) - u.y * s), (u.z * u.y * (1-c) + u.x * s), (c + u.z * u.z * (1-c))      )
+
+	-- Return rotated vector
+	return new( m1 * self, m2 * self, m3 * self )
+end
+
+function vector:rotate_inplace(phi, axis)
+	self = self:rotated(phi, axis)
+end
+
+function vector:perpendicular()
+	return new(-self.y, self.x, 0)
+end
+
+function vector:projectOn(v)
+	assert(isvector(v), "invalid argument: cannot project vector on " .. type(v))
+	-- (self * v) * v / v:len2()
+	local s = (self.x * v.x + self.y * v.y + self.z * v.z) / (v.x * v.x + v.y * v.y + v.z * v.z)
+	return new(s * v.x, s * v.y, s * v.z)
+end
+
+function vector:projectFrom(v)
+	assert(isvector(v), "invalid argument: cannot project vector on " .. type(v))
+	-- Does the reverse of projectOn.
+	local s = (v.x * v.x + v.y * v.y + v.z * v.z) / (self.x * v.x + self.y * v.y + self.z * v.z)
+	return new(s * v.x, s * v.y, s * v.z)
+end
+
+function vector:mirrorOn(v)
+	assert(isvector(v), "invalid argument: cannot mirror vector on " .. type(v))
+	-- 2 * self:projectOn(v) - self
+	local s = 2 * (self.x * v.x + self.y * v.y + self.z * v.z) / (v.x * v.x + v.y * v.y + v.z * v.z)
+	return new(s * v.x - self.x, s * v.y - self.y, s * v.z - self.z)
+end
+
+function vector:cross(v)
+	-- Cross product.
+	assert(isvector(v), "cross: wrong argument types (<vector> expected)")
+	return new(self.y*v.z - self.z*v.y, self.z*v.x - self.x*v.z, self.x*v.y - self.y*v.x)
+	--return self.x * v.y - self.y * v.x
+end
+
+-- ref.: http://blog.signalsondisplay.com/?p=336
+function vector:trim_inplace(maxLen)
+	local s = maxLen * maxLen / self:len2()
+	s = (s > 1 and 1) or math.sqrt(s)
+	self.x, self.y, self.z = self.x * s, self.y * s, self.z * s
+	return self
+end
+
+function vector:angleTo(other)
+	-- Only makes sense in 2D.
+	if other then
+		return atan2(self.y, self.x) - atan2(other.y, other.x)
+	end
+	return atan2(self.y, self.x)
+end
+
+function vector:angleBetween(other)
+	if other then
+		return acos(self*other / (self:len() * other:len()))
+	end
+	return 0
+end
+
+function vector:trimmed(maxLen)
+	return self:clone():trim_inplace(maxLen)
+end
+
+
+-- the module
+return setmetatable({new = new, isvector = isvector, zero = zero},
+{__call = function(_, ...) return new(...) end})

init.lua

@@ -0,0 +1,52 @@
+--[[
+                  .'@@@@@@@@@@@@@@#:                  
+              ,@@@@#;            .'@@@@+              
+           ,@@@'                      .@@@#           
+         +@@+            ....            .@@@         
+       ;@@;         '@@@@@@@@@@@@.          @@@       
+      @@#         @@@@@@@@++@@@@@@@;         `@@;     
+    .@@`         @@@@@#        #@@@@@          @@@    
+   `@@          @@@@@` Cirno's  `@@@@#          +@@   
+   @@          `@@@@@  Perfect   @@@@@           @@+  
+  @@+          ;@@@@+   Math     +@@@@+           @@  
+  @@           `@@@@@  Library   @@@@@@           #@' 
+ `@@            @@@@@@          @@@@@@@           `@@ 
+ :@@             #@@@@@@.    .@@@@@@@@@            @@ 
+ .@@               #@@@@@@@@@@@@;;@@@@@            @@ 
+  @@                  .;+@@#'.   ;@@@@@           :@@ 
+  @@`                            +@@@@+           @@. 
+  ,@@                            @@@@@           .@@  
+   @@#          ;;;;;.          `@@@@@           @@   
+    @@+         .@@@@@          @@@@@           @@`   
+     #@@         '@@@@@#`    ;@@@@@@          ;@@     
+      .@@'         @@@@@@@@@@@@@@@           @@#      
+        +@@'          '@@@@@@@;            @@@        
+          '@@@`                         '@@@          
+             #@@@;                  .@@@@:            
+                :@@@@@@@++;;;+#@@@@@@+`               
+                      .;'+++++;.                      
+--]]
+local current_folder = (...):gsub('%.init$', '') .. "."
+
+local cpml = {
+	_LICENSE = "CPML is distributed under the terms of the MIT license. See LICENSE.md.",
+	_URL = "https://github.com/shakesoda/cpml",
+	_VERSION = "0.0.1",
+	_DESCRIPTION = "Cirno's Perfect Math Library: Just about everything you need for 3D games. Hopefully."
+}
+
+local mat4 = require(current_folder .. "cpml.mat4")
+local vec3 = require(current_folder .. "cpml.vec3")
+local quat = require(current_folder .. "cpml.quat")
+local simplex = require(current_folder .. "cpml.simplex")
+local intersect = require(current_folder .. "cpml.intersect")
+local constants = require(current_folder .. "cpml.constants")
+
+cpml.mat4 = mat4
+cpml.vec3 = vec3
+cpml.quat = quat
+cpml.simplex = simplex
+cpml.intersect = intersect
+cpml.constants = constants
+
+return cpml