return function(R)
	require('./sys/math')(R)

	local S=math.sin
	local C=math.cos

	R.V4={
		EMPTY={0,0,0,0},
		ID={0,0,0,1},
	}
	
	R.M44={
		EMPTY={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}},
		ID={{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}},
	}
	

	-- V4 --

	-- e:
	-- 	{1,2,3,4}
	-- 	V4.EMPTY
	-- 	V4.IDENTITY
	function R.V4:new(e)
		local o={e={}}
		for i=1,4 do
			o.e[i]=e[i]
		end
		setmetatable(o,self)
		self.__index=self
		return o
	end

	function R.V4:clone()
		return R.V4:new(self.e)
	end

	function R.V4:add(v4)
		for x=1,4 do
			self.e[x]=self.e[x]+v4.e[x]
		end
	end

	function R.V4:addC(v4)
		local o=self:clone()
		o:add(v4)
		return o
	end

	function R.V4:sub(v4)
		for x=1,4 do
			self.e[x]=self.e[x]-v4.e[x]
		end
	end

	function R.V4:subC(v4)
		local o=self:clone()
		o:sub(v4)
		return o
	end

	function R.V4:mul(s)
		for x=1,4 do
			self.e[x]=self.e[x]*s
		end
	end

	function R.V4:mulC(s)
		local o=self:clone()
		o:mul(s)
		return o
	end

	function R.V4:normal()
		local l=self:length()
		-- Catch l==0
		l=(l<.001) and .001 or l
		self.e[1]=self.e[1]/l
		self.e[2]=self.e[2]/l
		self.e[3]=self.e[3]/l
	end

	function R.V4:normalC()
		local o=self:clone()
		o:normal()
		return o
	end

	function R.V4:length()
		return math.sqrt(
			self.e[1]*self.e[1]
			+self.e[2]*self.e[2]
			+self.e[3]*self.e[3]
		)
	end

	function R.V4:dot(v4)
		return
			self.e[1]*v4.e[1]
			+self.e[2]*v4.e[2]
			+self.e[3]*v4.e[3]
	end

	function R.V4:cross(v4)
		return R.V4:new({
			self.e[2]*v4.e[3]-self.e[3]*v4.e[2],
			self.e[3]*v4.e[1]-self.e[1]*v4.e[3],
			self.e[1]*v4.e[2]-self.e[2]*v4.e[1],
			1
		})
	end

	function R.V4:lerp(v4,t)
		return R.V4:new({
			R.lerp(t,self.e[1],v4.e[1]),
			R.lerp(t,self.e[2],v4.e[2]),
			R.lerp(t,self.e[3],v4.e[3]),
			1
		})
	end

	-- when working with V3 vectors for normals etc. and 
	-- we just want to rotate them
	function R.V4:mulM44RotE(m44)
		local out={}
		for y=1,3 do
			local acc=0
			for x=1,3 do
				acc=acc+self.e[x]*m44.e[y][x]
			end
			out[y]=acc
		end
		return out
	end

	function R.V4:mulM44Rot(m44,normalise)
		self.e=self:mulM44RotE(m44,normalise)
	end

	function R.V4:mulM44RotC(m44,normalise)
		return R.V4:new(self:mulM44RotE(m44,normalise))
	end

	function R.V4:mulM44E(m44,normalise)
		local out={}
		for y=1,4 do
			local acc=0
			for x=1,4 do
				acc=acc+self.e[x]*m44.e[y][x]
			end
			out[y]=acc
		end

		if normalise then
			if out[4]==0 then
				-- not sure if this is the right thing to do
				out[4]=0.00001
			end

			out[1]=out[1]/out[4]
			out[2]=out[2]/out[4]
			out[3]=out[3]/out[4]
		end

		return out
	end

	function R.V4:mulM44(m44,normalise)
		self.e=self:mulM44E(m44,normalise)
	end

	function R.V4:mulM44C(m44,normalise)
		return R.V4:new(self:mulM44E(m44,normalise))
	end

	function R.V4Normal(v1,v2,v3)
		local u=v2:subC(v1)
		local v=v3:subC(v1)
		local cross=u:cross(v)
		local l=cross:length()
		l=(l<.001) and .001 or l -- Catch l==0
		return R.V4:new({cross.e[1]/l,cross.e[2]/l,cross.e[3]/l,1})
	end

	-- M44 --

	function R.M44:new(e)
		local o={e={{},{},{},{}}}
		for y=1,4 do
			for x=1,4 do
				o.e[y][x]=e[y][x]
			end
		end
		setmetatable(o,self)
		self.__index=self
		return o
	end

	function R.M44:clone()
		return R.M44:new(self.e)
	end

	function R.M44:newRotX(a)
		local s,c=S(a),C(a)
		return R.M44:new({{1,0,0,0},{0,c,-s,0},{0,s,c,0},{0,0,0,1}})
	end

	function R.M44:newRotY(a)
		local s,c=S(a),C(a)
		return R.M44:new({
		{c,0,-s,0},
		{0,1,0,0},
		{s,0,c,0},
		{0,0,0,1}})
	end

	function R.M44:newRotZ(a)
		local s,c=S(a),C(a)
		return R.M44:new({{c,-s,0,0},{s,c,0,0},{0,0,1,0},{0,0,0,1}})
	end

	function R.M44:newScale(x,y,z)
		return R.M44:new({{x,0,0,0},{0,y,0,0},{0,0,z,0},{0,0,0,1}})
	end

	function R.M44:newTrans(x,y,z)
		return R.M44:new({{1,0,0,x},{0,1,0,y},{0,0,1,z},{0,0,0,1}})
	end

	-- Returns e[][]
	function R.M44:mulM44E(m44)
	local out={{},{},{},{}}
		for y=1,4 do
			for x=1,4 do
				local acc=0
				for i=1,4 do
					acc=acc+self.e[y][i]*m44.e[i][x]
				end
				out[y][x]=acc
			end
		end
		return out
	end

	function R.M44:mulM44(m44)
		self.e=self:mulM44E(m44)
	end

	function R.M44:mulM44C(m44)
		return R.M44:new(self:mulM44E(m44))
	end

	function R.NewM44LookAt(from,to,up)
		--[[
		-- https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/lookat-function/framing-lookat-function.html
		local forward=from:subC(to)
		forward:normal()
		local right=up:cross(forward)
		right:normal()
		local newUp=forward:cross(right)
		return R.M44:new({
			{right.e[1],newUp.e[1],forward.e[1],from.e[1]},
			{right.e[2],newUp.e[2],forward.e[2],from.e[2]},
			{right.e[3],newUp.e[3],forward.e[3],from.e[3]},
			{0,0,0,1}
		})
		--]]

		local z=(to:subC(from)):normalC()
		local x=up:cross(z):normalC()
		local y=z:cross(x):normalC()
		return R.M44:new({
			{x.e[1],x.e[2],x.e[3],-x:dot(from)},
			{y.e[1],y.e[2],y.e[3],-y:dot(from)},
			{z.e[1],z.e[2],z.e[3],-z:dot(from)},
			{0,0,0,1},
		})
	end

	function R.NewM44OrthProj(width,height,near,far)
		return R.M44:new({
			{1/width,0,0,0},
			{0,1/height,0,0},
			{0,0,-(2/(far-near)),-((far+near)/(far-near))},
			{0,0,0,1},
		})
	end

	-- viewAngle: 90
	-- near: .1
	-- far: 100
	-- https://www.scratchapixel.com/lessons/3d-basic-rendering/perspective-and-orthographic-projection-matrix/building-basic-perspective-projection-matrix.html
	function R.NewM44PersProj(viewAngle,near,far)
		local scale=1/math.tan(viewAngle * 0.5 * math.pi / 180);
		local p1=-(far/(near-far))
		local p2=-((far*near)/(far-near))

		return R.M44:new({
			{scale,0,0,0},
			{0,scale,0,0},
			{0,0,p1,-1},
			{0,0,p2,0},
		})
	end

	function R.M44:trace(name)
		if name==null then name="" end
		trace(string.format(name.."\n|%.2f, %.2f, %.2f, %.2f|\n|%.2f, %.2f, %.2f, %.2f|\n|%.2f, %.2f, %.2f, %.2f|\n|%.2f, %.2f, %.2f, %.2f|\n",
			self.e[1][1], self.e[1][2], self.e[1][3], self.e[1][4],
			self.e[2][1], self.e[2][2], self.e[2][3], self.e[2][4],
			self.e[3][1], self.e[3][2], self.e[3][3], self.e[3][4],
			self.e[4][1], self.e[4][2], self.e[4][3], self.e[4][4]))
	end
end
